!--------------------------------------------------------------------------
   !---------------------------- 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 ----------------------------------
   !--------------------------------------------------------------------------