Index: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 656)
+++ codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 657)
@@ -118,5 +118,5 @@
     REAL(rstd) :: D_il(iim*jjm,llm+1)    ! Thomas algorithm
     REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij
-    REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik
+    REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik, vreff
 
     INTEGER    :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext
@@ -282,5 +282,5 @@
     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
+    REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd, Rd_preff
     INTEGER    :: ij, l
 
@@ -292,8 +292,11 @@
 
 !$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)
-#define PHI_BOT_VAR phis
 #include "../kernels_hex/caldyn_solver.k90"
-#undef PHI_BOT_VAR
 #undef PHI_BOT
 !$OMP BARRIER
Index: codes/icosagcm/devel/src/kernels_hex/caldyn_mil.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_hex/caldyn_mil.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_hex/caldyn_mil.k90	(revision 657)
@@ -0,0 +1,22 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_mil ----------------------------------
+   IF (ll_begin==1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         m_il(ij,1) = .5*rhodz(ij,1)
+      END DO
+   END IF
+   DO l = ll_beginp1, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         m_il(ij,l) = .5*(rhodz(ij,l)+rhodz(ij,l-1))
+      END DO
+   END DO
+   IF(ll_endp1==llm+1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         m_il(ij,llm+1) = .5*rhodz(ij,llm+1 -1)
+      END DO
+   END IF
+   !---------------------------- caldyn_mil ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_hex/caldyn_solver.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_hex/caldyn_solver.k90	(revision 656)
+++ codes/icosagcm/devel/src/kernels_hex/caldyn_solver.k90	(revision 657)
@@ -1,26 +1,4 @@
    !--------------------------------------------------------------------------
    !---------------------------- caldyn_solver ----------------------------------
-   IF (ll_begin==1) THEN
-      !DIR$ SIMD
-      DO ij=ij_begin_ext, ij_end_ext
-         m_il(ij,1) = .5*rhodz(ij,1)
-      END DO
-   END IF
-   DO l = ll_beginp1, ll_end
-      !DIR$ SIMD
-      DO ij=ij_begin_ext, ij_end_ext
-         m_il(ij,l) = .5*(rhodz(ij,l)+rhodz(ij,l-1))
-      END DO
-   END DO
-   IF(ll_endp1==llm+1) THEN
-      !DIR$ SIMD
-      DO ij=ij_begin_ext, ij_end_ext
-         m_il(ij,llm+1) = .5*rhodz(ij,llm+1 -1)
-      END DO
-   END IF
-   !
-   IF(tau>0) THEN ! solve implicit problem for geopotential
-      CALL compute_NH_geopot(tau,PHI_BOT_VAR, rhodz, m_il, theta, W, geopot)
-   END IF
    !
    ! Compute pressure (pres) and Exner function (pk)
@@ -30,5 +8,6 @@
    gamma = 1./(1.-kappa)
    vreff = Rd*Treff/preff ! reference specific volume
-   Cvd = cpp-Rd
+   Cvd = 1./(cpp-Rd)
+   Rd_preff = kappa*cpp/preff
    DO l = ll_begin, ll_end
       !DIR$ SIMD
@@ -36,6 +15,7 @@
          SELECT CASE(caldyn_thermo)
          CASE(thermo_theta)
-            rho_ij = g*rhodz(ij,l)/(geopot(ij,l+1)-geopot(ij,l))
-            X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij
+            rho_ij = 1./(geopot(ij,l+1)-geopot(ij,l))
+            rho_ij = rho_ij*g*rhodz(ij,l)
+            X_ij = Rd_preff*theta(ij,l,1)*rho_ij
             ! kappa.theta.rho = p/exner
             ! => X = (p/p0)/(exner/Cp)
@@ -46,6 +26,7 @@
             pk(ij,l) = cpp*((pres(ij,l)/preff)**kappa) ! Exner
          CASE(thermo_entropy)
-            rho_ij = g*rhodz(ij,l)/(geopot(ij,l+1)-geopot(ij,l))
-            T_ij = Treff*exp( (theta(ij,l,1)+Rd*log(vreff*rho_ij))/Cvd )
+            rho_ij = 1./(geopot(ij,l+1)-geopot(ij,l))
+            rho_ij = rho_ij*g*rhodz(ij,l)
+            T_ij = Treff*exp( (theta(ij,l,1)+Rd*log(vreff*rho_ij))*Cvd )
             pres(ij,l) = rho_ij*Rd*T_ij
             pk(ij,l) = T_ij
Index: codes/icosagcm/devel/src/kernels_hex/compute_NH_geopot.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_hex/compute_NH_geopot.k90	(revision 656)
+++ codes/icosagcm/devel/src/kernels_hex/compute_NH_geopot.k90	(revision 657)
@@ -4,4 +4,5 @@
    g2=g*g
    gm2 = 1./g2
+   vreff = Treff*cpp/preff*kappa
    gamma = 1./(1.-kappa)
    !$OMP BARRIER
@@ -30,7 +31,6 @@
             DO ij=ij_omp_begin_ext,ij_omp_end_ext
                rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l))
-               X_ij = Treff*exp(theta(ij,l)/cpp) ! theta = Tref.exp(s/Cp)
-               X_ij = (cpp/preff)*kappa*X_ij*rho_ij
-               p_ik(ij,l) = preff*(X_ij**gamma)
+               X_ij = log(vreff*rho_ij) + theta(ij,l)/cpp
+               p_ik(ij,l) = preff*exp(X_ij*gamma)
                c2_mik = gamma*p_ik(ij,l)/(rho_ij*m_ik(ij,l)) ! c^2 = gamma*R*T = gamma*p/rho
                A_ik(ij,l) = c2_mik*(tau/g*rho_ij)**2
