Index: codes/icosagcm/devel/src/base/init_grid_param.f90
===================================================================
--- codes/icosagcm/devel/src/base/init_grid_param.f90	(revision 929)
+++ codes/icosagcm/devel/src/base/init_grid_param.f90	(revision 935)
@@ -7,4 +7,8 @@
   USE compute_caldyn_mod
   USE compute_pvort_only_mod
+  USE compute_theta_mod
+  USE compute_geopot_mod
+  USE compute_caldyn_fast_mod
+  USE compute_caldyn_slow_hydro_mod
 
   IMPLICIT NONE
@@ -64,5 +68,9 @@
     compute_pression_mid => compute_pression_mid_hex
     ! dynamics
-    compute_pvort_only   => compute_pvort_only_hex
+    compute_pvort_only        => compute_pvort_only_hex
+    compute_theta             => compute_theta_hex
+    compute_geopot            => compute_geopot_hex
+    compute_caldyn_fast       => compute_caldyn_fast_hex
+    compute_caldyn_slow_hydro => compute_caldyn_slow_hydro_hex
   END SUBROUTINE select_compute_hex
 
@@ -73,7 +81,10 @@
     compute_pression_mid => compute_pression_mid_unst
     ! dynamics
-    compute_pvort_only   => compute_pvort_only_unst
+    compute_pvort_only         => compute_pvort_only_unst
+    compute_theta              => compute_theta_unst
+!    compute_geopot           => compute_geopot_unst
+    compute_caldyn_fast        => compute_caldyn_fast_unst
+    compute_caldyn_slow_hydro  => compute_caldyn_slow_hydro_unst
   END SUBROUTINE select_compute_unst
-
   
 END MODULE init_grid_param_mod
Index: codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 929)
+++ codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 935)
@@ -5,6 +5,4 @@
   USE compute_caldyn_vert_mod, ONLY : compute_caldyn_vert => compute_caldyn_vert_manual
   USE compute_caldyn_vert_NH_mod, ONLY : compute_caldyn_vert_NH => compute_caldyn_vert_NH_manual
-  USE compute_theta_mod, ONLY : compute_theta => compute_theta_manual
-  USE compute_geopot_mod, ONLY : compute_geopot => compute_geopot_manual
   USE compute_caldyn_kv_mod, ONLY : compute_caldyn_kv
   USE compute_caldyn_Coriolis_mod, ONLY : compute_caldyn_Coriolis => compute_caldyn_Coriolis_manual
@@ -12,5 +10,4 @@
   USE compute_caldyn_slow_NH_mod, ONLY : compute_caldyn_slow_NH
   USE compute_caldyn_solver_mod, ONLY : compute_caldyn_solver
-  USE compute_caldyn_fast_mod, ONLY : compute_caldyn_fast => compute_caldyn_fast_manual
   USE compute_NH_geopot_mod, ONLY : compute_NH_geopot
   IMPLICIT NONE
@@ -35,5 +32,6 @@
     USE output_field_mod
     USE checksum_mod
-    USE compute_caldyn_mod, ONLY : compute_pvort_only
+    USE compute_caldyn_mod, ONLY : compute_pvort_only, compute_theta, &
+         compute_geopot, compute_caldyn_fast
     IMPLICIT NONE
     LOGICAL,INTENT(IN)    :: write_out
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn.f90	(revision 929)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn.f90	(revision 935)
@@ -1,20 +1,68 @@
 MODULE compute_caldyn_mod
+  USE prec, ONLY : rstd
   IMPLICIT NONE
   SAVE
+  
+  ! fake array dimensions, for interfaces
+  INTEGER, PARAMETER :: iim_jjm_i=1, iim_jjm_u=1, iim_jjm_v =1, llm_=1, llm1=1, nqdyn_=1
 
   INTERFACE
+     
+    SUBROUTINE comp_pvort_only(u,rhodz,qu,qv, hv)
+      IMPORT
+      REAL(rstd),INTENT(IN)  :: u(iim_jjm_i, llm_)
+      REAL(rstd),INTENT(INOUT) :: rhodz(iim_jjm_i, llm_)
+      REAL(rstd),INTENT(OUT) :: qu(iim_jjm_u, llm_)
+      REAL(rstd),INTENT(OUT) :: qv(iim_jjm_v, llm_)
+      REAL(rstd),INTENT(OUT) :: hv(iim_jjm_v, llm_)
+    END SUBROUTINE comp_pvort_only
 
-    SUBROUTINE comp_pvort_only(u,rhodz,qu,qv, hv_)
-      USE prec, ONLY : rstd
-      REAL(rstd),INTENT(IN)  :: u(1,1)
-      REAL(rstd),INTENT(INOUT) :: rhodz(1,1)
-      REAL(rstd),INTENT(OUT) :: qu(1,1)
-      REAL(rstd),INTENT(OUT) :: qv(1,1)
-      REAL(rstd),INTENT(OUT) :: hv_(1,1)
-    END SUBROUTINE comp_pvort_only
+    SUBROUTINE comp_theta(mass_col,theta_rhodz, rhodz,theta)
+      IMPORT
+      REAL(rstd),INTENT(IN)    :: mass_col(iim_jjm_i)
+      REAL(rstd),INTENT(IN)    :: theta_rhodz(iim_jjm_i, llm_, nqdyn_)
+      REAL(rstd),INTENT(INOUT) :: rhodz(iim_jjm_i, llm_)
+      REAL(rstd),INTENT(OUT)   :: theta(iim_jjm_i, llm_, nqdyn_)
+    END SUBROUTINE comp_theta
+
+    SUBROUTINE comp_geopot(rhodz,theta, ps,pk,geopot) 
+      IMPORT
+      REAL(rstd),INTENT(IN)    :: rhodz(iim_jjm_i, llm_)         ! rho*dz = mass per unit surface in each full model level
+      REAL(rstd),INTENT(IN)    :: theta(iim_jjm_i, llm_, nqdyn_) ! active scalars : theta/entropy, moisture, ...
+      REAL(rstd),INTENT(INOUT) :: ps(iim_jjm_i)                  ! surface pressure
+      REAL(rstd),INTENT(OUT)   :: pk(iim_jjm_i, llm_)            ! Exner function (compressible) /Lagrange multiplier (Boussinesq)
+      REAL(rstd),INTENT(INOUT) :: geopot(iim_jjm_i, llm1)        ! geopotential
+    END SUBROUTINE comp_geopot
+
+    SUBROUTINE comp_caldyn_fast(tau,theta,geopot, pk,berni,du,u)
+      IMPORT
+      REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
+      REAL(rstd),INTENT(IN)    :: theta(iim_jjm_i, llm_, nqdyn_)
+      REAL(rstd),INTENT(IN)    :: geopot(iim_jjm_i, llm1)
+      REAL(rstd),INTENT(INOUT) :: pk(iim_jjm_i, llm_)
+      REAL(rstd),INTENT(INOUT) :: berni(iim_jjm_i, llm_)  ! partial Bernoulli function
+      REAL(rstd),INTENT(INOUT) :: du(iim_jjm_u, llm_)
+      REAL(rstd),INTENT(INOUT) :: u(iim_jjm_u, llm_)   ! INOUT if tau>0
+    END SUBROUTINE comp_caldyn_fast
+
+  SUBROUTINE comp_caldyn_slow_hydro(zero, u,rhodz,hv,Kv, berni, hflux,du)
+    IMPORT
+    LOGICAL, INTENT(IN) :: zero
+    REAL(rstd),INTENT(IN)  :: u(iim_jjm_u, llm_)    ! prognostic "velocity"
+    REAL(rstd),INTENT(IN)  :: rhodz(iim_jjm_i, llm_)
+    REAL(rstd),INTENT(IN)  :: hv(iim_jjm_v, llm_)   ! height/mass averaged to vertices
+    REAL(rstd),INTENT(IN)  :: Kv(iim_jjm_v, llm_)   ! kinetic energy at vertices
+    REAL(rstd), INTENT(OUT) :: berni(iim_jjm_i, llm_)  ! Bernoulli function
+    REAL(rstd),INTENT(OUT) :: hflux(iim_jjm_u, llm_) ! hflux in kg/s 
+    REAL(rstd),INTENT(INOUT) :: du(iim_jjm_u, llm_)
+  END SUBROUTINE comp_caldyn_slow_hydro
 
   END INTERFACE
 
-  PROCEDURE(comp_pvort_only), POINTER :: compute_pvort_only => NULL()
+  PROCEDURE(comp_pvort_only),        POINTER :: compute_pvort_only        => NULL()
+  PROCEDURE(comp_theta),             POINTER :: compute_theta             => NULL()
+  PROCEDURE(comp_geopot),            POINTER :: compute_geopot            => NULL()
+  PROCEDURE(comp_caldyn_fast),       POINTER :: compute_caldyn_fast       => NULL()
+  PROCEDURE(comp_caldyn_slow_hydro), POINTER :: compute_caldyn_slow_hydro => NULL()
 
 END MODULE compute_caldyn_mod
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90	(revision 929)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90	(revision 935)
@@ -17,40 +17,38 @@
 #ifdef BEGIN_DYSL
 
-KERNEL(caldyn_fast)
-!
+{% macro case_caldyn_fast(name) %}
+  CASE({{name}})
+    FORALL_CELLS()
+      ON_PRIMAL
+        Phi_ik = .5*(geopot(CELL)+geopot(UP(CELL)))
+        {{ caller() }}
+      END_BLOCK
+    END_BLOCK
+{%- endmacro %}
+
+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)
+     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
+  {% call case_caldyn_fast('thermo_theta') %}
+     berni(CELL) = Phi_ik
+  {%- endcall %}
+  {% call case_caldyn_fast('thermo_entropy') %}
+     berni(CELL) = Phi_ik + pk(CELL)*(cpp-theta(CELL,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
+  {%- endcall %}
+  {% call case_caldyn_fast('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
+     cp_ik = cpp*(pk(ij,l)/Treff)**nu
+     berni(CELL) = Phi_ik + pk(CELL)*(cp_ik/(nu+1.)-theta(CELL,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
+  {%- endcall %}
   CASE DEFAULT
     PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo  ! FIXME
@@ -83,5 +81,5 @@
     DECLARE_INDICES
     DECLARE_EDGES
-    NUM          :: due, cp_ik
+    NUM          :: due, cp_ik, Phi_ik
     START_TRACE(id_fast, 4,0,2) ! primal, dual, edge
 #include "../kernels_unst/caldyn_fast.k90"
@@ -101,9 +99,10 @@
 
     INTEGER :: ij,l
-    REAL(rstd) :: cp_ik, qv, temp, chi, log_p_preff, due, due_right, due_lup, due_ldown
+    REAL(rstd) :: due, cp_ik, Phi_ik
 
     CALL trace_start("compute_caldyn_fast")
 #include "../kernels_hex/caldyn_fast.k90"
     CALL trace_end("compute_caldyn_fast")
+    
   END SUBROUTINE compute_caldyn_fast_hex
 
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_hydro.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_hydro.F90	(revision 929)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_hydro.F90	(revision 935)
@@ -73,11 +73,10 @@
     LOGICAL, INTENT(IN) :: zero
     REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
+    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: hv(2*iim*jjm,llm)   ! height/mass averaged to vertices
     REAL(rstd),INTENT(IN)  :: Kv(2*iim*jjm,llm)   ! kinetic energy at vertices
-    REAL(rstd),INTENT(IN)  :: hv(2*iim*jjm,llm)   ! height/mass averaged to vertices
-    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
+    REAL(rstd), INTENT(OUT) :: berni(iim*jjm,llm)  ! Bernoulli function
     REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
     REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
-    
-    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
     REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
     REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
Index: codes/icosagcm/devel/src/kernels_hex/caldyn_fast.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_hex/caldyn_fast.k90	(revision 929)
+++ codes/icosagcm/devel/src/kernels_hex/caldyn_fast.k90	(revision 935)
@@ -16,5 +16,6 @@
          !DIR$ SIMD
          DO ij=ij_begin, ij_end
-            berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+            Phi_ik = .5*(geopot(ij,l)+geopot(ij,l+1))
+            berni(ij,l) = Phi_ik
          END DO
       END DO
@@ -23,18 +24,18 @@
          !DIR$ SIMD
          DO ij=ij_begin, ij_end
-            berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-            berni(ij,l) = berni(ij,l) + pk(ij,l)*(cpp-theta(ij,l,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
+            Phi_ik = .5*(geopot(ij,l)+geopot(ij,l+1))
+            berni(ij,l) = Phi_ik + pk(ij,l)*(cpp-theta(ij,l,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
          END DO
       END DO
    CASE(thermo_variable_Cp)
-      ! thermodynamics with variable Cp
-      ! Cp(T) = Cp0 * (T/T0)^nu
-      ! => h = Cp(T).T/(nu+1)
       DO l = ll_begin, ll_end
          !DIR$ SIMD
          DO ij=ij_begin, ij_end
-            berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+            Phi_ik = .5*(geopot(ij,l)+geopot(ij,l+1))
+            ! thermodynamics with variable Cp
+            ! Cp(T) = Cp0 * (T/T0)^nu
+            ! => h = Cp(T).T/(nu+1)
             cp_ik = cpp*(pk(ij,l)/Treff)**nu
-            berni(ij,l) = berni(ij,l) + pk(ij,l)*(cp_ik/(nu+1.)-theta(ij,l,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
+            berni(ij,l) = Phi_ik + pk(ij,l)*(cp_ik/(nu+1.)-theta(ij,l,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
          END DO
       END DO
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90	(revision 929)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90	(revision 935)
@@ -19,5 +19,6 @@
          !DIR$ SIMD
          DO l = 1, llm
-            berni(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij))
+            Phi_ik = .5*(geopot(l,ij)+geopot(l+1,ij))
+            berni(l,ij) = Phi_ik
          END DO
       END DO
@@ -28,20 +29,20 @@
          !DIR$ SIMD
          DO l = 1, llm
-            berni(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij))
-            berni(l,ij) = berni(l,ij) + pk(l,ij)*(cpp-theta(l,ij,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
+            Phi_ik = .5*(geopot(l,ij)+geopot(l+1,ij))
+            berni(l,ij) = Phi_ik + pk(l,ij)*(cpp-theta(l,ij,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
          END DO
       END DO
       !$OMP END DO
    CASE(thermo_variable_Cp)
-      ! thermodynamics with variable Cp
-      ! Cp(T) = Cp0 * (T/T0)^nu
-      ! => h = Cp(T).T/(nu+1)
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
          !DIR$ SIMD
          DO l = 1, llm
-            berni(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij))
-            cp_ik = cpp*(pk(l,ij)/Treff)**nu
-            berni(l,ij) = berni(l,ij) + pk(l,ij)*(cp_ik/(nu+1.)-theta(l,ij,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
+            Phi_ik = .5*(geopot(l,ij)+geopot(l+1,ij))
+            ! thermodynamics with variable Cp
+            ! Cp(T) = Cp0 * (T/T0)^nu
+            ! => h = Cp(T).T/(nu+1)
+            cp_ik = cpp*(pk(ij,l)/Treff)**nu
+            berni(l,ij) = Phi_ik + pk(l,ij)*(cp_ik/(nu+1.)-theta(l,ij,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
          END DO
       END DO
Index: codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 929)
+++ codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 935)
@@ -203,5 +203,5 @@
   DECLARE_INDICES
   DECLARE_EDGES
-  NUM          :: due, cp_ik
+  NUM          :: due, cp_ik, Phi_ik
   START_TRACE(id_fast, 4,0,2) ! primal, dual, edge
 #include "../kernels_unst/caldyn_fast.k90"
