Index: codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 852)
+++ codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 853)
@@ -10,4 +10,5 @@
   USE compute_caldyn_slow_hydro_mod, ONLY : compute_caldyn_slow_hydro
   USE compute_caldyn_slow_NH_mod, ONLY : compute_caldyn_slow_NH
+  USE compute_caldyn_solver_mod, ONLY : compute_caldyn_solver
   IMPLICIT NONE
   PRIVATE
Index: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 852)
+++ codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 853)
@@ -13,5 +13,5 @@
   LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
 
-  PUBLIC :: compute_caldyn_solver, compute_caldyn_fast
+  PUBLIC :: compute_caldyn_fast,compute_NH_geopot
 
 CONTAINS
@@ -183,131 +183,4 @@
   END SUBROUTINE compute_NH_geopot
 
-  SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
-    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
-    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
-    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
-    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
-    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
-    REAL(rstd),INTENT(OUT)   :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces
-    REAL(rstd),INTENT(OUT)   :: pres(iim*jjm,llm)          ! pressure
-    REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1)
-    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1)
-    REAL(rstd),INTENT(OUT)   :: du(3*iim*jjm,llm)
-
-    REAL(rstd) :: berni(iim*jjm,llm)         ! (W/m_il)^2
-    REAL(rstd) :: berni1(iim*jjm)         ! (W/m_il)^2
-    REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd, Rd_preff
-    INTEGER    :: ij, l
-
-    CALL trace_start("compute_caldyn_solver")
-
-    Rd=cpp*kappa
-
-    IF(dysl) THEN
-
-!$OMP BARRIER
-
-#include "../kernels_hex/caldyn_mil.k90"
-  IF(tau>0) THEN ! solve implicit problem for geopotential
-    CALL compute_NH_geopot(tau,phis, rhodz, m_il, theta, W, geopot)
-  END IF
-#define PHI_BOT(ij) phis(ij)
-#include "../kernels_hex/caldyn_solver.k90"
-#undef PHI_BOT
-!$OMP BARRIER
-
-    ELSE
-
-#define BERNI(ij) berni1(ij)
-    ! FIXME : vertical OpenMP parallelism will not work
-
-    ! average m_ik to interfaces => m_il
-    !DIR$ SIMD
-    DO ij=ij_begin_ext,ij_end_ext
-       m_il(ij,1) = .5*rhodz(ij,1)
-    ENDDO
-    DO l=2,llm
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l))
-       ENDDO
-    ENDDO
-    !DIR$ SIMD
-    DO ij=ij_begin_ext,ij_end_ext
-       m_il(ij,llm+1) = .5*rhodz(ij,llm)
-    ENDDO
-
-    IF(tau>0) THEN ! solve implicit problem for geopotential
-       CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot)
-    END IF
-
-    ! Compute pressure, stored temporarily in pk
-    ! kappa = R/Cp
-    ! 1-kappa = Cv/Cp
-    ! Cp/Cv = 1/(1-kappa)
-    gamma = 1./(1.-kappa)
-    DO l=1,llm
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
-          X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij
-          ! kappa.theta.rho = p/exner 
-          ! => X = (p/p0)/(exner/Cp)
-          !      = (p/p0)^(1-kappa)
-          pk(ij,l) = preff*(X_ij**gamma)
-       ENDDO
-    ENDDO
-
-    ! Update W, compute tendencies
-    DO l=2,llm
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          dW(ij,l)   = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l)
-          W(ij,l)    = W(ij,l)+tau*dW(ij,l) ! update W
-          dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l)
-       ENDDO
-       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
-       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l))) 
-    ENDDO
-    ! Lower BC (FIXME : no orography yet !)
-    DO ij=ij_begin,ij_end         
-       dPhi(ij,1)=0
-       W(ij,1)=0
-       dW(ij,1)=0
-       dPhi(ij,llm+1)=0 ! rigid lid
-       W(ij,llm+1)=0
-       dW(ij,llm+1)=0
-    ENDDO
-    ! Upper BC p=ptop
-    !    DO ij=ij_omp_begin_ext,ij_omp_end_ext
-    !       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
-    !       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)
-    !    ENDDO
-    
-    ! Compute Exner function (needed by compute_caldyn_fast) and du=-g^2.grad(w^2)
-    DO l=1,llm
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
-          BERNI(ij) = (-.25*g*g)*(              &
-                 (W(ij,l)/m_il(ij,l))**2       &
-               + (W(ij,l+1)/m_il(ij,l+1))**2 )
-       ENDDO
-       DO ij=ij_begin,ij_end         
-          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
-          du(ij+u_lup,l)   = ne_lup  *(BERNI(ij)-BERNI(ij+t_lup))
-          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
-       ENDDO
-    ENDDO
-#undef BERNI
-
-    END IF ! dysl
-
-    CALL trace_end("compute_caldyn_solver")
-    
-  END SUBROUTINE compute_caldyn_solver
-  
   SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
     REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90	(revision 853)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90	(revision 853)
@@ -0,0 +1,145 @@
+MODULE compute_caldyn_solver_mod
+  USE grid_param, ONLY : llm
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC :: compute_caldyn_solver
+
+CONTAINS
+
+  SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
+    USE icosa
+    USE caldyn_vars_mod
+    USE trace
+    USE omp_para, ONLY : ll_begin, ll_end,ll_beginp1,ll_endp1
+    USE disvert_mod, ONLY : ptop
+    USE caldyn_kernels_hevi_mod
+    REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6
+    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
+    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
+    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
+    REAL(rstd),INTENT(OUT)   :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces
+    REAL(rstd),INTENT(OUT)   :: pres(iim*jjm,llm)          ! pressure
+    REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT)   :: du(3*iim*jjm,llm)
+
+    REAL(rstd) :: berni(iim*jjm,llm)         ! (W/m_il)^2
+    REAL(rstd) :: berni1(iim*jjm)         ! (W/m_il)^2
+    !REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd, Rd_preff
+    REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Cvd, Rd_preff
+    INTEGER    :: ij, l
+
+    CALL trace_start("compute_caldyn_solver")
+
+    !Rd=cpp*kappa
+
+    IF(dysl) THEN
+
+!$OMP BARRIER
+
+#include "../kernels_hex/caldyn_mil.k90"
+  IF(tau>0) THEN ! solve implicit problem for geopotential
+    CALL compute_NH_geopot(tau,phis, rhodz, m_il, theta, W, geopot)
+  END IF
+#define PHI_BOT(ij) phis(ij)
+#include "../kernels_hex/caldyn_solver.k90"
+#undef PHI_BOT
+!$OMP BARRIER
+
+    ELSE
+
+#define BERNI(ij) berni1(ij)
+    ! FIXME : vertical OpenMP parallelism will not work
+
+    ! average m_ik to interfaces => m_il
+    !DIR$ SIMD
+    DO ij=ij_begin_ext,ij_end_ext
+       m_il(ij,1) = .5*rhodz(ij,1)
+    ENDDO
+    DO l=2,llm
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l))
+       ENDDO
+    ENDDO
+    !DIR$ SIMD
+    DO ij=ij_begin_ext,ij_end_ext
+       m_il(ij,llm+1) = .5*rhodz(ij,llm)
+    ENDDO
+
+    IF(tau>0) THEN ! solve implicit problem for geopotential
+       CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot)
+    END IF
+
+    ! Compute pressure, stored temporarily in pk
+    ! kappa = R/Cp
+    ! 1-kappa = Cv/Cp
+    ! Cp/Cv = 1/(1-kappa)
+    gamma = 1./(1.-kappa)
+    DO l=1,llm
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
+          X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij
+          ! kappa.theta.rho = p/exner 
+          ! => X = (p/p0)/(exner/Cp)
+          !      = (p/p0)^(1-kappa)
+          pk(ij,l) = preff*(X_ij**gamma)
+       ENDDO
+    ENDDO
+
+    ! Update W, compute tendencies
+    DO l=2,llm
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          dW(ij,l)   = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l)
+          W(ij,l)    = W(ij,l)+tau*dW(ij,l) ! update W
+          dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l)
+       ENDDO
+       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
+       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l))) 
+    ENDDO
+    ! Lower BC (FIXME : no orography yet !)
+    DO ij=ij_begin,ij_end         
+       dPhi(ij,1)=0
+       W(ij,1)=0
+       dW(ij,1)=0
+       dPhi(ij,llm+1)=0 ! rigid lid
+       W(ij,llm+1)=0
+       dW(ij,llm+1)=0
+    ENDDO
+    ! Upper BC p=ptop
+    !    DO ij=ij_omp_begin_ext,ij_omp_end_ext
+    !       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
+    !       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)
+    !    ENDDO
+    
+    ! Compute Exner function (needed by compute_caldyn_fast) and du=-g^2.grad(w^2)
+    DO l=1,llm
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
+          BERNI(ij) = (-.25*g*g)*(              &
+                 (W(ij,l)/m_il(ij,l))**2       &
+               + (W(ij,l+1)/m_il(ij,l+1))**2 )
+       ENDDO
+       DO ij=ij_begin,ij_end         
+          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
+          du(ij+u_lup,l)   = ne_lup  *(BERNI(ij)-BERNI(ij+t_lup))
+          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
+       ENDDO
+    ENDDO
+#undef BERNI
+
+    END IF ! dysl
+
+    CALL trace_end("compute_caldyn_solver")
+    
+  END SUBROUTINE compute_caldyn_solver
+  
+END MODULE compute_caldyn_solver_mod
