Index: codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin
===================================================================
--- codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin	(revision 927)
+++ codes/icosagcm/devel/Python/src/kernels_caldyn_hevi.jin	(revision 928)
@@ -1,79 +1,2 @@
-KERNEL(caldyn_wflux)
-  SEQUENCE_C0
-    BODY('llm-1,1,-1')
-      ! cumulate mass flux convergence from top to bottom
-      convm(CELL) = convm(CELL) + convm(UP(CELL))
-    END_BLOCK
-    EPILOGUE(1)
-      dmass_col(HIDX(CELL)) = convm(CELL)
-    END_BLOCK
-    BODY('2,llm')
-      ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
-      wflux(CELL) = mass_bl(CELL) * dmass_col(HIDX(CELL)) - convm(CELL) 
-    END_BLOCK
-  END_BLOCK
-  ! make sure wflux is up to date
-  BARRIER
-END_BLOCK
-
-KERNEL(caldyn_dmass)
-  FORALL_CELLS()
-    ON_PRIMAL
-      convm(CELL) = mass_dbk(CELL) * dmass_col(HIDX(CELL))
-    END_BLOCK
-  END_BLOCK
-END_BLOCK
-
-KERNEL(caldyn_vert)
-
-  DO iq=1,nqdyn
-    FORALL_CELLS('2', 'llm')
-      ON_PRIMAL
-        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) + 0.5*(theta(CELL,iq)+theta(DOWN(CELL),iq))*wflux(CELL)  
-      END_BLOCK
-    END_BLOCK
-    FORALL_CELLS('1', 'llm-1')
-      ON_PRIMAL
-        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) - 0.5*(theta(CELL,iq)+theta(UP(CELL),iq))*wflux(UP(CELL))  
-      END_BLOCK
-    END_BLOCK
-  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)
-    FORALL_CELLS('2','llm')
-      ON_EDGES
-        wwuu(EDGE) = .25*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)+u(DOWN(EDGE))) ! Fu
-      END_BLOCK
-    END_BLOCK
-    ! make sure wwuu is up to date
-    BARRIER
-
-    FORALL_CELLS()
-      ON_EDGES
-        dFu_deta = wwuu(UP(EDGE))-wwuu(EDGE)                                          ! d/deta (F*u)
-	dF_deta  = .5*(wflux(UP(CELL1))+wflux(UP(CELL2))-(wflux(CELL1)+wflux(CELL2))) ! d/deta(F)
-        du(EDGE) = du(EDGE) - (dFu_deta-u(EDGE)*dF_deta) / (.5*(rhodz(CELL1)+rhodz(CELL2)))  ! (F/m)du/deta
-      END_BLOCK
-    END_BLOCK
-  ELSE  
-    FORALL_CELLS('2','llm')
-      ON_EDGES
-        wwuu(EDGE) = .5*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)-u(DOWN(EDGE)))
-      END_BLOCK
-    END_BLOCK
-
-    ! make sure wwuu is up to date
-    BARRIER
-
-    FORALL_CELLS()
-      ON_EDGES
-        du(EDGE) = du(EDGE) - (wwuu(EDGE)+wwuu(UP(EDGE))) / (rhodz(CELL1)+rhodz(CELL2))
-      END_BLOCK
-    END_BLOCK
-  END IF
-
-END_BLOCK
-
 KERNEL(gradient)
   FORALL_CELLS_EXT()
@@ -95,4 +18,2 @@
   END_BLOCK
 END_BLOCK
-
-
Index: codes/icosagcm/devel/XCodeML/parse_order.txt
===================================================================
--- codes/icosagcm/devel/XCodeML/parse_order.txt	(revision 927)
+++ codes/icosagcm/devel/XCodeML/parse_order.txt	(revision 928)
@@ -5,4 +5,5 @@
 ../src/base prec
 ../src/base profiling
+../src/base signature
 ../src/dcmip dcmip2016_baroclinic_wave
 ../src/dcmip dcmip2016_cyclone
@@ -12,4 +13,6 @@
 ../src/dcmip dcmip_initial_conditions_test_1_2_3_v5
 ../src/dcmip terminator
+../src/diagnostics compute_diagnostics
+../src/dynamics compute_caldyn
 ../src/output free_unit_mod
 ../ppsrc/base trace
@@ -30,15 +33,10 @@
 ../ppsrc/parallel omp_para
 ../ppsrc/parallel mpipara
-../ppsrc/unstructured data_unstructured
-../ppsrc/unstructured transfer_unstructured
-../ppsrc/unstructured transport_unstructured
 ../src/base getin
-../ppsrc/unstructured caldyn_unstructured
-../ppsrc/unstructured diags_unstructured
-../ppsrc/unstructured timestep_unstructured
 ../src/base earth_const
 ../src/base genmod
 ../src/sphere vector
 ../src/time time
+../src/output write_field_vars
 ../src/sphere spherical_geom
 ../src/sphere metric
@@ -46,18 +44,26 @@
 ../src/base dimensions
 ../src/base field
+../src/sphere geometry
+../src/output set_bounds
+../src/output write_field_mpi
+../ppsrc/physics physics_external
+../ppsrc/unstructured data_unstructured
+../ppsrc/unstructured transfer_unstructured
+../ppsrc/unstructured transport_unstructured
+../src/parallel transfert_mpi
 ../src/unstructured init_unstructured
-../ppsrc/physics physics_external
-../src/parallel transfert_mpi
+../src/output write_field
 ../ppsrc/parallel transfert
-../src/sphere geometry
-../src/output write_field
-../ppsrc/output xios_mod
+../ppsrc/unstructured caldyn_unstructured
+../ppsrc/unstructured diags_unstructured
+../ppsrc/unstructured timestep_unstructured
 ../src icosa_mod
 ../src/parallel checksum
 ../src/diagnostics kinetic
+../src/sphere compute_geometry
 ../src/physics dimphys
 ../src/physics physics_interface
-../src/dissip nudging_mod
 ../src/dynamics caldyn_vars
+../src/initial etat0_isothermal
 ../src/initial etat0_temperature
 ../src/initial q_sat
@@ -69,21 +75,22 @@
 ../src/vertical disvert_std
 ../src/vertical disvert_strato
-../src/output output_field
 ../ppsrc/diagnostics wind
+../ppsrc/dynamics compute_caldyn_kv
+../ppsrc/dynamics compute_caldyn_slow_NH
 ../ppsrc/transport advect
-../src/dcmip physics_dcmip
+../ppsrc/dynamics compute_pvort_only
 ../src/vertical disvert
 ../src/transport advect_tracer
-../src/output restart
-../src/output write_etat0
+../ppsrc/diagnostics compute_rhodz
 ../ppsrc/diagnostics diagflux
 ../ppsrc/dynamics caldyn_kernels_base
 ../ppsrc/dynamics caldyn_kernels_hevi
+../ppsrc/dynamics compute_NH_geopot
+../ppsrc/dynamics compute_theta
+../ppsrc/output xios_mod
 ../src/dcmip guided_ncar_mod
-../src/dcmip physics_dcmip2016
 ../src/diagnostics omega
 ../src/diagnostics pression
 ../src/diagnostics vorticity
-../src/dissip guided_mod
 ../src/dynamics caldyn_kernels
 ../src/initial etat0_bubble
@@ -97,24 +104,44 @@
 ../src/initial etat0_dcmip5
 ../src/initial etat0_jablonowsky06
-../src/initial etat0_start_file
 ../src/initial etat0_williamson
 ../src/time euler_scheme
 ../src/vertical vertical_interp
-../src/vertical vertical_remap
+../src/output output_field
+../src/output restart
+../src/output write_etat0
+../ppsrc/dynamics compute_caldyn_solver
+../src/dcmip physics_dcmip2016
+../src/dcmip physics_dcmip
 ../src/diagnostics exner
 ../src/diagnostics theta_rhodz
+../src/dissip nudging_mod
 ../src/dynamics caldyn_adv
-../src/initial etat0_database
+../src/initial etat0_collocated
 ../src/initial etat0_heldsz
+../src/initial etat0_start_file
 ../src/initial etat0_venus
 ../src/diagnostics geopotential_mod
+../src/dissip guided_mod
+../src/dissip sponge
+../src/initial etat0_academic
+../ppsrc/physics physics_lmdz_generic
+../ppsrc/diagnostics compute_pression
+../ppsrc/diagnostics compute_temperature
+../ppsrc/diagnostics compute_velocity
+../src/base init_grid_param
+../src/vertical vertical_remap
+../src/initial etat0_database
+../src/initial etat0
+../ppsrc/diagnostics compute_omega
+../ppsrc/dynamics compute_caldyn_vert
+../ppsrc/dynamics compute_caldyn_vert_NH
+../ppsrc/dynamics compute_geopot
+../ppsrc/dynamics compute_caldyn_Coriolis
+../ppsrc/dynamics compute_caldyn_fast
+../ppsrc/dynamics compute_caldyn_slow_hydro
 ../src/diagnostics observable
 ../src/dissip dissip_gcm
-../src/dissip sponge
 ../src/dynamics caldyn_hevi
-../src/initial etat0_academic
-../src/initial etat0
 ../src/time hevi_scheme
-../ppsrc/physics physics_lmdz_generic
 ../ppsrc/dynamics caldyn_gcm
 ../src/diagnostics check_conserve
@@ -123,4 +150,2 @@
 ../src/time explicit_scheme
 ../src/time timeloop_gcm
-../src icosa_init
-../src icosa_gcm
Index: codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90	(revision 927)
+++ codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90	(revision 928)
@@ -3,6 +3,5 @@
   USE transfert_mod
   USE caldyn_vars_mod
-  USE caldyn_kernels_hevi_mod
-  USE caldyn_kernels_base_mod
+  USE compute_caldyn_vert_mod, ONLY : compute_caldyn_vert => compute_caldyn_vert_manual
   USE caldyn_kernels_mod
   USE omp_para
@@ -171,4 +170,6 @@
     CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
     CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
+    CALL allocate_field(f_Ftheta,field_u,type_real,llm, name='Ftheta')
+    CALL allocate_field(f_berni, field_t,type_real,llm, name='berni')
     IF(.NOT.hydrostatic) THEN
        CALL allocate_field(f_Fel,      field_u,type_real,llm+1,name='F_el')
Index: codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 927)
+++ codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 928)
@@ -3,10 +3,10 @@
   USE transfert_mod
   USE caldyn_vars_mod
-  USE caldyn_kernels_hevi_mod
-  USE caldyn_kernels_base_mod
+  USE compute_caldyn_vert_mod, ONLY : compute_caldyn_vert => compute_caldyn_vert_manual
+  USE compute_caldyn_vert_NH_mod, ONLY : compute_caldyn_vert_NH => compute_caldyn_vert_NH_manual
   USE compute_theta_mod, ONLY : compute_theta => compute_theta_manual
   USE compute_geopot_mod, ONLY : compute_geopot => compute_geopot_manual
   USE compute_caldyn_kv_mod, ONLY : compute_caldyn_kv
-  USE compute_caldyn_Coriolis_mod, ONLY : compute_caldyn_Coriolis
+  USE compute_caldyn_Coriolis_mod, ONLY : compute_caldyn_Coriolis => compute_caldyn_Coriolis_manual
   USE compute_caldyn_slow_hydro_mod, ONLY : compute_caldyn_slow_hydro => compute_caldyn_slow_hydro_manual
   USE compute_caldyn_slow_NH_mod, ONLY : compute_caldyn_slow_NH
@@ -62,5 +62,5 @@
     REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
     REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:), Kv(:,:), hv(:,:)
-
+    REAL(rstd),POINTER :: berni(:,:), Ftheta(:,:)
 ! temporary shared variable
     REAL(rstd),POINTER  :: theta(:,:,:)  
@@ -136,6 +136,6 @@
        END IF
        u=f_u(ind)
-       Kv=f_Kv(ind) ! buffer for partial Bernoulli function
-       CALL compute_caldyn_fast(tau,theta,geopot, pk,Kv,du,u) ! computes du_fast and updates u
+       berni=f_berni(ind) ! buffer for partial Bernoulli function
+       CALL compute_caldyn_fast(tau,theta,geopot, pk,berni,du,u) ! computes du_fast and updates u
     ENDDO
     
@@ -186,5 +186,6 @@
           hv=f_hv(ind)
           Kv=f_Kv(ind)
-          CALL compute_caldyn_slow_hydro(.TRUE., u,mass,hv,Kv, Kv, hflux,du)
+          berni=f_berni(ind)
+          CALL compute_caldyn_slow_hydro(.TRUE., u,mass,hv,Kv, berni, hflux,du)
        ELSE
           W = f_W(ind)
@@ -197,5 +198,6 @@
           CALL compute_caldyn_slow_NH(u,mass,geopot,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
        END IF
-       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
+       Ftheta=f_Ftheta(ind)
+       CALL compute_caldyn_Coriolis(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
        IF(caldyn_eta==eta_mass) THEN
           wflux=f_wflux(ind)
Index: codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90	(revision 927)
+++ 	(revision )
@@ -1,203 +1,0 @@
-MODULE caldyn_kernels_base_mod
-  USE icosa
-  USE transfert_mod
-  USE disvert_mod
-  USE caldyn_vars_mod
-  USE omp_para
-  USE trace
-  IMPLICIT NONE
-  PRIVATE
-  SAVE
-
-  PUBLIC :: compute_caldyn_vert, compute_caldyn_vert_nh
-
-CONTAINS
-
-
-  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
-    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn)
-    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
-    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
-    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
-    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
-    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
-
-    ! temporary variable    
-    INTEGER :: i,j,ij,l,iq
-    REAL(rstd) :: p_ik, exner_ik, dF_deta, dFu_deta
-    INTEGER    :: ij_omp_begin, ij_omp_end
-
-    CALL trace_start("compute_caldyn_vert")
-
-!$OMP BARRIER
-
-    CALL distrib_level(ij_begin,ij_end, ij_omp_begin,ij_omp_end)
-
-    IF(dysl_caldyn_vert) THEN
-#define mass_bl(ij,l) bp(l)
-#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
-
-!!! cumulate mass flux convergence from top to bottom
-    DO  l = llm-1, 1, -1
-       !DIR$ SIMD
-       DO ij=ij_omp_begin,ij_omp_end         
-          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
-       ENDDO
-    ENDDO
-    !  ENDIF
-
-    !$OMP BARRIER
-    ! FLUSH on convm
-    ! compute dmass_col
-    IF (is_omp_first_level) THEN
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          ! dps/dt = -int(div flux)dz
-          dps(ij) = convm(ij,1)
-       ENDDO
-    ENDIF
-
-!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
-    DO l=ll_beginp1,ll_end
-       !    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
-          ! => w>0 for upward transport
-          wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
-       ENDDO
-    ENDDO
-
-    !--> flush wflux  
-    !$OMP BARRIER
-
-    DO iq=1,nqdyn
-       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 * &
-                  ( wflux(ij,l+1) * (theta(ij,l,iq) + theta(ij,l+1,iq)))
-          END DO
-       END DO
-       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 * &
-                  ( wflux(ij,l) * (theta(ij,l-1,iq) + theta(ij,l,iq) ) )
-          END DO
-       END DO
-    END DO
-
-    ! Compute vertical transport
-    DO l=ll_beginp1,ll_end
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          wwuu(ij+u_right,l) = 0.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) = 0.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) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
-       ENDDO
-    ENDDO
-
-    !--> flush wwuu  
-    !$OMP BARRIER
-
-    ! Add vertical transport 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)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
-          du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
-          du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
-       ENDDO
-    ENDDO
-
-    END IF ! dysl
-
-    CALL trace_end("compute_caldyn_vert")
-
-  END SUBROUTINE compute_caldyn_vert
-
-  SUBROUTINE compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
-    REAL(rstd),INTENT(IN) :: mass(iim*jjm,llm)
-    REAL(rstd),INTENT(IN) :: geopot(iim*jjm,llm+1)
-    REAL(rstd),INTENT(IN) :: W(iim*jjm,llm+1)
-    REAL(rstd),INTENT(IN) :: wflux(iim*jjm,llm+1)
-    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: dPhi(iim*jjm,llm+1)
-    REAL(rstd),INTENT(INOUT) :: dW(iim*jjm,llm+1)
-    REAL(rstd) :: W_etadot(iim*jjm,llm) ! vertical flux of vertical momentum
-    ! local arrays
-    REAL(rstd) :: eta_dot(iim*jjm, llm) ! eta_dot in full layers
-    REAL(rstd) :: wcov(iim*jjm,llm) ! covariant vertical momentum in full layers
-    ! indices and temporary values
-    INTEGER    :: ij, l
-    REAL(rstd) :: wflux_ij, w_ij
-
-    CALL trace_start("compute_caldyn_vert_nh")
-
-    IF(dysl) THEN
-!$OMP BARRIER
-#include "../kernels_hex/caldyn_vert_NH.k90"
-!$OMP BARRIER
-    ELSE
-#define ETA_DOT(ij) eta_dot(ij,1)
-#define WCOV(ij) wcov(ij,1)
-
-    DO l=ll_begin,ll_end
-       ! compute the local arrays
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          wflux_ij = .5*(wflux(ij,l)+wflux(ij,l+1))
-          w_ij = .5*(W(ij,l)+W(ij,l+1))/mass(ij,l)
-          W_etadot(ij,l) = wflux_ij*w_ij
-          ETA_DOT(ij) = wflux_ij / mass(ij,l)
-          WCOV(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l))
-       ENDDO
-       ! add NH term to du
-      !DIR$ SIMD
-      DO ij=ij_begin,ij_end
-          du(ij+u_right,l) = du(ij+u_right,l) &
-               - .5*(WCOV(ij+t_right)+WCOV(ij)) &
-               *ne_right*(ETA_DOT(ij+t_right)-ETA_DOT(ij))
-          du(ij+u_lup,l) = du(ij+u_lup,l) &
-               - .5*(WCOV(ij+t_lup)+WCOV(ij)) &
-               *ne_lup*(ETA_DOT(ij+t_lup)-ETA_DOT(ij))
-          du(ij+u_ldown,l) = du(ij+u_ldown,l) &
-               - .5*(WCOV(ij+t_ldown)+WCOV(ij)) &
-               *ne_ldown*(ETA_DOT(ij+t_ldown)-ETA_DOT(ij))
-       END DO
-    ENDDO
-    ! add NH terms to dW, dPhi
-    ! FIXME : TODO top and bottom
-    DO l=ll_beginp1,ll_end ! inner interfaces only
-       !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
-    DO l=ll_begin,ll_end
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end
-          dW(ij,l+1) = dW(ij,l+1) + W_etadot(ij,l) ! update inner+top interfaces
-          dW(ij,l)   = dW(ij,l)   - W_etadot(ij,l) ! update bottom+inner interfaces
-       END DO
-    END DO
-
-#undef ETA_DOT
-#undef WCOV
-
-    END IF ! dysl
-    CALL trace_end("compute_caldyn_vert_nh")
-
-  END SUBROUTINE compute_caldyn_vert_NH
-
-END MODULE caldyn_kernels_base_mod
Index: codes/icosagcm/devel/src/dynamics/caldyn_vars.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_vars.f90	(revision 927)
+++ codes/icosagcm/devel/src/dynamics/caldyn_vars.f90	(revision 928)
@@ -22,4 +22,5 @@
   TYPE(t_field),POINTER, PUBLIC :: f_qu(:), f_qv(:), f_Kv(:), f_hv(:), &
                            f_pk(:),f_wwuu(:),f_planetvel(:), &
+                           f_berni(:), f_Ftheta(:), &
                            f_Fel(:), f_gradPhi2(:), f_wil(:), f_Wetadot(:)
 
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 927)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 928)
@@ -1,4 +1,9 @@
 MODULE compute_caldyn_Coriolis_mod
-  USE grid_param, ONLY : llm
+  USE prec, ONLY : rstd
+  USE grid_param
+  USE earth_const
+  USE disvert_mod
+  USE omp_para
+  USE trace
   IMPLICIT NONE
   PRIVATE
@@ -6,5 +11,6 @@
 #include "../unstructured/unstructured.h90"
 
-  PUBLIC :: compute_caldyn_Coriolis
+  PUBLIC :: compute_caldyn_Coriolis_manual, &
+       compute_caldyn_Coriolis_unst, compute_caldyn_Coriolis_hex
 
 CONTAINS
@@ -55,13 +61,12 @@
 #endif END_DYSL
 
-  SUBROUTINE compute_coriolis_unst(hflux,theta,qu,Ftheta, convm,dtheta_rhodz,du)
+  SUBROUTINE compute_caldyn_coriolis_unst(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
     USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
-    USE grid_param, ONLY : nqdyn
     USE data_unstructured_mod, ONLY : id_coriolis, primal_num, dual_num, edge_num, &
           left, right,primal_deg,primal_edge,primal_ne,trisk_deg,wee,trisk,Ai, &
           enter_trace, exit_trace
-    FIELD_U     :: hflux, Ftheta, qu, du
-    FIELD_MASS  :: convm
-    FIELD_THETA :: theta, dtheta_rhodz
+    FIELD_U     :: hflux, Ftheta, qu, du ! IN, BUF, IN, INOUT
+    FIELD_MASS  :: convm                 ! BUF
+    FIELD_THETA :: theta, dtheta_rhodz   ! IN, OUT
     DECLARE_INDICES
     DECLARE_EDGES
@@ -70,30 +75,42 @@
 #include "../kernels_unst/coriolis.k90"
     STOP_TRACE
-  END SUBROUTINE compute_coriolis_unst
-
-  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
+  END SUBROUTINE compute_caldyn_coriolis_unst
+
+  SUBROUTINE compute_caldyn_Coriolis_hex(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
     USE icosa
-    USE trace
     USE caldyn_vars_mod
-    USE omp_para, ONLY : ll_begin, ll_end
-    REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) ! active scalars
-    REAL(rstd),INTENT(IN)  :: qu(3*iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
-    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(IN)    :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars
+    REAL(rstd),INTENT(IN)    :: qu(3*iim*jjm,llm)
+    REAL(rstd), INTENT(OUT)  :: Ftheta(3*iim*jjm,llm)  ! potential temperature flux
+    REAL(rstd),INTENT(OUT)   :: convm(iim*jjm,llm)  ! mass flux convergence
+    REAL(rstd),INTENT(OUT)   :: dtheta_rhodz(iim*jjm,llm,nqdyn)
     REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
-    
-    REAL(rstd) :: Ftheta(3*iim*jjm,llm)  ! potential temperature flux
     REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
     INTEGER :: ij,iq,l,kdown
-
     CALL trace_start("compute_caldyn_Coriolis")
 
-    IF(dysl_caldyn_coriolis) THEN
-
 #include "../kernels_hex/coriolis.k90"
 
-    ELSE
-#define FTHETA(ij) Ftheta(ij,1)
+    CALL trace_end("compute_caldyn_Coriolis")
+  END SUBROUTINE compute_caldyn_Coriolis_hex
+
+  SUBROUTINE compute_caldyn_Coriolis_manual(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
+    USE icosa
+    USE caldyn_vars_mod
+    REAL(rstd),INTENT(IN)    :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars
+    REAL(rstd),INTENT(IN)    :: qu(3*iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: Ftheta(3*iim*jjm,llm)  ! ignored in favor of local buffer
+    REAL(rstd),INTENT(OUT)   :: convm(iim*jjm,llm)  ! mass flux convergence
+    REAL(rstd),INTENT(OUT)   :: dtheta_rhodz(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
+    REAL(rstd)  :: buf_Ftheta(3*iim*jjm)  ! local buffer for potential temperature flux
+    REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
+    INTEGER :: ij,iq,l,kdown
+
+    CALL trace_start("compute_caldyn_Coriolis")
+
+#define FTHETA(ij) buf_Ftheta(ij)
 
     DO l=ll_begin, ll_end
@@ -275,9 +292,7 @@
 #undef FTHETA
 
-    END IF ! dysl
-
     CALL trace_end("compute_caldyn_Coriolis")
 
-  END SUBROUTINE compute_caldyn_Coriolis
+  END SUBROUTINE compute_caldyn_Coriolis_manual
   
 END MODULE compute_caldyn_Coriolis_mod
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90	(revision 927)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90	(revision 928)
@@ -121,5 +121,4 @@
     USE omp_para, ONLY : ll_begin, ll_end,ll_beginp1,ll_endp1
     USE disvert_mod, ONLY : ptop
-    USE caldyn_kernels_hevi_mod
     USE compute_NH_geopot_mod
     REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_vert.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_vert.F90	(revision 928)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_vert.F90	(revision 928)
@@ -0,0 +1,228 @@
+MODULE compute_caldyn_vert_mod
+  USE prec, ONLY : rstd
+  USE caldyn_vars_mod
+  USE grid_param
+  USE disvert_mod
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+  PRIVATE
+  SAVE
+
+  PUBLIC :: compute_caldyn_vert_manual, &
+       compute_caldyn_vert_hex
+CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(caldyn_wflux)
+  SEQUENCE_C0
+    BODY('llm-1,1,-1')
+      ! cumulate mass flux convergence from top to bottom
+      convm(CELL) = convm(CELL) + convm(UP(CELL))
+    END_BLOCK
+    EPILOGUE(1)
+      dmass_col(HIDX(CELL)) = convm(CELL)
+    END_BLOCK
+    BODY('2,llm')
+      ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
+      wflux(CELL) = mass_bl(CELL) * dmass_col(HIDX(CELL)) - convm(CELL) 
+    END_BLOCK
+  END_BLOCK
+  ! make sure wflux is up to date
+  BARRIER
+END_BLOCK
+
+KERNEL(caldyn_dmass)
+  FORALL_CELLS()
+    ON_PRIMAL
+      convm(CELL) = mass_dbk(CELL) * dmass_col(HIDX(CELL))
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+KERNEL(caldyn_vert)
+  DO iq=1,nqdyn
+    FORALL_CELLS('2', 'llm')
+      ON_PRIMAL
+        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) + 0.5*(theta(CELL,iq)+theta(DOWN(CELL),iq))*wflux(CELL)  
+      END_BLOCK
+    END_BLOCK
+    FORALL_CELLS('1', 'llm-1')
+      ON_PRIMAL
+        dtheta_rhodz(CELL,iq) = dtheta_rhodz(CELL,iq) - 0.5*(theta(CELL,iq)+theta(UP(CELL),iq))*wflux(UP(CELL))  
+      END_BLOCK
+    END_BLOCK
+  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)
+    FORALL_CELLS('2','llm')
+      ON_EDGES
+        wwuu(EDGE) = .25*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)+u(DOWN(EDGE))) ! Fu
+      END_BLOCK
+    END_BLOCK
+    ! make sure wwuu is up to date
+    BARRIER
+
+    FORALL_CELLS()
+      ON_EDGES
+        dFu_deta = wwuu(UP(EDGE))-wwuu(EDGE)                                          ! d/deta (F*u)
+	dF_deta  = .5*(wflux(UP(CELL1))+wflux(UP(CELL2))-(wflux(CELL1)+wflux(CELL2))) ! d/deta(F)
+        du(EDGE) = du(EDGE) - (dFu_deta-u(EDGE)*dF_deta) / (.5*(rhodz(CELL1)+rhodz(CELL2)))  ! (F/m)du/deta
+      END_BLOCK
+    END_BLOCK
+  ELSE  
+    FORALL_CELLS('2','llm')
+      ON_EDGES
+        wwuu(EDGE) = .5*(wflux(CELL1)+wflux(CELL2))*(u(EDGE)-u(DOWN(EDGE)))
+      END_BLOCK
+    END_BLOCK
+
+    ! make sure wwuu is up to date
+    BARRIER
+
+    FORALL_CELLS()
+      ON_EDGES
+        du(EDGE) = du(EDGE) - (wwuu(EDGE)+wwuu(UP(EDGE))) / (rhodz(CELL1)+rhodz(CELL2))
+      END_BLOCK
+    END_BLOCK
+  END IF
+END_BLOCK
+
+#endif END_DYSL
+
+  SUBROUTINE compute_caldyn_vert_hex(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
+    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
+    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
+    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
+
+    ! temporary variable    
+    INTEGER :: i,j,ij,l,iq
+    REAL(rstd) :: p_ik, exner_ik, dF_deta, dFu_deta
+    INTEGER    :: ij_omp_begin, ij_omp_end
+
+    CALL trace_start("compute_caldyn_vert")
+
+!$OMP BARRIER
+
+    CALL distrib_level(ij_begin,ij_end, ij_omp_begin,ij_omp_end)
+
+#define mass_bl(ij,l) bp(l)
+#define dmass_col(ij) dps(ij)
+#include "../kernels_hex/caldyn_wflux.k90"
+#include "../kernels_hex/caldyn_vert.k90"
+#undef mass_bl
+#undef dmass_col
+
+    CALL trace_end("compute_caldyn_vert")
+
+  END SUBROUTINE compute_caldyn_vert_hex
+
+  SUBROUTINE compute_caldyn_vert_manual(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
+    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
+    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
+    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
+
+    ! temporary variable    
+    INTEGER :: i,j,ij,l,iq
+    REAL(rstd) :: p_ik, exner_ik, dF_deta, dFu_deta
+    INTEGER    :: ij_omp_begin, ij_omp_end
+
+    CALL trace_start("compute_caldyn_vert")
+
+!$OMP BARRIER
+
+    CALL distrib_level(ij_begin,ij_end, ij_omp_begin,ij_omp_end)
+
+!!! cumulate mass flux convergence from top to bottom
+    DO  l = llm-1, 1, -1
+       !DIR$ SIMD
+       DO ij=ij_omp_begin,ij_omp_end         
+          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
+       ENDDO
+    ENDDO
+    !  ENDIF
+
+    !$OMP BARRIER
+    ! FLUSH on convm
+    ! compute dmass_col
+    IF (is_omp_first_level) THEN
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end         
+          ! dps/dt = -int(div flux)dz
+          dps(ij) = convm(ij,1)
+       ENDDO
+    ENDIF
+!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
+    DO l=ll_beginp1,ll_end
+       !    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end         
+          ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
+          ! => w>0 for upward transport
+          wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
+       ENDDO
+    ENDDO
+
+    !--> flush wflux  
+    !$OMP BARRIER
+
+    DO iq=1,nqdyn
+       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 * &
+                  ( wflux(ij,l+1) * (theta(ij,l,iq) + theta(ij,l+1,iq)))
+          END DO
+       END DO
+       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 * &
+                  ( wflux(ij,l) * (theta(ij,l-1,iq) + theta(ij,l,iq) ) )
+          END DO
+
+       END DO
+    END DO
+
+    ! Compute vertical transport
+    DO l=ll_beginp1,ll_end
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end         
+          wwuu(ij+u_right,l) = 0.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) = 0.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) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
+       ENDDO
+    ENDDO
+
+    !--> flush wwuu  
+    !$OMP BARRIER
+
+    ! Add vertical transport 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)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
+          du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
+          du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
+       ENDDO
+    ENDDO
+
+    CALL trace_end("compute_caldyn_vert")
+
+  END SUBROUTINE compute_caldyn_vert_manual
+
+END MODULE compute_caldyn_vert_mod
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_vert_NH.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_vert_NH.F90	(revision 928)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_vert_NH.F90	(revision 928)
@@ -0,0 +1,109 @@
+MODULE compute_caldyn_vert_NH_mod
+  USE prec, ONLY : rstd
+  USE grid_param
+  USE disvert_mod
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+  PRIVATE
+  SAVE
+
+  PUBLIC :: compute_caldyn_vert_nh_manual, &
+       compute_caldyn_vert_nh_hex
+
+CONTAINS
+
+  SUBROUTINE compute_caldyn_vert_NH_hex(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
+    REAL(rstd),INTENT(IN) :: mass(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN) :: W(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN) :: wflux(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: dPhi(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: dW(iim*jjm,llm+1)
+    REAL(rstd) :: W_etadot(iim*jjm,llm) ! vertical flux of vertical momentum
+    ! local arrays
+    REAL(rstd) :: eta_dot(iim*jjm, llm) ! eta_dot in full layers
+    REAL(rstd) :: wcov(iim*jjm,llm) ! covariant vertical momentum in full layers
+    ! indices and temporary values
+    INTEGER    :: ij, l
+    REAL(rstd) :: wflux_ij, w_ij
+
+    CALL trace_start("compute_caldyn_vert_nh")
+
+!$OMP BARRIER
+#include "../kernels_hex/caldyn_vert_NH.k90"
+!$OMP BARRIER
+    CALL trace_end("compute_caldyn_vert_nh")
+  END SUBROUTINE compute_caldyn_vert_NH_hex
+
+  SUBROUTINE compute_caldyn_vert_NH_manual(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
+    REAL(rstd),INTENT(IN) :: mass(iim*jjm,llm)
+    REAL(rstd),INTENT(IN) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN) :: W(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN) :: wflux(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: dPhi(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: dW(iim*jjm,llm+1)
+    REAL(rstd) :: W_etadot(iim*jjm,llm) ! vertical flux of vertical momentum
+    ! local arrays
+    REAL(rstd) :: eta_dot(iim*jjm, llm) ! eta_dot in full layers
+    REAL(rstd) :: wcov(iim*jjm,llm) ! covariant vertical momentum in full layers
+    ! indices and temporary values
+    INTEGER    :: ij, l
+    REAL(rstd) :: wflux_ij, w_ij
+
+    CALL trace_start("compute_caldyn_vert_nh")
+
+#define ETA_DOT(ij) eta_dot(ij,1)
+#define WCOV(ij) wcov(ij,1)
+
+    DO l=ll_begin,ll_end
+       ! compute the local arrays
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          wflux_ij = .5*(wflux(ij,l)+wflux(ij,l+1))
+          w_ij = .5*(W(ij,l)+W(ij,l+1))/mass(ij,l)
+          W_etadot(ij,l) = wflux_ij*w_ij
+          ETA_DOT(ij) = wflux_ij / mass(ij,l)
+          WCOV(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l))
+       ENDDO
+       ! add NH term to du
+      !DIR$ SIMD
+      DO ij=ij_begin,ij_end
+          du(ij+u_right,l) = du(ij+u_right,l) &
+               - .5*(WCOV(ij+t_right)+WCOV(ij)) &
+               *ne_right*(ETA_DOT(ij+t_right)-ETA_DOT(ij))
+          du(ij+u_lup,l) = du(ij+u_lup,l) &
+               - .5*(WCOV(ij+t_lup)+WCOV(ij)) &
+               *ne_lup*(ETA_DOT(ij+t_lup)-ETA_DOT(ij))
+          du(ij+u_ldown,l) = du(ij+u_ldown,l) &
+               - .5*(WCOV(ij+t_ldown)+WCOV(ij)) &
+               *ne_ldown*(ETA_DOT(ij+t_ldown)-ETA_DOT(ij))
+       END DO
+    ENDDO
+    ! add NH terms to dW, dPhi
+    ! FIXME : TODO top and bottom
+    DO l=ll_beginp1,ll_end ! inner interfaces only
+       !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
+    DO l=ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end
+          dW(ij,l+1) = dW(ij,l+1) + W_etadot(ij,l) ! update inner+top interfaces
+          dW(ij,l)   = dW(ij,l)   - W_etadot(ij,l) ! update bottom+inner interfaces
+       END DO
+    END DO
+
+#undef ETA_DOT
+#undef WCOV
+
+    CALL trace_end("compute_caldyn_vert_nh")
+
+  END SUBROUTINE compute_caldyn_vert_NH_manual
+
+END MODULE compute_caldyn_vert_NH_mod
