Index: codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin
===================================================================
--- codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin	(revision 874)
+++ codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin	(revision 876)
@@ -250,65 +250,2 @@
 END_BLOCK
 
-KERNEL(caldyn_slow_NH)
-  FORALL_CELLS_EXT('1', 'llm+1')
-     ON_PRIMAL
-        CST_IF(IS_INNER_INTERFACE, w_il(CELL) = 2.*W(CELL)/(rhodz(KDOWN(CELL))+rhodz(KUP(CELL))) )
-        CST_IF(IS_BOTTOM_LEVEL,    w_il(CELL) = 2.*W(CELL)/rhodz(KUP(CELL)) )
-        CST_IF(IS_TOP_INTERFACE,   w_il(CELL) = 2.*W(CELL)/rhodz(KDOWN(CELL)) )
-     END_BLOCK
-  END_BLOCK
-  FORALL_CELLS_EXT('1', 'llm+1')
-     ON_EDGES
-       ! compute DePhi, v_el, G_el, F_el
-       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
-       ! while DePhil, W_el and F_el do not
-       W_el         = .5*( W(CELL2)+W(CELL1) )
-       DePhil(EDGE) = SIGN*(Phi(CELL2)-Phi(CELL1))
-       F_el(EDGE)   = DePhil(EDGE)*W_el
-       W2_el        = .5*LE_DE * ( W(CELL1)*w_il(CELL1) + W(CELL2)*w_il(CELL2) )
-       v_el(EDGE)   = .5*LE_DE*(u(KUP(EDGE))+u(KDOWN(EDGE))) ! checked
-       G_el(EDGE)   = v_el(EDGE)*W_el - DePhil(EDGE)*W2_el
-     END_BLOCK
-  END_BLOCK
-
-  FORALL_CELLS_EXT('1', 'llm+1')
-    ! compute GradPhi2, dPhi, dW
-    ON_PRIMAL
-       gPhi2=0.
-       dP=0.
-       divG=0
-       FORALL_EDGES       
-         gPhi2 = gPhi2 + LE_DE*DePhil(EDGE)**2
-         dP = dP + LE_DE*DePhil(EDGE)*v_el(EDGE)
-         divG = divG + SIGN*G_el(EDGE) ! -div(G_el), G_el already has le_de
-       END_BLOCK
-       gradPhi2(CELL) = 1./(2.*AI) * gPhi2
-       dPhi(CELL) = gradPhi2(CELL)*w_il(CELL) - 1./(2.*AI)*dP
-       dW(CELL) = (-1./AI)*divG
-    END_BLOCK
-  END_BLOCK
-
-  ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
-  BARRIER
-
-  FORALL_CELLS_EXT()
-    ! Compute berni at scalar points
-    ON_PRIMAL
-       u2=0.
-       FORALL_EDGES       
-         u2 = u2 + LE_DE*u(EDGE)**2
-       END_BLOCK
-       berni(CELL) = 1./(4.*AI) * u2 - .25*( gradPhi2(CELL)*w_il(CELL)**2 + gradPhi2(UP(CELL))*w_il(UP(CELL))**2 )
-    END_BLOCK
-  END_BLOCK
-
-  FORALL_CELLS_EXT()
-    ON_EDGES
-    ! Compute mass flux and grad(berni)
-      uu = .5*(rhodz(CELL1)+rhodz(CELL2))*u(EDGE) - .5*( F_el(EDGE)+F_el(UP(EDGE)) )
-      hflux(EDGE) = LE_DE*uu
-      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
-    END_BLOCK
-  END_BLOCK
-
-END_BLOCK
Index: codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin
===================================================================
--- codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin	(revision 874)
+++ codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin	(revision 876)
@@ -127,42 +127,2 @@
 END_BLOCK
 
-KERNEL(coriolis)
-!
-  DO iq=1,nqdyn
-    FORALL_CELLS_EXT()
-      ON_EDGES
-        Ftheta(EDGE) = .5*(theta(CELL1,iq)+theta(CELL2,iq))*hflux(EDGE)
-      END_BLOCK
-    END_BLOCK
-    FORALL_CELLS()
-      ON_PRIMAL
-        divF=0.
-	FORALL_EDGES
-	  divF = divF + Ftheta(EDGE)*SIGN
-	END_BLOCK
-        dtheta_rhodz(CELL,iq) = -divF / AI
-      END_BLOCK
-    END_BLOCK
-  END DO ! iq
-!
-  FORALL_CELLS()
-    ON_PRIMAL
-      divF=0.
-      FORALL_EDGES
-        divF = divF + hflux(EDGE)*SIGN
-      END_BLOCK
-      convm(CELL) = -divF / AI
-    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
-
-END_BLOCK
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 874)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 876)
@@ -4,7 +4,71 @@
   PRIVATE
 
+#include "../unstructured/unstructured.h90"
+
   PUBLIC :: compute_caldyn_Coriolis
 
 CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(coriolis)
+!
+  DO iq=1,nqdyn
+    FORALL_CELLS_EXT()
+      ON_EDGES
+        Ftheta(EDGE) = .5*(theta(CELL1,iq)+theta(CELL2,iq))*hflux(EDGE)
+      END_BLOCK
+    END_BLOCK
+    FORALL_CELLS()
+      ON_PRIMAL
+        divF=0.
+        FORALL_EDGES
+          divF = divF + Ftheta(EDGE)*SIGN
+        END_BLOCK
+        dtheta_rhodz(CELL,iq) = -divF / AI
+      END_BLOCK
+    END_BLOCK
+  END DO ! iq
+!
+  FORALL_CELLS()
+    ON_PRIMAL
+      divF=0.
+      FORALL_EDGES
+        divF = divF + hflux(EDGE)*SIGN
+      END_BLOCK
+      convm(CELL) = -divF / AI
+    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
+
+END_BLOCK
+
+#endif END_DYSL
+
+  SUBROUTINE compute_coriolis_unst(hflux,theta,qu,Ftheta, convm,dtheta_rhodz,du)
+    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+    USE grid_param, ONLY : nqdyn
+    USE data_unstructured_mod, ONLY : id_coriolis, primal_num, dual_num, edge_num, &
+          left, right,primal_deg,primal_edge,primal_ne,trisk_deg,wee,trisk,Ai, &
+          enter_trace, exit_trace
+    FIELD_U     :: hflux, Ftheta, qu, du
+    FIELD_MASS  :: convm
+    FIELD_THETA :: theta, dtheta_rhodz
+    DECLARE_INDICES
+    DECLARE_EDGES
+    NUM :: divF, du_trisk
+    START_TRACE(id_coriolis, 3,4,0) ! primal, dual, edge
+#include "../kernels_unst/coriolis.k90"
+    STOP_TRACE
+  END SUBROUTINE compute_coriolis_unst
 
   SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_NH.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_NH.F90	(revision 874)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_NH.F90	(revision 876)
@@ -4,7 +4,91 @@
   PRIVATE
 
+#include "../unstructured/unstructured.h90"
+
   PUBLIC :: compute_caldyn_slow_NH
 
 CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(caldyn_slow_NH)
+  FORALL_CELLS_EXT('1', 'llm+1')
+     ON_PRIMAL
+        CST_IF(IS_INNER_INTERFACE, w_il(CELL) = 2.*W(CELL)/(rhodz(KDOWN(CELL))+rhodz(KUP(CELL))) )
+        CST_IF(IS_BOTTOM_LEVEL,    w_il(CELL) = 2.*W(CELL)/rhodz(KUP(CELL)) )
+        CST_IF(IS_TOP_INTERFACE,   w_il(CELL) = 2.*W(CELL)/rhodz(KDOWN(CELL)) )
+     END_BLOCK
+  END_BLOCK
+  FORALL_CELLS_EXT('1', 'llm+1')
+     ON_EDGES
+       ! compute DePhi, v_el, G_el, F_el
+       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+       ! while DePhil, W_el and F_el do not 
+       W_el         = .5*( W(CELL2)+W(CELL1) )
+       DePhil(EDGE) = SIGN*(Phi(CELL2)-Phi(CELL1))
+       F_el(EDGE)   = DePhil(EDGE)*W_el
+       W2_el        = .5*LE_DE * ( W(CELL1)*w_il(CELL1) + W(CELL2)*w_il(CELL2) )
+       v_el(EDGE)   = .5*LE_DE*(u(KUP(EDGE))+u(KDOWN(EDGE))) ! checked
+       G_el(EDGE)   = v_el(EDGE)*W_el - DePhil(EDGE)*W2_el
+     END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT('1', 'llm+1')
+    ! compute GradPhi2, dPhi, dW
+    ON_PRIMAL
+       gPhi2=0.
+       dP=0.
+       divG=0
+       FORALL_EDGES    
+         gPhi2 = gPhi2 + LE_DE*DePhil(EDGE)**2
+         dP = dP + LE_DE*DePhil(EDGE)*v_el(EDGE)
+         divG = divG + SIGN*G_el(EDGE) ! -div(G_el), G_el already has le_de
+       END_BLOCK
+       gradPhi2(CELL) = 1./(2.*AI) * gPhi2
+       dPhi(CELL) = gradPhi2(CELL)*w_il(CELL) - 1./(2.*AI)*dP
+       dW(CELL) = (-1./AI)*divG
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
+  BARRIER
+
+  FORALL_CELLS_EXT()
+    ! Compute berni at scalar points
+    ON_PRIMAL
+       u2=0.
+       FORALL_EDGES    
+         u2 = u2 + LE_DE*u(EDGE)**2
+       END_BLOCK
+       berni(CELL) = 1./(4.*AI) * u2 - .25*( gradPhi2(CELL)*w_il(CELL)**2 + gradPhi2(UP(CELL))*w_il(UP(CELL))**2 )
+    END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT()
+    ON_EDGES
+    ! Compute mass flux and grad(berni)
+      uu = .5*(rhodz(CELL1)+rhodz(CELL2))*u(EDGE) - .5*( F_el(EDGE)+F_el(UP(EDGE)) )
+      hflux(EDGE) = LE_DE*uu
+      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
+
+#endif END_DYSL
+
+  SUBROUTINE compute_caldyn_vert_NH_unst(mass,geopot,W,wflux, eta_dot,wcov,W_etadot, du,dPhi,dW)
+    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+    USE data_unstructured_mod, ONLY : left,right,edge_num,primal_num,dual_num,id_vert_NH, &
+        enter_trace, exit_trace
+    FIELD_MASS   :: mass, eta_dot, wcov, W_etadot ! IN, BUF*3
+    FIELD_GEOPOT :: geopot,W,wflux,dPhi,dW ! IN*3, INOUT*2
+    FIELD_U      :: du ! INOUT
+    DECLARE_INDICES
+    NUM :: w_ij, wflux_ij
+    START_TRACE(id_vert_NH, 6,0,1)
+#include "../kernels_unst/caldyn_vert_NH.k90"
+    STOP_TRACE
+  END SUBROUTINE compute_caldyn_vert_NH_unst
 
   SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
