Index: codes/icosagcm/trunk/src/base/earth_const.f90
===================================================================
--- codes/icosagcm/trunk/src/base/earth_const.f90	(revision 579)
+++ codes/icosagcm/trunk/src/base/earth_const.f90	(revision 580)
@@ -28,4 +28,6 @@
   !$OMP THREADPRIVATE(boussinesq)
   !$OMP THREADPRIVATE(hydrostatic) 
+  LOGICAL :: dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro
+  !$OMP THREADPRIVATE(dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro)
 
 CONTAINS
Index: codes/icosagcm/trunk/src/dynamics/caldyn.f90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn.f90	(revision 579)
+++ codes/icosagcm/trunk/src/dynamics/caldyn.f90	(revision 580)
@@ -2,8 +2,9 @@
   USE icosa
   PRIVATE
+  SAVE
   CHARACTER(LEN=255),SAVE :: caldyn_type
 !$OMP THREADPRIVATE(caldyn_type)
-  
-  PUBLIC init_caldyn, caldyn, caldyn_BC
+ 
+  PUBLIC init_caldyn, caldyn, caldyn_BC, dysl
   
 CONTAINS
@@ -28,5 +29,4 @@
     END SELECT
         
-      
   END SUBROUTINE init_caldyn
 
Index: codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90	(revision 579)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90	(revision 580)
@@ -22,11 +22,20 @@
     REAL(rstd),POINTER :: planetvel(:)
 
-#ifdef CPP_DYSL
-       IF(is_master) PRINT *,'CPP_DYSL : Using macro-generated compute kernels'
-#endif
-
     hydrostatic=.TRUE.
     CALL getin("hydrostatic",hydrostatic)
   
+    dysl=.NOT.hydrostatic ! dysl defaults to .FALSE. if hydrostatic, .TRUE. if NH                                              
+    CALL getin("dysl",dysl)
+    dysl_geopot=dysl
+    CALL getin("dysl_geopot",dysl_geopot)
+    dysl_caldyn_fast=dysl
+    CALL getin("dysl_caldyn_fast",dysl_caldyn_fast)
+    dysl_slow_hydro=dysl
+    CALL getin("dysl_slow_hydro",dysl_slow_hydro)
+    dysl_pvort_only=dysl
+    CALL getin("dysl_pvort_only",dysl_pvort_only)
+    dysl_caldyn_coriolis=dysl
+    CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis)
+
     def='energy'
     CALL getin('caldyn_conserv',def)
Index: codes/icosagcm/trunk/src/dynamics/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_hevi.f90	(revision 579)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_hevi.f90	(revision 580)
@@ -48,5 +48,5 @@
     TYPE(t_field) :: f_dPhi_fast(:)
     
-    REAL(rstd),POINTER :: ps(:), dps(:)
+    REAL(rstd),POINTER :: ps(:), dps(:), phis(:)
     REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
     REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
@@ -115,4 +115,5 @@
           CALL compute_geopot(mass,theta, ps,pk,geopot)
        ELSE
+          phis = f_phis(ind)
           W = f_W(ind)
           dW = f_dW_fast(ind)
@@ -121,5 +122,5 @@
           m_il = f_wil(ind)
           pres = f_gradPhi2(ind)
-          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
+          CALL compute_caldyn_solver(tau,phis, mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
        END IF
        u=f_u(ind)
Index: codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90	(revision 579)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90	(revision 580)
@@ -49,9 +49,7 @@
     Rd = kappa*cpp
 
-#ifdef CPP_DYSL
-!#if 0
+    IF(dysl_geopot) THEN
 #include "../kernels/compute_geopot.k90"
-#else
-
+    ELSE
     ! Pressure is computed first top-down (temporarily stored in pk)
     ! Then Exner pressure and geopotential are computed bottom-up
@@ -164,5 +162,5 @@
     END IF
 
-#endif
+    END IF ! dysl
 
     !ym flush geopot
@@ -314,9 +312,9 @@
     CALL trace_start("compute_caldyn_vert_nh")
 
-#ifdef CPP_DYSL
+    IF(dysl) THEN
 !$OMP BARRIER
 #include "../kernels/caldyn_vert_NH.k90"
 !$OMP BARRIER
-#else
+    ELSE
 #define ETA_DOT(ij) eta_dot(ij,1)
 #define WCOV(ij) wcov(ij,1)
@@ -365,6 +363,6 @@
 #undef ETA_DOT
 #undef WCOV
-#endif
-
+
+    END IF ! dysl
     CALL trace_end("compute_caldyn_vert_nh")
 
Index: codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90	(revision 579)
+++ codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90	(revision 580)
@@ -9,8 +9,7 @@
   PRIVATE
 
-  REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME
+  REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6
 
   LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
-  LOGICAL, PARAMETER :: rigid=.TRUE.
 
   PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, &
@@ -62,7 +61,8 @@
     CALL trace_start("compute_pvort_only")  
 !!! Compute shallow-water potential vorticity
-#ifdef CPP_DYSL
+    IF(dysl_pvort_only) THEN
 #include "../kernels/pvort_only.k90"
-#else
+    ELSE
+
     radius_m2=radius**(-2)
     DO l = ll_begin,ll_end
@@ -94,11 +94,13 @@
 
     ENDDO
-#endif
+   
+    END IF ! dysl
     CALL trace_end("compute_pvort_only")
 
   END SUBROUTINE compute_pvort_only
 
-  SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il)
+  SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il)
     REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs
+    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
     REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm)
     REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1)
@@ -124,7 +126,9 @@
     ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
 
-#ifdef CPP_DYSL
+    IF(dysl) THEN
+#define PHI_BOT(ij) phis(ij)
 #include "../kernels/compute_NH_geopot.k90"
-#else
+#undef PHI_BOT
+    ELSE
 !    FIXME : vertical OpenMP parallelism will not work
    
@@ -164,5 +168,5 @@
           B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot
           R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  &
-               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-Phi_bot) )
+               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-phis(ij)) )
        ENDDO
        ! inner interfaces
@@ -260,10 +264,11 @@
     END DO ! Newton-Raphson
 
-#endif
+    END IF ! dysl
     
   END SUBROUTINE compute_NH_geopot
 
-  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
+  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)
@@ -278,4 +283,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
     INTEGER    :: ij, l
@@ -285,10 +291,17 @@
     Rd=cpp*kappa
 
-#ifdef CPP_DYSL
+    IF(dysl) THEN
+
 !$OMP BARRIER
+#define PHI_BOT(ij) phis(ij)
+#define PHI_BOT_VAR phis
 #include "../kernels/caldyn_solver.k90"
+#undef PHI_BOT_VAR
+#undef PHI_BOT
 !$OMP BARRIER
-#else
-#define BERNI(ij) berni(ij,1)
+
+    ELSE
+
+#define BERNI(ij) berni1(ij)
     ! FIXME : vertical OpenMP parallelism will not work
 
@@ -310,5 +323,5 @@
 
     IF(tau>0) THEN ! solve implicit problem for geopotential
-       CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot)
+       CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot)
     END IF
 
@@ -372,5 +385,6 @@
     ENDDO
 #undef BERNI
-#endif
+
+    END IF ! dysl
 
     CALL trace_end("compute_caldyn_solver")
@@ -396,7 +410,8 @@
     Rd=cpp*kappa
 
-#ifdef CPP_DYSL
+    IF(dysl_caldyn_fast) THEN
 #include "../kernels/caldyn_fast.k90"
-#else
+    ELSE
+
     ! Compute Bernoulli term
     IF(boussinesq) THEN
@@ -496,5 +511,6 @@
        END IF
     END DO
-#endif       
+
+    END IF ! dysl
     CALL trace_end("compute_caldyn_fast")
 
@@ -515,7 +531,9 @@
     CALL trace_start("compute_caldyn_Coriolis")
 
-#ifdef CPP_DYSL
+    IF(dysl_caldyn_coriolis) THEN
+
 #include "../kernels/coriolis.k90"
-#else
+
+    ELSE
 #define FTHETA(ij) Ftheta(ij,1)
 
@@ -653,5 +671,6 @@
     END SELECT
 #undef FTHETA
-#endif
+
+    END IF ! dysl
 
     CALL trace_end("compute_caldyn_Coriolis")
@@ -667,4 +686,5 @@
     
     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
     INTEGER :: ij,l
@@ -672,10 +692,13 @@
     CALL trace_start("compute_caldyn_slow_hydro")
 
-#ifdef CPP_DYSL
+    IF(dysl_slow_hydro) THEN
+
 #define BERNI(ij,l) berni(ij,l)
 #include "../kernels/caldyn_slow_hydro.k90"
 #undef BERNI
-#else
-#define BERNI(ij) berni(ij,1)
+
+     ELSE
+
+#define BERNI(ij) berni1(ij)
 
     DO l = ll_begin, ll_end
@@ -725,6 +748,8 @@
        END IF
     END DO
+
 #undef BERNI
-#endif
+
+    END IF ! dysl
     CALL trace_end("compute_caldyn_slow_hydro")    
   END SUBROUTINE compute_caldyn_slow_hydro
@@ -749,21 +774,26 @@
     REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu
 
-#ifdef CPP_DYSL
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
     REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W
     REAL(rstd) :: v_el(3*iim*jjm,llm+1)
-#else
-    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function
-    REAL(rstd) :: G_el(3*iim*jjm) ! horizontal flux of W
-    REAL(rstd) :: v_el(3*iim*jjm)
-#endif
+
+    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
+    REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W
+    REAL(rstd) :: v_el1(3*iim*jjm)
 
     CALL trace_start("compute_caldyn_slow_NH")
 
-#ifdef CPP_DYSL
+    IF(dysl) THEN
+
 !$OMP BARRIER
 #include "../kernels/caldyn_slow_NH.k90"
 !$OMP BARRIER
-#else
+     
+     ELSE
+
+#define BERNI(ij) berni1(ij)
+#define G_EL(ij) G_el1(ij)
+#define V_EL(ij) v_el1(ij)
+
     DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
        IF(l==1) THEN
@@ -792,6 +822,6 @@
           W2_el = .5*le_de(ij+u_right) * &
                ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
-          v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
-          G_el(ij+u_right) = v_el(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
+          V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
+          G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
           ! Compute on edge 'lup'
           W_el  = .5*( W(ij,l)+W(ij+t_lup,l) )
@@ -800,6 +830,6 @@
           W2_el = .5*le_de(ij+u_lup) * &
                ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
-          v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
-          G_el(ij+u_lup) = v_el(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
+          V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
+          G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
           ! Compute on edge 'ldown'
           W_el  = .5*( W(ij,l)+W(ij+t_ldown,l) )
@@ -808,6 +838,6 @@
           W2_el = .5*le_de(ij+u_ldown) * &
                ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
-          v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
-          G_el(ij+u_ldown) = v_el(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
+          V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
+          G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
        END DO
        ! compute GradPhi2, dPhi, dW
@@ -820,21 +850,20 @@
                le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  &
                le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 )
-!          gradPhi2(ij,l) = 0. ! FIXME !!
 
           dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & 
-               ( DePhil(ij+u_right,l)*v_el(ij+u_right) + & ! -v.gradPhi, 
-                 DePhil(ij+u_rup,l)*v_el(ij+u_rup) +     & ! v_el already has le_de
-                 DePhil(ij+u_lup,l)*v_el(ij+u_lup) +     &
-                 DePhil(ij+u_left,l)*v_el(ij+u_left) +   &
-                 DePhil(ij+u_ldown,l)*v_el(ij+u_ldown) + &
-                 DePhil(ij+u_rdown,l)*v_el(ij+u_rdown) )
+               ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi, 
+                 DePhil(ij+u_rup,l)*V_EL(ij+u_rup) +     & ! v_el already has le_de
+                 DePhil(ij+u_lup,l)*V_EL(ij+u_lup) +     &
+                 DePhil(ij+u_left,l)*V_EL(ij+u_left) +   &
+                 DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + &
+                 DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) )
 
           dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el), 
-               ne_right*G_el(ij+u_right) +  & ! G_el already has le_de
-               ne_rup*G_el(ij+u_rup) +      &
-               ne_lup*G_el(ij+u_lup) +      &  
-               ne_left*G_el(ij+u_left) +    &
-               ne_ldown*G_el(ij+u_ldown) +  &
-               ne_rdown*G_el(ij+u_rdown))
+               ne_right*G_EL(ij+u_right) +  & ! G_el already has le_de
+               ne_rup*G_EL(ij+u_rup) +      &
+               ne_lup*G_EL(ij+u_lup) +      &  
+               ne_left*G_EL(ij+u_left) +    &
+               ne_ldown*G_EL(ij+u_ldown) +  &
+               ne_rdown*G_EL(ij+u_rdown))
        END DO
     END DO
@@ -843,5 +872,5 @@
        ! Compute berni at scalar points
        DO ij=ij_begin_ext, ij_end_ext
-          berni(ij) = &
+          BERNI(ij) = &
                1/(4*Ai(ij))*( &
                     le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
@@ -860,18 +889,23 @@
                     -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))
           hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)
-          du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
+          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
           ! Compute on edge 'lup'
           uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) &
                   -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))
           hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)
-          du(ij+u_lup,l) = ne_lup*(berni(ij)-berni(ij+t_lup))
+          du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
           ! Compute on edge 'ldown'
           uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) &
                     -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))
           hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)
-          du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
+          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
        END DO
     END DO
-#endif
+
+#undef V_EL
+#undef G_EL
+#undef BERNI
+
+    END IF ! dysl
 
     CALL trace_end("compute_caldyn_slow_NH")
Index: codes/icosagcm/trunk/src/kernels/caldyn_fast.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/caldyn_fast.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/caldyn_fast.k90	(revision 580)
@@ -0,0 +1,50 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_fast ----------------------------------
+   !
+   SELECT CASE(caldyn_thermo)
+   CASE(thermo_boussinesq)
+      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
+   CASE(thermo_theta)
+      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))
+         END DO
+      END DO
+   CASE(thermo_entropy)
+      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))
+            berni(ij,l) = berni(ij,l) + pk(ij,l)*(cpp-theta(ij,l,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
+         END DO
+      END DO
+   CASE DEFAULT
+      PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo ! FIXME
+      STOP
+   END SELECT
+   !
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         due = .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)
+         du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due
+         u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+         due = .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)
+         du(ij+u_lup,l) = du(ij+u_lup,l) - ne_lup*due
+         u(ij+u_lup,l) = u(ij+u_lup,l) + tau*du(ij+u_lup,l)
+         due = .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_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due
+         u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+      END DO
+   END DO
+   !
+   !---------------------------- caldyn_fast ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/caldyn_slow_NH.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/caldyn_slow_NH.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/caldyn_slow_NH.k90	(revision 580)
@@ -0,0 +1,181 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_slow_NH ----------------------------------
+   IF (ll_begin==1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         w_il(ij,1) = 2.*W(ij,1)/rhodz(ij,1)
+      END DO
+   END IF
+   DO l = ll_beginp1, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,l-1)+rhodz(ij,l))
+      END DO
+   END DO
+   IF(ll_endp1==llm+1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         w_il(ij,llm+1) = 2.*W(ij,llm+1)/rhodz(ij,llm)
+      END DO
+   END IF
+   IF (ll_begin==1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ! compute DePhi, v_el, G_el, F_el
+         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+         ! while DePhil, W_el and F_el do not
+         W_el = .5*( W(ij+t_right,1)+W(ij,1) )
+         DePhil(ij+u_right,1) = ne_right*(Phi(ij+t_right,1)-Phi(ij,1))
+         F_el(ij+u_right,1) = DePhil(ij+u_right,1)*W_el
+         W2_el = .5*le_de(ij+u_right) * ( W(ij,1)*w_il(ij,1) + W(ij+t_right,1)*w_il(ij+t_right,1) )
+         v_el(ij+u_right,1) = .5*le_de(ij+u_right)*(u(ij+u_right,1)+u(ij+u_right,1)) ! checked
+         G_el(ij+u_right,1) = v_el(ij+u_right,1)*W_el - DePhil(ij+u_right,1)*W2_el
+         ! compute DePhi, v_el, G_el, F_el
+         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+         ! while DePhil, W_el and F_el do not
+         W_el = .5*( W(ij+t_lup,1)+W(ij,1) )
+         DePhil(ij+u_lup,1) = ne_lup*(Phi(ij+t_lup,1)-Phi(ij,1))
+         F_el(ij+u_lup,1) = DePhil(ij+u_lup,1)*W_el
+         W2_el = .5*le_de(ij+u_lup) * ( W(ij,1)*w_il(ij,1) + W(ij+t_lup,1)*w_il(ij+t_lup,1) )
+         v_el(ij+u_lup,1) = .5*le_de(ij+u_lup)*(u(ij+u_lup,1)+u(ij+u_lup,1)) ! checked
+         G_el(ij+u_lup,1) = v_el(ij+u_lup,1)*W_el - DePhil(ij+u_lup,1)*W2_el
+         ! compute DePhi, v_el, G_el, F_el
+         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+         ! while DePhil, W_el and F_el do not
+         W_el = .5*( W(ij+t_ldown,1)+W(ij,1) )
+         DePhil(ij+u_ldown,1) = ne_ldown*(Phi(ij+t_ldown,1)-Phi(ij,1))
+         F_el(ij+u_ldown,1) = DePhil(ij+u_ldown,1)*W_el
+         W2_el = .5*le_de(ij+u_ldown) * ( W(ij,1)*w_il(ij,1) + W(ij+t_ldown,1)*w_il(ij+t_ldown,1) )
+         v_el(ij+u_ldown,1) = .5*le_de(ij+u_ldown)*(u(ij+u_ldown,1)+u(ij+u_ldown,1)) ! checked
+         G_el(ij+u_ldown,1) = v_el(ij+u_ldown,1)*W_el - DePhil(ij+u_ldown,1)*W2_el
+      END DO
+   END IF
+   DO l = ll_beginp1, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ! compute DePhi, v_el, G_el, F_el
+         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+         ! while DePhil, W_el and F_el do not
+         W_el = .5*( W(ij+t_right,l)+W(ij,l) )
+         DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
+         F_el(ij+u_right,l) = DePhil(ij+u_right,l)*W_el
+         W2_el = .5*le_de(ij+u_right) * ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
+         v_el(ij+u_right,l) = .5*le_de(ij+u_right)*(u(ij+u_right,l)+u(ij+u_right,l-1)) ! checked
+         G_el(ij+u_right,l) = v_el(ij+u_right,l)*W_el - DePhil(ij+u_right,l)*W2_el
+         ! compute DePhi, v_el, G_el, F_el
+         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+         ! while DePhil, W_el and F_el do not
+         W_el = .5*( W(ij+t_lup,l)+W(ij,l) )
+         DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
+         F_el(ij+u_lup,l) = DePhil(ij+u_lup,l)*W_el
+         W2_el = .5*le_de(ij+u_lup) * ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
+         v_el(ij+u_lup,l) = .5*le_de(ij+u_lup)*(u(ij+u_lup,l)+u(ij+u_lup,l-1)) ! checked
+         G_el(ij+u_lup,l) = v_el(ij+u_lup,l)*W_el - DePhil(ij+u_lup,l)*W2_el
+         ! compute DePhi, v_el, G_el, F_el
+         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+         ! while DePhil, W_el and F_el do not
+         W_el = .5*( W(ij+t_ldown,l)+W(ij,l) )
+         DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
+         F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el
+         W2_el = .5*le_de(ij+u_ldown) * ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
+         v_el(ij+u_ldown,l) = .5*le_de(ij+u_ldown)*(u(ij+u_ldown,l)+u(ij+u_ldown,l-1)) ! checked
+         G_el(ij+u_ldown,l) = v_el(ij+u_ldown,l)*W_el - DePhil(ij+u_ldown,l)*W2_el
+      END DO
+   END DO
+   IF(ll_endp1==llm+1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ! compute DePhi, v_el, G_el, F_el
+         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+         ! while DePhil, W_el and F_el do not
+         W_el = .5*( W(ij+t_right,llm+1)+W(ij,llm+1) )
+         DePhil(ij+u_right,llm+1) = ne_right*(Phi(ij+t_right,llm+1)-Phi(ij,llm+1))
+         F_el(ij+u_right,llm+1) = DePhil(ij+u_right,llm+1)*W_el
+         W2_el = .5*le_de(ij+u_right) * ( W(ij,llm+1)*w_il(ij,llm+1) + W(ij+t_right,llm+1)*w_il(ij+t_right,llm+1) )
+         v_el(ij+u_right,llm+1) = .5*le_de(ij+u_right)*(u(ij+u_right,llm)+u(ij+u_right,llm)) ! checked
+         G_el(ij+u_right,llm+1) = v_el(ij+u_right,llm+1)*W_el - DePhil(ij+u_right,llm+1)*W2_el
+         ! compute DePhi, v_el, G_el, F_el
+         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+         ! while DePhil, W_el and F_el do not
+         W_el = .5*( W(ij+t_lup,llm+1)+W(ij,llm+1) )
+         DePhil(ij+u_lup,llm+1) = ne_lup*(Phi(ij+t_lup,llm+1)-Phi(ij,llm+1))
+         F_el(ij+u_lup,llm+1) = DePhil(ij+u_lup,llm+1)*W_el
+         W2_el = .5*le_de(ij+u_lup) * ( W(ij,llm+1)*w_il(ij,llm+1) + W(ij+t_lup,llm+1)*w_il(ij+t_lup,llm+1) )
+         v_el(ij+u_lup,llm+1) = .5*le_de(ij+u_lup)*(u(ij+u_lup,llm)+u(ij+u_lup,llm)) ! checked
+         G_el(ij+u_lup,llm+1) = v_el(ij+u_lup,llm+1)*W_el - DePhil(ij+u_lup,llm+1)*W2_el
+         ! compute DePhi, v_el, G_el, F_el
+         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+         ! while DePhil, W_el and F_el do not
+         W_el = .5*( W(ij+t_ldown,llm+1)+W(ij,llm+1) )
+         DePhil(ij+u_ldown,llm+1) = ne_ldown*(Phi(ij+t_ldown,llm+1)-Phi(ij,llm+1))
+         F_el(ij+u_ldown,llm+1) = DePhil(ij+u_ldown,llm+1)*W_el
+         W2_el = .5*le_de(ij+u_ldown) * ( W(ij,llm+1)*w_il(ij,llm+1) + W(ij+t_ldown,llm+1)*w_il(ij+t_ldown,llm+1) )
+         v_el(ij+u_ldown,llm+1) = .5*le_de(ij+u_ldown)*(u(ij+u_ldown,llm)+u(ij+u_ldown,llm)) ! checked
+         G_el(ij+u_ldown,llm+1) = v_el(ij+u_ldown,llm+1)*W_el - DePhil(ij+u_ldown,llm+1)*W2_el
+      END DO
+   END IF
+   DO l = ll_begin, ll_endp1
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ! compute GradPhi2, dPhi, dW
+         gPhi2=0.
+         dP=0.
+         divG=0
+         gPhi2 = gPhi2 + le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2
+         dP = dP + le_de(ij+u_rup)*DePhil(ij+u_rup,l)*v_el(ij+u_rup,l)
+         divG = divG + ne_rup*G_el(ij+u_rup,l) ! -div(G_el), G_el already has le_de
+         gPhi2 = gPhi2 + le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2
+         dP = dP + le_de(ij+u_lup)*DePhil(ij+u_lup,l)*v_el(ij+u_lup,l)
+         divG = divG + ne_lup*G_el(ij+u_lup,l) ! -div(G_el), G_el already has le_de
+         gPhi2 = gPhi2 + le_de(ij+u_left)*DePhil(ij+u_left,l)**2
+         dP = dP + le_de(ij+u_left)*DePhil(ij+u_left,l)*v_el(ij+u_left,l)
+         divG = divG + ne_left*G_el(ij+u_left,l) ! -div(G_el), G_el already has le_de
+         gPhi2 = gPhi2 + le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2
+         dP = dP + le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)*v_el(ij+u_ldown,l)
+         divG = divG + ne_ldown*G_el(ij+u_ldown,l) ! -div(G_el), G_el already has le_de
+         gPhi2 = gPhi2 + le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2
+         dP = dP + le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)*v_el(ij+u_rdown,l)
+         divG = divG + ne_rdown*G_el(ij+u_rdown,l) ! -div(G_el), G_el already has le_de
+         gPhi2 = gPhi2 + le_de(ij+u_right)*DePhil(ij+u_right,l)**2
+         dP = dP + le_de(ij+u_right)*DePhil(ij+u_right,l)*v_el(ij+u_right,l)
+         divG = divG + ne_right*G_el(ij+u_right,l) ! -div(G_el), G_el already has le_de
+         gradPhi2(ij,l) = 1./(2.*Ai(ij)) * gPhi2
+         dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) - 1./(2.*Ai(ij))*dP
+         dW(ij,l) = (-1./Ai(ij))*divG
+      END DO
+   END DO
+   ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
+   !$OMP BARRIER
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ! Compute berni at scalar points
+         u2=0.
+         u2 = u2 + le_de(ij+u_rup)*u(ij+u_rup,l)**2
+         u2 = u2 + le_de(ij+u_lup)*u(ij+u_lup,l)**2
+         u2 = u2 + le_de(ij+u_left)*u(ij+u_left,l)**2
+         u2 = u2 + le_de(ij+u_ldown)*u(ij+u_ldown,l)**2
+         u2 = u2 + le_de(ij+u_rdown)*u(ij+u_rdown,l)**2
+         u2 = u2 + le_de(ij+u_right)*u(ij+u_right,l)**2
+         berni(ij,l) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(ij,l)*w_il(ij,l)**2 + gradPhi2(ij,l+1)*w_il(ij,l+1)**2 )
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ! Compute mass flux and grad(berni)
+         uu = .5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) - .5*( F_el(ij+u_right,l)+F_el(ij+u_right,l+1) )
+         hflux(ij+u_right,l) = le_de(ij+u_right)*uu
+         du(ij+u_right,l) = ne_right*(berni(ij,l)-berni(ij+t_right,l))
+         ! Compute mass flux and grad(berni)
+         uu = .5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) - .5*( F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1) )
+         hflux(ij+u_lup,l) = le_de(ij+u_lup)*uu
+         du(ij+u_lup,l) = ne_lup*(berni(ij,l)-berni(ij+t_lup,l))
+         ! Compute mass flux and grad(berni)
+         uu = .5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) - .5*( F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1) )
+         hflux(ij+u_ldown,l) = le_de(ij+u_ldown)*uu
+         du(ij+u_ldown,l) = ne_ldown*(berni(ij,l)-berni(ij+t_ldown,l))
+      END DO
+   END DO
+   !---------------------------- caldyn_slow_NH ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/caldyn_slow_hydro.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/caldyn_slow_hydro.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/caldyn_slow_hydro.k90	(revision 580)
@@ -0,0 +1,47 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_slow_hydro ----------------------------------
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         uu = .5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
+         hflux(ij+u_right,l) = uu*le_de(ij+u_right)
+         uu = .5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
+         hflux(ij+u_lup,l) = uu*le_de(ij+u_lup)
+         uu = .5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
+         hflux(ij+u_ldown,l) = uu*le_de(ij+u_ldown)
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         ke=0.d0
+         ke = ke + le_de(ij+u_rup)*u(ij+u_rup,l)**2
+         ke = ke + le_de(ij+u_lup)*u(ij+u_lup,l)**2
+         ke = ke + le_de(ij+u_left)*u(ij+u_left,l)**2
+         ke = ke + le_de(ij+u_ldown)*u(ij+u_ldown,l)**2
+         ke = ke + le_de(ij+u_rdown)*u(ij+u_rdown,l)**2
+         ke = ke + le_de(ij+u_right)*u(ij+u_right,l)**2
+         BERNI(ij,l)=ke*(.25/Ai(ij))
+      END DO
+   END DO
+   IF(zero) THEN
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            du(ij+u_right,l) = ne_right*(berni(ij,l)-berni(ij+t_right,l)) ! minus gradient
+            du(ij+u_lup,l) = ne_lup*(berni(ij,l)-berni(ij+t_lup,l)) ! minus gradient
+            du(ij+u_ldown,l) = ne_ldown*(berni(ij,l)-berni(ij+t_ldown,l)) ! minus gradient
+         END DO
+      END DO
+   ELSE
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            du(ij+u_right,l) = du(ij+u_right,l) + ne_right*(berni(ij,l)-berni(ij+t_right,l)) ! minus gradient
+            du(ij+u_lup,l) = du(ij+u_lup,l) + ne_lup*(berni(ij,l)-berni(ij+t_lup,l)) ! minus gradient
+            du(ij+u_ldown,l) = du(ij+u_ldown,l) + ne_ldown*(berni(ij,l)-berni(ij+t_ldown,l)) ! minus gradient
+         END DO
+      END DO
+   END IF
+   !---------------------------- caldyn_slow_hydro ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/caldyn_solver.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/caldyn_solver.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/caldyn_solver.k90	(revision 580)
@@ -0,0 +1,103 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- 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)
+   ! kappa = R/Cp
+   ! 1-kappa = Cv/Cp
+   ! Cp/Cv = 1/(1-kappa)
+   gamma = 1./(1.-kappa)
+   vreff = Rd*Treff/preff ! reference specific volume
+   Cvd = cpp-Rd
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         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
+            ! kappa.theta.rho = p/exner
+            ! => X = (p/p0)/(exner/Cp)
+            ! = (p/p0)^(1-kappa)
+            pres(ij,l) = preff*(X_ij**gamma) ! pressure
+            ! Compute Exner function (needed by compute_caldyn_fast)
+            ! other formulae possible if exponentiation is slow
+            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 )
+            pres(ij,l) = rho_ij*Rd*T_ij
+            pk(ij,l) = T_ij
+         CASE DEFAULT
+            STOP
+         END SELECT
+      END DO
+   END DO
+   ! We need a barrier here because we compute pres above and do a vertical difference below
+   !$OMP BARRIER
+   IF (ll_begin==1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ! Lower BC
+         dW(ij,1) = (1./g)*(pbot-rho_bot*(geopot(ij,1)-PHI_BOT(ij))-pres(ij,1)) - m_il(ij,1)
+         W(ij,1) = W(ij,1)+tau*dW(ij,1) ! update W
+         dPhi(ij,1) = g*g*W(ij,1)/m_il(ij,1)
+      END DO
+   END IF
+   DO l = ll_beginp1, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         dW(ij,l) = (1./g)*(pres(ij,l-1)-pres(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)
+      END DO
+   END DO
+   IF(ll_endp1==llm+1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ! Top BC
+         dW(ij,llm+1) = (1./g)*(pres(ij,llm+1 -1)-ptop) - m_il(ij,llm+1)
+         W(ij,llm+1) = W(ij,llm+1)+tau*dW(ij,llm+1) ! update W
+         dPhi(ij,llm+1) = g*g*W(ij,llm+1)/m_il(ij,llm+1)
+      END DO
+   END IF
+   ! We need a barrier here because we update W above and do a vertical average below
+   !$OMP BARRIER
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ! compute du = -0.5*g^2.grad(w^2)
+         berni(ij,l) = (-.25*g*g)*((W(ij,l)/m_il(ij,l))**2 + (W(ij,l+1)/m_il(ij,l+1))**2 )
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         du(ij+u_right,l) = ne_right*(berni(ij,l)-berni(ij+t_right,l))
+         du(ij+u_lup,l) = ne_lup*(berni(ij,l)-berni(ij+t_lup,l))
+         du(ij+u_ldown,l) = ne_ldown*(berni(ij,l)-berni(ij+t_ldown,l))
+      END DO
+   END DO
+   !---------------------------- caldyn_solver ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/caldyn_vert_NH.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/caldyn_vert_NH.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/caldyn_vert_NH.k90	(revision 580)
@@ -0,0 +1,71 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_vert_NH ----------------------------------
+   IF (ll_begin==1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         w_ij = ( W(ij,1)+.5*W(ij,1 +1) )/mass(ij,1)
+         wflux_ij = .5*(wflux(ij,1)+wflux(ij,1 +1))
+         W_etadot(ij,1) = wflux_ij*w_ij
+         eta_dot(ij,1) = wflux_ij / mass(ij,1)
+         wcov(ij,1) = w_ij*(geopot(ij,1 +1)-geopot(ij,1))
+      END DO
+   END IF
+   DO l = ll_beginp1, ll_endm1
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         w_ij = .5*( W(ij,l)+W(ij,l+1) )/mass(ij,l)
+         wflux_ij = .5*(wflux(ij,l)+wflux(ij,l+1))
+         W_etadot(ij,l) = wflux_ij*w_ij
+         eta_dot(ij,l) = wflux_ij / mass(ij,l)
+         wcov(ij,l) = w_ij*(geopot(ij,l+1)-geopot(ij,l))
+      END DO
+   END DO
+   IF(ll_end==llm) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         w_ij = ( .5*W(ij,llm)+W(ij,llm+1) )/mass(ij,llm)
+         wflux_ij = .5*(wflux(ij,llm)+wflux(ij,llm+1))
+         W_etadot(ij,llm) = wflux_ij*w_ij
+         eta_dot(ij,llm) = wflux_ij / mass(ij,llm)
+         wcov(ij,llm) = w_ij*(geopot(ij,llm+1)-geopot(ij,llm))
+      END DO
+   END IF
+   ! add NH term to du
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         du(ij+u_right,l) = du(ij+u_right,l) - .5*(wcov(ij,l)+wcov(ij+t_right,l))*ne_right*(eta_dot(ij+t_right,l)-eta_dot(ij,l))
+         du(ij+u_lup,l) = du(ij+u_lup,l) - .5*(wcov(ij,l)+wcov(ij+t_lup,l))*ne_lup*(eta_dot(ij+t_lup,l)-eta_dot(ij,l))
+         du(ij+u_ldown,l) = du(ij+u_ldown,l) - .5*(wcov(ij,l)+wcov(ij+t_ldown,l))*ne_ldown*(eta_dot(ij+t_ldown,l)-eta_dot(ij,l))
+      END DO
+   END DO
+   ! add NH terms to dW, dPhi
+   ! FIXME : TODO top and bottom
+   DO l = ll_beginp1, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         dPhi(ij,l)=dPhi(ij,l)-wflux(ij,l)*(geopot(ij,l+1)-geopot(ij,l-1))/(mass(ij,l-1)+mass(ij,l))
+      END DO
+   END DO
+   ! We need a barrier here because we compute W_etadot above and do a vertical difference below
+   !$OMP BARRIER
+   IF (ll_begin==1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         dW(ij,1) = dW(ij,1) - W_etadot(ij,1)
+      END DO
+   END IF
+   DO l = ll_beginp1, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         dW(ij,l) = dW(ij,l) + W_etadot(ij,l-1) - W_etadot(ij,l)
+      END DO
+   END DO
+   IF(ll_endp1==llm+1) THEN
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         dW(ij,llm+1) = dW(ij,llm+1) + W_etadot(ij,llm+1 -1)
+      END DO
+   END IF
+   !---------------------------- caldyn_vert_NH ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/compute_NH_geopot.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/compute_NH_geopot.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/compute_NH_geopot.k90	(revision 580)
@@ -0,0 +1,119 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- compute_NH_geopot ----------------------------------
+   tau2_g=tau*tau/g
+   g2=g*g
+   gm2 = 1./g2
+   gamma = 1./(1.-kappa)
+   !$OMP BARRIER
+   ! compute Phi_star
+   DO l = 1,llm+1
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         Phi_star_il(ij,l) = Phi_il(ij,l) + tau*g2*(W_il(ij,l)/m_il(ij,l)-tau)
+      END DO
+   END DO
+   ! Newton-Raphson iteration : Phi_il contains current guess value
+   DO iter=1,2 ! 2 iterations should be enough
+      ! Compute pressure, A_ik
+      SELECT CASE(caldyn_thermo)
+      CASE(thermo_theta)
+         DO l = 1,llm
+            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 = (cpp/preff)*kappa*theta(ij,l)*rho_ij
+               p_ik(ij,l) = preff*(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
+            END DO
+         END DO
+      CASE(thermo_entropy)
+         DO l = 1,llm
+            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)
+               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
+            END DO
+         END DO
+      CASE DEFAULT
+         PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo
+         STOP
+      END SELECT
+      ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom
+      ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
+      ! Compute residual R_il and B_il
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         ! bottom interface l=1
+         ml_g2 = gm2*m_il(ij,1)
+         B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot
+         R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1)) &
+         + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-PHI_BOT(ij)) )
+      END DO
+      DO l = 2,llm
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            ! inner interfaces
+            ml_g2 = gm2*m_il(ij,l)
+            B_il(ij,l) = A_ik(ij,l)+A_ik(ij,l-1) + ml_g2
+            R_il(ij,l) = ml_g2*( Phi_il(ij,l)-Phi_star_il(ij,l)) &
+            + tau2_g*(p_ik(ij,l)-p_ik(ij,l-1))
+            ! consistency check : if Wil=0 and initial state is in hydrostatic balance
+            ! then Phi_star_il(ij,l) = Phi_il(ij,l) - tau^2*g^2
+            ! and residual = tau^2*(ml+(1/g)dl_pi)=0
+         END DO
+      END DO
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         ! top interface l=llm+1
+         ml_g2 = gm2*m_il(ij,llm+1)
+         B_il(ij,llm+1) = A_ik(ij,llm+1 -1) + ml_g2
+         R_il(ij,llm+1) = ml_g2*( Phi_il(ij,llm+1)-Phi_star_il(ij,llm+1)) &
+         + tau2_g*( ptop-p_ik(ij,llm+1 -1) )
+      END DO
+      !
+      ! Forward sweep :
+      ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
+      ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         X_ij = 1./B_il(ij,1)
+         C_ik(ij,1) = -A_ik(ij,1) * X_ij
+         D_il(ij,1) = R_il(ij,1) * X_ij
+      END DO
+      DO l = 2,llm
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            X_ij = 1./( B_il(ij,l) + A_ik(ij,l-1)*C_ik(ij,l-1) )
+            C_ik(ij,l) = -A_ik(ij,l) * X_ij
+            D_il(ij,l) = (R_il(ij,l)+A_ik(ij,l-1)*D_il(ij,l-1)) * X_ij
+         END DO
+      END DO
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         X_ij = 1./( B_il(ij,llm+1) + A_ik(ij,llm+1 -1)*C_ik(ij,llm+1 -1) )
+         D_il(ij,llm+1) = (R_il(ij,llm+1)+A_ik(ij,llm+1 -1)*D_il(ij,llm+1 -1)) * X_ij
+         ! Back substitution :
+         ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
+         ! + Newton-Raphson update
+         ! top interface l=llm+1
+         x_il(ij,llm+1) = D_il(ij,llm+1)
+         Phi_il(ij,llm+1) = Phi_il(ij,llm+1) - x_il(ij,llm+1)
+      END DO
+      DO l = llm,1,-1
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            ! Back substitution at lower interfaces
+            x_il(ij,l) = D_il(ij,l) - C_ik(ij,l)*x_il(ij,l+1)
+            Phi_il(ij,l) = Phi_il(ij,l) - x_il(ij,l)
+         END DO
+      END DO
+      IF(debug_hevi_solver) THEN
+         PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
+         PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
+         DO l=1,llm+1
+            WRITE(*,'(A,I2.1,I3.2,E9.2)'), '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:)))
+         END DO
+         DO l=1,llm+1
+            WRITE(*,'(A,I2.1,I3.2,E9.2)'), '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:)))
+         END DO
+      END IF
+   END DO ! Newton-Raphson
+   !$OMP BARRIER
+   debug_hevi_solver=.FALSE.
+   !---------------------------- compute_NH_geopot ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/compute_geopot.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/compute_geopot.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/compute_geopot.k90	(revision 580)
@@ -0,0 +1,102 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- compute_geopot ----------------------------------
+   SELECT CASE(caldyn_thermo)
+   CASE(thermo_boussinesq)
+      ! use hydrostatic balance with theta*rhodz to find pk (=Lagrange multiplier=pressure)
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         pk(ij,llm) = ptop + .5*g* theta(ij,llm,1)*rhodz(ij,llm)
+      END DO
+      DO l = llm-1,1,-1
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            pk(ij,l) = pk(ij,l+1) + (.5*g)*( theta(ij,l,1)*rhodz(ij,l) + theta(ij,l+1,1)*rhodz(ij,l+1) )
+         END DO
+      END DO
+      IF(caldyn_eta == eta_lag) THEN
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            ps(ij) = pk(ij,1) + .5*g* theta(ij,1,1)*rhodz(ij,1)
+         END DO
+      END IF
+      ! now pk contains the Lagrange multiplier (pressure)
+      ! specific volume 1 = dphi/g/rhodz
+      DO l = 1,llm
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
+         END DO
+      END DO
+   CASE(thermo_theta)
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         pk(ij,llm) = ptop + .5*g* rhodz(ij,llm)
+      END DO
+      DO l = llm-1,1,-1
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            pk(ij,l) = pk(ij,l+1) + (.5*g)*( rhodz(ij,l) + rhodz(ij,l+1) )
+         END DO
+      END DO
+      IF(caldyn_eta == eta_lag) THEN
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            ps(ij) = pk(ij,1) + .5*g* rhodz(ij,1)
+         END DO
+      END IF
+      DO l = 1,llm
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            p_ik = pk(ij,l)
+            exner_ik = cpp * (p_ik/preff) ** kappa
+            gv = (g*kappa)*theta(ij,l,1)*exner_ik/p_ik
+            pk(ij,l) = exner_ik
+            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
+         END DO
+      END DO
+   CASE(thermo_entropy)
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         pk(ij,llm) = ptop + .5*g* rhodz(ij,llm)
+      END DO
+      DO l = llm-1,1,-1
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            pk(ij,l) = pk(ij,l+1) + (.5*g)*( rhodz(ij,l) + rhodz(ij,l+1) )
+         END DO
+      END DO
+      IF(caldyn_eta == eta_lag) THEN
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            ps(ij) = pk(ij,1) + .5*g* rhodz(ij,1)
+         END DO
+      END IF
+      DO l = 1,llm
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            p_ik = pk(ij,l)
+            temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp)
+            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
+            pk(ij,l) = temp_ik
+            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
+         END DO
+      END DO
+   CASE(thermo_moist)
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         pk(ij,llm) = ptop + .5*g* rhodz(ij,llm)*(1.+theta(ij,llm,2))
+      END DO
+      DO l = llm-1,1,-1
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            pk(ij,l) = pk(ij,l+1) + (.5*g)*( rhodz(ij,l)*(1.+theta(ij,l,2)) + rhodz(ij,l+1)*(1.+theta(ij,l+1,2)) )
+         END DO
+      END DO
+      IF(caldyn_eta == eta_lag) THEN
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            ps(ij) = pk(ij,1) + .5*g* rhodz(ij,1)*(1.+theta(ij,1,2))
+         END DO
+      END IF
+      DO l = 1,llm
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            p_ik = pk(ij,l)
+            qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md
+            Rmix = Rd+qv*Rv
+            chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
+            temp_ik = Treff*exp(chi)
+            ! specific volume v = R*T/p
+            ! R = (Rd + qv.Rv)/(1+qv)
+            gv = g*Rmix*temp_ik/(p_ik*(1+qv))
+            pk(ij,l) = temp_ik
+            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
+         END DO
+      END DO
+   END SELECT
+   !---------------------------- compute_geopot ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/coriolis.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/coriolis.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/coriolis.k90	(revision 580)
@@ -0,0 +1,84 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- coriolis ----------------------------------
+   !
+   DO iq=1,nqdyn
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin_ext, ij_end_ext
+            Ftheta(ij+u_right,l) = .5*(theta(ij,l,iq)+theta(ij+t_right,l,iq))*hflux(ij+u_right,l)
+            Ftheta(ij+u_lup,l) = .5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq))*hflux(ij+u_lup,l)
+            Ftheta(ij+u_ldown,l) = .5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq))*hflux(ij+u_ldown,l)
+         END DO
+      END DO
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            divF=0.
+            divF = divF + Ftheta(ij+u_rup,l)*ne_rup
+            divF = divF + Ftheta(ij+u_lup,l)*ne_lup
+            divF = divF + Ftheta(ij+u_left,l)*ne_left
+            divF = divF + Ftheta(ij+u_ldown,l)*ne_ldown
+            divF = divF + Ftheta(ij+u_rdown,l)*ne_rdown
+            divF = divF + Ftheta(ij+u_right,l)*ne_right
+            dtheta_rhodz(ij,l,iq) = -divF / Ai(ij)
+         END DO
+      END DO
+   END DO ! iq
+   !
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         divF=0.
+         divF = divF + hflux(ij+u_rup,l)*ne_rup
+         divF = divF + hflux(ij+u_lup,l)*ne_lup
+         divF = divF + hflux(ij+u_left,l)*ne_left
+         divF = divF + hflux(ij+u_ldown,l)*ne_ldown
+         divF = divF + hflux(ij+u_rdown,l)*ne_rdown
+         divF = divF + hflux(ij+u_right,l)*ne_right
+         convm(ij,l) = -divF / Ai(ij)
+      END DO
+   END DO
+   !
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         du_trisk=0.
+         du_trisk = du_trisk + wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))
+         du_trisk = du_trisk + wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))
+         du_trisk = du_trisk + wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))
+         du_trisk = du_trisk + wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))
+         du_trisk = du_trisk + wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))
+         du_trisk = du_trisk + wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))
+         du_trisk = du_trisk + wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))
+         du_trisk = du_trisk + wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))
+         du_trisk = du_trisk + wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))
+         du_trisk = du_trisk + wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
+         du(ij+u_right,l) = du(ij+u_right,l) + .5*du_trisk
+         du_trisk=0.
+         du_trisk = du_trisk + wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l))
+         du_trisk = du_trisk + wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l))
+         du_trisk = du_trisk + wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l))
+         du_trisk = du_trisk + wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l))
+         du_trisk = du_trisk + wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l))
+         du_trisk = du_trisk + wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l))
+         du_trisk = du_trisk + wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l))
+         du_trisk = du_trisk + wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l))
+         du_trisk = du_trisk + wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l))
+         du_trisk = du_trisk + wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
+         du(ij+u_lup,l) = du(ij+u_lup,l) + .5*du_trisk
+         du_trisk=0.
+         du_trisk = du_trisk + wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l))
+         du_trisk = du_trisk + wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l))
+         du_trisk = du_trisk + wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l))
+         du_trisk = du_trisk + wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l))
+         du_trisk = du_trisk + wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l))
+         du_trisk = du_trisk + wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l))
+         du_trisk = du_trisk + wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l))
+         du_trisk = du_trisk + wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l))
+         du_trisk = du_trisk + wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l))
+         du_trisk = du_trisk + wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
+         du(ij+u_ldown,l) = du(ij+u_ldown,l) + .5*du_trisk
+      END DO
+   END DO
+   !---------------------------- coriolis ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/div.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/div.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/div.k90	(revision 580)
@@ -0,0 +1,17 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- div ----------------------------------
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         div_ij=0.
+         div_ij = div_ij + ne_rup*le_de(ij+u_rup)*u(ij+u_rup,l)
+         div_ij = div_ij + ne_lup*le_de(ij+u_lup)*u(ij+u_lup,l)
+         div_ij = div_ij + ne_left*le_de(ij+u_left)*u(ij+u_left,l)
+         div_ij = div_ij + ne_ldown*le_de(ij+u_ldown)*u(ij+u_ldown,l)
+         div_ij = div_ij + ne_rdown*le_de(ij+u_rdown)*u(ij+u_rdown,l)
+         div_ij = div_ij + ne_right*le_de(ij+u_right)*u(ij+u_right,l)
+         divu(ij,l) = div_ij / Ai(ij)
+      END DO
+   END DO
+   !---------------------------- div ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/gradient.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/gradient.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/gradient.k90	(revision 580)
@@ -0,0 +1,12 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- gradient ----------------------------------
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         grad(ij+u_right,l) = ne_right*(b(ij+t_right,l)-b(ij,l))
+         grad(ij+u_lup,l) = ne_lup*(b(ij+t_lup,l)-b(ij,l))
+         grad(ij+u_ldown,l) = ne_ldown*(b(ij+t_ldown,l)-b(ij,l))
+      END DO
+   END DO
+   !---------------------------- gradient ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/kernels/pvort_only.k90
===================================================================
--- codes/icosagcm/trunk/src/kernels/pvort_only.k90	(revision 580)
+++ codes/icosagcm/trunk/src/kernels/pvort_only.k90	(revision 580)
@@ -0,0 +1,35 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- pvort_only ----------------------------------
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         etav = 0.d0
+         etav = etav + ne_rup*u(ij+u_rup,l)
+         etav = etav + ne_left*u(ij+t_rup+u_left,l)
+         etav = etav + (-ne_lup)*u(ij+u_lup,l)
+         hv=0.
+         hv = hv + Riv2(ij,vup)*rhodz(ij,l)
+         hv = hv + Riv2(ij+t_rup,vldown)*rhodz(ij+t_rup,l)
+         hv = hv + Riv2(ij+t_lup,vrdown)*rhodz(ij+t_lup,l)
+         qv(ij+z_up,l) = (etav + fv(ij+z_up)*Av(ij+z_up) )/(hv*Av(ij+z_up))
+         etav = 0.d0
+         etav = etav + (-ne_rdown)*u(ij+u_rdown,l)
+         etav = etav + ne_right*u(ij+t_ldown+u_right,l)
+         etav = etav + ne_ldown*u(ij+u_ldown,l)
+         hv=0.
+         hv = hv + Riv2(ij,vdown)*rhodz(ij,l)
+         hv = hv + Riv2(ij+t_ldown,vrup)*rhodz(ij+t_ldown,l)
+         hv = hv + Riv2(ij+t_rdown,vlup)*rhodz(ij+t_rdown,l)
+         qv(ij+z_down,l) = (etav + fv(ij+z_down)*Av(ij+z_down) )/(hv*Av(ij+z_down))
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin, ij_end
+         qu(ij+u_right,l)=0.5d0*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))
+         qu(ij+u_lup,l)=0.5d0*(qv(ij+z_up,l)+qv(ij+z_lup,l))
+         qu(ij+u_ldown,l)=0.5d0*(qv(ij+z_ldown,l)+qv(ij+z_down,l))
+      END DO
+   END DO
+   !---------------------------- pvort_only ----------------------------------
+   !--------------------------------------------------------------------------
