Index: codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90	(revision 614)
@@ -0,0 +1,50 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_fast ----------------------------------
+   !
+   SELECT CASE(caldyn_thermo)
+   CASE(thermo_boussinesq)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         DO l = 1, llm
+            berni(l,ij) = pk(l,ij)
+            ! from now on pk contains the vertically-averaged geopotential
+            pk(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij))
+         END DO
+      END DO
+      !$OMP END DO
+   CASE(thermo_theta)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         DO l = 1, llm
+            berni(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij))
+         END DO
+      END DO
+      !$OMP END DO
+   CASE(thermo_entropy)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         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)
+         END DO
+      END DO
+      !$OMP END DO
+   CASE DEFAULT
+      PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo ! FIXME
+      STOP
+   END SELECT
+   !
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      DO l = 1, llm
+         due = .5*(theta(l,ij_left,1)+theta(l,ij_right,1))*(pk(l,ij_right)-pk(l,ij_left)) + berni(l,ij_right)-berni(l,ij_left)
+         du(l,edge) = du(l,edge) - 1.*due
+         u(l,edge) = u(l,edge) + tau*du(l,edge)
+      END DO
+   END DO
+   !$OMP END DO
+   !
+   !---------------------------- caldyn_fast ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90	(revision 614)
@@ -0,0 +1,110 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_slow_NH ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      kdown = 1
+      kup = 1
+      l=1
+      w_il(l,ij) = 2.*W(l,ij)/rhodz(kup,ij)
+      DO l = 2, llm
+         kdown = l-1
+         kup = l
+         w_il(l,ij) = 2.*W(l,ij)/(rhodz(kdown,ij)+rhodz(kup,ij))
+      END DO
+      kdown = llm
+      kup = llm
+      l=llm+1
+      w_il(l,ij) = 2.*W(l,ij)/rhodz(kdown,ij)
+   END DO
+   !$OMP END DO
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      kdown = 1
+      kup = 1
+      l=1
+      ! 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(l,ij_right)+W(l,ij_left) )
+      DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left))
+      F_el(l,edge) = DePhil(l,edge)*W_el
+      W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) )
+      v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked
+      G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el
+      DO l = 2, llm
+         kdown = l-1
+         kup = l
+         ! 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(l,ij_right)+W(l,ij_left) )
+         DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left))
+         F_el(l,edge) = DePhil(l,edge)*W_el
+         W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) )
+         v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked
+         G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el
+      END DO
+      kdown = llm
+      kup = llm
+      l=llm+1
+      ! 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(l,ij_right)+W(l,ij_left) )
+      DePhil(l,edge) = 1.*(Phi(l,ij_right)-Phi(l,ij_left))
+      F_el(l,edge) = DePhil(l,edge)*W_el
+      W2_el = .5*le_de(edge) * ( W(l,ij_left)*w_il(l,ij_left) + W(l,ij_right)*w_il(l,ij_right) )
+      v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked
+      G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el
+   END DO
+   !$OMP END DO
+   ! compute GradPhi2, dPhi, dW
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      DO l = 1, llm+1
+         gPhi2=0.
+         dP=0.
+         divG=0
+         DO iedge = 1, primal_deg(ij)
+            edge = primal_edge(iedge,ij)
+            gPhi2 = gPhi2 + le_de(edge)*DePhil(l,edge)**2
+            dP = dP + le_de(edge)*DePhil(l,edge)*v_el(l,edge)
+            divG = divG + primal_ne(iedge,ij)*G_el(l,edge) ! -div(G_el), G_el already has le_de
+         END DO
+         gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
+         dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
+         dW(l,ij) = (-1./Ai(ij))*divG
+      END DO
+   END DO
+   !$OMP 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
+   ! Compute berni at scalar points
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      DO l = 1, llm
+         u2=0.
+         DO iedge = 1, primal_deg(ij)
+            edge = primal_edge(iedge,ij)
+            u2 = u2 + le_de(edge)*u(l,edge)**2
+         END DO
+         berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
+      END DO
+   END DO
+   !$OMP END DO
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      DO l = 1, llm
+         ! Compute mass flux and grad(berni)
+         uu = .5*(rhodz(l,ij_left)+rhodz(l,ij_right))*u(l,edge) - .5*( F_el(l,edge)+F_el(l+1,edge) )
+         hflux(l,edge) = le_de(edge)*uu
+         du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right))
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- caldyn_slow_NH ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_slow_hydro.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_slow_hydro.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_slow_hydro.k90	(revision 614)
@@ -0,0 +1,47 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_slow_hydro ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      DO l = 1, llm
+         uu = .5*(rhodz(l,ij_left)+rhodz(l,ij_right))*u(l,edge)
+         hflux(l,edge) = uu*le_de(edge)
+      END DO
+   END DO
+   !$OMP END DO
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      DO l = 1, llm
+         ke=0.d0
+         DO iedge = 1, primal_deg(ij)
+            edge = primal_edge(iedge,ij)
+            ke = ke + le_de(edge)*u(l,edge)**2
+         END DO
+         BERNI(l,ij)=ke*(.25/Ai(ij))
+      END DO
+   END DO
+   !$OMP END DO
+   IF(zero) THEN
+      !$OMP DO SCHEDULE(STATIC)
+      DO edge = 1, edge_num
+         ij_left = left(edge)
+         ij_right = right(edge)
+         DO l = 1, llm
+            du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right)) ! minus gradient
+         END DO
+      END DO
+      !$OMP END DO
+   ELSE
+      !$OMP DO SCHEDULE(STATIC)
+      DO edge = 1, edge_num
+         ij_left = left(edge)
+         ij_right = right(edge)
+         DO l = 1, llm
+            du(l,edge) = du(l,edge) + 1.*(berni(l,ij_left)-berni(l,ij_right)) ! minus gradient
+         END DO
+      END DO
+      !$OMP END DO
+   END IF
+   !---------------------------- caldyn_slow_hydro ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90	(revision 614)
@@ -0,0 +1,98 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_solver ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      l=1
+      m_il(l,ij) = .5*rhodz(l,ij)
+      DO l = 2, llm
+         m_il(l,ij) = .5*(rhodz(l,ij)+rhodz(l-1,ij))
+      END DO
+      l=llm+1
+      m_il(l,ij) = .5*rhodz(l-1,ij)
+   END DO
+   !$OMP END DO
+   !
+   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
+   SELECT CASE(caldyn_thermo)
+   CASE(thermo_theta)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         DO l = 1, llm
+            rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij))
+            X_ij = (cpp/preff)*kappa*theta(l,ij,1)*rho_ij
+            ! kappa.theta.rho = p/exner
+            ! => X = (p/p0)/(exner/Cp)
+            ! = (p/p0)^(1-kappa)
+            pres(l,ij) = preff*(X_ij**gamma) ! pressure
+            ! Compute Exner function (needed by compute_caldyn_fast)
+            ! other formulae possible if exponentiation is slow
+            pk(l,ij) = cpp*((pres(l,ij)/preff)**kappa) ! Exner
+         END DO
+      END DO
+      !$OMP END DO
+   CASE(thermo_entropy)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         DO l = 1, llm
+            rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij))
+            T_ij = Treff*exp( (theta(l,ij,1)+Rd*log(vreff*rho_ij))/Cvd )
+            pres(l,ij) = rho_ij*Rd*T_ij
+            pk(l,ij) = T_ij
+         END DO
+      END DO
+      !$OMP END DO
+   CASE DEFAULT
+      STOP
+   END SELECT
+   ! We need a barrier here because we compute pres above and do a vertical difference below
+   !$OMP BARRIER
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      l=1
+      ! Lower BC
+      dW(l,ij) = (1./g)*(pbot-rho_bot*(geopot(l,ij)-PHI_BOT(ij))-pres(l,ij)) - m_il(l,ij)
+      W(l,ij) = W(l,ij)+tau*dW(l,ij) ! update W
+      dPhi(l,ij) = g*g*W(l,ij)/m_il(l,ij)
+      DO l = 2, llm
+         dW(l,ij) = (1./g)*(pres(l-1,ij)-pres(l,ij)) - m_il(l,ij)
+         W(l,ij) = W(l,ij)+tau*dW(l,ij) ! update W
+         dPhi(l,ij) = g*g*W(l,ij)/m_il(l,ij)
+      END DO
+      l=llm+1
+      ! Top BC
+      dW(l,ij) = (1./g)*(pres(l-1,ij)-ptop) - m_il(l,ij)
+      W(l,ij) = W(l,ij)+tau*dW(l,ij) ! update W
+      dPhi(l,ij) = g*g*W(l,ij)/m_il(l,ij)
+   END DO
+   !$OMP END DO
+   ! We need a barrier here because we update W above and do a vertical average below
+   !$OMP BARRIER
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      DO l = 1, llm
+         ! compute du = -0.5*g^2.grad(w^2)
+         berni(l,ij) = (-.25*g*g)*((W(l,ij)/m_il(l,ij))**2 + (W(l+1,ij)/m_il(l+1,ij))**2 )
+      END DO
+   END DO
+   !$OMP END DO
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      DO l = 1, llm
+         du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right))
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- caldyn_solver ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90	(revision 614)
@@ -0,0 +1,63 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_vert ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij=1,primal_num
+      DO l = llm-1,1,-1
+         ! cumulate mass flux convergence from top to bottom
+         convm(l,ij) = convm(l,ij) + convm(l+1,ij)
+      END DO
+      l=1
+      dmass_col(ij) = convm(l,ij)
+      DO l = 2,llm
+         ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
+         wflux(l,ij) = mass_bl(l,ij) * dmass_col(ij) - convm(l,ij)
+      END DO
+   END DO
+   !$OMP END DO
+   ! make sure wflux is up to date
+   !$OMP BARRIER
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      DO l = 2, llm
+         wwuu(l,edge) = .5*(wflux(l,ij_left)+wflux(l,ij_right))*(u(l,edge)-u(l-1,edge))
+      END DO
+   END DO
+   !$OMP END DO
+   ! make sure wwuu is up to date
+   !$OMP BARRIER
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      DO l = 1, llm
+         convm(l,ij) = mass_dbk(l,ij) * dmass_col(ij) ! FIXME : we shoud prognose mass_col
+      END DO
+   END DO
+   !$OMP END DO
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      DO l = 1, llm
+         du(l,edge) = du(l,edge) - (wwuu(l,edge)+wwuu(l+1,edge)) / (rhodz(l,ij_left)+rhodz(l,ij_right))
+      END DO
+   END DO
+   !$OMP END DO
+   DO iq=1,nqdyn
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         DO l = 2, llm
+            dtheta_rhodz(l,ij,iq) = dtheta_rhodz(l,ij,iq) + 0.5*(theta(l,ij,iq)+theta(l-1,ij,iq))*wflux(l,ij)
+         END DO
+      END DO
+      !$OMP END DO
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         DO l = 1, llm-1
+            dtheta_rhodz(l,ij,iq) = dtheta_rhodz(l,ij,iq) - 0.5*(theta(l,ij,iq)+theta(l+1,ij,iq))*wflux(l+1,ij)
+         END DO
+      END DO
+      !$OMP END DO
+   END DO
+   !---------------------------- caldyn_vert ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_vert_NH.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_vert_NH.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_vert_NH.k90	(revision 614)
@@ -0,0 +1,59 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_vert_NH ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      l=1
+      w_ij = ( W(l,ij)+.5*W(l+1,ij) )/mass(l,ij)
+      wflux_ij = .5*(wflux(l,ij)+wflux(l+1,ij))
+      W_etadot(l,ij) = wflux_ij*w_ij
+      eta_dot(l,ij) = wflux_ij / mass(l,ij)
+      wcov(l,ij) = w_ij*(geopot(l+1,ij)-geopot(l,ij))
+      DO l = 2, llm-1
+         w_ij = .5*( W(l,ij)+W(l+1,ij) )/mass(l,ij)
+         wflux_ij = .5*(wflux(l,ij)+wflux(l+1,ij))
+         W_etadot(l,ij) = wflux_ij*w_ij
+         eta_dot(l,ij) = wflux_ij / mass(l,ij)
+         wcov(l,ij) = w_ij*(geopot(l+1,ij)-geopot(l,ij))
+      END DO
+      l=llm
+      w_ij = ( .5*W(l,ij)+W(l+1,ij) )/mass(l,ij)
+      wflux_ij = .5*(wflux(l,ij)+wflux(l+1,ij))
+      W_etadot(l,ij) = wflux_ij*w_ij
+      eta_dot(l,ij) = wflux_ij / mass(l,ij)
+      wcov(l,ij) = w_ij*(geopot(l+1,ij)-geopot(l,ij))
+   END DO
+   !$OMP END DO
+   ! add NH term to du
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      DO l = 1, llm
+         du(l,edge) = du(l,edge) - .5*(wcov(l,ij_left)+wcov(l,ij_right))*1.*(eta_dot(l,ij_right)-eta_dot(l,ij_left))
+      END DO
+   END DO
+   !$OMP END DO
+   ! add NH terms to dW, dPhi
+   ! FIXME : TODO top and bottom
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      DO l = 2, llm
+         dPhi(l,ij)=dPhi(l,ij)-wflux(l,ij)*(geopot(l+1,ij)-geopot(l-1,ij))/(mass(l-1,ij)+mass(l,ij))
+      END DO
+   END DO
+   !$OMP END DO
+   ! We need a barrier here because we compute W_etadot above and do a vertical difference below
+   !$OMP BARRIER
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      l=1
+      dW(l,ij) = dW(l,ij) - W_etadot(l,ij)
+      DO l = 2, llm
+         dW(l,ij) = dW(l,ij) + W_etadot(l-1,ij) - W_etadot(l,ij)
+      END DO
+      l=llm+1
+      dW(l,ij) = dW(l,ij) + W_etadot(l-1,ij)
+   END DO
+   !$OMP END DO
+   !---------------------------- caldyn_vert_NH ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90	(revision 614)
@@ -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
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij=1,primal_num
+      DO l = 1,llm+1
+         Phi_star_il(l,ij) = Phi_il(l,ij) + tau*g2*(W_il(l,ij)/m_il(l,ij)-tau)
+      END DO
+   END DO
+   !$OMP 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)
+         !$OMP DO SCHEDULE(STATIC)
+         DO ij=1,primal_num
+            DO l = 1,llm
+               rho_ij = (g*m_ik(l,ij))/(Phi_il(l+1,ij)-Phi_il(l,ij))
+               X_ij = (cpp/preff)*kappa*theta(l,ij)*rho_ij
+               p_ik(l,ij) = preff*(X_ij**gamma)
+               c2_mik = gamma*p_ik(l,ij)/(rho_ij*m_ik(l,ij)) ! c^2 = gamma*R*T = gamma*p/rho
+               A_ik(l,ij) = c2_mik*(tau/g*rho_ij)**2
+            END DO
+         END DO
+         !$OMP END DO
+      CASE(thermo_entropy)
+         !$OMP DO SCHEDULE(STATIC)
+         DO ij=1,primal_num
+            DO l = 1,llm
+               rho_ij = (g*m_ik(l,ij))/(Phi_il(l+1,ij)-Phi_il(l,ij))
+               X_ij = Treff*exp(theta(l,ij)/cpp) ! theta = Tref.exp(s/Cp)
+               X_ij = (cpp/preff)*kappa*X_ij*rho_ij
+               p_ik(l,ij) = preff*(X_ij**gamma)
+               c2_mik = gamma*p_ik(l,ij)/(rho_ij*m_ik(l,ij)) ! c^2 = gamma*R*T = gamma*p/rho
+               A_ik(l,ij) = c2_mik*(tau/g*rho_ij)**2
+            END DO
+         END DO
+         !$OMP 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
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         ! Compute residual R_il and B_il
+         l=1
+         ! bottom interface l=1
+         ml_g2 = gm2*m_il(l,ij)
+         B_il(l,ij) = A_ik(l,ij) + ml_g2 + tau2_g*rho_bot
+         R_il(l,ij) = ml_g2*( Phi_il(l,ij)-Phi_star_il(l,ij)) &
+         + tau2_g*( p_ik(l,ij)-pbot+rho_bot*(Phi_il(l,ij)-PHI_BOT(ij)) )
+         DO l = 2,llm
+            ! inner interfaces
+            ml_g2 = gm2*m_il(l,ij)
+            B_il(l,ij) = A_ik(l,ij)+A_ik(l-1,ij) + ml_g2
+            R_il(l,ij) = ml_g2*( Phi_il(l,ij)-Phi_star_il(l,ij)) &
+            + tau2_g*(p_ik(l,ij)-p_ik(l-1,ij))
+            ! consistency check : if Wil=0 and initial state is in hydrostatic balance
+            ! then Phi_star_il(l,ij) = Phi_il(l,ij) - tau^2*g^2
+            ! and residual = tau^2*(ml+(1/g)dl_pi)=0
+         END DO
+         l=llm+1
+         ! top interface l=llm+1
+         ml_g2 = gm2*m_il(l,ij)
+         B_il(l,ij) = A_ik(l-1,ij) + ml_g2
+         R_il(l,ij) = ml_g2*( Phi_il(l,ij)-Phi_star_il(l,ij)) &
+         + tau2_g*( ptop-p_ik(l-1,ij) )
+         !
+         ! 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))
+         l=1
+         X_ij = 1./B_il(l,ij)
+         C_ik(l,ij) = -A_ik(l,ij) * X_ij
+         D_il(l,ij) = R_il(l,ij) * X_ij
+         DO l = 2,llm
+            X_ij = 1./( B_il(l,ij) + A_ik(l-1,ij)*C_ik(l-1,ij) )
+            C_ik(l,ij) = -A_ik(l,ij) * X_ij
+            D_il(l,ij) = (R_il(l,ij)+A_ik(l-1,ij)*D_il(l-1,ij)) * X_ij
+         END DO
+         l=llm+1
+         X_ij = 1./( B_il(l,ij) + A_ik(l-1,ij)*C_ik(l-1,ij) )
+         D_il(l,ij) = (R_il(l,ij)+A_ik(l-1,ij)*D_il(l-1,ij)) * 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(l,ij) = D_il(l,ij)
+         Phi_il(l,ij) = Phi_il(l,ij) - x_il(l,ij)
+         DO l = llm,1,-1
+            ! Back substitution at lower interfaces
+            x_il(l,ij) = D_il(l,ij) - C_ik(l,ij)*x_il(l+1,ij)
+            Phi_il(l,ij) = Phi_il(l,ij) - x_il(l,ij)
+         END DO
+      END DO
+      !$OMP 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/devel/src/kernels_unst/compute_geopot.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/compute_geopot.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/compute_geopot.k90	(revision 614)
@@ -0,0 +1,110 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- compute_geopot ----------------------------------
+   SELECT CASE(caldyn_thermo)
+   CASE(thermo_boussinesq)
+      ! use hydrostatic balance with theta*rhodz to find pk (=Lagrange multiplier=pressure)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         l=llm
+         pk(l,ij) = ptop + .5*g* theta(l,ij,1)*rhodz(l,ij)
+         DO l = llm-1,1,-1
+            pk(l,ij) = pk(l+1,ij) + (.5*g)*( theta(l,ij,1)*rhodz(l,ij) + theta(l+1,ij,1)*rhodz(l+1,ij) )
+         END DO
+         IF(caldyn_eta == eta_lag) THEN
+            l=1
+            ps(ij) = pk(l,ij) + .5*g* theta(l,ij,1)*rhodz(l,ij)
+         END IF
+      END DO
+      !$OMP END DO
+      ! now pk contains the Lagrange multiplier (pressure)
+      ! specific volume 1 = dphi/g/rhodz
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         DO l = 1,llm
+            geopot(l+1,ij) = geopot(l,ij) + g*rhodz(l,ij)
+         END DO
+      END DO
+      !$OMP END DO
+   CASE(thermo_theta)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         l=llm
+         pk(l,ij) = ptop + .5*g* rhodz(l,ij)
+         DO l = llm-1,1,-1
+            pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
+         END DO
+         IF(caldyn_eta == eta_lag) THEN
+            l=1
+            ps(ij) = pk(l,ij) + .5*g* rhodz(l,ij)
+         END IF
+      END DO
+      !$OMP END DO
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         DO l = 1,llm
+            p_ik = pk(l,ij)
+            exner_ik = cpp * (p_ik/preff) ** kappa
+            gv = (g*kappa)*theta(l,ij,1)*exner_ik/p_ik
+            pk(l,ij) = exner_ik
+            geopot(l+1,ij) = geopot(l,ij) + gv*rhodz(l,ij)
+         END DO
+      END DO
+      !$OMP END DO
+   CASE(thermo_entropy)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         l=llm
+         pk(l,ij) = ptop + .5*g* rhodz(l,ij)
+         DO l = llm-1,1,-1
+            pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
+         END DO
+         IF(caldyn_eta == eta_lag) THEN
+            l=1
+            ps(ij) = pk(l,ij) + .5*g* rhodz(l,ij)
+         END IF
+      END DO
+      !$OMP END DO
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         DO l = 1,llm
+            p_ik = pk(l,ij)
+            temp_ik = Treff*exp((theta(l,ij,1) + Rd*log(p_ik/preff))/cpp)
+            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
+            pk(l,ij) = temp_ik
+            geopot(l+1,ij) = geopot(l,ij) + gv*rhodz(l,ij)
+         END DO
+      END DO
+      !$OMP END DO
+   CASE(thermo_moist)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         l=llm
+         pk(l,ij) = ptop + .5*g* rhodz(l,ij)*(1.+theta(l,ij,2))
+         DO l = llm-1,1,-1
+            pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij)*(1.+theta(l,ij,2)) + rhodz(l+1,ij)*(1.+theta(l+1,ij,2)) )
+         END DO
+         IF(caldyn_eta == eta_lag) THEN
+            l=1
+            ps(ij) = pk(l,ij) + .5*g* rhodz(l,ij)*(1.+theta(l,ij,2))
+         END IF
+      END DO
+      !$OMP END DO
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         DO l = 1,llm
+            p_ik = pk(l,ij)
+            qv = theta(l,ij,2) ! water vaper mixing ratio = mv/md
+            Rmix = Rd+qv*Rv
+            chi = ( theta(l,ij,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(l,ij) = temp_ik
+            geopot(l+1,ij) = geopot(l,ij) + gv*rhodz(l,ij)
+         END DO
+      END DO
+      !$OMP END DO
+   END SELECT
+   !---------------------------- compute_geopot ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/coriolis.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/coriolis.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/coriolis.k90	(revision 614)
@@ -0,0 +1,54 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- coriolis ----------------------------------
+   !
+   DO iq=1,nqdyn
+      !$OMP DO SCHEDULE(STATIC)
+      DO edge = 1, edge_num
+         ij_left = left(edge)
+         ij_right = right(edge)
+         DO l = 1, llm
+            Ftheta(l,edge) = .5*(theta(l,ij_left,iq)+theta(l,ij_right,iq))*hflux(l,edge)
+         END DO
+      END DO
+      !$OMP END DO
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         DO l = 1, llm
+            divF=0.
+            DO iedge = 1, primal_deg(ij)
+               edge = primal_edge(iedge,ij)
+               divF = divF + Ftheta(l,edge)*primal_ne(iedge,ij)
+            END DO
+            dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
+         END DO
+      END DO
+      !$OMP END DO
+   END DO ! iq
+   !
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      DO l = 1, llm
+         divF=0.
+         DO iedge = 1, primal_deg(ij)
+            edge = primal_edge(iedge,ij)
+            divF = divF + hflux(l,edge)*primal_ne(iedge,ij)
+         END DO
+         convm(l,ij) = -divF / Ai(ij)
+      END DO
+   END DO
+   !$OMP END DO
+   !
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      DO l = 1, llm
+         du_trisk=0.
+         DO itrisk = 1, trisk_deg(edge)
+            edge_trisk = trisk(itrisk,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+         END DO
+         du(l,edge) = du(l,edge) + .5*du_trisk
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- coriolis ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/div.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/div.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/div.k90	(revision 614)
@@ -0,0 +1,16 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- div ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      DO l = 1, llm
+         div_ij=0.
+         DO iedge = 1, primal_deg(ij)
+            edge = primal_edge(iedge,ij)
+            div_ij = div_ij + primal_ne(iedge,ij)*le_de(edge)*u(l,edge)
+         END DO
+         divu(l,ij) = div_ij / Ai(ij)
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- div ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/gradient.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/gradient.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/gradient.k90	(revision 614)
@@ -0,0 +1,13 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- gradient ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      DO l = 1, llm
+         grad(l,edge) = 1.*(b(l,ij_right)-b(l,ij_left))
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- gradient ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/pvort_only.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/pvort_only.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/pvort_only.k90	(revision 614)
@@ -0,0 +1,30 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- pvort_only ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, dual_num
+      DO l = 1, llm
+         etav = 0.d0
+         DO iedge = 1, dual_deg(ij)
+            edge = dual_edge(iedge,ij)
+            etav = etav + dual_ne(iedge,ij)*u(l,edge)
+         END DO
+         hv=0.
+         DO ivertex = 1, dual_deg(ij)
+            vertex = dual_vertex(ivertex,ij)
+            hv = hv + Riv2(ivertex,ij)*rhodz(l,vertex)
+         END DO
+         qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij))
+      END DO
+   END DO
+   !$OMP END DO
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge = 1, edge_num
+      ij_up = up(edge)
+      ij_down = down(edge)
+      DO l = 1, llm
+         qu(l,edge)=0.5d0*(qv(l,ij_down)+qv(l,ij_up))
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- pvort_only ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/theta.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/theta.k90	(revision 614)
+++ codes/icosagcm/devel/src/kernels_unst/theta.k90	(revision 614)
@@ -0,0 +1,37 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- theta ----------------------------------
+   IF(caldyn_eta==eta_mass) THEN ! Compute mass
+      ! FIXME : here mass_col is computed from rhodz
+      ! so that the DOFs are the same whatever caldyn_eta
+      ! in DYNAMICO mass_col is prognosed rather than rhodz
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         l=0
+         mass_col(ij)=0.
+         DO l = 1,llm
+            mass_col(ij) = mass_col(ij) + rhodz(l,ij)
+         END DO
+      END DO
+      !$OMP END DO
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         DO l = 1, llm
+            ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass
+            ! m = mass_dak(l,ij)+(mass_col(ij)*g+ptop)*mass_dbk(l,ij)
+            ! rhodz(l,ij) = m/g
+            rhodz(l,ij) = mass_dak(l,ij) + mass_col(ij)*mass_dbk(l,ij)
+         END DO
+      END DO
+      !$OMP END DO
+   END IF
+   DO iq=1,nqdyn
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij = 1, primal_num
+         DO l = 1, llm
+            theta(l,ij,iq) = theta_rhodz(l,ij,iq)/rhodz(l,ij)
+         END DO
+      END DO
+      !$OMP END DO
+   END DO
+   !---------------------------- theta ----------------------------------
+   !--------------------------------------------------------------------------
