Index: codes/icosagcm/devel/DySL/codegen
===================================================================
--- codes/icosagcm/devel/DySL/codegen	(revision 941)
+++ codes/icosagcm/devel/DySL/codegen	(revision 941)
@@ -0,0 +1,60 @@
+#!/bin/bash
+
+function filter()
+{
+    PASS=$1 ; shift
+    cpp -P -DPASS_$PASS $* -include ../macros.jin
+}
+
+function log()
+{
+    echo $*
+    $*
+}
+
+function extract_dysl()
+{
+    sed -n '/BEGIN_DYSL/,/END_DYSL/{//b;p}' compute_pvort_only.F90 
+}
+
+function main()
+{
+    VARIANT=$SRCDIR/$1
+    cd $VARIANT
+    rm -f *~ tmp/*
+    mkdir -p tmp
+    cd tmp
+    cp $SRCDIR/macros.jin . # jinja macros common to all variants
+    cpp -P -DPASS_JINJA ../macros.jin >> macros.jin # variant-specific macros
+
+    # first expand user-defined CPP and jinja macros
+    cat $SRCDIR/kernels_*.jin | cpp -P > caldyn_all.jin
+    $PYTHON -u $ROOT/jinja caldyn_all.jin > caldyn_all.dysl
+    # then expand the code-generating CPP and jinja macros
+    cat caldyn_all.dysl | filter PRE1 | filter PRE2 >> macros.jin
+    $PYTHON -u $ROOT/jinja macros.jin > caldyn_all.k90
+    # extract each kernel and apply post-jinja CPP macros
+    KERNELS=$(grep KERNEL caldyn_all.dysl | cpp -P -D'KERNEL(x)=x')
+    for NAME in $KERNELS ; do
+	cat caldyn_all.k90 | filter POST1 -DKERNEL_${NAME} | filter POST2 | sed '/^$/d' > ../${NAME}.k90
+	head ../${NAME}.k90
+    done
+
+    cd ..
+    grep ERROR *.k90
+
+    cd $VARIANT
+    wc *.k90
+    cd $ROOT
+}
+
+set +x
+    
+ROOT=$PWD
+SRCDIR=$ROOT
+grep FIXME *.jin
+
+for x in $* ; do
+    echo "Generating code for variant : $x"
+    main $x
+done
Index: codes/icosagcm/devel/DySL/hex_master/macros.jin
===================================================================
--- codes/icosagcm/devel/DySL/hex_master/macros.jin	(revision 941)
+++ codes/icosagcm/devel/DySL/hex_master/macros.jin	(revision 941)
@@ -0,0 +1,2 @@
+#define HEX_MASTER
+#include "../hexagonal/macros.jin"
Index: codes/icosagcm/devel/DySL/hexagonal/macros.jin
===================================================================
--- codes/icosagcm/devel/DySL/hexagonal/macros.jin	(revision 941)
+++ codes/icosagcm/devel/DySL/hexagonal/macros.jin	(revision 941)
@@ -0,0 +1,327 @@
+#ifdef PASS_PRE1
+
+#define _AND_ &&
+#define _OR_ ||
+#define _NOT_ !
+#define IS_INNER_LAYER (_NOT_ (IS_BOTTOM_LEVEL _OR_ IS_TOP_LAYER))
+#define IS_INNER_INTERFACE (_NOT_ (IS_BOTTOM_LEVEL _OR_ IS_TOP_INTERFACE))
+
+#define CELL ij,LL
+#define VERTEX1 ij+{{vertex1}},LL
+#define VERTEX2 ij+{{vertex2}},LL
+#define EDGE ij+{{edge}},LL
+#define EDGE_TRISK ij+{{edge_trisk}},LL
+#define DUAL_CELL ij+{{triangle}},LL
+#define CELL1 ij,LL
+#define CELL2 ij+{{neighbour}},LL
+#define VERTEX {{vertex}},LL
+#define SIGN {{sign}}
+#define LE_DE le_de(ij+{{edge}})
+#define WEE wee(ij+{{edge}},{{itrisk}})
+#define DE de(ij+{{edge}})
+#define LE le(ij+{{edge}})
+#define RIV2 Riv2({{vertex}},{{vertex_dir}})
+#define AI Ai(ij)
+#define FV fv(ij+{{triangle}})
+#define AV Av(ij+{{triangle}})
+
+#endif
+
+#ifdef PASS_PRE2
+
+#ifdef HEX_MASTER
+#define BARRIER
+#else
+#define BARRIER !$OMP BARRIER
+#endif
+
+#define KERNEL(name) {% call define_kernel(#name) %}
+#define FORALL_CELLS(...) {% call forall_cells(__VA_ARGS__) %}
+#define FORALL_CELLS_EXT(...) {% call forall_cells_ext(__VA_ARGS__) %}
+#define ON_PRIMAL {% call(forall_edges) on_hexagons() %}
+#define ON_DUAL {% call(triangle,forall_edges,forall_vertices) on_triangles() %}
+#define ON_EDGES {% call(edge,neighbour,sign,vertex1,vertex2,forall_trisk) on_edges() %}
+#define FORALL_EDGES {% call(edge,sign,vertex1,vertex2) forall_edges() %}
+#define FORALL_VERTICES {% call(vertex,vertex_dir) forall_vertices() %}
+#define FORALL_TRISK {% call(itrisk,edge_trisk) forall_trisk() %}
+#define END_BLOCK {% endcall %}
+#define SEQUENCE_C0 {% call(body,at_level,ext) sequence() %}
+#define SEQUENCE_C1 {% call(body,at_level,ext) sequence_ext() %}
+#define SEQUENCE_E0 {% call(body,at_level,ext) sequence_edge() %}
+#define PROLOGUE(level) {% call(edge,neighbour) at_level(level,ext) %}
+#define BODY(range){% call(edge,neighbour) body(range,ext) %}
+#define EPILOGUE(level){% call(edge,neighbour) at_level(level,ext) %}
+
+#define CST_IF(condition, action) {{ cst_if(#condition, #action) }}
+#define CST_IFTHEN(condition) {{flat}}if condition
+#define CST_ELSEIF(condition) {{flat}}elif condition
+#define CST_ELSE {{flat}}else
+#define CST_ENDIF {{flat}}endif
+
+#endif
+
+#if PASS_JINJA
+
+#ifdef HEX_MASTER
+
+{% set ll_begin, ll_beginp1 = '1', '2' %}  
+{% set ll_endm1, ll_end, ll_endp1 = 'llm-1', 'llm', 'llm+1' %}  
+{% macro ij_omp(ext) -%}
+ij_begin{{ext}},ij_end{{ext}}
+{%- endmacro %}
+
+#else
+
+{% set ll_begin, ll_beginp1 = 'll_begin', 'll_beginp1' %}  
+{% set ll_endm1, ll_end, ll_endp1 = 'll_endm1', 'll_end', 'll_endp1' %}  
+{% macro ij_omp(ext) -%}
+ij_omp_begin{{ext}},ij_omp_end{{ext}}
+{%- endmacro %}
+
+#endif
+
+{% macro sequence() -%}
+{{ caller(body_primal,at_level_primal,'') }}
+{%- endmacro %}
+
+{% macro sequence_ext() -%}
+{{ caller(body_primal,at_level_primal,'_ext') }}
+{%- endmacro %}
+
+{% macro at_level_primal(lev,ext) -%}
+{{ define('LL', lev) }}
+!DIR$ SIMD
+DO ij={{ ij_omp(ext) }}
+  {{ caller() }}
+END DO
+{%- endmacro %}
+
+{% macro body_primal(range,ext) -%}
+{{ define('LL', 'l') }}
+DO l = {{ range }}
+  !DIR$ SIMD
+  DO ij={{ ij_omp(ext) }}
+    {{ caller() }}
+  END DO
+END DO
+{%- endmacro %}
+
+{% macro sequence_edge() -%}
+{{ caller(body_edge,at_level_edge,'') }}
+{%- endmacro %}
+
+{% macro at_level_edge(lev,ext) -%}
+{{ define('LL', lev) }}
+!DIR$ SIMD
+DO ij={{ ij_omp(ext) }}
+      {{ caller('u_right','t_right') }}
+      {{ caller('u_lup','t_lup') }}
+      {{ caller('u_ldown','t_ldown') }}
+END DO
+{%- endmacro %}
+
+{% macro body_edge(range,ext) -%}
+{{ define('LL', 'l') }}
+DO l = {{ range }}
+   !DIR$ SIMD
+   DO ij={{ ij_omp(ext) }}
+      {{ caller('u_right','t_right') }}
+      {{ caller('u_lup','t_lup') }}
+      {{ caller('u_ldown','t_ldown') }}
+  END DO
+END DO
+{%- endmacro %}
+
+{#---------------------- big macro generating nested DO loops -------------------------#}
+{# specialized code is generated for l=1 and l=llm,llm+1 if necessary #}
+{# start/end values for l are replaced by corresponding OpenMP values #}
+
+{% macro define_levels(ll,kdown,kup,thecode) %}
+  {{ define('LL', ll) }}
+  {{ define_cppmacro('KDOWN','(ij,ll)', kdown) if 'KDOWN' in thecode }}
+  {{ define_cppmacro('KUP','(ij,ll)', kup) if 'KUP'   in thecode }}
+{%- endmacro %}
+
+{% macro forall(start,end,ij_range,thecode) %}
+{% set start_omp={'1':ll_begin , '2':ll_beginp1}[start] %}  
+{% set end_omp={'llm-1':ll_endm1 , 'llm':ll_end, 'llm+1':ll_endp1}[end] %}
+{% set is_top_layer, is_top_inter = False, False %}
+
+{{ define('IS_TOP_LAYER', '_FALSE_') }}
+{{ define('IS_TOP_INTERFACE', '_FALSE_') }}
+
+{% if 'IS_BOTTOM_LEVEL' in thecode or 'KDOWN' in thecode%}
+{{ 'ERROR : using IS_BOTTOM_LEVEL in a loop starting at l=2' if start=='2' }}
+{# the code in the loop checks whether l==1, for the sake of performance 
+we shall write special code for l=1 and start the loop at l=2 #}
+
+IF ({{ll_begin}}==1) THEN
+  {{ define('IS_BOTTOM_LEVEL', '_TRUE_') }}
+  {{ define_levels('1', 'ij,1', 'ij,1', thecode) }}
+    !DIR$ SIMD
+    DO ij={{ij_range}}
+      {{ thecode }}
+    END DO
+END IF
+{% set start_omp=ll_beginp1 %}
+{% endif %}
+
+{% if 'IS_TOP_LAYER' in thecode %}
+{{ 'ERROR : using _IS_TOP_LAYER_ in a loop ending at l=llm+1' if end=='llm+1' }}
+{# the code checks whether l==llm, write special code for l=llm and end the loop at l=llm-1 #}
+{% set end_omp, is_top_layer = ll_endm1, True %}
+{% endif %}
+
+{% if 'IS_TOP_INTERFACE' in thecode or 'KUP' in thecode %}
+{# the code checks whether l==llm+1, write special code for l=llm+1 and end the loop at l=llm #}
+{{ 'ERROR : using _IS_TOP_INTERFACE_ in a loop ending at l=llm' if end=='llm' }}
+{% set end_omp, is_top_inter = ll_end, True %}
+{% endif %}
+
+{{ define('IS_BOTTOM_LEVEL', '_FALSE_') }}
+{{ define_levels('l', 'ij,l-1', 'ij,l', thecode) }}
+DO l = {{start_omp}}, {{end_omp}}
+  !DIR$ SIMD
+  DO ij={{ij_range}}
+  {{ thecode }}
+  END DO
+END DO
+
+{% if is_top_layer %}
+IF({{ll_end}}==llm) THEN
+  {{ define('IS_TOP_LAYER', '_TRUE_') }}
+  {{ define_levels('llm', 'ij,llm-1', 'ij,llm', thecode) }}
+  !DIR$ SIMD
+  DO ij={{ij_range}}
+    {{ thecode }}
+  END DO
+END IF
+{% endif %}
+
+{% if is_top_inter %}
+IF({{ll_endp1}}==llm+1) THEN
+  {{ define('IS_TOP_INTERFACE', '_TRUE_') }}
+  {{ define_levels('llm+1', 'ij,llm', 'ij,llm', thecode) }}
+  !DIR$ SIMD
+  DO ij={{ij_range}}
+    {{ thecode }}
+  END DO
+END IF
+{% endif %}
+
+{% endmacro %}
+
+{#-------------------------------------------------------------------------------------#}
+
+{% macro forall_cells(start='1', end='llm') -%}
+{{ forall(start,end,'ij_begin, ij_end', caller()) }}
+{% endmacro %}
+
+{% macro forall_cells_ext(start='1', end='llm') -%}
+{{ forall(start,end,'ij_begin_ext, ij_end_ext', caller()) }}
+{% endmacro %}
+
+{% macro forall_edges_hex() -%}
+{{ caller('u_rup','ne_rup', 'z_rup','z_up') }}
+{{ caller('u_lup','ne_lup', 'z_lup','z_up') }}
+{{ caller('u_left','ne_left','z_lup','z_ldown') }}
+{{ caller('u_ldown','ne_ldown', 'z_ldown','z_down') }}
+{{ caller('u_rdown','ne_rdown', 'z_rdown','z_down') }}
+{{ caller('u_right','ne_right', 'z_rup','z_rdown') }}
+{%- endmacro %}
+
+{% macro on_hexagons() -%}
+{{ caller(forall_edges_hex) }}
+{%- endmacro %}
+
+{% macro forall_edges_zup() -%}
+{{ caller('u_rup','ne_rup', 'TODO','TODO') }}
+{{ caller('t_rup+u_left','ne_left', 'TODO','TODO') }}
+{{ caller('u_lup','(-ne_lup)', 'TODO','TODO') }}
+{%- endmacro %}
+
+{% macro forall_edges_zdown() -%}
+{{ caller('u_rdown','(-ne_rdown)', 'TODO', 'TODO') }}
+{{ caller('t_ldown+u_right','ne_right', 'TODO', 'TODO') }}
+{{ caller('u_ldown','ne_ldown', 'TODO', 'TODO') }}
+{%- endmacro %}
+
+{% macro forall_vertices_zup() -%}
+{{ caller('ij','vup') }}
+{{ caller('ij+t_rup','vldown') }}
+{{ caller('ij+t_lup','vrdown') }}
+{%- endmacro %}
+
+{% macro forall_vertices_zdown() -%}
+{{ caller('ij','vdown') }}
+{{ caller('ij+t_ldown','vrup') }}
+{{ caller('ij+t_rdown','vlup') }}
+{%- endmacro %}
+
+{% macro on_triangles() -%}
+{{ caller('z_up',forall_edges_zup,forall_vertices_zup) }}
+{{ caller('z_down',forall_edges_zdown,forall_vertices_zdown) }}
+{%- endmacro %}
+
+{% macro trisk_right() -%}
+{{ caller('1,1', 'u_rup') }}
+{{ caller('2,1', 'u_lup') }}
+{{ caller('3,1', 'u_left') }}
+{{ caller('4,1', 'u_ldown') }}
+{{ caller('5,1', 'u_rdown') }}
+{{ caller('1,2', 't_right+u_ldown') }}
+{{ caller('2,2', 't_right+u_rdown') }}
+{{ caller('3,2', 't_right+u_right') }}
+{{ caller('4,2', 't_right+u_rup') }}
+{{ caller('5,2', 't_right+u_lup') }}
+{%- endmacro %}
+
+{% macro trisk_lup() -%}
+{{ caller('1,1', 'u_left') }}
+{{ caller('2,1', 'u_ldown') }}
+{{ caller('3,1', 'u_rdown') }}
+{{ caller('4,1', 'u_right') }}
+{{ caller('5,1', 'u_rup') }}
+{{ caller('1,2', 't_lup+u_right') }}
+{{ caller('2,2', 't_lup+u_rup') }}
+{{ caller('3,2', 't_lup+u_lup') }}
+{{ caller('4,2', 't_lup+u_left') }}
+{{ caller('5,2', 't_lup+u_ldown') }}
+{%- endmacro %}
+
+{% macro trisk_ldown() -%}
+{{ caller('1,1', 'u_rdown') }}
+{{ caller('2,1', 'u_right') }}
+{{ caller('3,1', 'u_rup') }}
+{{ caller('4,1', 'u_lup') }}
+{{ caller('5,1', 'u_left') }}
+{{ caller('1,2', 't_ldown+u_lup') }}
+{{ caller('2,2', 't_ldown+u_left') }}
+{{ caller('3,2', 't_ldown+u_ldown') }}
+{{ caller('4,2', 't_ldown+u_rdown') }}
+{{ caller('5,2', 't_ldown+u_right') }}
+{%- endmacro %}
+
+{% macro on_edges() -%}
+{{ caller('u_right','t_right','ne_right','z_rdown','z_rup',trisk_right) }}
+{{ caller('u_lup','t_lup','ne_lup','z_up','z_lup', trisk_lup) }}
+{{ caller('u_ldown','t_ldown','ne_ldown','z_ldown','z_down', trisk_ldown) }}
+{%- endmacro %}
+
+{% set llm = 'llm' %}
+
+#endif
+
+#ifdef PASS_POST1
+#define _TRUE_ (0==0)
+#define _FALSE_ (0==1)
+#define CELL ij,LL
+#endif
+
+#ifdef PASS_POST2
+#define UP(ij,l) ij,l+1
+#define DOWN(ij,l) ij,l-1
+#define HIDX(ij,l) ij
+#define VIDX(ij,l) l
+#define AT_LEVEL(ij,l,lev) ij,lev
+#endif
Index: codes/icosagcm/devel/DySL/jinja
===================================================================
--- codes/icosagcm/devel/DySL/jinja	(revision 941)
+++ codes/icosagcm/devel/DySL/jinja	(revision 941)
@@ -0,0 +1,42 @@
+from jinja2 import Environment, FileSystemLoader
+from sys import argv
+from textwrap import dedent
+from re import search # regular expressions
+
+def pretty(line,indent):
+    line = dedent(line)
+    if search('^#', line): # check if line is a preprocessor directive
+       print line
+    else:
+       print (3*indent)*' ' + line
+    	
+env=Environment(loader=FileSystemLoader('.'))
+tpl=env.get_template(argv[1])
+source = tpl.render()
+source = [line for line in source.split('\n') if line.strip() != '']
+
+# regular expressions matching Fortran keywords
+# \b detects word boundaries
+# $ and ( must be escaped : \$, \(
+# THIS|THAT means THIS or THAT
+# raw Python strings r'blabla' avoid doubling '\'
+noop      = r'!\$OMP'
+increase  = r'\bTHEN\b|\bDO\b|\bSELECT CASE\b'
+inter     = r'\bELSE\b|\bCASE\(|\bCASE DEFAULT\b'
+decrease  = r'\bEND IF\b|\bENDIF\b|\bEND DO\b|\bENDDO\b|\bEND SELECT\b' 
+
+indent=1
+
+for line in source :
+    if search(noop, line) :
+        pretty(line,indent)
+    elif search(decrease, line) :
+        indent = indent-1
+        pretty(line,indent)
+    elif search(increase,line):
+        pretty(line,indent)
+        indent = indent+1
+    elif search(inter,line):
+        pretty(line,indent-1)
+    else:
+        pretty(line,indent)
Index: codes/icosagcm/devel/DySL/kernels_caldyn_NH.jin
===================================================================
--- codes/icosagcm/devel/DySL/kernels_caldyn_NH.jin	(revision 941)
+++ codes/icosagcm/devel/DySL/kernels_caldyn_NH.jin	(revision 941)
@@ -0,0 +1,63 @@
+{% macro compute_p_and_Aik() %}
+{% set p_and_c2_from_theta=caller %}
+    SEQUENCE_C1
+      BODY('1,llm')
+          rho_ij = (g*m_ik(CELL))/(Phi_il(UP(CELL))-Phi_il(CELL))
+          {{ p_and_c2_from_theta() }}
+          A_ik(CELL) = c2_mik*(tau/g*rho_ij)**2
+      END_BLOCK
+    END_BLOCK
+{%- endmacro %}
+
+KERNEL(caldyn_mil)
+  FORALL_CELLS_EXT('1', 'llm+1')
+    ON_PRIMAL
+       CST_IF(IS_BOTTOM_LEVEL,    m_il(CELL) = .5*rhodz(CELL)                     )
+       CST_IF(IS_INNER_INTERFACE, m_il(CELL) = .5*(rhodz(CELL)+rhodz(DOWN(CELL))) )
+       CST_IF(IS_TOP_INTERFACE,   m_il(CELL) = .5*rhodz(DOWN(CELL))               )
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+KERNEL(caldyn_vert_NH)
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      CST_IF(IS_BOTTOM_LEVEL,     w_ij = ( W(CELL)+.5*W(UP(CELL)) )/mass(CELL) )
+      CST_IF(IS_INNER_LAYER , w_ij = .5*( W(CELL)+W(UP(CELL)) )/mass(CELL) )
+      CST_IF(IS_TOP_LAYER,        w_ij = ( .5*W(CELL)+W(UP(CELL)) )/mass(CELL) )
+      wflux_ij = .5*(wflux(CELL)+wflux(UP(CELL)))
+      W_etadot(CELL) = wflux_ij*w_ij
+      eta_dot(CELL) = wflux_ij / mass(CELL)
+      wcov(CELL) = w_ij*(geopot(UP(CELL))-geopot(CELL))
+    END_BLOCK
+  END_BLOCK
+  ! add NH term to du
+  FORALL_CELLS()
+    ON_EDGES
+      du(EDGE) = du(EDGE) - .5*(wcov(CELL1)+wcov(CELL2))*SIGN*(eta_dot(CELL2)-eta_dot(CELL1))
+    END_BLOCK
+  END_BLOCK
+  ! add NH terms to dW, dPhi
+  ! FIXME : TODO top and bottom
+  FORALL_CELLS('2','llm')
+    ON_PRIMAL
+      dPhi(CELL)=dPhi(CELL)-wflux(CELL)*(geopot(UP(CELL))-geopot(DOWN(CELL)))/(mass(DOWN(CELL))+mass(CELL))
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we compute W_etadot above and do a vertical difference below
+  BARRIER
+
+  FORALL_CELLS('1','llm+1')
+    ON_PRIMAL
+      CST_IFTHEN(IS_BOTTOM_LEVEL)
+         dW(CELL) = dW(CELL) - W_etadot(CELL)
+      CST_ELSEIF(IS_TOP_INTERFACE)
+         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL))
+      CST_ELSE
+         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL)) - W_etadot(CELL)
+      CST_ENDIF
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
Index: codes/icosagcm/devel/DySL/kernels_transport.jin
===================================================================
--- codes/icosagcm/devel/DySL/kernels_transport.jin	(revision 941)
+++ codes/icosagcm/devel/DySL/kernels_transport.jin	(revision 941)
@@ -0,0 +1,349 @@
+KERNEL(advect_vert)
+
+  FORALL_CELLS_EXT('2','llm')
+    ON_PRIMAL
+      dzqw(CELL) = q(CELL)-q(DOWN(CELL))
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we compute dzqw above and do a vertical average below
+  BARRIER  
+
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      CST_IFTHEN(IS_BOTTOM_LEVEL || IS_TOP_LAYER)
+        dzq(CELL)=0.
+      CST_ELSE
+        dzq_down=dzqw(CELL)
+        dzq_up=dzqw(UP(CELL)
+        IF(dzq_up*dzq_down>0.) THEN
+           dzqm = .5*(dzqup+dzq_down)
+           dzqmax = max_slope * min( ABS(dzqw_down), ABS(dzq_up) )
+           dzq(CELL) = sign( min(abs(dzqm),dzqmax) , dzqm )  ! sign(a,b)=a*sign(b) 
+        ELSE
+          dzq(CELL)=0.
+        END IF
+      CST_ENDIF
+    END_BLOCK
+  END_BLOCK
+   
+  ! We need a barrier here because we compute dzq above and do a vertical average below
+  BARRIER
+
+  FORALL_CELLS_EXT('1','llm+1')
+    ON_PRIMAL
+      CST_IFTHEN(IS_BOTTOM_LEVEL || IS_TOP_INTERFACE)
+        wq(CELL)=0.
+      CST_ELSE
+        w = fac*wfluxt(CELL)
+	IF(w>0.) THEN  ! upward transport, upwind side is at level l
+	  sigw = w/mass(DOWN(CELL))
+	  qq   = q(DOWN(CELL) + 0.5*(1.-sigw)*dzq(DOWN(CELL)) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
+        ELSE           ! downward transport, upwind side is at level l+1
+          sigw = w/mass(CELL)
+          qq   = q(CELL)      -0.5*(1.+sigw)*dzq(ij,l) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0               
+	ENDIF
+	wq(ij,l) = w*qq																 
+      CST_ENDIF
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we compute wq above and do a vertical difference below
+  BARRIER
+
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      newmass = mass(CELL) + fac*(wflux(CELL)-wflux(UP(CELL)))
+      q(CELL) = ( q(CELL)*mass(CELL) + wq(CELL)-wq(UP(CELL)) ) / newmass
+      IF(update_mass) mass(CELL) = newmass
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
+
+KERNEL(advect_horiz)
+  ! evaluate tracer field at point cc using piecewise linear reconstruction
+  ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
+  ! sign*hfluxt>0 iff outgoing
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      IF(SIGN*hfluxt(EDGE)>0.) THEN
+        qe = qi(CELL1)
+        qe = qe + (cc(EDGE,1)-xyz_i(HIDX(CELL1),1))*gradq3d(CELL1,1) 
+        qe = qe + (cc(EDGE,2)-xyz_i(HIDX(CELL1),2))*gradq3d(CELL1,2) 
+        qe = qe + (cc(EDGE,3)-xyz_i(HIDX(CELL1),3))*gradq3d(CELL1,3) 
+      ELSE
+        qe = qi(CELL2)
+        qe = qe + (cc(EDGE,1)-xyz_i(HIDX(CELL2),1))*gradq3d(CELL2,1) 
+        qe = qe + (cc(EDGE,2)-xyz_i(HIDX(CELL2),2))*gradq3d(CELL2,2) 
+        qe = qe + (cc(EDGE,3)-xyz_i(HIDX(CELL2),3))*gradq3d(CELL2,3) 
+      END IF
+      qflux(EDGE) = hfluxt(EDGE)*qe
+      IF(diagflux_on) qfluxt(EDGE) = qfluxt(EDGE)+qflux(EDGE)
+    END_BLOCK
+  END_BLOCK
+
+  ! update q and, if update_mass, update mass
+  FORALL_CELLS()
+    ON_PRIMAL
+      dmass=0.
+      dq=0.
+      FORALL_EDGES
+        dmass = dmass + SIGN*hfluxt(EDGE)
+        dq = dq + SIGN*qflux(EDGE)
+      END_BLOCK
+      newmass = mass(CELL) - dmass/AI
+      qi(CELL) = (qi(CELL)*mass(CELL)-dq/AI) / newmass
+      IF(update_mass) mass(CELL)=newmass
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
+
+{% macro perot_reconstruction() %}
+  ! Perot reconstruction based on Gauss theorem
+  ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
+  FORALL_CELLS()
+    ON_PRIMAL
+      cx=centroid(HIDX(CELL),1)
+      cy=centroid(HIDX(CELL),2)
+      cz=centroid(HIDX(CELL),3)
+      ux=0. ; uy=0. ; uz=0.
+      FORALL_EDGES
+        {{ caller() }}
+	ux = ux + ue_le*(.5*(xyz_v(HIDX(VERTEX1),1)+xyz_v(HIDX(VERTEX2),1))-cx)
+	uy = uy + ue_le*(.5*(xyz_v(HIDX(VERTEX1),2)+xyz_v(HIDX(VERTEX2),2))-cy)
+	uz = uz + ue_le*(.5*(xyz_v(HIDX(VERTEX1),3)+xyz_v(HIDX(VERTEX2),3))-cz)
+      END_BLOCK
+      fac = scale*(1./AI)
+      ucenter(CELL,1)=ux*fac
+      ucenter(CELL,2)=uy*fac
+      ucenter(CELL,3)=uz*fac
+    END_BLOCK
+  END_BLOCK
+{% endmacro %}
+
+KERNEL(wind_centered)
+{% call() perot_reconstruction() %}
+ue_le = SIGN*ue(EDGE)*le(HIDX(EDGE))
+{% endcall %}
+END_BLOCK
+
+KERNEL(flux_centered)
+! NB : here the input data is a flux and has already the factor l_e in it
+! Input data is rescaled by factor scale
+{% call() perot_reconstruction() %}
+ue_le = SIGN*ue(EDGE)
+{% endcall %}
+END_BLOCK
+
+/*------------------------------------- Energetics --------------------------------*/
+
+{% macro energy_flux(energy) %}
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      {{ caller() }}
+      {{ energy }}(CELL) = {{ energy }}(CELL) + frac*rhodz(CELL)*energy
+      pk(CELL) = energy
+    END_BLOCK
+  END_BLOCK
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      {{ energy }}_flux(EDGE) = {{ energy }}_flux(EDGE) + .5*massflux(EDGE)*(pk(CELL1)+pk(CELL2))
+    END_BLOCK
+  END_BLOCK
+{% endmacro %}
+
+KERNEL(energy_fluxes)
+  ! First diagnose geopotential and temperature, column-wise
+  BARRIER
+  SEQUENCE_C1
+    PROLOGUE(llm)
+      pk(CELL) = ptop + .5*g*rhodz(CELL)
+    END_BLOCK
+    BODY('llm-1,1,-1')
+      pk(CELL) = pk(UP(CELL)) + (.5*g)*( rhodz(CELL)+rhodz(UP(CELL)) )
+    END_BLOCK
+  END_BLOCK
+  ! NB : at this point pressure is stored in array pk
+  !      pk then serves as buffer to store temperature
+  SELECT CASE(caldyn_thermo)
+  CASE(thermo_theta)
+    GEOPOT('temp_ik')
+      theta_ik = theta_rhodz(CELL,1)/rhodz(CELL)
+      theta(CELL) = theta_ik
+      temp_ik = theta_ik*(p_ik/preff)**kappa
+      gv = (g*Rd)*temp_ik/p_ik
+    END_GEOPOT
+  CASE(thermo_entropy)
+    GEOPOT('temp_ik')
+      theta_ik = theta_rhodz(CELL,1)/rhodz(CELL)
+      temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+      theta(CELL) = Treff*exp(theta_ik/cpp)
+      gv = (g*Rd)*temp_ik/p_ik   ! specific volume v = Rd*T/p
+    END_GEOPOT
+  END SELECT
+  BARRIER
+  ! Now accumulate energies and energy fluxes
+  ! NB : at this point temperature is stored in array pk
+  !      pk then serves as buffer to store energy
+  ! enthalpy
+  {% call energy_flux('enthalpy') %}
+    energy = cpp*pk(CELL)
+  {% endcall %}
+  ! potential energy
+  {% call energy_flux('epot') %}
+    energy = .5*(geopot(UP(CELL))+geopot(CELL))
+  {% endcall %}
+  ! theta
+  {% call energy_flux('thetat') %}
+    energy = theta(CELL)
+  {% endcall %}
+  ! kinetic energy
+  {% call energy_flux('ekin') %}
+    energy=0.d0
+    FORALL_EDGES
+      energy = energy + le(HIDX(EDGE))*de(HIDX(EDGE))*ue(EDGE)**2
+    END_BLOCK
+    energy = energy * (.25/AI)
+  {% endcall %}
+  ! ulon
+  {% call energy_flux('ulon') %}
+    cx=centroid(HIDX(CELL),1)
+    cy=centroid(HIDX(CELL),2)
+    cz=centroid(HIDX(CELL),3)
+    ux=0. ; uy=0. ; uz=0.
+    FORALL_EDGES
+      ue_le = SIGN*ue(EDGE)*le(HIDX(EDGE))
+      ux = ux + ue_le*(.5*(xyz_v(HIDX(VERTEX1),1)+xyz_v(HIDX(VERTEX2),1))-cx)
+      uy = uy + ue_le*(.5*(xyz_v(HIDX(VERTEX1),2)+xyz_v(HIDX(VERTEX2),2))-cy)
+      uz = uz + ue_le*(.5*(xyz_v(HIDX(VERTEX1),3)+xyz_v(HIDX(VERTEX2),3))-cz)
+    END_BLOCK
+    ulon_i = ux*elon_i(HIDX(CELL),1) + uy*elon_i(HIDX(CELL),2) + uz*elon_i(HIDX(CELL),3)
+    energy = ulon_i*(1./AI)
+  {% endcall %}
+
+END_BLOCK
+
+KERNEL(remap_eta)
+! IN  : rhodz, mass_dak, mass_dbk
+! TMP : mass_col, cur_lev, new_rhodz_cum
+! OUT : rhodz, old_rhodz, eta
+    SEQUENCE_C1
+      PROLOGUE('1')
+        rhodz_cum(CELL)=0.
+        cur_lev(HIDX(CELL))=1
+        eta(CELL)=1.
+        new_rhodz_cum(HIDX(CELL))=0.
+      END_BLOCK
+      BODY('1,llm')
+        rhodz_cum(UP(CELL)) = rhodz_cum(CELL) + rhodz(CELL)
+      END_BLOCK
+      EPILOGUE('llm+1')
+        mass_col(HIDX(CELL)) = rhodz_cum(CELL)
+      END_BLOCK
+
+      BODY('1,llm')
+        old_rhodz(CELL) = rhodz(CELL)
+        rhodz(CELL) = mass_dak(CELL) + mass_col(HIDX(CELL))*mass_dbk(CELL)
+        rhodz_cum_target = new_rhodz_cum(HIDX(CELL)) + rhodz(CELL)
+        DO level = cur_lev(HIDX(CELL)), llm
+           rhodz_cum_levp1 = rhodz_cum(AT_LEVEL(CELL,level+1))
+           IF(rhodz_cum_target<=rhodz_cum_levp1) EXIT
+        END DO
+        IF(level>llm) level=llm
+        rhodz_cum_lev = rhodz_cum(AT_LEVEL(CELL,level))
+        ! now rhodz_cum_lev <= rhodz_cum_target <= rhodz_cum_levp1  
+        cur_lev(HIDX(CELL)) = level
+        new_rhodz_cum(HIDX(CELL)) = rhodz_cum_target
+        eta(UP(CELL)) = level + (rhodz_cum_target-rhodz_cum_lev)/(rhodz_cum_levp1-rhodz_cum_lev)
+      END_BLOCK
+    END_BLOCK
+END_BLOCK
+
+KERNEL(remap_theta)
+! IN  : thetarhodz, eta
+! TMP : thetarhodz_cum, new_thetarhodz_cum
+! OUT : thetarhodz
+    SEQUENCE_C1
+      PROLOGUE('1')
+        thetarhodz_cum(CELL)=0.
+      END_BLOCK
+      BODY('1,llm')
+        thetarhodz_cum(UP(CELL)) = thetarhodz_cum(CELL) + thetarhodz(CELL)
+      END_BLOCK
+        
+      BODY('1,llm+1')
+        X = eta(CELL)
+        level = MIN(llm,FLOOR(X)) ! eta=llm+1 => level=llm, X=1
+        X = X-level
+        new_thetarhodz_cum(CELL) = thetarhodz_cum(AT_LEVEL(CELL,level))+X*thetarhodz(AT_LEVEL(CELL,level))
+      END_BLOCK
+      BODY('1,llm')
+        thetarhodz(CELL) = new_thetarhodz_cum(UP(CELL)) - new_thetarhodz_cum(CELL)
+      END_BLOCK
+    END_BLOCK
+END_BLOCK
+
+KERNEL(remap_u)
+! IN  : u, old_rhodz, rhodz, eta
+! TMP : urhodz_cum, new_urhodz_cum
+! OUT : u
+    SEQUENCE_E0
+      PROLOGUE('1')
+        urhodz_cum(EDGE)=0.
+      END_BLOCK
+      BODY('1,llm')
+        urhodz(EDGE) = u(EDGE)*(old_rhodz(CELL1)+old_rhodz(CELL2))
+        urhodz_cum(UP(EDGE)) = urhodz_cum(EDGE) + urhodz(EDGE)
+      END_BLOCK
+        
+      BODY('1,llm+1')
+        X = .5*(eta(CELL1)+eta(CELL2))
+        level = MIN(llm,FLOOR(X))
+        X = X-level
+        new_urhodz_cum(EDGE) = urhodz_cum(AT_LEVEL(EDGE,level))+X*urhodz(AT_LEVEL(EDGE,level))
+      END_BLOCK
+      BODY('1,llm')
+        u(EDGE) = (new_urhodz_cum(UP(EDGE)) - new_urhodz_cum(EDGE)) / (rhodz(CELL1)+rhodz(CELL2))
+      END_BLOCK
+    END_BLOCK
+END_BLOCK
+
+KERNEL(scalar_laplacian)
+  FORALL_CELLS_EXT()
+    ON_EDGES
+       grad(EDGE) = SIGN*(b(CELL2)-b(CELL1))
+    END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      div_ij=0.
+      FORALL_EDGES
+        div_ij = div_ij + SIGN*LE_DE*grad(EDGE)
+      END_BLOCK
+      divu(CELL) = div_ij / AI
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
+
+KERNEL(curl_laplacian)
+  FORALL_CELLS_EXT()
+    ON_DUAL
+      etav = 0.d0
+      FORALL_EDGES
+         etav = etav + SIGN*u(EDGE)
+      END_BLOCK
+      qv(DUAL_CELL) = etav/AV
+    END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT()
+    ON_EDGES
+       curlcurl(EDGE) = SIGN*(qv(VERTEX1)-qv(VERTEX2))*(1./LE_DE)
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
Index: codes/icosagcm/devel/DySL/macros.jin
===================================================================
--- codes/icosagcm/devel/DySL/macros.jin	(revision 941)
+++ codes/icosagcm/devel/DySL/macros.jin	(revision 941)
@@ -0,0 +1,50 @@
+{#  ---------------------------------- Macros common to all variants ---------------------------------------- #}
+
+{% set flat='#' %}
+
+{% macro cst_if(condition, action) %}
+{{ flat }}if {{ condition }}
+{{ action }}
+{{ flat }}endif
+{%- endmacro %}
+ 
+{% macro define(key,val) %} 
+{{ flat }}undef {{ key }}
+{{ flat }}define {{ key }} {{ val }} 
+{%- endmacro %}
+
+{% macro define_cppmacro(key,args,val) %} 
+{{ flat }}undef {{ key }} 
+{{ flat }}define {{ key+args }} {{ val }} 
+{%- endmacro %}
+
+{% macro undef(key) %} 
+{{ flat }}undef {{ key }}
+{{ flat }}define {{ key }} {{ key }}_UNDEFINED_HERE
+{%- endmacro %}
+
+{% macro cdef(code, key,val) %} 
+{{ define(key,val) if key in code }}
+{%- endmacro %}
+
+{% macro cundef(code, keys) %} 
+{% for key in keys %}
+{{ undef(key) if key in code }}
+{%- endfor %}
+{%- endmacro %}
+
+
+{% macro dummy() -%}
+  {{ caller() }}
+{%- endmacro %}
+
+{% macro define_kernel(name) %} 
+{{ flat }}ifdef KERNEL_{{name}}
+!--------------------------------------------------------------------------
+!---------------------------- {{ name }} ----------------------------------
+  {{ caller() }}
+!---------------------------- {{ name }} ----------------------------------
+!--------------------------------------------------------------------------
+
+{{ flat }}endif
+{%- endmacro %}
Index: codes/icosagcm/devel/DySL/unstructured/macros.jin
===================================================================
--- codes/icosagcm/devel/DySL/unstructured/macros.jin	(revision 941)
+++ codes/icosagcm/devel/DySL/unstructured/macros.jin	(revision 941)
@@ -0,0 +1,410 @@
+#ifdef PASS_PRE1
+#define _AND_ &&
+#define _OR_ ||
+#define _NOT_ !
+#define IS_INNER_LAYER (_NOT_ (IS_BOTTOM_LEVEL _OR_ IS_TOP_LAYER))
+#define IS_INNER_INTERFACE (_NOT_ (IS_BOTTOM_LEVEL _OR_ IS_TOP_INTERFACE))
+#endif
+
+#ifdef PASS_PRE2
+
+#define BARRIER !$OMP BARRIER
+#define IS_INNER_INTERFACE (.NOT. (IS_BOTTOM_LEVEL .OR. IS_TOP_INTERFACE))
+#define KERNEL(name) {% call() define_kernel(#name) %}
+#define FORALL_CELLS(...) {% call(lrange,flags) forall_cells(__VA_ARGS__) %}
+#define FORALL_CELLS_EXT(...) {% call(lrange,flags) forall_cells(__VA_ARGS__) %}
+#define ON_PRIMAL {% call(mesh, degree=0) on_mesh('primal',lrange,flags) %}
+#define ON_DUAL {% call(mesh, degree=0) on_mesh('dual',lrange,flags) %}
+#define ON_EDGES {% call(mesh, degree=0) on_edges('edge',lrange,flags) %}
+#define FORALL_EDGES {% call forall_edges(mesh,degree) %}
+#define FORALL_VERTICES {% call forall_vertices(mesh,degree) %}
+#define FORALL_TRISK {% call forall_trisk(degree) %}
+#define END_BLOCK {% endcall %}
+#define SEQUENCE_C0 {% call(body,at_level) sequence_primal() %} 
+#define SEQUENCE_C1 {% call(body,at_level) sequence_primal() %} 
+#define SEQUENCE_E0 {% call(body,at_level) sequence_edge() %} 
+#define PROLOGUE(level) {% call at_level(level) %}
+#define EPILOGUE(level) {% call at_level(level) %}
+#define BODY(range) {% call body(range) %}
+
+#define CST_IF(condition, action) {{ cst_if(#condition, #action) }}
+#define CST_IFTHEN(condition) {{flat}}if condition
+#define CST_ELSEIF(condition) {{flat}}elif condition
+#define CST_ELSE {{flat}}else
+#define CST_ENDIF {{flat}}endif
+
+#endif
+
+#ifdef PASS_JINJA
+
+{#  -------------------------- Design notes ------------------------------
+
+The syntax 
+{% call(ARGNAMES) macro(ARGS) %}
+...
+{% endcall %}
+passes ARGS to macro ; this macro "calls" the body using the syntax 'caller(arg1,...)' ; these arguments
+are available in the body under the names ARGNAMES, as if it were a macro caller(ARGNAMES).
+
+Using its optional argument 'options', forall_cells() sets variable 'flags' 
+and passes it on to on_mesh() / on_edges()
+
+Similarly variable 'lrange' is passed from forall_cells to on_mesh() / on_egdes()  
+
+-------------------------------------------------------------------------- #}
+
+{# NB : degree = number of elements in stencil, for instance number of edges of primal cell #}
+{# define the degrees for which a special version of code is generated #}
+{% set degree_primal, degree_dual, degree_trisk = (4,5,6), (3,4), (4,10) %}
+
+{# ----------------------- SEQUENCE, FORALL ----------------------- #}
+
+{% macro sequence_primal() -%}
+!$OMP DO SCHEDULE(STATIC)                                                                            
+  DO ij=1,primal_num
+    {{ caller(body_primal, at_level_primal) }}
+  END DO
+!$OMP END DO
+{%- endmacro %}
+
+{% macro sequence_edge() -%}
+{% set thecode = caller(body_edge, at_level_edge) %}
+!$OMP DO SCHEDULE(STATIC)                                                                            
+  DO edge=1,edge_num
+{{ 'ij_left = left(edge)'   if 'CELL1'    in thecode }}
+{{ 'ij_right = right(edge)' if 'CELL2'   in thecode }}
+{{ 'ij_up = up(edge)'       if 'VERTEX1' in thecode }}
+{{ 'ij_down = down(edge)'   if 'VERTEX2' in thecode }}
+    {{ thecode }}
+  END DO
+!$OMP END DO
+{%- endmacro %}
+
+{% macro body_primal(range) -%}
+  {{ define('CELL','l,ij') }}
+    DO l = {{ range }}
+      {{ caller() }}
+    END DO
+  {{ undef('CELL') }}
+{%- endmacro %}
+
+{% macro body_edge(range) -%}
+{% set thecode = caller() %}
+{{ define('EDGE', 'l,edge') }}
+{{ cdef(thecode, 'LE_DE', 'le_de(edge)') }}
+{{ cdef(thecode, 'LE', 'le(edge)') }}
+{{ cdef(thecode, 'DE', 'de(edge)') }}
+{{ cdef(thecode, 'SIGN', '1.') }}
+{{ cdef(thecode, 'CELL1', 'l,ij_left') }}
+{{ cdef(thecode, 'CELL2',  'l,ij_right') }}
+{{ cdef(thecode, 'VERTEX1', 'l,ij_down') }}
+{{ cdef(thecode, 'VERTEX2', 'l,ij_up') }}
+    DO l = {{ range }}
+      {{ thecode }}
+    END DO
+{{ cundef(thecode, ('EDGE','LE_DE','SIGN','CELL1','CELL2','VERTEX1','VERTEX2') ) }}
+{%- endmacro %}
+
+{% macro at_level_primal(lev) -%}
+  {{ define('CELL','%s,ij'%lev) }}
+  {{ caller() }}
+  {{ undef('CELL') }}
+{%- endmacro %}
+
+{% macro at_level_edge(lev) -%}
+  {{ define('EDGE','%s,edge'%lev) }}
+  {{ caller() }}
+  {{ undef('EDGE') }}
+{%- endmacro %}
+
+{% macro forall_cells(start='1', end='llm', options=[]) -%}
+  {% set flags=[options] if options is string else options %}
+  {{ caller((start,end),flags) }}
+{%- endmacro %}
+
+{# ---------------------- PRELOAD STENCIL INDICES ---------------------- #}
+
+{% macro get_edges(thecode, mesh,degree) %}
+{% for iedge in range(1,degree+1) %} 
+{{ 'edge%d = %s_edge(%d,ij)'%(iedge,mesh,iedge) if 'EDGE' in thecode }} 
+{% endfor %}
+{% for iedge in range(1,degree+1) %} 
+{{ 'le_de%d = le_de(edge%d)'%(iedge,iedge) if 'LE_DE' in thecode }}      
+{% endfor %}
+{% for iedge in range(1,degree+1) %} 
+{{ 'sign%d = %s_ne(%d,ij)'%(iedge,mesh,iedge) if 'SIGN' in thecode }}   
+{% endfor %}
+{% for iedge in range(1,degree+1) %} 
+{{ 'ij_up%d = up(edge%d)'%(iedge,iedge) if 'VERTEX1' in thecode }}
+{% endfor %}
+{% for iedge in range(1,degree+1) %} 
+{{ 'ij_down%d = down(edge%d)'%(iedge,iedge) if 'VERTEX2' in thecode }}
+{% endfor %}
+{% endmacro %}
+
+{% macro get_vertices(thecode, mesh,degree) %}
+{% for ivertex in range(1,degree+1) %} 
+{{ 'vertex%d = %s_vertex(%d,ij)'%(ivertex,mesh,ivertex) if 'VERTEX' in thecode }} 
+{% endfor %}
+{% endmacro %}
+
+{# ------------------------------ MESHES --------------------------- #}
+
+{# argument 'code' is the body enclosed by ON_XXX ... END_BLOCK ; it takes arguments : mesh,degree=0 #}
+
+{% macro vloop_unroll(mesh, lrange, flags, code) %}
+
+{% set thecode, has_primal, has_dual, has_trisk, has_none = code(mesh), False, False, False, True %}
+
+{% if 'primal_deg' in thecode %}
+! this VLOOP iterates over primal cell edges
+{% set has_primal, has_none = True,False %}
+{% endif %}
+{% if 'dual_deg' in thecode %}
+! this VLOOP iterates over dual cell edges
+{% set has_dual,has_none = True,False %}
+{% endif %}
+{% if 'trisk_deg' in thecode %}
+! this VLOOP iterates over the TRISK stencil
+{% set has_trisk,has_none = True,False %}
+{% endif %}
+
+{% if has_none %}
+{{ vloop(mesh, lrange, flags, code) }}
+{% endif %}
+
+{% if has_primal %}
+SELECT CASE(primal_deg(ij))
+{% for degree in degree_primal %}
+CASE({{ degree }})
+{{ get_edges(thecode, mesh, degree) }}
+{{ vloop(mesh, lrange, flags, code, degree) }}
+{% endfor %}
+CASE DEFAULT
+{{ vloop(mesh, lrange, flags, code) }}
+END SELECT
+{% endif %}
+
+{% if has_dual %}
+SELECT CASE(dual_deg(ij))
+{% for degree in degree_dual %}
+CASE({{ degree }})
+{{ get_edges(thecode, mesh, degree) }}
+{{ get_vertices(thecode,mesh,degree) }}
+{{ vloop(mesh, lrange, flags, code, degree) }}
+{% endfor %}
+CASE DEFAULT
+{{ vloop(mesh, lrange, flags, code) }}
+END SELECT
+{% endif %}
+
+{% if has_trisk %}
+SELECT CASE(trisk_deg(edge))
+{% for degree in degree_trisk %}
+CASE({{ degree }})
+{{ vloop(mesh, lrange, flags, code, degree) }}
+{% endfor %}
+CASE DEFAULT
+{{ vloop(mesh, lrange, flags, code) }}
+END SELECT
+{% endif %}
+
+{%- endmacro %}
+
+{% macro vloop(mesh, lrange, flags, code, degree=0) %}
+
+{% set start,end = lrange %}
+{% set thecode, is_top_layer, is_top_inter = code(mesh,degree), False, False %}
+
+{{ define('IS_TOP_LAYER', '_FALSE_') }}
+{{ define('IS_TOP_INTERFACE', '_FALSE_') }}
+
+{% if 'IS_BOTTOM_LEVEL' in thecode or 'KDOWN' in thecode%}
+{# the code in the loop checks whether l==1, for the sake of performance 
+we shall write special code for l=1 and start the loop at l=2 #}
+{{ define('IS_BOTTOM_LEVEL', '_TRUE_') }}
+{{ 'kdown = 1' if 'KDOWN' in thecode }}
+{{ 'kup   = 1' if 'KUP'   in thecode }}
+l=1
+{{ thecode }}
+{% set start='2' %}
+{% endif %}
+
+{% if 'IS_TOP_LAYER' in thecode %}
+{{ 'ERROR : using IS_TOP_LAYER in a loop ending at l=llm+1' if end=='llm+1' }}
+{# the code checks whether l==llm, write special code for l=llm and end the loop at l=llm-1 #}
+{% set end, is_top_layer = 'llm-1', True %}
+{% endif %}
+
+{% if 'IS_TOP_INTERFACE' in thecode or 'KUP' in thecode %}
+{# the code checks whether l==llm+1, write special code for l=llm+1 and end the loop at l=llm #}
+{{ 'ERROR : using IS_TOP_INTERFACE in a loop ending at l=llm' if end=='llm' }}
+{% set end, is_top_inter ='llm', True %}
+{% endif %}
+
+{{ define('IS_BOTTOM_LEVEL', '_FALSE_') }}
+!DIR$ SIMD
+  DO l = {{start}}, {{end}}
+    {{ 'kdown = l-1' if 'KDOWN' in thecode }}
+    {{ 'kup   = l' if 'KUP'   in thecode }}
+    {{ thecode }}
+  END DO
+
+{% if is_top_layer %}
+{{ define('IS_TOP_LAYER', '_TRUE_') }}
+{{ 'kdown = llm-1' if 'KDOWN' in thecode }}
+{{ 'kup   = llm' if 'KUP'   in thecode }}
+l=llm
+{{ thecode }}
+{% endif %}
+
+{% if is_top_inter %}
+{{ define('IS_TOP_INTERFACE', '_TRUE_') }}
+{{ 'kdown = llm' if 'KDOWN' in thecode }}
+{{ 'kup   = llm' if 'KUP'   in thecode }}
+l=llm+1
+{{ thecode }}
+{% endif %}
+
+{% endmacro %}
+
+{% macro on_mesh(mesh,lrange,flags) -%}
+{{ define('CELL','l,ij') if mesh=='primal' }}
+{{ define('DUAL_CELL', 'l,ij') if mesh=='dual'}}
+!$OMP DO SCHEDULE(STATIC)                                                                            
+DO ij = 1, {{ mesh }}_num
+{{ vloop_unroll(mesh, lrange, flags, caller) }}
+END DO
+!$OMP END DO
+{{ undef('CELL') }}
+{{ undef('DUAL_CELL') }}
+{%- endmacro %}
+
+{# ------------------------------ STENCILS --------------------------- #}
+
+{% macro on_edges(mesh,lrange,flags) -%}
+{% set thecode = caller(mesh) %}
+{{ define('EDGE', 'l,edge') }}
+{{ cdef(thecode, 'LE_DE', 'le_de(edge)') }}
+{{ cdef(thecode, 'DE', 'de(edge)') }}
+{{ cdef(thecode, 'LE', 'le(edge)') }}
+{{ cdef(thecode, 'SIGN', '1.') }}
+{{ cdef(thecode, 'CELL1', 'l,ij_left') }}
+{{ cdef(thecode, 'CELL2',  'l,ij_right') }}
+{{ cdef(thecode, 'VERTEX1', 'l,ij_down') }}
+{{ cdef(thecode, 'VERTEX2', 'l,ij_up') }}
+!$OMP DO SCHEDULE(STATIC)
+DO edge = 1, edge_num
+{{ 'ij_left = left(edge)'   if 'CELL1'    in thecode }}
+{{ 'ij_right = right(edge)' if 'CELL2'   in thecode }}
+{{ 'ij_up = up(edge)'       if 'VERTEX1' in thecode }}
+{{ 'ij_down = down(edge)'   if 'VERTEX2' in thecode }}
+{{ vloop_unroll(mesh, lrange, flags, caller) }}
+END DO
+!$OMP END DO
+{{ cundef(thecode, ('EDGE','LE_DE','SIGN','CELL1','CELL2','VERTEX1','VERTEX2') ) }}
+{%- endmacro %}
+
+{% macro forall_edges(mesh,degree) -%}
+{% set thecode = caller() %}
+
+{% if degree>1 %}
+{% for iedge in range(1,degree+1) %}
+{{ cdef(thecode, 'EDGE', 'l,edge%d'%iedge) }}
+{{ cdef(thecode, 'SIGN', 'sign%d'%iedge) }}
+{{ cdef(thecode, 'LE_DE', 'le_de%d'%iedge) }}
+{{ cdef(thecode, 'VERTEX1', 'l,ij_up%d'%iedge) }}
+{{ cdef(thecode, 'VERTEX2', 'l,ij_down%d'%iedge) }}
+  {{ thecode }}
+{% endfor %}
+{% else %}
+{{ cdef(thecode, 'EDGE', 'l,edge') }}
+{{ cdef(thecode, 'SIGN', mesh + '_ne(iedge,ij)') }}
+{{ cdef(thecode, 'LE_DE', 'le_de(edge)') }}
+{{ cdef(thecode, 'VERTEX1', 'l,ij_up') }}
+{{ cdef(thecode, 'VERTEX2', 'l,ij_down') }}
+DO iedge = 1, {{ mesh }}_deg(ij)
+  edge = {{ mesh }}_edge(iedge,ij)
+  {{ 'ij_up = up(edge)'       if 'VERTEX1' in thecode }}
+  {{ 'ij_down = down(edge)'   if 'VERTEX2' in thecode }}
+  {{ thecode }}
+END DO
+{% endif %}
+{{ cundef(thecode, ('EDGE', 'SIGN', 'LE_DE', 'VERTEX1', 'VERTEX2') ) }}
+{%- endmacro %}
+
+{% macro forall_trisk(degree) -%}
+{% set thecode = caller() %}
+{{ define('EDGE_TRISK', 'l,edge_trisk') }}
+{% if degree>1 %}
+{% for itrisk in range(1,degree+1) %}
+  itrisk = {{ itrisk }}
+  edge_trisk = trisk({{ itrisk }},edge)
+  {{ thecode }}
+{% endfor %}
+{% else %}
+DO itrisk = 1, trisk_deg(edge)
+  edge_trisk = trisk(itrisk,edge)
+  {{ thecode }}
+END DO
+{% endif %}
+{{ undef('EDGE_TRISK') }}
+{%- endmacro %}
+
+{% macro forall_vertices(mesh,degree) -%}
+{% set thecode = caller() %}
+{% if degree>1 %}
+{% for ivertex in range(1,degree+1) %}
+{{ define('RIV2', 'Riv2(%d,ij)'%ivertex ) }}
+{{ define('VERTEX', 'l,vertex%d'%ivertex ) }}
+  {{ thecode }}
+{% endfor %}
+{% else %}
+{{ define('RIV2', 'Riv2(ivertex,ij)') }}
+{{ define('VERTEX', 'l,vertex') }}
+DO ivertex = 1, {{ mesh }}_deg(ij)
+  vertex = {{ mesh }}_vertex(ivertex,ij)
+  {{ thecode }}
+END DO
+{% endif %}
+{{ undef('VERTEX') }}
+{{ undef('RIV2') }}
+{%- endmacro %}
+
+{# --------------------------------------------------------- #}
+
+{% set llm='llm' %}
+
+{{ undef('SIGN') }}
+{{ undef('CELL') }}
+{{ undef('CELL1') }}
+{{ undef('CELL2') }}
+{{ undef('VERTEX1') }}
+{{ undef('VERTEX2') }}
+{{ undef('EDGE_TRISK') }}
+
+{# --------------------- END JINJA ------------------------- #}
+
+#endif
+
+#ifdef PASS_POST1
+
+#define _TRUE_ (0==0)
+#define _FALSE_ (0==1)
+
+#define AI Ai(ij)
+#define AV Av(ij)
+#define FV fv(ij)
+#define WEE wee(itrisk,edge,1)
+#define edge_ne(iedge,ij) 1.
+#endif
+
+#ifdef PASS_POST2
+#define AT_LEVEL(l,ij,lev) lev,ij
+#define KUP(l,ij) kup,ij
+#define KDOWN(l,ij) kdown,ij
+#define DOWN(l,ij) l-1,ij
+#define UP(l,ij) l+1,ij
+#define HIDX(l,ij) ij
+#define VIDX(l,ij) l
+#endif
Index: codes/icosagcm/devel/arch/arch-X64_IRENE.dev
===================================================================
--- codes/icosagcm/devel/arch/arch-X64_IRENE.dev	(revision 941)
+++ codes/icosagcm/devel/arch/arch-X64_IRENE.dev	(revision 941)
@@ -0,0 +1,7 @@
+module load python/2.7.14
+. $CCCWORKDIR/IRENE/virtualenv/GrandChallenge/bin/activate
+
+. linux.dev
+
+# Needed for Intel OpenMP
+KERNELS_LDFLAGS="$NETCDF_LIBDIR $NETCDF_LIB -liomp5 $MKL_LDFLAGS $PARMETIS_LDFLAGS"
Index: codes/icosagcm/devel/arch/arch-X64_OSX.dev
===================================================================
--- codes/icosagcm/devel/arch/arch-X64_OSX.dev	(revision 941)
+++ codes/icosagcm/devel/arch/arch-X64_OSX.dev	(revision 941)
@@ -0,0 +1,43 @@
+export CODEGEN_ENV_SET=yes
+
+#----------- Useful directories ---------
+
+BASEDIR=$(dirname ${BASH_SOURCE[0]})
+export BASEDIR=$(cd $BASEDIR/../.. ; pwd)
+
+#--------- To create shared libraries -------
+
+export F90=mpif90
+export LINK=mpicxx
+export SOFLAGS="-dynamiclib -undefined dynamic_lookup -Wl,-headerpad_max_install_names -lgfortran -lparmetis $NETCDF_LIBDIR -fopenmp -llapack -lcblas -lf77blas -latlas"
+export WHOLE_ARCHIVE='-Wl,-all_load'
+export NOWHOLE_ARCHIVE=''
+export LDD="otool -L"
+
+#----------- To compile kernels ------------
+
+export C99=/usr/bin/clang
+export CFLAGS="-fPIC $PARMETIS_CFLAGS"
+export C99FLAGS="-std=c99 -I /opt/local/include/openmpi-gcc48 -I /opt/local/include $CLFAGS"
+export FFLAGS_BASE="-I $XIOS_ROOT/inc -fPIC -g -openmp"
+export FFLAGS_PROD="-O3"
+export FFLAGS_DBG="-O3 -fcheck=bounds"
+export FFLAGS="$FFLAGS_BASE $FFLAGS_PROD"
+export LDFLAGS="-fopenmp -lxios"
+
+#---------------- To execute ----------------
+
+# Python will look for modules (*.py) and extensions (*.so) in PYTHONPATH
+# Shared libs they depend on should be in a standard place or a path included in LD_LIBRARY_PATH
+export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:$DYNAMICO_ROOT/lib"
+export PYTHONPATH="$DYNAMICO_ROOT/lib"
+export PYTHON=$(which python)
+
+#----------------- Garbage kept for the record -------------
+# export PYTHON_ROOT=/opt/local/Library/Frameworks/Python.framework/Versions/2.7/
+
+# export LIBS="-L /opt/local/lib/openmpi-gcc48/ -L /opt/local/lib -L/opt/local/lib/gcc48 /opt/local/lib/gcc48/libgfortran.a -lgcc_ext.10.5 -lmpi -lparmetis"
+# export MPIRUN=mpirun-openmpi-gcc48
+# export NETCDF_LDFLAGS="-L/opt/local/lib -lnetcdff -lnetcdf -lstdc++"
+#export CC=/usr/bin/clang
+#export CFLAGS="-fno-strict-aliasing -fno-common -dynamic -pipe -Os -fwrapv -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/include -I/opt/local/Library/Frameworks/Python.framework/Versions/2.7/include/python2.7"
Index: codes/icosagcm/devel/arch/arch-ifort_CICLAD.dev
===================================================================
--- codes/icosagcm/devel/arch/arch-ifort_CICLAD.dev	(revision 941)
+++ codes/icosagcm/devel/arch/arch-ifort_CICLAD.dev	(revision 941)
@@ -0,0 +1,40 @@
+# arch.env and arch.path should be sourced before sourcing this script
+module load python/2.7-anaconda
+
+# If not done yet, certain python packages must be installed in a dynamico-only conda environment of yours :
+#   conda create -n dynamico python
+#   source activate dynamico
+#   pip install scipy mpi4py jinja2 matplotlib cython cmake
+
+# netCDF4 must be installed manually in order to use the netCDF library from the modules (requires Cython, see above)
+# for unpack Python/extern/netCDF4-1.2.9.tar.gz then
+#   LDFLAGS='-lirc' python setup.py install
+# NOTE : libirc.so is the intel runtime library and is needed at runtime.
+# NOTE : setup.py installs a *.egg file which contains the netCDF4 and netcdftime module. Deactivate then activate again your conda environment if you want to test that they are working.
+
+# CICLAD does not have a parmetis module
+# compiling parmetis uses cmake, installed above using pip because the system-wide cmake might not work with the modules
+# for this unpack Python/extern/parmetis-4.0.3.tar.gz then
+#   make config shared=1 prefix=$HOME/.conda/envs/dynamico
+#   make install
+# manually install libmetis.so and metis.h ...
+#   cp $(find . -name 'libmetis.so') $HOME/.conda/envs/dynamico/lib/
+#   cp $(find . -name 'metis.h') $HOME/.conda/envs/dynamico/include/
+# NOTE : the tar.gz is derived from http://glaros.dtc.umn.edu/gkhome/fetch/sw/parmetis/parmetis-4.0.3.tar.gz 
+#        with modifications along the lines of
+#        https://stackoverflow.com/questions/35058709/building-shared-parmetis-4-0-3
+
+source activate dynamico
+PARMETIS_CFLAGS="-I$HOME/.conda/envs/dynamico/include"
+PARMETIS_LDFLAGS="-L$HOME/.conda/envs/dynamico/lib -lparmetis"
+
+export NETCDF_LDFLAGS="$(nf-config --flibs)"
+
+KERNELS_LDFLAGS="$NETCDF_LIBDIR $NETCDF_LIB -liomp5 $MKL_LIBDIR $MKL_LIB $PARMETIS_LDFLAGS"
+# check that packages work
+for x in netCDF4 mpi4py numpy jinja2 ; do python -c "import $x" ; done
+python -c "import ctypes ; parmetis=ctypes.cdll.LoadLibrary('libparmetis.so')"
+
+export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/.conda/envs/dynamico/lib
+
+. linux.dev
Index: codes/icosagcm/devel/arch/linux.dev
===================================================================
--- codes/icosagcm/devel/arch/linux.dev	(revision 941)
+++ codes/icosagcm/devel/arch/linux.dev	(revision 941)
@@ -0,0 +1,39 @@
+export CODEGEN_ENV_SET=yes
+
+#----------- Useful directories ---------
+
+BASEDIR=$(dirname ${BASH_SOURCE[0]})
+export BASEDIR=$(cd $BASEDIR/../.. ; pwd)
+
+#--------- To create shared libraries -------
+
+export F90=mpif90
+export LINK=$F90
+export SOFLAGS="-shared"
+export WHOLE_ARCHIVE='-Wl,--whole-archive'
+export NO_WHOLE_ARCHIVE='-Wl,--no-whole-archive'
+export LDD="ldd"
+
+#----------- To compile kernels ------------
+
+export C99=mpicc
+export CFLAGS="-fPIC $PARMETIS_CFLAGS"
+export C99FLAGS="-std=c99 $CFLAGS"
+export FFLAGS_BASE="-I $XIOS_INCDIR -fPIC -g -openmp"
+export FFLAGS_PROD="-O3"
+export FFLAGS_DBG="-O3 -fcheck=bounds"
+export FFLAGS="$FFLAGS_BASE $FFLAGS_PROD"
+
+#---------------- To execute ----------------
+
+# Python will look for modules (*.py) and extensions (*.so) in PYTHONPATH
+# Shared libs they depend on should be in a standard place or a path included in LD_LIBRARY_PATH
+export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:$DYNAMICO_ROOT/lib"
+export PYTHONPATH="$DYNAMICO_ROOT/lib"
+export PYTHON=$(which python)
+
+# disable Intel OpenMP in MKL (loaded by numpy) to avoid conflict with our openMP
+export MKL_THREADING_LAYER=SEQUENTIAL
+
+# select a Matplotlib backend that works without an X11 connection
+export MPLBACKEND="agg"
Index: codes/icosagcm/devel/make_DySL
===================================================================
--- codes/icosagcm/devel/make_DySL	(revision 941)
+++ codes/icosagcm/devel/make_DySL	(revision 941)
@@ -0,0 +1,65 @@
+#!/bin/bash
+
+function log()
+{
+    echo $*
+    $*
+}
+
+#----------- Extract DySL from Fortran, invoked by : './make_DySL extract' ------------
+
+function cmd_extract()
+{
+    XT=$KERNELS/kernels_extracted.jin
+    rm -f $XT
+    touch $XT
+    cd $DYNAMICO_ROOT/src
+    for F90 in */*.F90 ; do
+	sed -n '/BEGIN_DYSL/,/END_DYSL/{//b;p}' $F90 >> $XT
+    done
+    cd $DYNAMICO_ROOT
+}
+
+#--------------- Generate kernels, invoked by : './make_DySL' ---------------
+
+function update()
+{
+    diff $2 $1 || log cp -pr $1 $2
+}
+
+function update_kernels()
+{
+    cd $KERNELS/$1
+    for x in *.k90 ; do
+	update "$x" "$DYNAMICO_ROOT/src/$2/$x"
+    done
+}
+
+function cmd_() 
+{
+    cmd_clean
+    cmd_extract
+    cd $KERNELS
+    ./codegen hexagonal hex_master unstructured
+    update_kernels hexagonal kernels_hex
+    update_kernels hex_master kernels_hex_master
+    update_kernels unstructured kernels_unst
+# make sure that ./rebuild recompiles dynamics
+    cd $DYNAMICO_ROOT
+#    touch src/dynamics/*.F90 src/diagnostics/*.F90 src/transport/*.F90 src/unstructured/*.F90
+}
+
+#------------------------- script starts here ------------------------
+
+DYNAMICO_ROOT=$PWD
+KERNELS=$DYNAMICO_ROOT/DySL
+
+. arch.env
+. arch.path
+echo "Selected architecture : $CURRENT_ARCH"
+if [[ "$CURRENT_ARCH" == "" ]] ; then echo "Current architecture undefined or unsupported" ; exit ; fi
+
+cd arch
+. "arch-${CURRENT_ARCH}.dev"
+cd ..
+time cmd_$1
Index: codes/icosagcm/devel/make_python
===================================================================
--- codes/icosagcm/devel/make_python	(revision 940)
+++ codes/icosagcm/devel/make_python	(revision 941)
@@ -7,55 +7,9 @@
 }
 
-function full_dir()
-{
-    (cd -P $1 ; pwd)
-}
-
 function cmd_clean() # cleanup, invoked by : './make_python clean'
 {
     rm -rf lib/*.so lib/dynamico
-    cd $KERNELS/src
+    cd $DYNAMICO_ROOT/Python/src
     rm -rf *.o *.so build/
-}
-
-#--------------- Generate kernels, invoked by : './make_python kernels' ---------------
-
-function update()
-{
-    diff $2 $1 || log cp -pr $1 $2
-}
-
-function cmd_extract_dysl()
-{
-    XT=$DYNAMICO_ROOT/Python/src/kernels_extracted.jin
-    rm -f $XT
-    touch $XT
-    cd $DYNAMICO_ROOT/src
-    for F90 in */*.F90 ; do
-	sed -n '/BEGIN_DYSL/,/END_DYSL/{//b;p}' $F90 >> $XT
-    done
-    cd $DYNAMICO_ROOT
-}
-
-function update_kernels()
-{
-    cd $KERNELS/src/$1
-    for x in *.k90 ; do
-	update "$x" "$DYNAMICO_ROOT/src/$2/$x"
-    done
-}
-
-function cmd_kernels() # this function is invoked by : './make_python kernels'
-{
-    cmd_clean
-    cmd_extract_dysl
-    cd $KERNELS
-    ./codegen hexagonal hex_master unstructured
-    update_kernels hexagonal kernels_hex
-    update_kernels hex_master kernels_hex_master
-    update_kernels unstructured kernels_unst
-# make sure that ./rebuild recompiles dynamics
-    cd $DYNAMICO_ROOT
-#    touch src/dynamics/*.F90 src/diagnostics/*.F90 src/transport/*.F90 src/unstructured/*.F90
 }
 
@@ -120,5 +74,4 @@
 DYNAMICO_ROOT=$PWD
 LIBDIR=$DYNAMICO_ROOT/lib
-KERNELS=$DYNAMICO_ROOT/Python
 
 . arch.env
@@ -126,8 +79,9 @@
 
 echo "Selected architecture : $CURRENT_ARCH"
-
 if [[ "$CURRENT_ARCH" == "" ]] ; then echo "Current architecture undefined or unsupported" ; exit ; fi
 
-cd Python/env ;  . ${CURRENT_ARCH}.env ; cd $DYNAMICO_ROOT
+cd arch
+. arch-${CURRENT_ARCH}.dev
+cd $DYNAMICO_ROOT
 
 time cmd_$1
