SUBMODULE(compute_caldyn_mod) compute_caldyn_slow_NH_mod
   USE prec, ONLY: rstd
   IMPLICIT NONE

CONTAINS

   MODULE SUBROUTINE compute_caldyn_slow_NH(u, rhodz, Phi, W, F_el, gradPhi2, w_il, hflux, du, dPhi, dW)
      REAL(rstd), INTENT(IN)  :: u(:, :)   ! prognostic "velocity"
      REAL(rstd), INTENT(IN)  :: rhodz(:, :)  ! rho*dz
      REAL(rstd), INTENT(IN)  :: Phi(:, :)  ! prognostic geopotential
      REAL(rstd), INTENT(IN)  :: W(:, :)    ! prognostic vertical momentum

      REAL(rstd), INTENT(OUT) :: F_el(:, :) ! NH mass flux
      REAL(rstd), INTENT(OUT) :: gradPhi2(:, :) ! grad_Phi**2
      REAL(rstd), INTENT(OUT) :: w_il(:, :) ! Wil/mil

      REAL(rstd), INTENT(OUT) :: hflux(:, :) ! hflux in kg/s
      REAL(rstd), INTENT(OUT) :: du(:, :)
      REAL(rstd), INTENT(OUT) :: dW(:, :)
      REAL(rstd), INTENT(OUT) :: dPhi(:, :)

      CALL compute_caldyn_slow_NH_hex(u, rhodz, Phi, W, F_el, gradPhi2, w_il, hflux, du, dPhi, dW)
   END SUBROUTINE

   SUBROUTINE compute_caldyn_slow_NH_hex(u, rhodz, Phi, W, F_el, gradPhi2, w_il, hflux, du, dPhi, dW)
      USE caldyn_vars_mod, ONLY: dysl
      USE icosa
      USE trace
      USE omp_para, ONLY: ll_begin, ll_end, ll_beginp1, ll_endp1

      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), INTENT(OUT) :: w_il(iim*jjm, llm + 1) ! Wil/mil
      REAL(rstd), INTENT(OUT) :: F_el(3*iim*jjm, llm + 1) ! NH mass flux
      REAL(rstd), INTENT(OUT) :: gradPhi2(iim*jjm, llm + 1) ! grad_Phi**2

      REAL(rstd) :: DePhil(3*iim*jjm, llm + 1) ! grad(Phi)

      INTEGER :: ij, l, kdown, kup
      REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu

      REAL(rstd) :: berni(iim*jjm, llm)  ! Bernoulli function
      REAL(rstd) :: G_el(3*iim*jjm, llm + 1) ! horizontal flux of W
      REAL(rstd) :: v_el(3*iim*jjm, llm + 1)

      REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
      REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W
      REAL(rstd) :: v_el1(3*iim*jjm)

      CALL trace_start("compute_caldyn_slow_NH")

      IF (dysl) THEN

!$OMP BARRIER

         !--------------------------------------------------------------------------
         !---------------------------- caldyn_slow_NH ----------------------------------
         IF (ll_begin == 1) THEN
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               w_il(ij, 1) = 2.*W(ij, 1)/rhodz(ij, 1)
            END DO
         END IF
         DO l = ll_beginp1, ll_end
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               w_il(ij, l) = 2.*W(ij, l)/(rhodz(ij, l - 1) + rhodz(ij, l))
            END DO
         END DO
         IF (ll_endp1 == llm + 1) THEN
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               w_il(ij, llm + 1) = 2.*W(ij, llm + 1)/rhodz(ij, llm)
            END DO
         END IF
         IF (ll_begin == 1) THEN
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               ! 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 do not
               W_el = .5*(W(ij + t_right, 1) + W(ij, 1))
               DePhil(ij + u_right, 1) = ne_right*(Phi(ij + t_right, 1) - Phi(ij, 1))
               F_el(ij + u_right, 1) = DePhil(ij + u_right, 1)*W_el
               W2_el = .5*le_de(ij + u_right)*(W(ij, 1)*w_il(ij, 1) + W(ij + t_right, 1)*w_il(ij + t_right, 1))
               v_el(ij + u_right, 1) = .5*le_de(ij + u_right)*(u(ij + u_right, 1) + u(ij + u_right, 1)) ! checked
               G_el(ij + u_right, 1) = v_el(ij + u_right, 1)*W_el - DePhil(ij + u_right, 1)*W2_el
               ! 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 do not
               W_el = .5*(W(ij + t_lup, 1) + W(ij, 1))
               DePhil(ij + u_lup, 1) = ne_lup*(Phi(ij + t_lup, 1) - Phi(ij, 1))
               F_el(ij + u_lup, 1) = DePhil(ij + u_lup, 1)*W_el
               W2_el = .5*le_de(ij + u_lup)*(W(ij, 1)*w_il(ij, 1) + W(ij + t_lup, 1)*w_il(ij + t_lup, 1))
               v_el(ij + u_lup, 1) = .5*le_de(ij + u_lup)*(u(ij + u_lup, 1) + u(ij + u_lup, 1)) ! checked
               G_el(ij + u_lup, 1) = v_el(ij + u_lup, 1)*W_el - DePhil(ij + u_lup, 1)*W2_el
               ! 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 do not
               W_el = .5*(W(ij + t_ldown, 1) + W(ij, 1))
               DePhil(ij + u_ldown, 1) = ne_ldown*(Phi(ij + t_ldown, 1) - Phi(ij, 1))
               F_el(ij + u_ldown, 1) = DePhil(ij + u_ldown, 1)*W_el
               W2_el = .5*le_de(ij + u_ldown)*(W(ij, 1)*w_il(ij, 1) + W(ij + t_ldown, 1)*w_il(ij + t_ldown, 1))
               v_el(ij + u_ldown, 1) = .5*le_de(ij + u_ldown)*(u(ij + u_ldown, 1) + u(ij + u_ldown, 1)) ! checked
               G_el(ij + u_ldown, 1) = v_el(ij + u_ldown, 1)*W_el - DePhil(ij + u_ldown, 1)*W2_el
            END DO
         END IF
         DO l = ll_beginp1, ll_end
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               ! 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 do not
               W_el = .5*(W(ij + t_right, l) + W(ij, 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, l) = .5*le_de(ij + u_right)*(u(ij + u_right, l) + u(ij + u_right, l - 1)) ! checked
               G_el(ij + u_right, l) = v_el(ij + u_right, l)*W_el - DePhil(ij + u_right, l)*W2_el
               ! 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 do not
               W_el = .5*(W(ij + t_lup, l) + W(ij, 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, l) = .5*le_de(ij + u_lup)*(u(ij + u_lup, l) + u(ij + u_lup, l - 1)) ! checked
               G_el(ij + u_lup, l) = v_el(ij + u_lup, l)*W_el - DePhil(ij + u_lup, l)*W2_el
               ! 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 do not
               W_el = .5*(W(ij + t_ldown, l) + W(ij, 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, l) = .5*le_de(ij + u_ldown)*(u(ij + u_ldown, l) + u(ij + u_ldown, l - 1)) ! checked
               G_el(ij + u_ldown, l) = v_el(ij + u_ldown, l)*W_el - DePhil(ij + u_ldown, l)*W2_el
            END DO
         END DO
         IF (ll_endp1 == llm + 1) THEN
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               ! 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 do not
               W_el = .5*(W(ij + t_right, llm + 1) + W(ij, llm + 1))
               DePhil(ij + u_right, llm + 1) = ne_right*(Phi(ij + t_right, llm + 1) - Phi(ij, llm + 1))
               F_el(ij + u_right, llm + 1) = DePhil(ij + u_right, llm + 1)*W_el
               W2_el = .5*le_de(ij + u_right)* &
               &       (W(ij, llm + 1)*w_il(ij, llm + 1) + W(ij + t_right, llm + 1)*w_il(ij + t_right, llm + 1))
               v_el(ij + u_right, llm + 1) = .5*le_de(ij + u_right)*(u(ij + u_right, llm) + u(ij + u_right, llm)) ! checked
               G_el(ij + u_right, llm + 1) = v_el(ij + u_right, llm + 1)*W_el - DePhil(ij + u_right, llm + 1)*W2_el
               ! 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 do not
               W_el = .5*(W(ij + t_lup, llm + 1) + W(ij, llm + 1))
               DePhil(ij + u_lup, llm + 1) = ne_lup*(Phi(ij + t_lup, llm + 1) - Phi(ij, llm + 1))
               F_el(ij + u_lup, llm + 1) = DePhil(ij + u_lup, llm + 1)*W_el
               W2_el = .5*le_de(ij + u_lup)* &
               &       (W(ij, llm + 1)*w_il(ij, llm + 1) + W(ij + t_lup, llm + 1)*w_il(ij + t_lup, llm + 1))
               v_el(ij + u_lup, llm + 1) = .5*le_de(ij + u_lup)*(u(ij + u_lup, llm) + u(ij + u_lup, llm)) ! checked
               G_el(ij + u_lup, llm + 1) = v_el(ij + u_lup, llm + 1)*W_el - DePhil(ij + u_lup, llm + 1)*W2_el
               ! 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 do not
               W_el = .5*(W(ij + t_ldown, llm + 1) + W(ij, llm + 1))
               DePhil(ij + u_ldown, llm + 1) = ne_ldown*(Phi(ij + t_ldown, llm + 1) - Phi(ij, llm + 1))
               F_el(ij + u_ldown, llm + 1) = DePhil(ij + u_ldown, llm + 1)*W_el
               W2_el = .5*le_de(ij + u_ldown)* &
               &       (W(ij, llm + 1)*w_il(ij, llm + 1) + W(ij + t_ldown, llm + 1)*w_il(ij + t_ldown, llm + 1))
               v_el(ij + u_ldown, llm + 1) = .5*le_de(ij + u_ldown)*(u(ij + u_ldown, llm) + u(ij + u_ldown, llm)) ! checked
               G_el(ij + u_ldown, llm + 1) = v_el(ij + u_ldown, llm + 1)*W_el - DePhil(ij + u_ldown, llm + 1)*W2_el
            END DO
         END IF
         DO l = ll_begin, ll_endp1
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               ! compute GradPhi2, dPhi, dW
               gPhi2 = 0.
               dP = 0.
               divG = 0
               gPhi2 = gPhi2 + le_de(ij + u_rup)*DePhil(ij + u_rup, l)**2
               dP = dP + le_de(ij + u_rup)*DePhil(ij + u_rup, l)*v_el(ij + u_rup, l)
               divG = divG + ne_rup*G_el(ij + u_rup, l) ! -div(G_el), G_el already has le_de
               gPhi2 = gPhi2 + le_de(ij + u_lup)*DePhil(ij + u_lup, l)**2
               dP = dP + le_de(ij + u_lup)*DePhil(ij + u_lup, l)*v_el(ij + u_lup, l)
               divG = divG + ne_lup*G_el(ij + u_lup, l) ! -div(G_el), G_el already has le_de
               gPhi2 = gPhi2 + le_de(ij + u_left)*DePhil(ij + u_left, l)**2
               dP = dP + le_de(ij + u_left)*DePhil(ij + u_left, l)*v_el(ij + u_left, l)
               divG = divG + ne_left*G_el(ij + u_left, l) ! -div(G_el), G_el already has le_de
               gPhi2 = gPhi2 + le_de(ij + u_ldown)*DePhil(ij + u_ldown, l)**2
               dP = dP + le_de(ij + u_ldown)*DePhil(ij + u_ldown, l)*v_el(ij + u_ldown, l)
               divG = divG + ne_ldown*G_el(ij + u_ldown, l) ! -div(G_el), G_el already has le_de
               gPhi2 = gPhi2 + le_de(ij + u_rdown)*DePhil(ij + u_rdown, l)**2
               dP = dP + le_de(ij + u_rdown)*DePhil(ij + u_rdown, l)*v_el(ij + u_rdown, l)
               divG = divG + ne_rdown*G_el(ij + u_rdown, l) ! -div(G_el), G_el already has le_de
               gPhi2 = gPhi2 + le_de(ij + u_right)*DePhil(ij + u_right, l)**2
               dP = dP + le_de(ij + u_right)*DePhil(ij + u_right, l)*v_el(ij + u_right, l)
               divG = divG + ne_right*G_el(ij + u_right, l) ! -div(G_el), G_el already has le_de
               gradPhi2(ij, l) = 1./(2.*Ai(ij))*gPhi2
               dPhi(ij, l) = gradPhi2(ij, l)*w_il(ij, l) - 1./(2.*Ai(ij))*dP
               dW(ij, l) = (-1./Ai(ij))*divG
            END DO
         END DO
         ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
         !$OMP BARRIER
         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               ! Compute berni at scalar points
               u2 = 0.
               u2 = u2 + le_de(ij + u_rup)*u(ij + u_rup, l)**2
               u2 = u2 + le_de(ij + u_lup)*u(ij + u_lup, l)**2
               u2 = u2 + le_de(ij + u_left)*u(ij + u_left, l)**2
               u2 = u2 + le_de(ij + u_ldown)*u(ij + u_ldown, l)**2
               u2 = u2 + le_de(ij + u_rdown)*u(ij + u_rdown, l)**2
               u2 = u2 + le_de(ij + u_right)*u(ij + u_right, l)**2
               berni(ij, l) = 1./(4.*Ai(ij))*u2 - .25*(gradPhi2(ij, l)*w_il(ij, l)**2 + gradPhi2(ij, l + 1)*w_il(ij, l + 1)**2)
            END DO
         END DO
         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               ! Compute mass flux and grad(berni)
               uu = .5*(rhodz(ij, l) + rhodz(ij + t_right, l))*u(ij + u_right, l) - &
               &    .5*(F_el(ij + u_right, l) + F_el(ij + u_right, l + 1))
               hflux(ij + u_right, l) = le_de(ij + u_right)*uu
               du(ij + u_right, l) = ne_right*(berni(ij, l) - berni(ij + t_right, l))
               ! Compute mass flux and grad(berni)
               uu = .5*(rhodz(ij, l) + rhodz(ij + t_lup, l))*u(ij + u_lup, l) - &
               &    .5*(F_el(ij + u_lup, l) + F_el(ij + u_lup, l + 1))
               hflux(ij + u_lup, l) = le_de(ij + u_lup)*uu
               du(ij + u_lup, l) = ne_lup*(berni(ij, l) - berni(ij + t_lup, l))
               ! Compute mass flux and grad(berni)
               uu = .5*(rhodz(ij, l) + rhodz(ij + t_ldown, l))*u(ij + u_ldown, l) - &
               &    .5*(F_el(ij + u_ldown, l) + F_el(ij + u_ldown, l + 1))
               hflux(ij + u_ldown, l) = le_de(ij + u_ldown)*uu
               du(ij + u_ldown, l) = ne_ldown*(berni(ij, l) - berni(ij + t_ldown, l))
            END DO
         END DO
         !---------------------------- caldyn_slow_NH ----------------------------------
         !--------------------------------------------------------------------------
!$OMP BARRIER

      ELSE

         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_el1(ij + u_right) = .5*le_de(ij + u_right)*(u(ij + u_right, kup) + u(ij + u_right, kdown)) ! checked
               G_el1(ij + u_right) = v_el1(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_el1(ij + u_lup) = .5*le_de(ij + u_lup)*(u(ij + u_lup, kup) + u(ij + u_lup, kdown)) ! checked
               G_el1(ij + u_lup) = v_el1(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_el1(ij + u_ldown) = .5*le_de(ij + u_ldown)*(u(ij + u_ldown, kup) + u(ij + u_ldown, kdown)) ! checked
               G_el1(ij + u_ldown) = v_el1(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_el1(ij + u_right) + & ! -v.gradPhi,
                              DePhil(ij + u_rup, l)*v_el1(ij + u_rup) + & ! v_el already has le_de
                              DePhil(ij + u_lup, l)*v_el1(ij + u_lup) + &
                              DePhil(ij + u_left, l)*v_el1(ij + u_left) + &
                              DePhil(ij + u_ldown, l)*v_el1(ij + u_ldown) + &
                              DePhil(ij + u_rdown, l)*v_el1(ij + u_rdown))

               dW(ij, l) = -1./Ai(ij)*( & ! -div(G_el),
                           ne_right*G_el1(ij + u_right) + & ! G_el already has le_de
                           ne_rup*G_el1(ij + u_rup) + &
                           ne_lup*G_el1(ij + u_lup) + &
                           ne_left*G_el1(ij + u_left) + &
                           ne_ldown*G_el1(ij + u_ldown) + &
                           ne_rdown*G_el1(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
               berni1(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*(berni1(ij) - berni1(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*(berni1(ij) - berni1(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*(berni1(ij) - berni1(ij + t_ldown))
            END DO
         END DO

      END IF ! dysl

      CALL trace_end("compute_caldyn_slow_NH")

   END SUBROUTINE compute_caldyn_slow_NH_hex

END SUBMODULE compute_caldyn_slow_NH_mod
