Index: codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 853)
+++ codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 854)
@@ -11,4 +11,5 @@
   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
   IMPLICIT NONE
   PRIVATE
Index: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 853)
+++ codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 854)
@@ -13,5 +13,5 @@
   LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
 
-  PUBLIC :: compute_caldyn_fast,compute_NH_geopot
+  PUBLIC :: compute_NH_geopot
 
 CONTAINS
@@ -183,126 +183,3 @@
   END SUBROUTINE compute_NH_geopot
 
-  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
-    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
-    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0
-    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
-    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
-    REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm)
-    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
-    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
-
-    INTEGER :: i,j,ij,l
-    REAL(rstd) :: cp_ik, qv, temp, chi, nu, due, due_right, due_lup, due_ldown
-
-    CALL trace_start("compute_caldyn_fast")
-
-    IF(dysl_caldyn_fast) THEN
-#include "../kernels_hex/caldyn_fast.k90"
-    ELSE
-
-    ! Compute Bernoulli term
-    IF(boussinesq) THEN
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             berni(ij,l) = pk(ij,l)
-             ! from now on pk contains the vertically-averaged geopotential
-             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-          END DO
-       END DO
-    ELSE ! compressible
-
-       DO l=ll_begin,ll_end
-          SELECT CASE(caldyn_thermo)
-          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-             END DO
-          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s) 
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
-             END DO
-          CASE(thermo_moist) 
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
-                ! Bd = Phi + gibbs_d
-                ! Bv = Phi + gibbs_v
-                ! pk=temperature, theta=entropy
-                qv = theta(ij,l,2)
-                temp = pk(ij,l)
-                chi = log(temp/Treff)
-                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + temp*(cpp*(1.-chi)+Rd*nu)
-                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + temp*(cppv*(1.-chi)+Rv*nu)
-            END DO
-          END SELECT
-       END DO
-
-    END IF ! Boussinesq/compressible
-
-!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
-    DO l=ll_begin,ll_end
-       IF(caldyn_thermo == thermo_moist) THEN
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
-                       *(pk(ij+t_right,l)-pk(ij,l))             &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
-                       *(berniv(ij+t_right,l)-berniv(ij,l))
-             
-             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
-                       *(pk(ij+t_lup,l)-pk(ij,l))               &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
-                       *(berniv(ij+t_lup,l)-berniv(ij,l))
-             
-             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
-                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
-                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
-             
-             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
-             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
-             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
-             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
-          END DO
-       ELSE
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
-                  *(pk(ij+t_right,l)-pk(ij,l))        &
-                  +  berni(ij+t_right,l)-berni(ij,l)
-             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
-                  *(pk(ij+t_lup,l)-pk(ij,l))          &
-                  +  berni(ij+t_lup,l)-berni(ij,l)
-             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
-                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
-                  +   berni(ij+t_ldown,l)-berni(ij,l)
-             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
-             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
-             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
-             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
-          END DO
-       END IF
-    END DO
-
-    END IF ! dysl
-    CALL trace_end("compute_caldyn_fast")
-
-  END SUBROUTINE compute_caldyn_fast
-
 END MODULE caldyn_kernels_hevi_mod
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90	(revision 854)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_fast.F90	(revision 854)
@@ -0,0 +1,138 @@
+MODULE compute_caldyn_fast_mod
+  USE grid_param, ONLY : llm
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC :: compute_caldyn_fast
+
+CONTAINS
+
+  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
+    USE icosa
+    USE trace
+    USE caldyn_vars_mod
+    USE omp_para, ONLY : ll_begin, ll_end
+    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
+    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0
+    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm)
+    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
+    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
+
+    INTEGER :: i,j,ij,l
+    !REAL(rstd) :: cp_ik, qv, temp, chi, nu, due, due_right, due_lup, due_ldown
+    REAL(rstd) :: cp_ik, qv, temp, chi, due, due_right, due_lup, due_ldown
+
+    CALL trace_start("compute_caldyn_fast")
+
+    IF(dysl_caldyn_fast) THEN
+#include "../kernels_hex/caldyn_fast.k90"
+    ELSE
+
+    ! Compute Bernoulli term
+    IF(boussinesq) THEN
+       DO l=ll_begin,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             berni(ij,l) = pk(ij,l)
+             ! from now on pk contains the vertically-averaged geopotential
+             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+          END DO
+       END DO
+    ELSE ! compressible
+
+       DO l=ll_begin,ll_end
+          SELECT CASE(caldyn_thermo)
+          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+             END DO
+          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s) 
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
+             END DO
+          CASE(thermo_moist) 
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
+                ! Bd = Phi + gibbs_d
+                ! Bv = Phi + gibbs_v
+                ! pk=temperature, theta=entropy
+                qv = theta(ij,l,2)
+                temp = pk(ij,l)
+                chi = log(temp/Treff)
+                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + temp*(cpp*(1.-chi)+Rd*nu)
+                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + temp*(cppv*(1.-chi)+Rv*nu)
+            END DO
+          END SELECT
+       END DO
+
+    END IF ! Boussinesq/compressible
+
+!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
+    DO l=ll_begin,ll_end
+       IF(caldyn_thermo == thermo_moist) THEN
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
+                       *(pk(ij+t_right,l)-pk(ij,l))             &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
+                       *(berniv(ij+t_right,l)-berniv(ij,l))
+             
+             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
+                       *(pk(ij+t_lup,l)-pk(ij,l))               &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
+                       *(berniv(ij+t_lup,l)-berniv(ij,l))
+             
+             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
+                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
+                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
+             
+             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+          END DO
+       ELSE
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
+                  *(pk(ij+t_right,l)-pk(ij,l))        &
+                  +  berni(ij+t_right,l)-berni(ij,l)
+             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
+                  *(pk(ij+t_lup,l)-pk(ij,l))          &
+                  +  berni(ij+t_lup,l)-berni(ij,l)
+             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
+                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
+                  +   berni(ij+t_ldown,l)-berni(ij,l)
+             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+          END DO
+       END IF
+    END DO
+
+    END IF ! dysl
+    CALL trace_end("compute_caldyn_fast")
+
+  END SUBROUTINE compute_caldyn_fast
+
+END MODULE compute_caldyn_fast_mod
