Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 360)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 361)
@@ -2,21 +2,7 @@
   USE icosa
   USE transfert_mod
+  USE caldyn_kernels_mod
   IMPLICIT NONE
   PRIVATE
-
-  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
 
   PUBLIC init_caldyn, caldyn_BC, caldyn, req_ps, req_mass
@@ -163,6 +149,4 @@
 !$OMP THREADPRIVATE(first)
     
-    ! MPI messages need to be sent at first call to caldyn
-    ! This is needed only once : the next ones will be sent by timeloop
     IF (first) THEN 
       first=.FALSE.
@@ -175,4 +159,6 @@
       CALL init_message(f_u,req_e1_vect,req_u)
       CALL init_message(f_qu,req_e1_scal,req_qu)
+      ! Overlapping com/compute (deactivated) : MPI messages need to be sent at first call to caldyn
+      ! This is needed only once : the next ones will be sent by timeloop
 !      IF(caldyn_eta==eta_mass) THEN
 !         CALL send_message(f_ps,req_ps) 
@@ -186,12 +172,12 @@
     CALL trace_start("caldyn")
 
-      IF(caldyn_eta==eta_mass) THEN
-         CALL send_message(f_ps,req_ps) 
-      ELSE
-         CALL send_message(f_mass,req_mass) 
-      END IF
-
-    CALL send_message(f_theta_rhodz,req_theta_rhodz) 
-    CALL send_message(f_u,req_u)
+    IF(caldyn_eta==eta_mass) THEN
+       CALL send_message(f_ps,req_ps) ! COM00
+    ELSE
+       CALL send_message(f_mass,req_mass) ! COM00
+    END IF
+    
+    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
+    CALL send_message(f_u,req_u) ! COM02
 
     SELECT CASE(caldyn_conserv)
@@ -208,10 +194,10 @@
           qu=f_qu(ind)
           qv=f_qv(ind)
-          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
+          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv) ! COM00 COM01 COM02
        ENDDO
 !       CALL checksum(f_mass)
 !       CALL checksum(f_theta)
 
-       CALL send_message(f_qu,req_qu)
+       CALL send_message(f_qu,req_qu) ! COM03 wait_message in caldyn_horiz
 !       CALL wait_message(req_qu)
 
@@ -287,588 +273,4 @@
 END SUBROUTINE caldyn
 
-SUBROUTINE compute_planetvel(planetvel)
-  USE wind_mod
-  REAL(rstd),INTENT(OUT)  :: planetvel(iim*3*jjm)
-  REAL(rstd) :: ulon(iim*3*jjm)
-  REAL(rstd) :: ulat(iim*3*jjm)
-  REAL(rstd) :: lon,lat
-  INTEGER :: ij
-  DO ij=ij_begin_ext,ij_end_ext
-     ulon(ij+u_right)=radius*omega*cos(lat_e(ij+u_right))
-     ulat(ij+u_right)=0
-
-     ulon(ij+u_lup)=radius*omega*cos(lat_e(ij+u_lup))
-     ulat(ij+u_lup)=0
-
-     ulon(ij+u_ldown)=radius*omega*cos(lat_e(ij+u_ldown))
-     ulat(ij+u_ldown)=0
-  END DO
-  CALL compute_wind2D_perp_from_lonlat_compound(ulon, ulat, planetvel)
-END SUBROUTINE compute_planetvel
-    
-SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
-  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)  :: u(iim*3*jjm,llm)
-  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)
-  REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
-  
-  INTEGER :: i,j,ij,l
-  REAL(rstd) :: etav,hv, m
-!  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity  
-  
-  CALL trace_start("compute_pvort")  
-
-  IF(caldyn_eta==eta_mass) THEN
-     CALL wait_message(req_ps)  
-  ELSE
-     CALL wait_message(req_mass)
-  END IF
-  CALL wait_message(req_theta_rhodz) 
-
-  IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
-     DO l = ll_begin,ll_end
-        CALL test_message(req_u) 
-!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
-        CALL test_message(req_u) 
-!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 wait_message(req_u)   
-  
-!!! 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")
-                                                       
-  END SUBROUTINE compute_pvort
-  
-  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
-
-  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, 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(INOUT) :: pk(iim*jjm,llm) ! Exner function
-    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
-
-    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_horiz")
-
-!    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 + 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
-          DO ij=ij_begin,ij_end         
-                
-                berni(ij,l) = pk(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 )  
-                ! 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)) &
-                     + 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) * (             &
-                               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) )
-        
-        
-             du(ij+u_lup,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) )
-                
-             du(ij+u_ldown,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) )
-
-       ENDDO
-    ENDDO
- 
-    CALL trace_end("compute_caldyn_horiz")
-
-END SUBROUTINE compute_caldyn_horiz
-
-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_geopot")
-
-    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
-
 !-------------------------------- Diagnostics ----------------------------
 
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 361)
@@ -0,0 +1,152 @@
+MODULE caldyn_hevi_mod
+  USE icosa
+  USE transfert_mod
+  USE caldyn_kernels_mod
+  USE caldyn_gcm_mod
+
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC caldyn_hevi
+
+CONTAINS
+  
+  SUBROUTINE caldyn_hevi(write_out,tau, f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
+       f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du_slow, f_du_fast)
+    USE icosa
+    USE observable_mod
+    USE disvert_mod, ONLY : caldyn_eta, eta_mass
+    USE vorticity_mod
+    USE kinetic_mod
+    USE theta2theta_rhodz_mod
+    USE wind_mod
+    USE mpipara
+    USE trace
+    USE omp_para
+    USE output_field_mod
+    USE checksum_mod
+    IMPLICIT NONE
+    LOGICAL,INTENT(IN)    :: write_out
+    REAL(rstd), INTENT(IN) :: tau
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_mass(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_q(:)
+    TYPE(t_field),POINTER :: f_geopot(:)
+    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
+    TYPE(t_field) :: f_dps(:)
+    TYPE(t_field) :: f_dmass(:)
+    TYPE(t_field) :: f_dtheta_rhodz(:)
+    TYPE(t_field) :: f_du_slow(:)
+    TYPE(t_field) :: f_du_fast(:)
+    
+    REAL(rstd),POINTER :: ps(:), dps(:)
+    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
+    REAL(rstd),POINTER :: du(:,:), hflux(:,:), wflux(:,:)
+    REAL(rstd),POINTER :: u(:,:), qu(:,:), qv(:,:)
+
+! temporary shared variable
+    REAL(rstd),POINTER  :: theta(:,:)  
+    REAL(rstd),POINTER  :: pk(:,:)
+    REAL(rstd),POINTER  :: geopot(:,:)
+    REAL(rstd),POINTER  :: convm(:,:) 
+    REAL(rstd),POINTER  :: wwuu(:,:)
+        
+    INTEGER :: ind
+    LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+    
+    IF (first) THEN 
+      first=.FALSE.
+      IF(caldyn_eta==eta_mass) THEN
+         CALL init_message(f_ps,req_i1,req_ps)
+      ELSE
+         CALL init_message(f_mass,req_i1,req_mass)
+      END IF
+      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
+      CALL init_message(f_u,req_e1_vect,req_u)
+      CALL init_message(f_qu,req_e1_scal,req_qu)
+    ENDIF
+    
+    CALL trace_start("caldyn")
+    
+    IF(caldyn_eta==eta_mass) THEN
+       CALL send_message(f_ps,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
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps=f_ps(ind)
+       u=f_u(ind)
+       du=f_du_fast(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       mass=f_mass(ind)
+       theta = f_theta(ind)
+       pk = f_pk(ind)
+       geopot = f_geopot(ind)  
+       CALL compute_theta(ps,theta_rhodz, mass,theta)
+       CALL compute_geopot(ps,mass,theta, pk,geopot)
+       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot, du) ! computes du_fast and updates u
+    ENDDO
+    
+    CALL send_message(f_u,req_u) ! COM02
+    CALL wait_message(req_u)   ! COM02
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       mass=f_mass(ind)
+       qu=f_qu(ind)
+       qv=f_qv(ind)
+       CALL compute_pvort_only(u,mass,qu,qv)
+    ENDDO
+    
+    CALL send_message(f_qu,req_qu) ! COM03
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       mass=f_mass(ind)
+       theta = f_theta(ind)
+       qu=f_qu(ind)
+       hflux=f_hflux(ind)
+       convm = f_dmass(ind)
+       dtheta_rhodz=f_dtheta_rhodz(ind)
+       du=f_du_slow(ind)
+       CALL compute_caldyn_slow(u,mass,qu,theta, hflux,convm,dtheta_rhodz,du) ! COM03
+       IF(caldyn_eta==eta_mass) THEN
+          wflux=f_wflux(ind)
+          wwuu=f_wwuu(ind)
+          dps=f_dps(ind)
+          CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
+       END IF
+    ENDDO
+    
+!$OMP BARRIER
+    !    CALL check_mass_conservation(f_ps,f_dps)
+    CALL trace_end("caldyn_hevi")
+!!$OMP BARRIER
+    
+  END SUBROUTINE caldyn_hevi
+
+END MODULE caldyn_hevi_mod
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 361)
@@ -0,0 +1,1011 @@
+MODULE caldyn_kernels_mod
+  USE icosa
+  USE transfert_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
+
+CONTAINS
+
+  SUBROUTINE compute_planetvel(planetvel)
+    USE wind_mod
+    REAL(rstd),INTENT(OUT)  :: planetvel(iim*3*jjm)
+    REAL(rstd) :: ulon(iim*3*jjm)
+    REAL(rstd) :: ulat(iim*3*jjm)
+    REAL(rstd) :: lon,lat
+    INTEGER :: ij
+    DO ij=ij_begin_ext,ij_end_ext
+       ulon(ij+u_right)=radius*omega*cos(lat_e(ij+u_right))
+       ulat(ij+u_right)=0
+
+       ulon(ij+u_lup)=radius*omega*cos(lat_e(ij+u_lup))
+       ulat(ij+u_lup)=0
+
+       ulon(ij+u_ldown)=radius*omega*cos(lat_e(ij+u_ldown))
+       ulat(ij+u_ldown)=0
+    END DO
+    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)
+    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)  :: u(iim*3*jjm,llm)
+    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)
+    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
+
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: etav,hv, m
+    !  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity  
+
+    CALL trace_start("compute_pvort")  
+
+    IF(caldyn_eta==eta_mass) THEN
+       CALL wait_message(req_ps)  ! COM00
+    ELSE
+       CALL wait_message(req_mass) ! COM00
+    END IF
+    CALL wait_message(req_theta_rhodz) ! COM01
+
+    IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
+       DO l = ll_begin,ll_end
+          CALL test_message(req_u) 
+          !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
+          CALL test_message(req_u) 
+          !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 wait_message(req_u) ! COM02
+
+!!! 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")
+  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 **********************
+
+  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, 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(INOUT) :: pk(iim*jjm,llm) ! Exner function
+    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
+
+    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_horiz")
+
+    !    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) ! COM03
+
+       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 + 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
+          DO ij=ij_begin,ij_end         
+
+             berni(ij,l) = pk(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 )  
+             ! 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)) &
+                  + 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) * (             &
+               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) )
+
+
+          du(ij+u_lup,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) )
+
+          du(ij+u_ldown,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) )
+
+       ENDDO
+    ENDDO
+
+    CALL trace_end("compute_caldyn_horiz")
+
+  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/hevi_scheme.f90
===================================================================
--- codes/icosagcm/trunk/src/hevi_scheme.f90	(revision 360)
+++ codes/icosagcm/trunk/src/hevi_scheme.f90	(revision 361)
@@ -48,5 +48,5 @@
     USE time_mod
     USE disvert_mod
-    USE caldyn_mod
+    USE caldyn_hevi_mod
     LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
     INTEGER :: it,j,l, ind
@@ -54,18 +54,10 @@
 
     DO j=1,nb_stage
-       !       CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), &
-       !            taujj(j), f_phis, f_geopot, &
-       !            f_ps,f_mass,f_theta_rhodz,f_u, &
-       !            f_hflux, f_wflux, &
-       !            f_dps_slow(:,j), f_dmass_slow(:,j), f_dtheta_rhodz_slow(:,j), &
-       !            f_du_slow(:,j), f_du_fast(:,j) )
-
-       CALL caldyn((j==1) .AND. (MOD(it,itau_out)==0), &
+       CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), &
             f_phis, f_ps,f_mass,f_theta_rhodz,f_u,f_q, &
             f_geopot, f_hflux, f_wflux, &
             f_dps_slow(:,j), f_dmass_slow(:,j), f_dtheta_rhodz_slow(:,j), &
-            f_du_slow(:,j))
-
-       ! cumulate mass fluxes for transport scheme
+            f_du_slow(:,j), f_du_fast(:,j) )
+       ! accumulate mass fluxes for transport scheme
        DO ind=1,ndomain
           IF (.NOT. assigned_domain(ind)) CYCLE
@@ -75,5 +67,4 @@
           CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, wj(j), fluxt_zero(ind))
        END DO
-
        ! update model state
        DO l=1,j
@@ -85,5 +76,5 @@
           CALL update(bjl(l,j), f_theta_rhodz, f_dtheta_rhodz_slow(:,l))
           CALL update(bjl(l,j), f_u, f_du_slow(:,l))
-!          CALL update(cjl(l,j), f_u, f_du_fast(:,l)) ! FIXME - only slow terms right now
+          CALL update(cjl(l,j), f_u, f_du_fast(:,l))
        END DO
     END DO
