Index: codes/icosagcm/devel/src/base/earth_const.f90
===================================================================
--- codes/icosagcm/devel/src/base/earth_const.f90	(revision 623)
+++ codes/icosagcm/devel/src/base/earth_const.f90	(revision 624)
@@ -28,6 +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)
+  LOGICAL :: dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro, dysl_caldyn_vert
+  !$OMP THREADPRIVATE(dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro, dysl_caldyn_vert)
 
 CONTAINS
Index: codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90	(revision 623)
+++ codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90	(revision 624)
@@ -34,4 +34,6 @@
     dysl_caldyn_coriolis=dysl
     CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis)
+    dysl_caldyn_vert=dysl
+    CALL getin("dysl_caldyn_vert",dysl_caldyn_vert)
 
     def='energy'
Index: codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90	(revision 623)
+++ codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90	(revision 624)
@@ -180,8 +180,9 @@
     REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
+    INTEGER, PARAMETER :: caldyn_vert_cons=1, caldyn_vert_variant=caldyn_vert_cons
 
     ! temporary variable    
     INTEGER :: i,j,ij,l,iq
-    REAL(rstd) :: p_ik, exner_ik
+    REAL(rstd) :: p_ik, exner_ik, dF_deta, dFu_deta
     INTEGER    :: ij_omp_begin, ij_omp_end
 
@@ -190,17 +191,16 @@
     CALL distrib_level(ij_begin,ij_end, ij_omp_begin,ij_omp_end)
 
-    !    REAL(rstd) :: wwuu(iim*3*jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp
-    ! need to be understood
-
-    !    wwuu=wwuu_out
-    CALL trace_start("compute_caldyn_vert")
+    IF(dysl_caldyn_vert) THEN
+#define mass_bl(ij,l) bp(ij)
+#define dmass_col(ij) dps(ij)
+#include "../kernels_hex/caldyn_wflux.k90"
+#include "../kernels_hex/caldyn_vert.k90"
+#undef mass_bl
+#undef dmass_col
+    ELSE
 
     !$OMP BARRIER   
 !!! cumulate mass flux convergence from top to bottom
-    !  IF (is_omp_level_master) THEN
     DO  l = llm-1, 1, -1
-       !    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
-
-!!$OMP DO SCHEDULE(STATIC) 
        !DIR$ SIMD
        DO ij=ij_omp_begin,ij_omp_end         
@@ -212,7 +212,5 @@
     !$OMP BARRIER
     ! FLUSH on convm
-!!!!!!!!!!!!!!!!!!!!!!!!!
-
-    ! compute dps
+    ! compute dmass_col
     IF (is_omp_first_level) THEN
        !DIR$ SIMD
@@ -277,12 +275,5 @@
     ENDDO
 
-    !  DO l=ll_beginp1,ll_end
-    !!DIR$ SIMD
-    !      DO ij=ij_begin,ij_end         
-    !        wwuu_out(ij+u_right,l) = wwuu(ij+u_right,l)
-    !        wwuu_out(ij+u_lup,l)   = wwuu(ij+u_lup,l)
-    !        wwuu_out(ij+u_ldown,l) = wwuu(ij+u_ldown,l)
-    !     ENDDO
-    !   ENDDO
+    END IF ! dysl
 
     CALL trace_end("compute_caldyn_vert")
Index: codes/icosagcm/devel/src/kernels_hex/caldyn_vert.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_hex/caldyn_vert.k90	(revision 624)
+++ codes/icosagcm/devel/src/kernels_hex/caldyn_vert.k90	(revision 624)
@@ -0,0 +1,64 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_vert ----------------------------------
+   DO iq=1,nqdyn
+      DO l = ll_beginp1, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            dtheta_rhodz(ij,l,iq) = dtheta_rhodz(ij,l,iq) + 0.5*(theta(ij,l,iq)+theta(ij,l-1,iq))*wflux(ij,l)
+         END DO
+      END DO
+      DO l = ll_begin, ll_endm1
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            dtheta_rhodz(ij,l,iq) = dtheta_rhodz(ij,l,iq) - 0.5*(theta(ij,l,iq)+theta(ij,l+1,iq))*wflux(ij,l+1)
+         END DO
+      END DO
+   END DO
+   IF(caldyn_vert_variant == caldyn_vert_cons) THEN
+      ! conservative vertical transport of momentum : (F/m)du/deta = 1/m (d/deta(Fu)-u.dF/deta)
+      DO l = ll_beginp1, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            wwuu(ij+u_right,l) = .25*(wflux(ij,l)+wflux(ij+t_right,l))*(u(ij+u_right,l)+u(ij+u_right,l-1)) ! Fu
+            wwuu(ij+u_lup,l) = .25*(wflux(ij,l)+wflux(ij+t_lup,l))*(u(ij+u_lup,l)+u(ij+u_lup,l-1)) ! Fu
+            wwuu(ij+u_ldown,l) = .25*(wflux(ij,l)+wflux(ij+t_ldown,l))*(u(ij+u_ldown,l)+u(ij+u_ldown,l-1)) ! Fu
+         END DO
+      END DO
+      ! make sure wwuu is up to date
+      !$OMP BARRIER
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            dFu_deta = wwuu(ij+u_right,l+1)-wwuu(ij+u_right,l) ! d/deta (F*u)
+            dF_deta = .5*(wflux(ij,l+1)+wflux(ij+t_right,l+1)-(wflux(ij,l)+wflux(ij+t_right,l))) ! d/deta(F)
+            du(ij+u_right,l) = du(ij+u_right,l) - (dFu_deta-u(ij+u_right,l)*dF_deta) / (.5*(rhodz(ij,l)+rhodz(ij+t_right,l))) ! (F/m)du/deta
+            dFu_deta = wwuu(ij+u_lup,l+1)-wwuu(ij+u_lup,l) ! d/deta (F*u)
+            dF_deta = .5*(wflux(ij,l+1)+wflux(ij+t_lup,l+1)-(wflux(ij,l)+wflux(ij+t_lup,l))) ! d/deta(F)
+            du(ij+u_lup,l) = du(ij+u_lup,l) - (dFu_deta-u(ij+u_lup,l)*dF_deta) / (.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))) ! (F/m)du/deta
+            dFu_deta = wwuu(ij+u_ldown,l+1)-wwuu(ij+u_ldown,l) ! d/deta (F*u)
+            dF_deta = .5*(wflux(ij,l+1)+wflux(ij+t_ldown,l+1)-(wflux(ij,l)+wflux(ij+t_ldown,l))) ! d/deta(F)
+            du(ij+u_ldown,l) = du(ij+u_ldown,l) - (dFu_deta-u(ij+u_ldown,l)*dF_deta) / (.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))) ! (F/m)du/deta
+         END DO
+      END DO
+   ELSE
+      DO l = ll_beginp1, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            wwuu(ij+u_right,l) = .5*(wflux(ij,l)+wflux(ij+t_right,l))*(u(ij+u_right,l)-u(ij+u_right,l-1))
+            wwuu(ij+u_lup,l) = .5*(wflux(ij,l)+wflux(ij+t_lup,l))*(u(ij+u_lup,l)-u(ij+u_lup,l-1))
+            wwuu(ij+u_ldown,l) = .5*(wflux(ij,l)+wflux(ij+t_ldown,l))*(u(ij+u_ldown,l)-u(ij+u_ldown,l-1))
+         END DO
+      END DO
+      ! make sure wwuu is up to date
+      !$OMP BARRIER
+      DO l = ll_begin, ll_end
+         !DIR$ SIMD
+         DO ij=ij_begin, ij_end
+            du(ij+u_right,l) = du(ij+u_right,l) - (wwuu(ij+u_right,l)+wwuu(ij+u_right,l+1)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
+            du(ij+u_lup,l) = du(ij+u_lup,l) - (wwuu(ij+u_lup,l)+wwuu(ij+u_lup,l+1)) / (rhodz(ij,l)+rhodz(ij+t_lup,l))
+            du(ij+u_ldown,l) = du(ij+u_ldown,l) - (wwuu(ij+u_ldown,l)+wwuu(ij+u_ldown,l+1)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
+         END DO
+      END DO
+   END IF
+   !---------------------------- caldyn_vert ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_hex/caldyn_wflux.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_hex/caldyn_wflux.k90	(revision 624)
+++ codes/icosagcm/devel/src/kernels_hex/caldyn_wflux.k90	(revision 624)
@@ -0,0 +1,21 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_wflux ----------------------------------
+   DO l = llm-1,1,-1
+      DO ij=ij_omp_begin,ij_omp_end
+         ! cumulate mass flux convergence from top to bottom
+         convm(ij,l) = convm(ij,l) + convm(ij,l+1)
+      END DO
+   END DO
+   DO ij=ij_omp_begin,ij_omp_end
+      dmass_col(ij) = convm(ij,1)
+   END DO
+   DO l = 2,llm
+      DO ij=ij_omp_begin,ij_omp_end
+         ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
+         wflux(ij,l) = mass_bl(ij,l) * dmass_col(ij) - convm(ij,l)
+      END DO
+   END DO
+   ! make sure wflux is up to date
+   !$OMP BARRIER
+   !---------------------------- caldyn_wflux ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_dmass.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_dmass.k90	(revision 624)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_dmass.k90	(revision 624)
@@ -0,0 +1,11 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_dmass ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      DO l = 1, llm
+         convm(l,ij) = mass_dbk(l,ij) * dmass_col(ij)
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- caldyn_dmass ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90	(revision 623)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90	(revision 624)
@@ -1,47 +1,4 @@
    !--------------------------------------------------------------------------
    !---------------------------- 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)
@@ -60,4 +17,50 @@
       !$OMP END DO
    END DO
+   IF(caldyn_vert_variant == caldyn_vert_cons) THEN
+      ! conservative vertical transport of momentum : (F/m)du/deta = 1/m (d/deta(Fu)-u.dF/deta)
+      !$OMP DO SCHEDULE(STATIC)
+      DO edge = 1, edge_num
+         ij_left = left(edge)
+         ij_right = right(edge)
+         DO l = 2, llm
+            wwuu(l,edge) = .25*(wflux(l,ij_left)+wflux(l,ij_right))*(u(l,edge)+u(l-1,edge)) ! Fu
+         END DO
+      END DO
+      !$OMP END DO
+      ! make sure wwuu 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 = 1, llm
+            dFu_deta = wwuu(l+1,edge)-wwuu(l,edge) ! d/deta (F*u)
+            dF_deta = .5*(wflux(l+1,ij_left)+wflux(l+1,ij_right)-(wflux(l,ij_left)+wflux(l,ij_right))) ! d/deta(F)
+            du(l,edge) = du(l,edge) - (dFu_deta-u(l,edge)*dF_deta) / (.5*(rhodz(l,ij_left)+rhodz(l,ij_right))) ! (F/m)du/deta
+         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 = 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 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
+   END IF
    !---------------------------- caldyn_vert ----------------------------------
    !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_wflux.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_wflux.k90	(revision 624)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_wflux.k90	(revision 624)
@@ -0,0 +1,20 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_wflux ----------------------------------
+   !$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
+   !---------------------------- caldyn_wflux ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 623)
+++ codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 624)
@@ -2,10 +2,13 @@
   USE ISO_C_BINDING
   USE OMP_LIB
+#ifdef CPP_USING_XIOS
   USE xios
+#endif
   IMPLICIT NONE
   PRIVATE
   SAVE
 
-#define BINDC(thename) BIND(C, name='dynamico_' // #thename)
+#define BINDC_(thename) BIND(C, name=#thename)
+#define BINDC(thename) BINDC_(dynamico_ ## thename)
 !
 #define DBL REAL(C_DOUBLE)
@@ -15,8 +18,11 @@
 #define INDEX INTEGER(C_INT)
 #define PHI_BOT(ij) Phi_bot
-#define PHI_BOT_VAR 0.d0
-!
-  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2, thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4
-  INTEGER(C_INT),  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, nb_threads=1
+#define PHI_BOT_VAR 0.
+!
+  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2, &
+       thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4, &
+       caldyn_vert_cons=1
+  INTEGER(C_INT),  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, &
+       caldyn_vert_variant=caldyn_vert_cons, nb_threads=1
   LOGICAL(C_BOOL), BIND(C) :: hydrostatic=.TRUE., debug_hevi_solver=.TRUE., rigid=.TRUE.
 !
@@ -40,5 +46,8 @@
 !
   INTEGER(C_INT), BIND(C) :: comm_icosa
+
+#ifdef CPP_USING_XIOS
   TYPE(xios_context) :: ctx_hdl
+#endif
 
 CONTAINS
@@ -199,5 +208,8 @@
   DOUBLE2(llm+1, edge_num) :: wwuu
   DECLARE_INDICES
+  DBL :: dF_deta, dFu_deta
   wwuu=0.
+#include "../kernels_unst/caldyn_wflux.k90"
+#include "../kernels_unst/caldyn_dmass.k90"
 #include "../kernels_unst/caldyn_vert.k90"
 END SUBROUTINE caldyn_vert
@@ -392,4 +404,6 @@
 !---------------------------------------------- XIOS -----------------------------------------
 
+#ifdef CPP_USING_XIOS
+
   SUBROUTINE setup_xios() BINDC(setup_xios)
     ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
@@ -443,3 +457,5 @@
   END SUBROUTINE call_xios_update_calendar
 
+#endif
+
 END MODULE caldyn_unstructured_mod
