Index: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 536)
+++ codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 536)
@@ -0,0 +1,840 @@
+MODULE caldyn_kernels_hevi_mod
+  USE icosa
+  USE trace
+  USE omp_para
+  USE disvert_mod
+  USE transfert_mod
+  USE caldyn_kernels_base_mod
+  IMPLICIT NONE
+  PRIVATE
+
+  REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME
+
+  LOGICAL, PARAMETER :: debug_hevi_solver = .FALSE.
+
+  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, &
+       compute_caldyn_slow_hydro, compute_caldyn_slow_NH, &
+       compute_caldyn_solver, compute_caldyn_fast
+
+CONTAINS
+
+  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
+    REAL(rstd),INTENT(IN)    :: ps(iim*jjm)
+    REAL(rstd),INTENT(IN)    :: theta_rhodz(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT)   :: theta(iim*jjm,llm,nqdyn)
+    INTEGER :: ij,l,iq
+    REAL(rstd) :: m
+    CALL trace_start("compute_theta")
+
+    IF(caldyn_eta==eta_mass) THEN ! Compute mass
+       DO l = ll_begin,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l) ! ps is actually Ms
+             rhodz(ij,l) = m/g
+          END DO
+       END DO
+    END IF
+
+    DO l = ll_begin,ll_end
+       DO iq=1,nqdyn
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             theta(ij,l,iq) = theta_rhodz(ij,l,iq)/rhodz(ij,l)
+          END DO
+       END DO
+    END DO
+
+    CALL trace_end("compute_theta")
+  END SUBROUTINE compute_theta
+
+  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
+    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
+
+    INTEGER :: ij,l
+    REAL(rstd) :: etav,hv,radius_m2
+
+    CALL trace_start("compute_pvort_only")  
+!!! Compute shallow-water potential vorticity
+#ifdef CPP_DYSL
+#include "../kernels/pvort_only.k90"
+#else
+    radius_m2=radius**(-2)
+    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)   &
+                  + ne_left * u(ij+t_rup+u_left,l)          &
+                  - ne_lup  * u(ij+u_lup,l) )                               
+          hv =   Riv2(ij,vup)          * rhodz(ij,l)        &
+               + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)  &
+               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
+          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
+          
+          etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)   &
+               + ne_right * u(ij+t_ldown+u_right,l)               &
+               - ne_rdown * u(ij+u_rdown,l) )
+          hv =   Riv2(ij,vdown)        * rhodz(ij,l)              &
+               + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)      &
+               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
+          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
+       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
+#endif
+    CALL trace_end("compute_pvort_only")
+
+  END SUBROUTINE compute_pvort_only
+
+  SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il)
+    REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs
+    REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: W_il(iim*jjm,llm+1)   ! vertical momentum
+    REAL(rstd),INTENT(INOUT) :: Phi_il(iim*jjm,llm+1) ! geopotential
+
+    REAL(rstd) :: Phi_star_il(iim*jjm,llm+1) 
+    REAL(rstd) :: p_ik(iim*jjm,llm)      ! pressure
+    REAL(rstd) :: R_il(iim*jjm,llm+1)    ! rhs of tridiag problem
+    REAL(rstd) :: x_il(iim*jjm,llm+1)    ! solution of tridiag problem
+    REAL(rstd) :: A_ik(iim*jjm,llm)      ! off-diagonal coefficients of tridiag problem
+    REAL(rstd) :: B_il(iim*jjm,llm+1)    ! diagonal coefficients of tridiag problem
+    REAL(rstd) :: C_ik(iim*jjm,llm)      ! Thomas algorithm
+    REAL(rstd) :: D_il(iim*jjm,llm+1)    ! Thomas algorithm
+    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij
+    REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik
+
+    INTEGER    :: iter, ij, l
+
+!    FIXME : vertical OpenMP parallelism will not work
+   
+    tau2_g=tau*tau/g
+    g2=g*g
+    gm2 = g**-2
+    gamma = 1./(1.-kappa)
+    
+    ! compute Phi_star
+    DO l=1,llm+1
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          Phi_star_il(ij,l) = Phi_il(ij,l) + tau*g2*(W_il(ij,l)/m_il(ij,l)-tau)
+       ENDDO
+    ENDDO
+
+    ! Newton-Raphson iteration : Phi_il contains current guess value
+    DO iter=1,5 ! 2 iterations should be enough
+
+       ! Compute pressure, A_ik
+       DO l=1,llm
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l))
+             X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
+             p_ik(ij,l) = preff*(X_ij**gamma)
+             c2_mik = gamma*p_ik(ij,l)/(rho_ij*m_ik(ij,l)) ! c^2 = gamma*R*T = gamma*p/rho
+             A_ik(ij,l) = c2_mik*(tau/g*rho_ij)**2
+          ENDDO
+       ENDDO
+
+       ! Compute residual, B_il
+       ! bottom interface l=1
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          ml_g2 = gm2*m_il(ij,1) 
+          B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot
+          R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  &
+               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-Phi_bot) )
+       ENDDO
+       ! inner interfaces
+       DO l=2,llm
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             ml_g2 = gm2*m_il(ij,l) 
+             B_il(ij,l) = A_ik(ij,l)+A_ik(ij,l-1) + ml_g2
+             R_il(ij,l) = ml_g2*( Phi_il(ij,l)-Phi_star_il(ij,l))  &
+                     + tau2_g*(p_ik(ij,l)-p_ik(ij,l-1))
+             ! consistency check : if Wil=0 and initial state is in hydrostatic balance
+             ! then Phi_star_il(ij,l) = Phi_il(ij,l) - tau^2*g^2
+             ! and residual = tau^2*(ml+(1/g)dl_pi)=0
+          ENDDO
+       ENDDO
+       ! top interface l=llm+1
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          ml_g2 = gm2*m_il(ij,llm+1) 
+          B_il(ij,llm+1) = A_ik(ij,llm) + ml_g2
+          R_il(ij,llm+1) = ml_g2*( Phi_il(ij,llm+1)-Phi_star_il(ij,llm+1))  &
+               + tau2_g*( ptop-p_ik(ij,llm) )
+       ENDDO
+
+       ! FIXME later
+       ! the lines below modify the tridiag problem 
+       ! for flat, rigid boundary conditions at top and bottom :
+       ! zero out A(1), A(llm), R(1), R(llm+1)
+       ! => x(l)=0 at l=1,llm+1
+       DO ij=ij_begin_ext,ij_end_ext
+          A_ik(ij,1) = 0.
+          A_ik(ij,llm) = 0.
+          R_il(ij,1) = 0.
+          R_il(ij,llm+1) = 0.
+       ENDDO
+
+       IF(debug_hevi_solver) THEN ! print Linf(residual)
+          PRINT *, '[hevi_solver] R,p', iter, MAXVAL(ABS(R_il)), MAXVAL(p_ik)
+       END IF
+
+       ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
+       ! Forward sweep :
+       ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)), 
+       ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
+       ! bottom interface l=1
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          X_ij = 1./B_il(ij,1)
+          C_ik(ij,1) = -A_ik(ij,1) * X_ij 
+          D_il(ij,1) =  R_il(ij,1) * X_ij
+       ENDDO
+       ! inner interfaces/layers
+       DO l=2,llm
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             X_ij = 1./(B_il(ij,l) + A_ik(ij,l-1)*C_ik(ij,l-1))
+             C_ik(ij,l) = -A_ik(ij,l) * X_ij 
+             D_il(ij,l) =  (R_il(ij,l)+A_ik(ij,l-1)*D_il(ij,l-1)) * X_ij
+          ENDDO
+       ENDDO
+       ! top interface l=llm+1
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          X_ij = 1./(B_il(ij,llm+1) + A_ik(ij,llm)*C_ik(ij,llm))
+          D_il(ij,llm+1) =  (R_il(ij,llm+1)+A_ik(ij,llm)*D_il(ij,llm)) * X_ij
+       ENDDO
+       
+       ! Back substitution :
+       ! x(i) = D(i)-C(i)x(i+1), x(N+1)=0
+       ! + Newton-Raphson update
+       x_il=0. ! FIXME
+       ! top interface l=llm+1
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          x_il(ij,llm+1) = D_il(ij,llm+1)
+          Phi_il(ij,llm+1) = Phi_il(ij,llm+1) - x_il(ij,llm+1)
+       ENDDO
+       ! lower interfaces
+       DO l=llm,1,-1
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             x_il(ij,l) = D_il(ij,l) - C_ik(ij,l)*x_il(ij,l+1)
+             Phi_il(ij,l) = Phi_il(ij,l) - x_il(ij,l)
+          ENDDO
+       ENDDO
+
+       IF(debug_hevi_solver) THEN
+          PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
+          PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
+          DO l=1,llm+1
+             WRITE(*,'(A,I2.1,I3.2,E9.2)'), '[hevi_solver] x', iter,l, MAXVAL(ABS(x_il(:,l)))
+          END DO
+       END IF
+
+    END DO ! Newton-Raphson
+    
+  END SUBROUTINE compute_NH_geopot
+
+  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)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
+    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
+
+    CALL trace_start("compute_caldyn_solver")
+
+    ! FIXME : vertical OpenMP parallelism will not work
+
+    ! average m_ik to interfaces => m_il
+    !DIR$ SIMD
+    DO ij=ij_begin_ext,ij_end_ext
+       m_il(ij,1) = .5*rhodz(ij,1)
+    ENDDO
+    DO l=2,llm
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l))
+       ENDDO
+    ENDDO
+    !DIR$ SIMD
+    DO ij=ij_begin_ext,ij_end_ext
+       m_il(ij,llm+1) = .5*rhodz(ij,llm)
+    ENDDO
+
+    IF(tau>0) THEN ! solve implicit problem for geopotential
+       CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot)
+    END IF
+
+    ! Compute pressure, stored temporarily in pk
+    ! kappa = R/Cp
+    ! 1-kappa = Cv/Cp
+    ! Cp/Cv = 1/(1-kappa)
+    gamma = 1./(1.-kappa)
+    DO l=1,llm
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
+          X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
+          ! kappa.theta.rho = p/exner 
+          ! => X = (p/p0)/(exner/Cp)
+          !      = (p/p0)^(1-kappa)
+          pk(ij,l) = preff*(X_ij**gamma)
+       ENDDO
+    ENDDO
+
+    ! Update W, compute tendencies
+    DO l=2,llm
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          dW(ij,l)   = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l)
+          W(ij,l)    = W(ij,l)+tau*dW(ij,l) ! update W
+          dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l)
+       ENDDO
+       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
+       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l))) 
+    ENDDO
+    ! Lower BC (FIXME : no orography yet !)
+    DO ij=ij_begin,ij_end         
+       dPhi(ij,1)=0
+       W(ij,1)=0
+       dW(ij,1)=0
+       dPhi(ij,llm+1)=0 ! rigid lid
+       W(ij,llm+1)=0
+       dW(ij,llm+1)=0
+    ENDDO
+    ! Upper BC p=ptop
+    !    DO ij=ij_omp_begin_ext,ij_omp_end_ext
+    !       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
+    !       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)
+    !    ENDDO
+    
+    ! 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) = (-.25*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)-berni(ij+t_right))
+          du(ij+u_lup,l)   = ne_lup  *(berni(ij)-berni(ij+t_lup))
+          du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
+       ENDDO
+    ENDDO
+    
+    CALL trace_end("compute_caldyn_solver")
+    
+  END SUBROUTINE compute_caldyn_solver
+  
+  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
+    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
+    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0
+    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm)
+    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
+    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
+
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: Rd, qv, temp, chi, nu, due, due_right, due_lup, due_ldown
+
+    CALL trace_start("compute_caldyn_fast")
+
+    Rd=cpp*kappa
+
+#ifdef CPP_DYSL
+#include "../kernels/caldyn_fast.k90"
+#else
+    ! Compute Bernoulli term
+    IF(boussinesq) THEN
+       DO l=ll_begin,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             berni(ij,l) = pk(ij,l)
+             ! from now on pk contains the vertically-averaged geopotential
+             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+          END DO
+       END DO
+    ELSE ! compressible
+
+       DO l=ll_begin,ll_end
+          SELECT CASE(caldyn_thermo)
+          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+             END DO
+          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s) 
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
+             END DO
+          CASE(thermo_moist) 
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
+                ! Bd = Phi + gibbs_d
+                ! Bv = Phi + gibbs_v
+                ! pk=temperature, theta=entropy
+                qv = theta(ij,l,2)
+                temp = pk(ij,l)
+                chi = log(temp/Treff)
+                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + temp*(cpp*(1.-chi)+Rd*nu)
+                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + temp*(cppv*(1.-chi)+Rv*nu)
+            END DO
+          END SELECT
+       END DO
+
+    END IF ! Boussinesq/compressible
+
+!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
+    DO l=ll_begin,ll_end
+       IF(caldyn_thermo == thermo_moist) THEN
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
+                       *(pk(ij+t_right,l)-pk(ij,l))             &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
+                       *(berniv(ij+t_right,l)-berniv(ij,l))
+             
+             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
+                       *(pk(ij+t_lup,l)-pk(ij,l))               &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
+                       *(berniv(ij+t_lup,l)-berniv(ij,l))
+             
+             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
+                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
+                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
+                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
+                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
+             
+             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*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)
+          END DO
+       ELSE
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
+                  *(pk(ij+t_right,l)-pk(ij,l))        &
+                  +  berni(ij+t_right,l)-berni(ij,l)
+             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
+                  *(pk(ij+t_lup,l)-pk(ij,l))          &
+                  +  berni(ij+t_lup,l)-berni(ij,l)
+             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
+                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
+                  +   berni(ij+t_ldown,l)-berni(ij,l)
+             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*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)
+          END DO
+       END IF
+    END DO
+#endif       
+    CALL trace_end("compute_caldyn_fast")
+
+  END SUBROUTINE compute_caldyn_fast
+
+  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,nqdyn) ! active scalars
+    REAL(rstd),INTENT(IN)  :: qu(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
+    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
+    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
+    
+    REAL(rstd) :: Ftheta(3*iim*jjm)  ! potential temperature flux
+    REAL(rstd) :: uu_right, uu_lup, uu_ldown
+    INTEGER :: ij,iq,l,kdown
+
+    CALL trace_start("compute_caldyn_Coriolis")
+
+    DO l=ll_begin, ll_end
+       ! compute theta flux
+       DO iq=1,nqdyn
+       !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext      
+             Ftheta(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
+                                  * hflux(ij+u_right,l)
+             Ftheta(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
+                  * hflux(ij+u_lup,l)
+             Ftheta(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
+                  * hflux(ij+u_ldown,l)
+          END DO
+          ! horizontal divergence of fluxes
+       !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             ! dtheta_rhodz =  -div(flux.theta)
+             dtheta_rhodz(ij,l,iq)= &
+                  -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) )
+          END DO
+       END DO
+
+       !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))
+       END DO ! ij
+    END DO ! llm
+
+!!! Compute potential vorticity (Coriolis) contribution to du
+    SELECT CASE(caldyn_conserv)
+
+    CASE(energy) ! energy-conserving TRiSK
+
+       DO l=ll_begin,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             uu_right = &
+                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
+                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
+                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
+                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
+                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
+                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
+                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
+                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
+                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
+                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
+             uu_lup = &
+                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
+                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
+                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
+                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
+                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
+                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
+                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
+                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
+                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
+                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
+             uu_ldown = &
+                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
+                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
+                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
+                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
+                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
+                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
+                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
+                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
+                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
+                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
+             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
+
+    CASE(enstrophy) ! enstrophy-conserving TRiSK
+
+       DO l=ll_begin,ll_end
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             uu_right = &
+                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
+                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
+                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
+                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
+                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
+                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
+                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
+                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
+                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
+                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
+             uu_lup = &
+                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
+                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
+                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
+                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
+                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
+                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
+                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
+                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
+                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
+                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
+             uu_ldown = &
+                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
+                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
+                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
+                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
+                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
+                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
+                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
+                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
+                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
+                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
+
+             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
+
+    CALL trace_end("compute_caldyn_Coriolis")
+
+  END SUBROUTINE compute_caldyn_Coriolis
+
+  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero)
+    LOGICAL, INTENT(IN) :: zero
+    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(INOUT) :: 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")
+
+    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         
+          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)
+       IF(zero) THEN
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end
+             du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
+             du(ij+u_lup,l)   = ne_lup*(berni(ij)-berni(ij+t_lup))
+             du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
+          END DO
+       ELSE
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end
+             du(ij+u_right,l) = du(ij+u_right,l) + &
+                  ne_right*(berni(ij)-berni(ij+t_right))
+             du(ij+u_lup,l)   = du(ij+u_lup,l) + &
+                  ne_lup*(berni(ij)-berni(ij+t_lup))
+             du(ij+u_ldown,l) = du(ij+u_ldown,l) + &
+                  ne_ldown*(berni(ij)-berni(ij+t_ldown))
+          END DO
+       END IF
+    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
+       ! below : "checked" means "formula also valid when kup=kdown (top/bottom)"
+       ! 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)) ! checked
+       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)) ! checked
+          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)) ! checked
+          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)) ! checked
+          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 )
+!          gradPhi2(ij,l) = 0. ! FIXME !!
+
+          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
+    ! FIXME !!
+ !   F_el(:,:)=0. 
+ !   dPhi(:,:)=0.
+ !   dW(:,:)=0.
+
+    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 ) &
+               - .25*( 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)-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)-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)-berni(ij+t_ldown))
+       END DO
+    END DO
+    ! FIXME !!
+    ! du(:,:)=0.
+    ! hflux(:,:)=0.
+
+    CALL trace_end("compute_caldyn_slow_NH")
+
+  END SUBROUTINE compute_caldyn_slow_NH
+
+END MODULE caldyn_kernels_hevi_mod
Index: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.f90	(revision 535)
+++ 	(revision )
@@ -1,837 +1,0 @@
-MODULE caldyn_kernels_hevi_mod
-  USE icosa
-  USE trace
-  USE omp_para
-  USE disvert_mod
-  USE transfert_mod
-  USE caldyn_kernels_base_mod
-  IMPLICIT NONE
-  PRIVATE
-
-  REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME
-
-  LOGICAL, PARAMETER :: debug_hevi_solver = .FALSE.
-
-  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, &
-       compute_caldyn_slow_hydro, compute_caldyn_slow_NH, &
-       compute_caldyn_solver, compute_caldyn_fast
-
-CONTAINS
-
-  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
-    REAL(rstd),INTENT(IN)    :: ps(iim*jjm)
-    REAL(rstd),INTENT(IN)    :: theta_rhodz(iim*jjm,llm,nqdyn)
-    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT)   :: theta(iim*jjm,llm,nqdyn)
-    INTEGER :: ij,l,iq
-    REAL(rstd) :: m
-    CALL trace_start("compute_theta")
-
-    IF(caldyn_eta==eta_mass) THEN ! Compute mass
-       DO l = ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext
-             m = mass_dak(l)+(ps(ij)*g+ptop)*mass_dbk(l) ! ps is actually Ms
-             rhodz(ij,l) = m/g
-          END DO
-       END DO
-    END IF
-
-    DO l = ll_begin,ll_end
-       DO iq=1,nqdyn
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext
-             theta(ij,l,iq) = theta_rhodz(ij,l,iq)/rhodz(ij,l)
-          END DO
-       END DO
-    END DO
-
-    CALL trace_end("compute_theta")
-  END SUBROUTINE compute_theta
-
-  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
-    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
-
-    INTEGER :: ij,l
-    REAL(rstd) :: etav,hv,radius_m2
-
-    CALL trace_start("compute_pvort_only")  
-!!! Compute shallow-water potential vorticity
-    radius_m2=radius**(-2)
-    DO l = ll_begin,ll_end
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          etav= 1./Av(ij+z_up)*(  ne_rup  * u(ij+u_rup,l)   &
-                  + ne_left * u(ij+t_rup+u_left,l)          &
-                  - ne_lup  * u(ij+u_lup,l) )                               
-          hv =   Riv2(ij,vup)          * rhodz(ij,l)        &
-               + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)  &
-               + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
-          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
-          
-          etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)   &
-               + ne_right * u(ij+t_ldown+u_right,l)               &
-               - ne_rdown * u(ij+u_rdown,l) )
-          hv =   Riv2(ij,vdown)        * rhodz(ij,l)              &
-               + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)      &
-               + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
-          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
-       ENDDO
-
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end
-          qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))  
-          qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))  
-          qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))  
-       END DO
-
-    ENDDO
-
-    CALL trace_end("compute_pvort_only")
-
-  END SUBROUTINE compute_pvort_only
-
-  SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il)
-    REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs
-    REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1)
-    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: W_il(iim*jjm,llm+1)   ! vertical momentum
-    REAL(rstd),INTENT(INOUT) :: Phi_il(iim*jjm,llm+1) ! geopotential
-
-    REAL(rstd) :: Phi_star_il(iim*jjm,llm+1) 
-    REAL(rstd) :: p_ik(iim*jjm,llm)      ! pressure
-    REAL(rstd) :: R_il(iim*jjm,llm+1)    ! rhs of tridiag problem
-    REAL(rstd) :: x_il(iim*jjm,llm+1)    ! solution of tridiag problem
-    REAL(rstd) :: A_ik(iim*jjm,llm)      ! off-diagonal coefficients of tridiag problem
-    REAL(rstd) :: B_il(iim*jjm,llm+1)    ! diagonal coefficients of tridiag problem
-    REAL(rstd) :: C_ik(iim*jjm,llm)      ! Thomas algorithm
-    REAL(rstd) :: D_il(iim*jjm,llm+1)    ! Thomas algorithm
-    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij
-    REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik
-
-    INTEGER    :: iter, ij, l
-
-!    FIXME : vertical OpenMP parallelism will not work
-   
-    tau2_g=tau*tau/g
-    g2=g*g
-    gm2 = g**-2
-    gamma = 1./(1.-kappa)
-    
-    ! compute Phi_star
-    DO l=1,llm+1
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          Phi_star_il(ij,l) = Phi_il(ij,l) + tau*g2*(W_il(ij,l)/m_il(ij,l)-tau)
-       ENDDO
-    ENDDO
-
-    ! Newton-Raphson iteration : Phi_il contains current guess value
-    DO iter=1,5 ! 2 iterations should be enough
-
-       ! Compute pressure, A_ik
-       DO l=1,llm
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext
-             rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l))
-             X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
-             p_ik(ij,l) = preff*(X_ij**gamma)
-             c2_mik = gamma*p_ik(ij,l)/(rho_ij*m_ik(ij,l)) ! c^2 = gamma*R*T = gamma*p/rho
-             A_ik(ij,l) = c2_mik*(tau/g*rho_ij)**2
-          ENDDO
-       ENDDO
-
-       ! Compute residual, B_il
-       ! bottom interface l=1
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          ml_g2 = gm2*m_il(ij,1) 
-          B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot
-          R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  &
-               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-Phi_bot) )
-       ENDDO
-       ! inner interfaces
-       DO l=2,llm
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext
-             ml_g2 = gm2*m_il(ij,l) 
-             B_il(ij,l) = A_ik(ij,l)+A_ik(ij,l-1) + ml_g2
-             R_il(ij,l) = ml_g2*( Phi_il(ij,l)-Phi_star_il(ij,l))  &
-                     + tau2_g*(p_ik(ij,l)-p_ik(ij,l-1))
-             ! consistency check : if Wil=0 and initial state is in hydrostatic balance
-             ! then Phi_star_il(ij,l) = Phi_il(ij,l) - tau^2*g^2
-             ! and residual = tau^2*(ml+(1/g)dl_pi)=0
-          ENDDO
-       ENDDO
-       ! top interface l=llm+1
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          ml_g2 = gm2*m_il(ij,llm+1) 
-          B_il(ij,llm+1) = A_ik(ij,llm) + ml_g2
-          R_il(ij,llm+1) = ml_g2*( Phi_il(ij,llm+1)-Phi_star_il(ij,llm+1))  &
-               + tau2_g*( ptop-p_ik(ij,llm) )
-       ENDDO
-
-       ! FIXME later
-       ! the lines below modify the tridiag problem 
-       ! for flat, rigid boundary conditions at top and bottom :
-       ! zero out A(1), A(llm), R(1), R(llm+1)
-       ! => x(l)=0 at l=1,llm+1
-       DO ij=ij_begin_ext,ij_end_ext
-          A_ik(ij,1) = 0.
-          A_ik(ij,llm) = 0.
-          R_il(ij,1) = 0.
-          R_il(ij,llm+1) = 0.
-       ENDDO
-
-       IF(debug_hevi_solver) THEN ! print Linf(residual)
-          PRINT *, '[hevi_solver] R,p', iter, MAXVAL(ABS(R_il)), MAXVAL(p_ik)
-       END IF
-
-       ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
-       ! Forward sweep :
-       ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)), 
-       ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
-       ! bottom interface l=1
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          X_ij = 1./B_il(ij,1)
-          C_ik(ij,1) = -A_ik(ij,1) * X_ij 
-          D_il(ij,1) =  R_il(ij,1) * X_ij
-       ENDDO
-       ! inner interfaces/layers
-       DO l=2,llm
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext
-             X_ij = 1./(B_il(ij,l) + A_ik(ij,l-1)*C_ik(ij,l-1))
-             C_ik(ij,l) = -A_ik(ij,l) * X_ij 
-             D_il(ij,l) =  (R_il(ij,l)+A_ik(ij,l-1)*D_il(ij,l-1)) * X_ij
-          ENDDO
-       ENDDO
-       ! top interface l=llm+1
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          X_ij = 1./(B_il(ij,llm+1) + A_ik(ij,llm)*C_ik(ij,llm))
-          D_il(ij,llm+1) =  (R_il(ij,llm+1)+A_ik(ij,llm)*D_il(ij,llm)) * X_ij
-       ENDDO
-       
-       ! Back substitution :
-       ! x(i) = D(i)-C(i)x(i+1), x(N+1)=0
-       ! + Newton-Raphson update
-       x_il=0. ! FIXME
-       ! top interface l=llm+1
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          x_il(ij,llm+1) = D_il(ij,llm+1)
-          Phi_il(ij,llm+1) = Phi_il(ij,llm+1) - x_il(ij,llm+1)
-       ENDDO
-       ! lower interfaces
-       DO l=llm,1,-1
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext
-             x_il(ij,l) = D_il(ij,l) - C_ik(ij,l)*x_il(ij,l+1)
-             Phi_il(ij,l) = Phi_il(ij,l) - x_il(ij,l)
-          ENDDO
-       ENDDO
-
-       IF(debug_hevi_solver) THEN
-          PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
-          PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
-          DO l=1,llm+1
-             WRITE(*,'(A,I2.1,I3.2,E9.2)'), '[hevi_solver] x', iter,l, MAXVAL(ABS(x_il(:,l)))
-          END DO
-       END IF
-
-    END DO ! Newton-Raphson
-    
-  END SUBROUTINE compute_NH_geopot
-
-  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)
-    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
-    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
-    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
-
-    CALL trace_start("compute_caldyn_solver")
-
-    ! FIXME : vertical OpenMP parallelism will not work
-
-    ! average m_ik to interfaces => m_il
-    !DIR$ SIMD
-    DO ij=ij_begin_ext,ij_end_ext
-       m_il(ij,1) = .5*rhodz(ij,1)
-    ENDDO
-    DO l=2,llm
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l))
-       ENDDO
-    ENDDO
-    !DIR$ SIMD
-    DO ij=ij_begin_ext,ij_end_ext
-       m_il(ij,llm+1) = .5*rhodz(ij,llm)
-    ENDDO
-
-    IF(tau>0) THEN ! solve implicit problem for geopotential
-       CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot)
-    END IF
-
-    ! Compute pressure, stored temporarily in pk
-    ! kappa = R/Cp
-    ! 1-kappa = Cv/Cp
-    ! Cp/Cv = 1/(1-kappa)
-    gamma = 1./(1.-kappa)
-    DO l=1,llm
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
-          X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
-          ! kappa.theta.rho = p/exner 
-          ! => X = (p/p0)/(exner/Cp)
-          !      = (p/p0)^(1-kappa)
-          pk(ij,l) = preff*(X_ij**gamma)
-       ENDDO
-    ENDDO
-
-    ! Update W, compute tendencies
-    DO l=2,llm
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          dW(ij,l)   = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l)
-          W(ij,l)    = W(ij,l)+tau*dW(ij,l) ! update W
-          dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l)
-       ENDDO
-       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
-       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l))) 
-    ENDDO
-    ! Lower BC (FIXME : no orography yet !)
-    DO ij=ij_begin,ij_end         
-       dPhi(ij,1)=0
-       W(ij,1)=0
-       dW(ij,1)=0
-       dPhi(ij,llm+1)=0 ! rigid lid
-       W(ij,llm+1)=0
-       dW(ij,llm+1)=0
-    ENDDO
-    ! Upper BC p=ptop
-    !    DO ij=ij_omp_begin_ext,ij_omp_end_ext
-    !       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
-    !       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)
-    !    ENDDO
-    
-    ! 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) = (-.25*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)-berni(ij+t_right))
-          du(ij+u_lup,l)   = ne_lup  *(berni(ij)-berni(ij+t_lup))
-          du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
-       ENDDO
-    ENDDO
-    
-    CALL trace_end("compute_caldyn_solver")
-    
-  END SUBROUTINE compute_caldyn_solver
-  
-  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
-    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
-    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0
-    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
-    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
-    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
-    REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm)
-    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
-    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
-
-    INTEGER :: i,j,ij,l
-    REAL(rstd) :: Rd, qv, temp, chi, nu, due_right, due_lup, due_ldown
-
-    CALL trace_start("compute_caldyn_fast")
-
-    Rd=cpp*kappa
-
-    ! Compute Bernoulli term
-    IF(boussinesq) THEN
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             berni(ij,l) = pk(ij,l)
-             ! from now on pk contains the vertically-averaged geopotential
-             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-          END DO
-       END DO
-    ELSE ! compressible
-
-       DO l=ll_begin,ll_end
-          SELECT CASE(caldyn_thermo)
-          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-             END DO
-          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s) 
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
-             END DO
-          CASE(thermo_moist) 
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
-                ! Bd = Phi + gibbs_d
-                ! Bv = Phi + gibbs_v
-                ! pk=temperature, theta=entropy
-                qv = theta(ij,l,2)
-                temp = pk(ij,l)
-                chi = log(temp/Treff)
-                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + temp*(cpp*(1.-chi)+Rd*nu)
-                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + temp*(cppv*(1.-chi)+Rv*nu)
-            END DO
-          END SELECT
-       END DO
-
-    END IF ! Boussinesq/compressible
-
-!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
-    DO l=ll_begin,ll_end
-       IF(caldyn_thermo == thermo_moist) THEN
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
-                       *(pk(ij+t_right,l)-pk(ij,l))             &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
-                       *(berniv(ij+t_right,l)-berniv(ij,l))
-             
-             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
-                       *(pk(ij+t_lup,l)-pk(ij,l))               &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
-                       *(berniv(ij+t_lup,l)-berniv(ij,l))
-             
-             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
-                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
-                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
-             
-             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*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)
-          END DO
-       ELSE
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
-                  *(pk(ij+t_right,l)-pk(ij,l))        &
-                  +  berni(ij+t_right,l)-berni(ij,l)
-             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
-                  *(pk(ij+t_lup,l)-pk(ij,l))          &
-                  +  berni(ij+t_lup,l)-berni(ij,l)
-             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
-                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
-                  +   berni(ij+t_ldown,l)-berni(ij,l)
-             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*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)
-          END DO
-       END IF
-    END DO
-       
-    CALL trace_end("compute_caldyn_fast")
-
-  END SUBROUTINE compute_caldyn_fast
-
-  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,nqdyn) ! active scalars
-    REAL(rstd),INTENT(IN)  :: qu(3*iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
-    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
-    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
-    
-    REAL(rstd) :: Ftheta(3*iim*jjm)  ! potential temperature flux
-    REAL(rstd) :: uu_right, uu_lup, uu_ldown
-    INTEGER :: ij,iq,l,kdown
-
-    CALL trace_start("compute_caldyn_Coriolis")
-
-    DO l=ll_begin, ll_end
-       ! compute theta flux
-       DO iq=1,nqdyn
-       !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext      
-             Ftheta(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
-                                  * hflux(ij+u_right,l)
-             Ftheta(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
-                  * hflux(ij+u_lup,l)
-             Ftheta(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
-                  * hflux(ij+u_ldown,l)
-          END DO
-          ! horizontal divergence of fluxes
-       !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             ! dtheta_rhodz =  -div(flux.theta)
-             dtheta_rhodz(ij,l,iq)= &
-                  -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) )
-          END DO
-       END DO
-
-       !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))
-       END DO ! ij
-    END DO ! llm
-
-!!! Compute potential vorticity (Coriolis) contribution to du
-    SELECT CASE(caldyn_conserv)
-
-    CASE(energy) ! energy-conserving TRiSK
-
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             uu_right = &
-                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
-                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
-                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
-                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
-                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
-                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
-                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
-                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
-                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
-                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
-             uu_lup = &
-                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
-                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
-                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
-                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
-                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
-                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
-                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
-                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
-                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
-                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
-             uu_ldown = &
-                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
-                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
-                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
-                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
-                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
-                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
-                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
-                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
-                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
-                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
-             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
-
-    CASE(enstrophy) ! enstrophy-conserving TRiSK
-
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             uu_right = &
-                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
-                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
-                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
-                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
-                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
-                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
-                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
-                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
-                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
-                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
-             uu_lup = &
-                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
-                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
-                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
-                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
-                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
-                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
-                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
-                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
-                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
-                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
-             uu_ldown = &
-                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
-                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
-                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
-                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
-                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
-                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
-                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
-                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
-                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
-                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
-
-             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
-
-    CALL trace_end("compute_caldyn_Coriolis")
-
-  END SUBROUTINE compute_caldyn_Coriolis
-
-  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero)
-    LOGICAL, INTENT(IN) :: zero
-    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(INOUT) :: 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
-       !  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         
-          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)
-       IF(zero) THEN
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end
-             du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
-             du(ij+u_lup,l)   = ne_lup*(berni(ij)-berni(ij+t_lup))
-             du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
-          END DO
-       ELSE
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end
-             du(ij+u_right,l) = du(ij+u_right,l) + &
-                  ne_right*(berni(ij)-berni(ij+t_right))
-             du(ij+u_lup,l)   = du(ij+u_lup,l) + &
-                  ne_lup*(berni(ij)-berni(ij+t_lup))
-             du(ij+u_ldown,l) = du(ij+u_ldown,l) + &
-                  ne_ldown*(berni(ij)-berni(ij+t_ldown))
-          END DO
-       END IF
-    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
-       ! below : "checked" means "formula also valid when kup=kdown (top/bottom)"
-       ! 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)) ! checked
-       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)) ! checked
-          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)) ! checked
-          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)) ! checked
-          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 )
-!          gradPhi2(ij,l) = 0. ! FIXME !!
-
-          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
-    ! FIXME !!
- !   F_el(:,:)=0. 
- !   dPhi(:,:)=0.
- !   dW(:,:)=0.
-
-    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 ) &
-               - .25*( 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)-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)-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)-berni(ij+t_ldown))
-       END DO
-    END DO
-    ! FIXME !!
-    ! du(:,:)=0.
-    ! hflux(:,:)=0.
-
-    CALL trace_end("compute_caldyn_slow_NH")
-
-  END SUBROUTINE compute_caldyn_slow_NH
-
-END MODULE caldyn_kernels_hevi_mod
