Index: codes/icosagcm/trunk/src/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 368)
+++ codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 369)
@@ -110,5 +110,7 @@
        pk = f_pk(ind)
        geopot = f_geopot(ind)
+       du=f_du_fast(ind)
        IF(hydrostatic) THEN
+          du(:,:)=0.
           CALL compute_geopot(ps,mass,theta, pk,geopot)
        ELSE
@@ -116,8 +118,7 @@
           dW = f_dW_fast(ind)
           dPhi = f_dPhi_fast(ind)
-          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W,dPhi,dW) ! computes d(Phi,W)_fast and updates Phi,W
+          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W,dPhi,dW,du) ! computes d(Phi,W)_fast and updates Phi,W
        END IF
        u=f_u(ind)
-       du=f_du_fast(ind)
        CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot,du) ! computes du_fast and updates du
     ENDDO
@@ -138,5 +139,6 @@
     
     CALL send_message(f_qu,req_qu) ! COM03
-    
+    CALL wait_message(req_qu) ! COM03
+       
     DO ind=1,ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
@@ -151,10 +153,24 @@
        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(hydrostatic) THEN
+          CALL compute_caldyn_slow_hydro(u,mass,hflux,du)
+       ELSE
+          W = f_W(ind)
+          dW = f_dW_slow(ind)
+          geopot = f_geopot(ind)
+          dPhi = f_dPhi_slow(ind)
+          CALL compute_caldyn_slow_NH(u,mass,geopot,W, hflux,du,dPhi,dW)
+       END IF
+       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
        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)
+          IF(hydrostatic) THEN
+             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
+          ELSE
+             PRINT *, 'NH dynamics limited to Lagrangian vertical coordinate so far !' 
+             STOP
+          END IF
        END IF
     ENDDO
Index: codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 368)
+++ codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 369)
@@ -1,4 +1,7 @@
 MODULE caldyn_kernels_hevi_mod
   USE icosa
+  USE trace
+  USE omp_para
+  USE disvert_mod
   USE transfert_mod
   USE caldyn_kernels_base_mod
@@ -10,5 +13,6 @@
   LOGICAL, PARAMETER :: debug_hevi_solver = .FALSE.
 
-  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_slow, &
+  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, &
+       compute_caldyn_slow_hydro, compute_caldyn_slow_NH, &
        compute_caldyn_solver, compute_caldyn_fast
 
@@ -16,9 +20,4 @@
 
   SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
-    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)
@@ -55,8 +54,4 @@
 
   SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
-    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)
@@ -124,7 +119,4 @@
 
   SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il)
-    USE omp_para
-    USE disvert_mod
-    IMPLICIT NONE
     REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs
     REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm)
@@ -282,10 +274,5 @@
   END SUBROUTINE compute_NH_geopot
 
-  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, dPhi,dW)
-    USE disvert_mod
-    USE exner_mod
-    USE trace
-    USE omp_para
-    IMPLICIT NONE
+  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, dPhi,dW,du)
     REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
@@ -296,6 +283,8 @@
     REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1)
     REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT)   :: du(3*iim*jjm,llm)
 
     REAL(rstd) :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces
+    REAL(rstd) :: berni(iim*jjm)             ! (W/m_il)^2
     REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij
     INTEGER    :: ij, l
@@ -342,5 +331,5 @@
     ENDDO
 
-    ! Update W, compute tendencies for geopotential and vertical momentum
+    ! Update W, compute tendencies
     DO l=2,llm
        !DIR$ SIMD
@@ -368,9 +357,16 @@
     !    ENDDO
     
-    ! Compute Exner function (needed by compute_caldyn_fast)
+    ! Compute Exner function (needed by compute_caldyn_fast) and du=g^-2.grad(w^2)
     DO l=1,llm
        !DIR$ SIMD
        DO ij=ij_begin_ext,ij_end_ext
           pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
+          berni(ij) = (-.5*g*g)*( (W(ij,l)/m_il(ij,l))**2 &
+               + (W(ij,l+1)/m_il(ij,l+1))**2 )
+       ENDDO
+       DO ij=ij_begin,ij_end         
+          du(ij+u_right,l) = ne_right*berni(ij)+ne_left *berni(ij+t_right)
+          du(ij+u_lup,l)   = ne_lup  *berni(ij)+ne_rdown*berni(ij+t_lup)
+          du(ij+u_ldown,l) = ne_ldown*berni(ij)+ne_rup  *berni(ij+t_ldown)
        ENDDO
     ENDDO
@@ -381,10 +377,4 @@
   
   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)   ! OUT if tau>0
@@ -393,5 +383,5 @@
     REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
-    REAL(rstd),INTENT(OUT)   :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm)
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
 
@@ -424,5 +414,5 @@
     END IF ! Boussinesq/compressible
 
-!!! u:=u+tau*du, du = gradients of Bernoulli and Exner functions
+!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
     DO l=ll_begin,ll_end
        !DIR$ SIMD
@@ -437,15 +427,10 @@
                          *(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
+          du(ij+u_right,l) = du(ij+u_right,l) + due_right
+          du(ij+u_lup,l)   = du(ij+u_lup,l)   + due_lup
+          du(ij+u_ldown,l) = du(ij+u_ldown,l) + due_ldown
+          u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+          u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+          u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
        ENDDO
     ENDDO
@@ -455,58 +440,30 @@
   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)
+  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
+    REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
     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(IN)  :: qu(3*iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
     REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
-    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),INTENT(INOUT) :: du(3*iim*jjm,llm)
+    
+    REAL(rstd) :: Ftheta(3*iim*jjm)  ! potential temperature flux
     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
+    INTEGER :: ij,l,kdown
+
+    CALL trace_start("compute_caldyn_Coriolis")
+
+    DO l=ll_begin, ll_end
+       ! compute theta flux
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext      
+          Ftheta(ij+u_right) = 0.5*(theta(ij,l)+theta(ij+t_right,l)) &
+                                  * hflux(ij+u_right,l)
+          Ftheta(ij+u_lup)   = 0.5*(theta(ij,l)+theta(ij+t_lup,l)) &
+                                  * hflux(ij+u_lup,l)
+          Ftheta(ij+u_ldown) = 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
+                                  * hflux(ij+u_ldown,l)
+       ENDDO
+       ! compute horizontal divergence of fluxes
        DO ij=ij_begin,ij_end         
           ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
@@ -517,23 +474,18 @@
                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
-
+          dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right)    +  &
+               ne_rup*Ftheta(ij+u_rup)        +  &  
+               ne_lup*Ftheta(ij+u_lup)        +  &  
+               ne_left*Ftheta(ij+u_left)      +  &  
+               ne_ldown*Ftheta(ij+u_ldown)    +  &  
+               ne_rdown*Ftheta(ij+u_rdown))
+       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
@@ -573,13 +525,7 @@
                   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
+             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown
           ENDDO
        ENDDO
@@ -623,70 +569,203 @@
                   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
-
+
+             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown
+          END DO
+       END DO
     CASE DEFAULT
        STOP
     END SELECT
 
-    ! Compute bernouilli term = Kinetic Energy
+    CALL trace_end("compute_caldyn_Coriolis")
+
+  END SUBROUTINE compute_caldyn_Coriolis
+
+  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du)
+    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
+    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
+    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
+    
+    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function
+    REAL(rstd) :: uu_right, uu_lup, uu_ldown
+    INTEGER :: ij,l
+
+    CALL trace_start("compute_caldyn_slow_hydro")
+
     le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect
-    DO l=ll_begin,ll_end
+
+    DO l = ll_begin, ll_end
+       !  Compute mass 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)
+          uu_right= uu_right*le_de(ij+u_right)
+          uu_lup  = uu_lup  *le_de(ij+u_lup)
+          uu_ldown= uu_ldown*le_de(ij+u_ldown)
+          hflux(ij+u_right,l)=uu_right
+          hflux(ij+u_lup,l)  =uu_lup
+          hflux(ij+u_ldown,l)=uu_ldown
+       ENDDO
+       ! Compute Bernoulli=kinetic energy
        !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
+          berni(ij) = &
+               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 )  
+       ENDDO
+       ! Compute du=grad(Bernoulli)
        !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
+             du(ij+u_right,l) = ne_right*berni(ij)+ne_left*berni(ij+t_right)
+             du(ij+u_lup,l)   = ne_lup*berni(ij)+ne_rdown*berni(ij+t_lup)
+             du(ij+u_ldown,l) = ne_ldown*berni(ij)+ne_rup*berni(ij+t_ldown)
        END DO
-    ENDDO
-
-    CALL trace_end("compute_caldyn_slow")
-
-  END SUBROUTINE compute_caldyn_slow
+    END DO
+
+    CALL trace_end("compute_caldyn_slow_hydro")    
+  END SUBROUTINE compute_caldyn_slow_hydro
+
+  SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, hflux,du,dPhi,dW)
+    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
+    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)  ! rho*dz
+    REAL(rstd),INTENT(IN)  :: Phi(iim*jjm,llm+1)  ! prognostic geopotential
+    REAL(rstd),INTENT(IN)  :: W(iim*jjm,llm+1)    ! prognostic vertical momentum
+
+    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
+    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)
+
+    REAL(rstd) :: w_il(3*iim*jjm,llm+1) ! Wil/mil
+    REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux
+    REAL(rstd) :: GradPhi2(3*iim*jjm,llm+1) ! grad_Phi**2
+    REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi)
+    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function
+    REAL(rstd) :: G_el(3*iim*jjm) ! horizontal flux of W
+    REAL(rstd) :: v_el(3*iim*jjm)
+    
+    INTEGER :: ij,l,kdown,kup
+    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, W2_el
+
+    CALL trace_start("compute_caldyn_slow_NH")
+
+    le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect
+
+    DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
+       IF(l==1) THEN
+          kdown=1
+       ELSE
+          kdown=l-1
+       END IF
+       IF(l==llm+1) THEN
+          kup=llm
+       ELSE
+          kup=l
+       END IF
+       ! compute mil, wil=Wil/mil
+       DO ij=ij_begin_ext, ij_end_ext
+          w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup))
+       END DO
+       ! compute DePhi, v_el, G_el, F_el
+       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+       ! while DePhil, W_el and F_el don't
+       DO ij=ij_begin_ext, ij_end_ext
+          ! Compute on edge 'right'
+          W_el  = .5*( W(ij,l)+W(ij+t_right,l) )
+          DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
+          F_el(ij+u_right,l)   = DePhil(ij+u_right,l)*W_el
+          W2_el = .5*le_de(ij+u_right) * &
+               ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
+          v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown))
+          G_el(ij+u_right) = v_el(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
+          ! Compute on edge 'lup'
+          W_el  = .5*( W(ij,l)+W(ij+t_lup,l) )
+          DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
+          F_el(ij+u_lup,l)   = DePhil(ij+u_lup,l)*W_el
+          W2_el = .5*le_de(ij+u_lup) * &
+               ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
+          v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown) )
+          G_el(ij+u_lup) = v_el(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
+          ! Compute on edge 'ldown'
+          W_el  = .5*( W(ij,l)+W(ij+t_ldown,l) )
+          DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
+          F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el
+          W2_el = .5*le_de(ij+u_ldown) * &
+               ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
+          v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown) )
+          G_el(ij+u_ldown) = v_el(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
+       END DO
+       ! compute GradPhi2, dPhi, dW
+       DO ij=ij_begin_ext, ij_end_ext
+          gradPhi2(ij,l) = &
+               1/(2*Ai(ij))*(le_de(ij+u_right)*DePhil(ij+u_right,l)**2 + &
+               le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2 +      &
+               le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2 +      &
+               le_de(ij+u_left)*DePhil(ij+u_left,l)**2 +    &
+               le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  &
+               le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 )
+          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l)-1/(2*Ai(ij))* & 
+               ( DePhil(ij+u_right,l)*v_el(ij+u_right) + & ! -v.gradPhi, 
+                 DePhil(ij+u_rup,l)*v_el(ij+u_rup) +     & ! v_el already has le_de
+                 DePhil(ij+u_lup,l)*v_el(ij+u_lup) +     &
+                 DePhil(ij+u_left,l)*v_el(ij+u_left) +   &
+                 DePhil(ij+u_ldown,l)*v_el(ij+u_ldown) + &
+                 DePhil(ij+u_rdown,l)*v_el(ij+u_rdown) )
+          dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el), 
+               ne_right*G_el(ij+u_right) +  & ! G_el already has le_de
+               ne_rup*G_el(ij+u_rup) +      &
+               ne_lup*G_el(ij+u_lup) +      &  
+               ne_left*G_el(ij+u_left) +    &
+               ne_ldown*G_el(ij+u_ldown) +  &
+               ne_rdown*G_el(ij+u_rdown))
+       END DO
+    END DO
+    
+    DO l=ll_begin, ll_end ! compute on k levels (layers)
+       ! Compute berni at scalar points
+       DO ij=ij_begin_ext, ij_end_ext
+          berni(ij) = &
+               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 ) &
+               - .5*( gradPhi2(ij,l)  *w_il(ij,l)**2 +   &
+                      gradPhi2(ij,l+1)*w_il(ij,l+1)**2 )
+       END DO
+       ! Compute mass flux and grad(berni) at edges
+       DO ij=ij_begin_ext, ij_end_ext
+          ! Compute on edge 'right'
+          uu_right = 0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) &
+                    -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))
+          hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)
+          du(ij+u_right,l) = ne_right*berni(ij)+ne_left*berni(ij+t_right)
+          ! Compute on edge 'lup'
+          uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) &
+                  -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))
+          hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)
+          du(ij+u_lup,l) = ne_lup*berni(ij)+ne_rdown*berni(ij+t_lup)
+          ! Compute on edge 'ldown'
+          uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) &
+                    -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))
+          hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)
+          du(ij+u_ldown,l) = ne_ldown*berni(ij)+ne_rup*berni(ij+t_ldown)
+       END DO
+    END DO
+
+    CALL trace_end("compute_caldyn_slow_NH")
+
+  END SUBROUTINE compute_caldyn_slow_NH
 
 END MODULE caldyn_kernels_hevi_mod
Index: codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90
===================================================================
--- codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90	(revision 368)
+++ codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90	(revision 369)
@@ -1230,6 +1230,6 @@
 			    	Om      = 0.d0,			&	! Rotation Rate of Earth
                             	as      = a/X,			&	! New Radius of small Earth     
-!			    	u0      = 20.d0,		&	! Reference Velocity 
-			    	u0      = 0.d0,		        &	! FIXME : no zonal wind for NH tests
+			    	u0      = 20.d0,		&	! Reference Velocity 
+!			    	u0      = 0.d0,		        &	! FIXME : no zonal wind for NH tests
                             	Teq     = 300.d0,		&	! Temperature at Equator    
 			    	Peq     = 100000.d0,		&	! Reference PS at Equator
