Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 939)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 940)
@@ -1,4 +1,5 @@
 MODULE compute_caldyn_Coriolis_mod
   USE prec, ONLY : rstd
+  USE caldyn_vars_mod, ONLY : caldyn_conserv, conserv_energy, conserv_enstrophy, conserv_gassmann
   USE grid_param
   USE earth_const
@@ -46,15 +47,27 @@
     END_BLOCK
   END_BLOCK
-!
-  FORALL_CELLS()
-    ON_EDGES
-      du_trisk=0.
-      FORALL_TRISK
-        du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK))
-      END_BLOCK
-      du(EDGE) = du(EDGE) + .5*du_trisk
-    END_BLOCK
-  END_BLOCK
-
+  !
+  SELECT CASE(caldyn_conserv)
+  CASE(conserv_energy) ! energy-conserving TRiSK
+     FORALL_CELLS()
+       ON_EDGES
+         du_trisk=0.
+         FORALL_TRISK
+           du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK))
+         END_BLOCK
+         du(EDGE) = du(EDGE) + .5*du_trisk
+       END_BLOCK
+     END_BLOCK
+  CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK
+     FORALL_CELLS()
+       ON_EDGES
+         du_trisk=0.
+         FORALL_TRISK
+           du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)
+         END_BLOCK
+         du(EDGE) = du(EDGE) + du_trisk*qu(EDGE)
+       END_BLOCK
+     END_BLOCK
+  END SELECT
 END_BLOCK
 
@@ -65,5 +78,5 @@
     USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
          id_coriolis, left, right, primal_deg, primal_edge, primal_ne, &
-         trisk_deg, trisk, wee ! FIXME wee
+         trisk_deg, trisk
          
     FIELD_U     :: hflux, Ftheta, qu, du ! IN, BUF, IN, INOUT
@@ -80,5 +93,4 @@
   SUBROUTINE compute_caldyn_Coriolis_hex(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
     USE icosa
-    USE caldyn_vars_mod
     REAL(rstd),INTENT(IN)    :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
     REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars
Index: codes/icosagcm/devel/src/kernels_hex/coriolis.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_hex/coriolis.k90	(revision 939)
+++ codes/icosagcm/devel/src/kernels_hex/coriolis.k90	(revision 940)
@@ -40,45 +40,90 @@
    END DO
    !
-   DO l = ll_begin, ll_end
-      !DIR$ SIMD
-      DO ij=ij_begin, ij_end
-         du_trisk=0.
-         du_trisk = du_trisk + wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))
-         du_trisk = du_trisk + wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))
-         du_trisk = du_trisk + wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))
-         du_trisk = du_trisk + wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))
-         du_trisk = du_trisk + wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))
-         du_trisk = du_trisk + wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))
-         du_trisk = du_trisk + wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))
-         du_trisk = du_trisk + wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))
-         du_trisk = du_trisk + wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))
-         du_trisk = du_trisk + wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
-         du(ij+u_right,l) = du(ij+u_right,l) + .5*du_trisk
-         du_trisk=0.
-         du_trisk = du_trisk + wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l))
-         du_trisk = du_trisk + wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l))
-         du_trisk = du_trisk + wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l))
-         du_trisk = du_trisk + wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l))
-         du_trisk = du_trisk + wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l))
-         du_trisk = du_trisk + wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l))
-         du_trisk = du_trisk + wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l))
-         du_trisk = du_trisk + wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l))
-         du_trisk = du_trisk + wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l))
-         du_trisk = du_trisk + wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
-         du(ij+u_lup,l) = du(ij+u_lup,l) + .5*du_trisk
-         du_trisk=0.
-         du_trisk = du_trisk + wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l))
-         du_trisk = du_trisk + wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l))
-         du_trisk = du_trisk + wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l))
-         du_trisk = du_trisk + wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l))
-         du_trisk = du_trisk + wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l))
-         du_trisk = du_trisk + wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l))
-         du_trisk = du_trisk + wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l))
-         du_trisk = du_trisk + wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l))
-         du_trisk = du_trisk + wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l))
-         du_trisk = du_trisk + wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
-         du(ij+u_ldown,l) = du(ij+u_ldown,l) + .5*du_trisk
+   SELECT CASE(caldyn_conserv)
+   CASE(conserv_energy) ! energy-conserving TRiSK
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            du_trisk=0.
+            du_trisk = du_trisk + wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))
+            du_trisk = du_trisk + wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))
+            du_trisk = du_trisk + wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))
+            du_trisk = du_trisk + wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))
+            du_trisk = du_trisk + wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))
+            du_trisk = du_trisk + wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))
+            du_trisk = du_trisk + wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))
+            du_trisk = du_trisk + wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))
+            du_trisk = du_trisk + wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))
+            du_trisk = du_trisk + wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
+            du(ij+u_right,l) = du(ij+u_right,l) + .5*du_trisk
+            du_trisk=0.
+            du_trisk = du_trisk + wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l))
+            du_trisk = du_trisk + wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l))
+            du_trisk = du_trisk + wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l))
+            du_trisk = du_trisk + wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l))
+            du_trisk = du_trisk + wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l))
+            du_trisk = du_trisk + wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l))
+            du_trisk = du_trisk + wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l))
+            du_trisk = du_trisk + wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l))
+            du_trisk = du_trisk + wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l))
+            du_trisk = du_trisk + wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
+            du(ij+u_lup,l) = du(ij+u_lup,l) + .5*du_trisk
+            du_trisk=0.
+            du_trisk = du_trisk + wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l))
+            du_trisk = du_trisk + wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l))
+            du_trisk = du_trisk + wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l))
+            du_trisk = du_trisk + wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l))
+            du_trisk = du_trisk + wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l))
+            du_trisk = du_trisk + wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l))
+            du_trisk = du_trisk + wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l))
+            du_trisk = du_trisk + wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l))
+            du_trisk = du_trisk + wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l))
+            du_trisk = du_trisk + wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
+            du(ij+u_ldown,l) = du(ij+u_ldown,l) + .5*du_trisk
+         END DO
       END DO
-   END DO
+   CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            du_trisk=0.
+            du_trisk = du_trisk + wee(ij+u_right,1,1)*hflux(ij+u_rup,l)
+            du_trisk = du_trisk + wee(ij+u_right,2,1)*hflux(ij+u_lup,l)
+            du_trisk = du_trisk + wee(ij+u_right,3,1)*hflux(ij+u_left,l)
+            du_trisk = du_trisk + wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)
+            du_trisk = du_trisk + wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)
+            du_trisk = du_trisk + wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)
+            du_trisk = du_trisk + wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)
+            du_trisk = du_trisk + wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)
+            du_trisk = du_trisk + wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)
+            du_trisk = du_trisk + wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
+            du(ij+u_right,l) = du(ij+u_right,l) + du_trisk*qu(ij+u_right,l)
+            du_trisk=0.
+            du_trisk = du_trisk + wee(ij+u_lup,1,1)*hflux(ij+u_left,l)
+            du_trisk = du_trisk + wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)
+            du_trisk = du_trisk + wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)
+            du_trisk = du_trisk + wee(ij+u_lup,4,1)*hflux(ij+u_right,l)
+            du_trisk = du_trisk + wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)
+            du_trisk = du_trisk + wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)
+            du_trisk = du_trisk + wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)
+            du_trisk = du_trisk + wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)
+            du_trisk = du_trisk + wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)
+            du_trisk = du_trisk + wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
+            du(ij+u_lup,l) = du(ij+u_lup,l) + du_trisk*qu(ij+u_lup,l)
+            du_trisk=0.
+            du_trisk = du_trisk + wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)
+            du_trisk = du_trisk + wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)
+            du_trisk = du_trisk + wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)
+            du_trisk = du_trisk + wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)
+            du_trisk = du_trisk + wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)
+            du_trisk = du_trisk + wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)
+            du_trisk = du_trisk + wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)
+            du_trisk = du_trisk + wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)
+            du_trisk = du_trisk + wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)
+            du_trisk = du_trisk + wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
+            du(ij+u_ldown,l) = du(ij+u_ldown,l) + du_trisk*qu(ij+u_ldown,l)
+         END DO
+      END DO
+   END SELECT
    !---------------------------- coriolis ----------------------------------
    !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/coriolis.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/coriolis.k90	(revision 939)
+++ codes/icosagcm/devel/src/kernels_unst/coriolis.k90	(revision 940)
@@ -176,75 +176,150 @@
    !$OMP END DO
    !
-   !$OMP DO SCHEDULE(STATIC)
-   DO edge = 1, edge_num
-      ! this VLOOP iterates over the TRISK stencil
-      SELECT CASE(trisk_deg(edge))
-      CASE(4)
-         !DIR$ SIMD
-         DO l = 1, llm
-            du_trisk=0.
-            itrisk = 1
-            edge_trisk = trisk(1,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 2
-            edge_trisk = trisk(2,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 3
-            edge_trisk = trisk(3,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 4
-            edge_trisk = trisk(4,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            du(l,edge) = du(l,edge) + .5*du_trisk
-         END DO
-      CASE(10)
-         !DIR$ SIMD
-         DO l = 1, llm
-            du_trisk=0.
-            itrisk = 1
-            edge_trisk = trisk(1,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 2
-            edge_trisk = trisk(2,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 3
-            edge_trisk = trisk(3,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 4
-            edge_trisk = trisk(4,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 5
-            edge_trisk = trisk(5,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 6
-            edge_trisk = trisk(6,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 7
-            edge_trisk = trisk(7,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 8
-            edge_trisk = trisk(8,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 9
-            edge_trisk = trisk(9,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            itrisk = 10
-            edge_trisk = trisk(10,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            du(l,edge) = du(l,edge) + .5*du_trisk
-         END DO
-      CASE DEFAULT
-         !DIR$ SIMD
-         DO l = 1, llm
-            du_trisk=0.
-            DO itrisk = 1, trisk_deg(edge)
-               edge_trisk = trisk(itrisk,edge)
-               du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-            END DO
-            du(l,edge) = du(l,edge) + .5*du_trisk
-         END DO
-      END SELECT
-   END DO
-   !$OMP END DO
+   SELECT CASE(caldyn_conserv)
+   CASE(conserv_energy) ! energy-conserving TRiSK
+      !$OMP DO SCHEDULE(STATIC)
+      DO edge = 1, edge_num
+         ! this VLOOP iterates over the TRISK stencil
+         SELECT CASE(trisk_deg(edge))
+         CASE(4)
+            !DIR$ SIMD
+            DO l = 1, llm
+               du_trisk=0.
+               itrisk = 1
+               edge_trisk = trisk(1,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 2
+               edge_trisk = trisk(2,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 3
+               edge_trisk = trisk(3,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 4
+               edge_trisk = trisk(4,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               du(l,edge) = du(l,edge) + .5*du_trisk
+            END DO
+         CASE(10)
+            !DIR$ SIMD
+            DO l = 1, llm
+               du_trisk=0.
+               itrisk = 1
+               edge_trisk = trisk(1,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 2
+               edge_trisk = trisk(2,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 3
+               edge_trisk = trisk(3,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 4
+               edge_trisk = trisk(4,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 5
+               edge_trisk = trisk(5,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 6
+               edge_trisk = trisk(6,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 7
+               edge_trisk = trisk(7,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 8
+               edge_trisk = trisk(8,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 9
+               edge_trisk = trisk(9,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               itrisk = 10
+               edge_trisk = trisk(10,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               du(l,edge) = du(l,edge) + .5*du_trisk
+            END DO
+         CASE DEFAULT
+            !DIR$ SIMD
+            DO l = 1, llm
+               du_trisk=0.
+               DO itrisk = 1, trisk_deg(edge)
+                  edge_trisk = trisk(itrisk,edge)
+                  du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+               END DO
+               du(l,edge) = du(l,edge) + .5*du_trisk
+            END DO
+         END SELECT
+      END DO
+      !$OMP END DO
+   CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK
+      !$OMP DO SCHEDULE(STATIC)
+      DO edge = 1, edge_num
+         ! this VLOOP iterates over the TRISK stencil
+         SELECT CASE(trisk_deg(edge))
+         CASE(4)
+            !DIR$ SIMD
+            DO l = 1, llm
+               du_trisk=0.
+               itrisk = 1
+               edge_trisk = trisk(1,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 2
+               edge_trisk = trisk(2,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 3
+               edge_trisk = trisk(3,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 4
+               edge_trisk = trisk(4,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               du(l,edge) = du(l,edge) + du_trisk*qu(l,edge)
+            END DO
+         CASE(10)
+            !DIR$ SIMD
+            DO l = 1, llm
+               du_trisk=0.
+               itrisk = 1
+               edge_trisk = trisk(1,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 2
+               edge_trisk = trisk(2,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 3
+               edge_trisk = trisk(3,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 4
+               edge_trisk = trisk(4,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 5
+               edge_trisk = trisk(5,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 6
+               edge_trisk = trisk(6,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 7
+               edge_trisk = trisk(7,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 8
+               edge_trisk = trisk(8,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 9
+               edge_trisk = trisk(9,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               itrisk = 10
+               edge_trisk = trisk(10,edge)
+               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               du(l,edge) = du(l,edge) + du_trisk*qu(l,edge)
+            END DO
+         CASE DEFAULT
+            !DIR$ SIMD
+            DO l = 1, llm
+               du_trisk=0.
+               DO itrisk = 1, trisk_deg(edge)
+                  edge_trisk = trisk(itrisk,edge)
+                  du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)
+               END DO
+               du(l,edge) = du(l,edge) + du_trisk*qu(l,edge)
+            END DO
+         END SELECT
+      END DO
+      !$OMP END DO
+   END SELECT
    !---------------------------- coriolis ----------------------------------
    !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 939)
+++ codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 940)
@@ -163,4 +163,5 @@
   FIELD_MASS  :: convm
   FIELD_THETA :: theta, dtheta_rhodz
+  INTEGER, PARAMETER :: conserv_energy=1, conserv_enstrophy=2, caldyn_conserv = conserv_enstrophy ! FIXME
   DECLARE_INDICES
   DECLARE_EDGES
@@ -211,4 +212,26 @@
 !----------------------------- Unused -----------------------------
 
+#ifdef BEGIN_DYSL
+KERNEL(gradient)
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      grad(EDGE) = SIGN*(b(CELL2)-b(CELL1))
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+KERNEL(div)
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      div_ij=0.
+      FORALL_EDGES
+        div_ij = div_ij + SIGN*LE_DE*u(EDGE)
+      END_BLOCK
+      divu(CELL) = div_ij / AI
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+#endif END_DYSL
+
 SUBROUTINE gradient(b,grad) BINDC(gradient)
   FIELD_MASS :: b
Index: codes/icosagcm/devel/src/unstructured/data_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 939)
+++ codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 940)
@@ -43,6 +43,6 @@
   NUM1(max_nb_stage), BIND(C)              :: tauj       ! diagonal of fast Butcher tableau
   NUM2(max_nb_stage,max_nb_stage), BIND(C) :: cslj, cflj ! slow and fast modified Butcher tableaus
-  NUM2(:,:), POINTER          :: centroid, xyz_v, Riv2, wee, ap,bp, mass_bl, mass_dak, mass_dbk
-
+  NUM2(:,:), POINTER          :: centroid, xyz_v, Riv2, ap,bp, mass_bl, mass_dak, mass_dbk
+  NUM3(:,:,:), POINTER        :: wee
   INTEGER(C_INT), BIND(C) :: comm_icosa
 
@@ -147,4 +147,5 @@
 #define ALLOC1(v,n1) IF(ASSOCIATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1))
 #define ALLOC2(v,n1,n2) IF(ASSOCIATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1,n2))
+#define ALLOC3(v,n1,n2,n3) IF(ASSOCIATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1,n2,n3))
 
   SUBROUTINE init_mesh( & 
@@ -225,5 +226,5 @@
     ALLOC1(le_de,edge_num)
     ALLOC2(Riv2, max_dual_deg, dual_num)
-    ALLOC2(wee, max_trisk_deg, edge_num)
+    ALLOC3(wee, max_trisk_deg, edge_num, 1)
     Ai(:) = Ai_(:)
     Av(:) = Av_(:)
@@ -231,5 +232,5 @@
     le_de(:) = le_de_(:)
     Riv2(:,:)=Riv2_(:,:)
-    wee(:,:) = wee_(:,:)
+    wee(:,:,1) = wee_(:,:)
     IF(is_mpi_master) THEN
        PRINT *, 'Max Ai : ',    MAXVAL(ABS(Ai))
Index: codes/icosagcm/devel/src/unstructured/init_unstructured.f90
===================================================================
--- codes/icosagcm/devel/src/unstructured/init_unstructured.f90	(revision 939)
+++ codes/icosagcm/devel/src/unstructured/init_unstructured.f90	(revision 940)
@@ -282,5 +282,6 @@
     ALLOCATE(Riv2, source = Ddata_read2)
     CALL read_from_gridfile(id_nc, 'float', 'wee')
-    ALLOCATE(wee, source = Ddata_read2)
+    ALLOCATE(wee(SIZE(Ddata_read2,1), SIZE(Ddata_read2,2), 1))
+    wee(:,:,1) = Ddata_read2(:,:)
 
     ! read cell centers and bounds
