Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 361)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 362)
@@ -2,10 +2,11 @@
   USE icosa
   USE transfert_mod
+  USE caldyn_kernels_base_mod
   USE caldyn_kernels_mod
   IMPLICIT NONE
   PRIVATE
 
-  PUBLIC init_caldyn, caldyn_BC, caldyn, req_ps, req_mass
-
+  PUBLIC init_caldyn, caldyn_BC, caldyn
+  
 CONTAINS
   
@@ -13,5 +14,4 @@
     USE icosa
     USE observable_mod
-!    USE exner_mod
     USE mpipara
     USE omp_para
Index: codes/icosagcm/trunk/src/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 361)
+++ codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 362)
@@ -2,10 +2,9 @@
   USE icosa
   USE transfert_mod
-  USE caldyn_kernels_mod
+  USE caldyn_kernels_base_mod
+  USE caldyn_kernels_hevi_mod
   USE caldyn_gcm_mod
-
   IMPLICIT NONE
   PRIVATE
-
   PUBLIC caldyn_hevi
 
@@ -75,15 +74,9 @@
     IF(caldyn_eta==eta_mass) THEN
        CALL send_message(f_ps,req_ps) ! COM00
+       CALL wait_message(req_ps) ! COM00
     ELSE
        CALL send_message(f_mass,req_mass) ! COM00
-    END IF
-
-    ! Block below moved from caldyn_pvort when splitting it into caldyn_theta and caldyn_pvort_only
-    IF(caldyn_eta==eta_mass) THEN
-       CALL wait_message(req_ps) ! COM00
-    ELSE
        CALL wait_message(req_mass) ! COM00
     END IF
-
     CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
     CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
@@ -100,5 +93,5 @@
        theta = f_theta(ind)
        pk = f_pk(ind)
-       geopot = f_geopot(ind)  
+       geopot = f_geopot(ind)
        CALL compute_theta(ps,theta_rhodz, mass,theta)
        CALL compute_geopot(ps,mass,theta, pk,geopot)
Index: codes/icosagcm/trunk/src/caldyn_kernels.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels.f90	(revision 361)
+++ codes/icosagcm/trunk/src/caldyn_kernels.f90	(revision 362)
@@ -2,21 +2,10 @@
   USE icosa
   USE transfert_mod
+  USE caldyn_kernels_base_mod
   IMPLICIT NONE
-
-  INTEGER, PARAMETER :: energy=1, enstrophy=2
-  TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:)
-  REAL(rstd),SAVE,POINTER :: out_u(:,:), p(:,:), qu(:,:)
-  !$OMP THREADPRIVATE(out_u, p, qu)
-
-  ! temporary shared variable for caldyn
-  TYPE(t_field),POINTER :: f_pk(:)
-  TYPE(t_field),POINTER :: f_wwuu(:)
-  TYPE(t_field),POINTER :: f_planetvel(:)
-
-  INTEGER :: caldyn_conserv
-  !$OMP THREADPRIVATE(caldyn_conserv) 
-
-  TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu
-
+  PRIVATE
+
+  PUBLIC :: compute_planetvel, compute_pvort, compute_geopot, &
+       compute_caldyn_horiz, compute_caldyn_vert
 CONTAINS
 
@@ -40,6 +29,4 @@
     CALL compute_wind2D_perp_from_lonlat_compound(ulon, ulat, planetvel)
   END SUBROUTINE compute_planetvel
-
-  !********************** compute_pvort = compute_theta + compute_pvort_only ******************
 
   SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
@@ -60,6 +47,4 @@
     INTEGER :: i,j,ij,l
     REAL(rstd) :: etav,hv, m
-    !  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity  
-
     CALL trace_start("compute_pvort")  
 
@@ -130,190 +115,4 @@
     CALL trace_end("compute_pvort")
   END SUBROUTINE compute_pvort
-
-  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
-    USE icosa
-    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
-    USE exner_mod
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
-    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
-    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
-    INTEGER :: ij,l
-    REAL(rstd) :: m
-    CALL trace_start("compute_theta")  
-
-    IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
-       DO l = ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext
-             m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
-             rhodz(ij,l) = m
-             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-          ENDDO
-       ENDDO
-    ELSE ! Compute only theta
-       DO l = ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext
-             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-          ENDDO
-       ENDDO
-    END IF
-
-    CALL trace_end("compute_theta")
-  END SUBROUTINE compute_theta
-
-  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
-    USE icosa
-    USE exner_mod
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
-    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
-
-    INTEGER :: ij,l
-    REAL(rstd) :: etav,hv
-
-    CALL trace_start("compute_pvort_only")  
-!!! Compute shallow-water potential vorticity
-    DO l = ll_begin,ll_end
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
-               + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
-               - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
-
-          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
-               + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
-               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
-
-          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
-
-          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
-               + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
-               - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
-
-          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
-               + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
-               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
-
-          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
-
-       ENDDO
-
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end
-          qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))  
-          qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))  
-          qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))  
-       END DO
-
-    ENDDO
-
-    CALL trace_end("compute_pvort_only")
-
-  END SUBROUTINE compute_pvort_only
-
-  !**************************** Geopotential *****************************
-
-  SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot) 
-    USE icosa
-    USE disvert_mod
-    USE exner_mod
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
-    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
-    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature
-    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)       ! Exner function
-    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
-
-    INTEGER :: i,j,ij,l
-    REAL(rstd) :: p_ik, exner_ik
-    INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext
-
-
-    CALL trace_start("compute_geopot")
-
-    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
-    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1
-    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
-
-    IF(caldyn_eta==eta_mass) THEN
-
-!!! Compute exner function and geopotential
-       DO l = 1,llm
-          !DIR$ SIMD
-          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
-             !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
-             exner_ik = cpp * (p_ik/preff) ** kappa
-             pk(ij,l) = exner_ik
-             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
-          ENDDO
-       ENDDO
-       !       ENDIF
-    ELSE 
-       ! We are using a Lagrangian vertical coordinate
-       ! Pressure must be computed first top-down (temporarily stored in pk)
-       ! Then Exner pressure and geopotential are computed bottom-up
-       ! Notice that the computation below should work also when caldyn_eta=eta_mass          
-
-       IF(boussinesq) THEN ! compute only geopotential : pressure pk will be computed in compute_caldyn_horiz
-          ! specific volume 1 = dphi/g/rhodz
-          !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
-          DO l = 1,llm
-             !DIR$ SIMD
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-                geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
-             ENDDO
-          ENDDO
-       ELSE ! non-Boussinesq, compute geopotential and Exner pressure
-          ! uppermost layer
-
-          !DIR$ SIMD
-          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
-          END DO
-          ! other layers
-          DO l = llm-1, 1, -1
-             !DIR$ SIMD
-             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
-          ! surface pressure (for diagnostics)
-          DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
-          END DO
-
-          ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-          DO l = 1,llm
-             !DIR$ SIMD
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-                p_ik = pk(ij,l)
-                exner_ik = cpp * (p_ik/preff) ** kappa
-                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
-                pk(ij,l) = exner_ik
-             ENDDO
-          ENDDO
-       END IF
-
-    END IF
-
-    !ym flush geopot
-    !$OMP BARRIER
-
-    CALL trace_end("compute_geopot")
-
-  END SUBROUTINE compute_geopot
 
   !************************* caldyn_horiz = caldyn_fast + caldyn_slow **********************
@@ -490,24 +289,4 @@
 !!! Compute bernouilli term = Kinetic Energy + geopotential
     IF(boussinesq) THEN
-       ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) 
-       ! uppermost layer
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext         
-          pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
-       END DO
-       ! other layers
-       DO l = llm-1, 1, -1
-          !          !$OMP DO SCHEDULE(STATIC) 
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext         
-             pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
-          END DO
-       END DO
-       ! surface pressure (for diagnostics) FIXME
-       ! DO ij=ij_begin_ext,ij_end_ext         
-       !   ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1)
-       ! END DO
-       ! now pk contains the Lagrange multiplier (pressure)
-
        DO l=ll_begin,ll_end
           !DIR$ SIMD
@@ -572,440 +351,3 @@
   END SUBROUTINE compute_caldyn_horiz
 
-  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot, du)
-    USE icosa
-    USE disvert_mod
-    USE exner_mod
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
-    REAL(rstd), INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs
-    REAL(rstd),INTENT(INOUT)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
-    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
-    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
-    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
-    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
-    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
-
-    INTEGER :: i,j,ij,l
-
-    CALL trace_start("compute_caldyn_fast")
-
-    !    CALL wait_message(req_theta_rhodz) 
-    
-    ! Compute bernouilli term
-    IF(boussinesq) THEN
-       ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) 
-       ! uppermost layer
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext         
-          pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
-       END DO
-       ! other layers
-       DO l = llm-1, 1, -1
-          !          !$OMP DO SCHEDULE(STATIC) 
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext         
-             pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
-          END DO
-       END DO
-       ! now pk contains the Lagrange multiplier (pressure)
-       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))
-          ENDDO
-       ENDDO
-
-    ELSE ! compressible
-
-       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))
-          ENDDO
-       ENDDO
-
-    END IF ! Boussinesq/compressible
-
-!!! u:=u+tau*du, du = gradients of Bernoulli and Exner functions
-    DO l=ll_begin,ll_end
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-
-          du(ij+u_right,l) = 1/de(ij+u_right) * (             &
-               0.5*(theta(ij,l)+theta(ij+t_right,l))                &
-               *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
-               + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
-          u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
-
-
-          du(ij+u_lup,l) = 1/de(ij+u_lup) * (                  &
-               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
-               *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
-               + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
-          u(ij+u_lup,l) = u(ij+u_lup,l) + tau*du(ij+u_lup,l)
-
-          du(ij+u_ldown,l) = 1/de(ij+u_ldown) * (             &
-               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
-               *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
-               + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
-          u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
-
-       ENDDO
-    ENDDO
-
-    CALL trace_end("compute_caldyn_fast")
-
-  END SUBROUTINE compute_caldyn_fast
-
-  SUBROUTINE compute_caldyn_slow(u,rhodz,qu,theta, hflux,convm, dtheta_rhodz, du)
-    USE icosa
-    USE disvert_mod
-    USE exner_mod
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
-    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
-    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
-
-    REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s
-    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
-    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
-
-    REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi)
-    REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
-    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
-    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
-
-    INTEGER :: i,j,ij,l
-    REAL(rstd) :: ww,uu
-
-    CALL trace_start("compute_caldyn_slow")
-
-    !    CALL wait_message(req_theta_rhodz) 
-
-    DO l = ll_begin, ll_end
-!!!  Compute mass and theta fluxes
-       IF (caldyn_conserv==energy) CALL test_message(req_qu) 
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext         
-          hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
-          hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
-          hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
-
-          Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
-          Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
-          Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
-       ENDDO
-
-!!! compute horizontal divergence of fluxes
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
-          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
-               ne_rup*hflux(ij+u_rup,l)       +  &  
-               ne_lup*hflux(ij+u_lup,l)       +  &  
-               ne_left*hflux(ij+u_left,l)     +  &  
-               ne_ldown*hflux(ij+u_ldown,l)   +  &  
-               ne_rdown*hflux(ij+u_rdown,l))    
-
-          ! signe ? attention d (rho theta dz)
-          ! dtheta_rhodz =  -div(flux.theta)
-          dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
-               ne_rup*Ftheta(ij+u_rup,l)        +  &  
-               ne_lup*Ftheta(ij+u_lup,l)        +  &  
-               ne_left*Ftheta(ij+u_left,l)      +  &  
-               ne_ldown*Ftheta(ij+u_ldown,l)    +  &  
-               ne_rdown*Ftheta(ij+u_rdown,l))
-       ENDDO
-
-    END DO
-
-!!! Compute potential vorticity (Coriolis) contribution to du
-
-    SELECT CASE(caldyn_conserv)
-    CASE(energy) ! energy-conserving TRiSK
-
-       CALL wait_message(req_qu)
-
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-
-             uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
-                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
-                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
-                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
-                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
-                  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))+         &
-                  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))+         &
-                  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))+         &
-                  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))+             &
-                  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) = .5*uu/de(ij+u_right)
-
-             uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
-                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
-                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
-                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
-                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
-                  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)) +          &
-                  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)) +              &
-                  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)) +              &
-                  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)) +            &
-                  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) = .5*uu/de(ij+u_lup)
-
-
-             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
-                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
-                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
-                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
-                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
-                  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)) +          &
-                  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)) +        &
-                  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)) +      &
-                  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)) +      &
-                  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) = .5*uu/de(ij+u_ldown)
-
-          ENDDO
-       ENDDO
-
-    CASE(enstrophy) ! enstrophy-conserving TRiSK
-
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-
-             uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
-                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
-                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
-                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
-                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
-                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
-                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
-                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
-                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
-                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
-             du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
-
-
-             uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
-                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
-                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
-                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
-                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
-                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
-                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
-                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
-                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
-                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
-             du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
-
-             uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
-                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
-                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
-                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
-                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
-                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
-                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
-                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
-                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
-                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
-             du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
-
-          ENDDO
-       ENDDO
-
-    CASE DEFAULT
-       STOP
-    END SELECT
-
-    ! Compute bernouilli term = Kinetic Energy
-    IF(boussinesq) THEN
-
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-
-             berni(ij,l) = &
-                  1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
-                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
-                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
-                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
-                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
-                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
-          ENDDO
-       ENDDO
-
-    ELSE ! compressible
-
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-
-             berni(ij,l) =  &
-                  + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
-                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
-                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
-                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
-                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
-                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
-          ENDDO
-       ENDDO
-
-    END IF ! Boussinesq/compressible
-
-!!! Add gradients of Bernoulli and Exner functions 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) + 1/de(ij+u_right)            &
-               * ( ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
-          du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup)                 &
-               * ( ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
-          du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown)            &
-               * ( ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
-       ENDDO
-    ENDDO
-
-    CALL trace_end("compute_caldyn_slow")
-
-  END SUBROUTINE compute_caldyn_slow
-
-  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
-    USE icosa
-    USE disvert_mod
-    USE exner_mod
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
-    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
-    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)
-    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
-
-    ! temporary variable    
-    INTEGER :: i,j,ij,l
-    REAL(rstd) :: p_ik, exner_ik
-    INTEGER    :: ij_omp_begin, ij_omp_end
-
-
-    CALL trace_start("compute_caldyn_vert")
-
-    CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin,ij_omp_end)
-    ij_omp_begin=ij_omp_begin+ij_begin-1
-    ij_omp_end=ij_omp_end+ij_begin-1
-
-    !    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")
-
-    !$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         
-          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
-       ENDDO
-    ENDDO
-    !  ENDIF
-
-    !$OMP BARRIER
-    ! FLUSH on convm
-!!!!!!!!!!!!!!!!!!!!!!!!!
-
-    ! compute dps
-    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) * g 
-       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 l=ll_begin,ll_endm1
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
-       ENDDO
-    ENDDO
-
-    DO l=ll_beginp1,ll_end
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
-       ENDDO
-    ENDDO
-
-
-    ! 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
-
-    !  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
-
-    CALL trace_end("compute_caldyn_vert")
-
-  END SUBROUTINE compute_caldyn_vert
-
 END MODULE caldyn_kernels_mod
Index: codes/icosagcm/trunk/src/caldyn_kernels_base.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 362)
+++ codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 362)
@@ -0,0 +1,267 @@
+MODULE caldyn_kernels_base_mod
+  USE icosa
+  USE transfert_mod
+  IMPLICIT NONE
+  PRIVATE
+
+  LOGICAL, PARAMETER, PUBLIC :: DEC = .TRUE.
+
+  INTEGER, PARAMETER,PUBLIC :: energy=1, enstrophy=2
+  TYPE(t_field),POINTER,PUBLIC :: f_out_u(:), f_qu(:), f_qv(:)
+  REAL(rstd),SAVE,POINTER :: out_u(:,:), p(:,:), qu(:,:)
+  !$OMP THREADPRIVATE(out_u, p, qu)
+
+  ! temporary shared variables for caldyn
+  TYPE(t_field),POINTER,PUBLIC :: f_pk(:),f_wwuu(:),f_planetvel(:)
+
+  INTEGER, PUBLIC :: caldyn_conserv
+  !$OMP THREADPRIVATE(caldyn_conserv) 
+
+  TYPE(t_message),PUBLIC :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu
+
+  PUBLIC :: compute_geopot, compute_caldyn_vert
+
+CONTAINS
+
+  !**************************** Geopotential *****************************
+
+  SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot) 
+    USE icosa
+    USE disvert_mod
+    USE exner_mod
+    USE trace
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
+    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature
+    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)       ! Exner function (compressible) /Lagrange multiplier (Boussinesq)
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
+
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: p_ik, exner_ik
+    INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext
+
+
+    CALL trace_start("compute_geopot")
+
+    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
+    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1
+    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
+
+    IF(caldyn_eta==eta_mass .AND. .NOT. DEC) THEN
+
+!!! Compute exner function and geopotential
+       DO l = 1,llm
+          !DIR$ SIMD
+          DO ij=ij_omp_begin_ext,ij_omp_end_ext
+             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
+             exner_ik = cpp * (p_ik/preff) ** kappa
+             pk(ij,l) = exner_ik
+             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
+          ENDDO
+       ENDDO
+    ELSE 
+       ! We are using a Lagrangian vertical coordinate
+       ! Pressure must be computed first top-down (temporarily stored in pk)
+       ! Then Exner pressure and geopotential are computed bottom-up
+       ! Notice that the computation below should work also when caldyn_eta=eta_mass          
+
+       IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier
+          ! specific volume 1 = dphi/g/rhodz
+          !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
+          DO l = 1,llm
+             !DIR$ SIMD
+             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
+                geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
+             ENDDO
+          ENDDO
+          ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) 
+          ! uppermost layer
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext         
+             pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
+          END DO
+          ! other layers
+          DO l = llm-1, 1, -1
+             !          !$OMP DO SCHEDULE(STATIC) 
+             !DIR$ SIMD
+             DO ij=ij_begin_ext,ij_end_ext         
+                pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
+             END DO
+          END DO
+          ! now pk contains the Lagrange multiplier (pressure)
+       ELSE ! non-Boussinesq, compute geopotential and Exner pressure
+          ! uppermost layer
+
+          !DIR$ SIMD
+          DO ij=ij_omp_begin_ext,ij_omp_end_ext
+             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
+          END DO
+          ! other layers
+          DO l = llm-1, 1, -1
+             !DIR$ SIMD
+             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
+          ! surface pressure (for diagnostics)
+          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
+          ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+          DO l = 1,llm
+             !DIR$ SIMD
+             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
+                p_ik = pk(ij,l)
+                exner_ik = cpp * (p_ik/preff) ** kappa
+                pk(ij,l) = exner_ik
+                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
+             ENDDO
+          ENDDO
+       END IF
+
+    END IF
+
+    !ym flush geopot
+    !$OMP BARRIER
+
+    CALL trace_end("compute_geopot")
+
+  END SUBROUTINE compute_geopot
+
+  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
+    USE icosa
+    USE disvert_mod
+    USE exner_mod
+    USE trace
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
+    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)
+    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
+
+    ! temporary variable    
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: p_ik, exner_ik
+    INTEGER    :: ij_omp_begin, ij_omp_end
+
+    CALL trace_start("compute_caldyn_vert")
+
+    CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin,ij_omp_end)
+    ij_omp_begin=ij_omp_begin+ij_begin-1
+    ij_omp_end=ij_omp_end+ij_begin-1
+
+    !    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")
+
+    !$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         
+          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
+       ENDDO
+    ENDDO
+    !  ENDIF
+
+    !$OMP BARRIER
+    ! FLUSH on convm
+!!!!!!!!!!!!!!!!!!!!!!!!!
+
+    ! compute dps
+    IF (is_omp_first_level) THEN
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end         
+          ! dps/dt = -int(div flux)dz
+          IF(DEC) THEN
+             dps(ij) = convm(ij,1)
+          ELSE
+             dps(ij) = convm(ij,1) * g 
+          END IF
+       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 l=ll_begin,ll_endm1
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end         
+          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
+       ENDDO
+    ENDDO
+
+    DO l=ll_beginp1,ll_end
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end         
+          dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
+       ENDDO
+    ENDDO
+
+
+    ! 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
+
+    !  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
+
+    CALL trace_end("compute_caldyn_vert")
+
+  END SUBROUTINE compute_caldyn_vert
+
+END MODULE caldyn_kernels_base_mod
Index: codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 362)
+++ codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 362)
@@ -0,0 +1,431 @@
+MODULE caldyn_kernels_hevi_mod
+  USE icosa
+  USE transfert_mod
+  USE caldyn_kernels_base_mod
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_fast, compute_caldyn_slow
+
+CONTAINS
+
+  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
+    USE icosa
+    USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop
+    USE exner_mod
+    USE trace
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
+    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
+    INTEGER :: ij,l
+    REAL(rstd) :: m
+    CALL trace_start("compute_theta")  
+
+    IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
+       DO l = ll_begin,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             IF(DEC) THEN  ! ps is actually Ms
+                m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l)
+             ELSE
+                m = mass_dak(l)+ps(ij)*mass_dbk(l)
+             END IF
+             rhodz(ij,l) = m/g
+             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+          ENDDO
+       ENDDO
+    ELSE ! Compute only theta
+       DO l = ll_begin,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+          ENDDO
+       ENDDO
+    END IF
+
+    CALL trace_end("compute_theta")
+  END SUBROUTINE compute_theta
+
+  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
+    USE icosa
+    USE exner_mod
+    USE trace
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
+
+    INTEGER :: ij,l
+    REAL(rstd) :: etav,hv,radius_m2
+
+    CALL trace_start("compute_pvort_only")  
+!!! Compute shallow-water potential vorticity
+    radius_m2=radius**(-2)
+    DO l = ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          IF(DEC) THEN
+             etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)         &
+                                   + ne_left * u(ij+t_rup+u_left,l)  &
+                                   - ne_lup  * u(ij+u_lup,l) )                               
+
+             hv =   Riv2(ij,vup)          * rhodz(ij,l)           &
+                  + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
+                  + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
+             qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
+
+             etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          &
+                                      + ne_right * u(ij+t_ldown+u_right,l)  &
+                                      - ne_rdown * u(ij+u_rdown,l) )
+             hv =   Riv2(ij,vdown)        * rhodz(ij,l)          &
+                  + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
+                  + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
+             qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
+          ELSE
+             etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)        * de(ij+u_rup)         &
+                                   + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
+                                   - ne_lup  * u(ij+u_lup,l)        * de(ij+u_lup) )                               
+
+             hv =   Riv2(ij,vup)          * rhodz(ij,l)           &
+                  + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
+                  + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
+             qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
+
+             etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
+                                      + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
+                                      - ne_rdown * u(ij+u_rdown,l)          * de(ij+u_rdown) )
+             hv =   Riv2(ij,vdown)        * rhodz(ij,l)          &
+                  + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
+                  + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
+             qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
+          END IF
+       ENDDO
+
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end
+          qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))  
+          qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))  
+          qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))  
+       END DO
+
+    ENDDO
+
+    CALL trace_end("compute_pvort_only")
+
+  END SUBROUTINE compute_pvort_only
+
+  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot, du)
+    USE icosa
+    USE disvert_mod
+    USE exner_mod
+    USE trace
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd), INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs
+    REAL(rstd),INTENT(INOUT)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
+    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
+    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
+    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
+    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
+
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: due_right, due_lup, due_ldown
+
+    CALL trace_start("compute_caldyn_fast")
+    
+    ! Compute bernouilli term
+    IF(boussinesq) THEN
+       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))
+          ENDDO
+       ENDDO
+
+    ELSE ! compressible
+
+       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))
+          ENDDO
+       ENDDO
+
+    END IF ! Boussinesq/compressible
+
+!!! u:=u+tau*du, du = gradients of Bernoulli and Exner functions
+    DO l=ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end         
+          due_right =  0.5*(theta(ij,l)+theta(ij+t_right,l))                  &
+                          *(ne_right*pk(ij,l)   +ne_left*pk(ij+t_right,l))    &
+                         +  ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l)
+          due_lup = 0.5*(theta(ij,l)+theta(ij+t_lup,l))                       &
+                       *(ne_lup*pk(ij,l)   +ne_rdown*pk(ij+t_lup,l))          &
+                      +  ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l)
+          due_ldown = 0.5*(theta(ij,l)+theta(ij+t_ldown,l))                   &
+                         *(ne_ldown*pk(ij,l)   +ne_rup*pk(ij+t_ldown,l))      &
+                        +  ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l)
+          IF(.NOT.DEC) THEN
+             due_right = due_right/de(ij+u_right)
+             due_lup   = due_lup/de(ij+u_lup)
+             due_ldown = due_ldown/de(ij+u_ldown)
+          END IF
+          du(ij+u_right,l) = due_right
+          du(ij+u_lup,l)   = due_lup
+          du(ij+u_ldown,l) = due_ldown
+          u(ij+u_right,l) = u(ij+u_right,l) + tau*due_right
+          u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*due_lup
+          u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*due_ldown
+       ENDDO
+    ENDDO
+
+    CALL trace_end("compute_caldyn_fast")
+
+  END SUBROUTINE compute_caldyn_fast
+
+  SUBROUTINE compute_caldyn_slow(u,rhodz,qu,theta, hflux,convm, dtheta_rhodz, du)
+    USE icosa
+    USE disvert_mod
+    USE exner_mod
+    USE trace
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
+    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
+
+    REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s
+    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
+    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
+
+    REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi)
+    REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
+    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
+    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
+
+    INTEGER :: ij,l
+    REAL(rstd) :: uu_right, uu_lup, uu_ldown
+
+    CALL trace_start("compute_caldyn_slow")
+
+    DO l = ll_begin, ll_end
+!!!  Compute mass and theta fluxes
+       IF (caldyn_conserv==energy) CALL test_message(req_qu) 
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
+          uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
+          uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
+          IF(DEC) THEN
+             uu_right= uu_right*le(ij+u_right)/de(ij+u_right)
+             uu_lup  = uu_lup  *le(ij+u_lup)/de(ij+u_lup)
+             uu_ldown= uu_ldown*le(ij+u_ldown)/de(ij+u_ldown)
+          ELSE
+             uu_right= uu_right*le(ij+u_right)
+             uu_lup  = uu_lup  *le(ij+u_lup)
+             uu_ldown= uu_ldown*le(ij+u_ldown)
+          END IF
+          hflux(ij+u_right,l)=uu_right
+          hflux(ij+u_lup,l)  =uu_lup
+          hflux(ij+u_ldown,l)=uu_ldown
+          Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*uu_right
+          Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*uu_lup
+          Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*uu_ldown
+       ENDDO
+
+!!! compute horizontal divergence of fluxes
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end         
+          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
+          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
+               ne_rup*hflux(ij+u_rup,l)       +  &  
+               ne_lup*hflux(ij+u_lup,l)       +  &  
+               ne_left*hflux(ij+u_left,l)     +  &  
+               ne_ldown*hflux(ij+u_ldown,l)   +  &  
+               ne_rdown*hflux(ij+u_rdown,l))    
+
+          ! signe ? attention d (rho theta dz)
+          ! dtheta_rhodz =  -div(flux.theta)
+          dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
+               ne_rup*Ftheta(ij+u_rup,l)        +  &  
+               ne_lup*Ftheta(ij+u_lup,l)        +  &  
+               ne_left*Ftheta(ij+u_left,l)      +  &  
+               ne_ldown*Ftheta(ij+u_ldown,l)    +  &  
+               ne_rdown*Ftheta(ij+u_rdown,l))
+       ENDDO
+
+    END DO
+
+!!! Compute potential vorticity (Coriolis) contribution to du
+
+    SELECT CASE(caldyn_conserv)
+    CASE(energy) ! energy-conserving TRiSK
+
+       CALL wait_message(req_qu)
+
+       DO l=ll_begin,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             uu_right = &
+                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
+                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
+                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
+                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
+                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
+                  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))+         &
+                  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))+         &
+                  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))+         &
+                  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))+             &
+                  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))
+             uu_lup = &
+                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
+                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
+                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
+                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
+                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
+                  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)) +          &
+                  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)) +              &
+                  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)) +              &
+                  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)) +            &
+                  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))
+             uu_ldown = &
+                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
+                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
+                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
+                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
+                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
+                  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)) +          &
+                  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)) +        &
+                  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)) +      &
+                  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)) +      &
+                  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))
+             IF(DEC) THEN
+                du(ij+u_right,l) = .5*uu_right
+                du(ij+u_lup,l)   = .5*uu_lup
+                du(ij+u_ldown,l)   = .5*uu_ldown
+             ELSE
+                du(ij+u_right,l) = .5*uu_right/de(ij+u_right)
+                du(ij+u_lup,l)   = .5*uu_lup  /de(ij+u_lup)
+                du(ij+u_ldown,l) = .5*uu_ldown/de(ij+u_ldown)
+             END IF
+          ENDDO
+       ENDDO
+
+    CASE(enstrophy) ! enstrophy-conserving TRiSK
+
+       DO l=ll_begin,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             uu_right = &
+                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
+                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
+                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
+                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
+                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
+                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
+                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
+                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
+                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
+                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
+             uu_lup = &
+                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
+                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
+                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
+                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
+                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
+                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
+                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
+                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
+                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
+                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
+             uu_ldown = &
+                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
+                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
+                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
+                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
+                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
+                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
+                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
+                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
+                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
+                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
+             IF(DEC) THEN
+                du(ij+u_right,l) = qu(ij+u_right,l)*uu_right
+                du(ij+u_lup,l)   = qu(ij+u_lup,l)  *uu_lup
+                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu_ldown
+             ELSE
+                du(ij+u_right,l) = qu(ij+u_right,l)*uu_right/de(ij+u_right)
+                du(ij+u_lup,l)   = qu(ij+u_lup,l)  *uu_lup  /de(ij+u_lup)
+                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu_ldown/de(ij+u_ldown)
+             END IF
+          ENDDO
+       ENDDO
+
+    CASE DEFAULT
+       STOP
+    END SELECT
+
+    ! Compute bernouilli term = Kinetic Energy
+    le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect
+    DO l=ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end         
+          IF(DEC) THEN
+             berni(ij,l) = &
+                  1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
+                  le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
+                  le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
+                  le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
+                  le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
+                  le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+          ELSE
+             berni(ij,l) = &
+                  1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
+                  le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
+                  le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
+                  le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
+                  le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
+                  le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+          END IF
+       ENDDO
+    ENDDO
+
+!!! Add gradients of Bernoulli and Exner functions to du
+    DO l=ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_begin,ij_end
+          IF(DEC) THEN
+             du(ij+u_right,l) = du(ij+u_right,l)                       &
+                  + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l)
+             du(ij+u_lup,l) = du(ij+u_lup,l)                           &
+                  + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l)
+             du(ij+u_ldown,l) = du(ij+u_ldown,l)                       &
+                  + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l)
+          ELSE
+             du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right)            &
+                  * ( ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
+             du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup)                 &
+                  * ( ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown)            &
+                  * ( ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
+          END IF
+       END DO
+    ENDDO
+
+    CALL trace_end("compute_caldyn_slow")
+
+  END SUBROUTINE compute_caldyn_slow
+
+END MODULE caldyn_kernels_hevi_mod
Index: codes/icosagcm/trunk/src/euler_scheme.f90
===================================================================
--- codes/icosagcm/trunk/src/euler_scheme.f90	(revision 361)
+++ codes/icosagcm/trunk/src/euler_scheme.f90	(revision 362)
@@ -19,6 +19,5 @@
   !$OMP THREADPRIVATE(nb_stage, matsuno_period, scheme, scheme_family)
 
-
-  PUBLIC :: euler_scheme, accumulate_fluxes
+  PUBLIC :: euler_scheme, accumulate_fluxes, legacy_to_DEC, DEC_to_legacy
 CONTAINS
 
@@ -143,3 +142,77 @@
   END SUBROUTINE accumulate_fluxes
 
+  SUBROUTINE legacy_to_DEC(f_ps, f_u)
+    USE icosa
+    USE disvert_mod
+    USE omp_para
+    USE trace
+    TYPE(t_field),POINTER :: f_ps(:), f_u(:)
+    REAL(rstd), POINTER :: ps(:), u(:,:)
+    INTEGER :: ind,ij,l
+    CALL trace_start("legacy_to_DEC")
+
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+
+       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN ! update ps
+          ps=f_ps(ind)
+          !$SIMD
+          DO ij=ij_begin,ij_end
+             ps(ij)=(ps(ij)-ptop)/g ! convert ps to column-integrated mass
+          ENDDO
+       END IF
+
+       u=f_u(ind)
+       DO l=ll_begin,ll_end
+          !$SIMD
+          DO ij=ij_begin,ij_end
+             u(ij+u_right,l)=u(ij+u_right,l)*de(ij+u_right)
+             u(ij+u_lup,l)=u(ij+u_lup,l)*de(ij+u_lup)
+             u(ij+u_ldown,l)=u(ij+u_ldown,l)*de(ij+u_ldown)
+          ENDDO
+       ENDDO
+    ENDDO
+
+    CALL trace_end("legacy_to_DEC")  
+  END SUBROUTINE Legacy_to_DEC
+
+  SUBROUTINE DEC_to_legacy(f_ps, f_u)
+    USE icosa
+    USE disvert_mod
+    USE omp_para
+    USE trace
+    TYPE(t_field),POINTER :: f_ps(:), f_u(:)
+    REAL(rstd), POINTER :: ps(:), u(:,:)
+    INTEGER :: ind,ij,l
+    CALL trace_start("legacy_to_DEC")
+
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+
+       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN
+          ps=f_ps(ind)
+          !$SIMD
+          DO ij=ij_begin,ij_end
+             ps(ij)=ptop+ps(ij)*g ! convert column-integrated mass to ps
+          ENDDO
+       ENDIF
+       
+       u=f_u(ind)
+       DO l=ll_begin,ll_end
+          !$SIMD
+          DO ij=ij_begin,ij_end
+             u(ij+u_right,l)=u(ij+u_right,l)/de(ij+u_right)
+             u(ij+u_lup,l)=u(ij+u_lup,l)/de(ij+u_lup)
+             u(ij+u_ldown,l)=u(ij+u_ldown,l)/de(ij+u_ldown)
+          ENDDO
+       ENDDO
+    ENDDO
+
+    CALL trace_end("DEC_to_legacy")  
+  END SUBROUTINE DEC_to_legacy
+
 END MODULE euler_scheme_mod
Index: codes/icosagcm/trunk/src/geometry.f90
===================================================================
--- codes/icosagcm/trunk/src/geometry.f90	(revision 361)
+++ codes/icosagcm/trunk/src/geometry.f90	(revision 362)
@@ -22,4 +22,5 @@
     TYPE(t_field),POINTER :: de(:)
     TYPE(t_field),POINTER :: le(:)
+    TYPE(t_field),POINTER :: le_de(:) ! le/de, 0. if de=0.
     TYPE(t_field),POINTER :: Riv(:)
     TYPE(t_field),POINTER :: Riv2(:)
@@ -69,4 +70,6 @@
   REAL(rstd),POINTER :: le(:)          ! lenght of a edge
 !$OMP THREADPRIVATE(le)
+  REAL(rstd),POINTER :: le_de(:)          ! le/de
+!$OMP THREADPRIVATE(le_de)
   REAL(rstd),POINTER :: Riv(:,:)       ! weight
 !$OMP THREADPRIVATE(Riv)
@@ -116,4 +119,5 @@
     CALL allocate_field(geom%de,field_u,type_real)
     CALL allocate_field(geom%le,field_u,type_real)
+    CALL allocate_field(geom%le_de,field_u,type_real)
     CALL allocate_field(geom%bi,field_t,type_real)
     CALL allocate_field(geom%Av,field_z,type_real)
@@ -150,4 +154,5 @@
     de=geom%de(ind)
     le=geom%le(ind)
+    le_de=geom%le_de(ind)
     Av=geom%Av(ind)
     Riv=geom%Riv(ind)
@@ -438,5 +443,9 @@
           CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup))
           CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown))
-         
+
+          le_de(n+u_right)=le(n+u_right)/de(n+u_right) ! NaN possible but should be harmless
+          le_de(n+u_lup)  =le(n+u_lup)  /de(n+u_lup)
+          le_de(n+u_ldown)=le(n+u_ldown)/de(n+u_ldown)
+
           Ai(n)=0
           DO k=0,5
@@ -528,4 +537,5 @@
             CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1)
             CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2)
+!            Riv(n,k+1)=0.5*(S1+S2)*(radius**2) ! Definition modified for DEC
             Riv(n,k+1)=0.5*(S1+S2)/Ai(n)
             Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1)
Index: codes/icosagcm/trunk/src/hevi_scheme.f90
===================================================================
--- codes/icosagcm/trunk/src/hevi_scheme.f90	(revision 361)
+++ codes/icosagcm/trunk/src/hevi_scheme.f90	(revision 362)
@@ -1,7 +1,8 @@
 MODULE hevi_scheme_mod
-  USE euler_scheme_mod
   USE prec
   USE domain_mod
   USE field_mod
+  USE euler_scheme_mod
+  USE caldyn_kernels_base_mod, ONLY : DEC
   IMPLICIT NONE
   PRIVATE
@@ -53,4 +54,5 @@
     REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
 
+    IF(DEC) CALL legacy_to_DEC(f_ps, f_u)
     DO j=1,nb_stage
        CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), &
@@ -79,4 +81,5 @@
        END DO
     END DO
+    IF(DEC) CALL DEC_to_legacy(f_ps, f_u)
   END SUBROUTINE HEVI_scheme
   
@@ -87,10 +90,12 @@
     REAL(rstd), POINTER :: y(:,:), dy(:,:)
     INTEGER :: ind
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       dy=f_dy(ind); y=f_y(ind)
-       CALL compute_update(w,y,dy)
-    ENDDO
+    IF(w /= 0.) THEN
+       DO ind=1,ndomain
+          IF (.NOT. assigned_domain(ind)) CYCLE
+          CALL swap_dimensions(ind)
+          dy=f_dy(ind); y=f_y(ind)
+          CALL compute_update(w,y,dy)
+       ENDDO
+    END IF
   END SUBROUTINE update
     
