Index: codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin
===================================================================
--- codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin	(revision 857)
+++ codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin	(revision 858)
@@ -127,55 +127,4 @@
 END_BLOCK
 
-KERNEL(caldyn_fast)
-!
-  SELECT CASE(caldyn_thermo)
-  CASE(thermo_boussinesq)
-    FORALL_CELLS()
-      ON_PRIMAL
-        berni(CELL) = pk(CELL)
-        ! from now on pk contains the vertically-averaged geopotential
-        pk(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
-      END_BLOCK
-    END_BLOCK
-  CASE(thermo_theta)
-    FORALL_CELLS()
-      ON_PRIMAL
-        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
-      END_BLOCK
-    END_BLOCK
-  CASE(thermo_entropy)
-    FORALL_CELLS()
-      ON_PRIMAL
-        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
-        berni(CELL) = berni(CELL) + pk(CELL)*(cpp-theta(CELL,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
-      END_BLOCK
-    END_BLOCK
-  CASE(thermo_variable_Cp)
-    ! thermodynamics with variable Cp
-    !           Cp(T) = Cp0 * (T/T0)^nu
-    ! =>            h = Cp(T).T/(nu+1)
-
-    FORALL_CELLS()
-      ON_PRIMAL
-        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
-        cp_ik = cpp*(pk(CELL)/Treff)**nu
-        berni(CELL) = berni(CELL) + pk(CELL)*(cp_ik/(nu+1.)-theta(CELL,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
-      END_BLOCK
-    END_BLOCK
-  CASE DEFAULT
-    PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo  ! FIXME
-    STOP
-  END SELECT
-!
-  FORALL_CELLS()
-    ON_EDGES
-      due = .5*(theta(CELL1,1)+theta(CELL2,1))*(pk(CELL2)-pk(CELL1)) + berni(CELL2)-berni(CELL1)
-      du(EDGE) = du(EDGE) - SIGN*due
-      u(EDGE) = u(EDGE) + tau*du(EDGE)
-    END_BLOCK
-  END_BLOCK  
-!
-END_BLOCK
-
 KERNEL(coriolis)
 !
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90	(revision 857)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90	(revision 858)
@@ -4,7 +4,85 @@
   PRIVATE
 
+#include "../unstructured/unstructured.h90"
+
   PUBLIC :: compute_caldyn_fast
 
 CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(caldyn_fast)
+!
+  SELECT CASE(caldyn_thermo)
+  CASE(thermo_boussinesq)
+    FORALL_CELLS()
+      ON_PRIMAL
+        berni(CELL) = pk(CELL)
+        ! from now on pk contains the vertically-averaged geopotential
+        pk(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
+      END_BLOCK
+    END_BLOCK
+  CASE(thermo_theta)
+    FORALL_CELLS()
+      ON_PRIMAL
+        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
+      END_BLOCK
+    END_BLOCK
+  CASE(thermo_entropy)
+    FORALL_CELLS()
+      ON_PRIMAL
+        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
+        berni(CELL) = berni(CELL) + pk(CELL)*(cpp-theta(CELL,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
+      END_BLOCK
+    END_BLOCK
+  CASE(thermo_variable_Cp)
+    ! thermodynamics with variable Cp
+    !           Cp(T) = Cp0 * (T/T0)^nu
+    ! =>            h = Cp(T).T/(nu+1)
+
+    FORALL_CELLS()
+      ON_PRIMAL
+        berni(CELL) = .5*(geopot(CELL)+geopot(UP(CELL)))
+        cp_ik = cpp*(pk(CELL)/Treff)**nu
+        berni(CELL) = berni(CELL) + pk(CELL)*(cp_ik/(nu+1.)-theta(CELL,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
+      END_BLOCK
+    END_BLOCK
+  CASE DEFAULT
+    PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo  ! FIXME
+    STOP
+  END SELECT
+!
+  FORALL_CELLS()
+    ON_EDGES
+      due = .5*(theta(CELL1,1)+theta(CELL2,1))*(pk(CELL2)-pk(CELL1)) + berni(CELL2)-berni(CELL1)
+      du(EDGE) = du(EDGE) - SIGN*due
+      u(EDGE) = u(EDGE) + tau*du(EDGE)
+    END_BLOCK
+  END_BLOCK  
+!
+END_BLOCK
+
+#endif END_DYSL
+
+  SUBROUTINE compute_caldyn_fast_unst(tau, pk,berni,theta,geopot, du,u)
+    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+    USE earth_const
+    USE grid_param, ONLY : nqdyn
+    USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
+         id_fast, primal_num, dual_num, edge_num, &
+         dual_deg, dual_edge, dual_ne, dual_vertex, &
+         up, down, left, right, Av, fv, Riv2
+    NUM, INTENT(IN) :: tau
+    FIELD_MASS   :: pk,berni  ! INOUT, OUT
+    FIELD_THETA  :: theta     ! IN
+    FIELD_GEOPOT :: geopot    ! IN
+    FIELD_U      :: u,du      ! INOUT,INOUT
+    DECLARE_INDICES
+    DECLARE_EDGES
+    NUM          :: due, cp_ik
+    START_TRACE(id_fast, 4,0,2) ! primal, dual, edge
+#include "../kernels_unst/caldyn_fast.k90"
+    STOP_TRACE
+  END SUBROUTINE compute_caldyn_fast_unst
 
   SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
