!--------------------------------------------------------------------------
   !---------------------------- caldyn_vert_NH ----------------------------------
   IF (ll_begin==1) THEN
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         w_ij = ( W(ij,1)+.5*W(ij,1 +1) )/mass(ij,1)
         wflux_ij = .5*(wflux(ij,1)+wflux(ij,1 +1))
         W_etadot(ij,1) = wflux_ij*w_ij
         eta_dot(ij,1) = wflux_ij / mass(ij,1)
         wcov(ij,1) = w_ij*(geopot(ij,1 +1)-geopot(ij,1))
      END DO
   END IF
   DO l = ll_beginp1, ll_endm1
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         w_ij = .5*( W(ij,l)+W(ij,l+1) )/mass(ij,l)
         wflux_ij = .5*(wflux(ij,l)+wflux(ij,l+1))
         W_etadot(ij,l) = wflux_ij*w_ij
         eta_dot(ij,l) = wflux_ij / mass(ij,l)
         wcov(ij,l) = w_ij*(geopot(ij,l+1)-geopot(ij,l))
      END DO
   END DO
   IF(ll_end==llm) THEN
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         w_ij = ( .5*W(ij,llm)+W(ij,llm+1) )/mass(ij,llm)
         wflux_ij = .5*(wflux(ij,llm)+wflux(ij,llm+1))
         W_etadot(ij,llm) = wflux_ij*w_ij
         eta_dot(ij,llm) = wflux_ij / mass(ij,llm)
         wcov(ij,llm) = w_ij*(geopot(ij,llm+1)-geopot(ij,llm))
      END DO
   END IF
   ! add NH term to du
   DO l = ll_begin, ll_end
      !DIR$ SIMD
      DO ij=ij_begin, ij_end
         du(ij+u_right,l) = du(ij+u_right,l) - .5*(wcov(ij,l)+wcov(ij+t_right,l))*ne_right*(eta_dot(ij+t_right,l)-eta_dot(ij,l))
         du(ij+u_lup,l) = du(ij+u_lup,l) - .5*(wcov(ij,l)+wcov(ij+t_lup,l))*ne_lup*(eta_dot(ij+t_lup,l)-eta_dot(ij,l))
         du(ij+u_ldown,l) = du(ij+u_ldown,l) - .5*(wcov(ij,l)+wcov(ij+t_ldown,l))*ne_ldown*(eta_dot(ij+t_ldown,l)-eta_dot(ij,l))
      END DO
   END DO
   ! add NH terms to dW, dPhi
   ! FIXME : TODO top and bottom
   DO l = ll_beginp1, ll_end
      !DIR$ SIMD
      DO ij=ij_begin, ij_end
         dPhi(ij,l)=dPhi(ij,l)-wflux(ij,l)*(geopot(ij,l+1)-geopot(ij,l-1))/(mass(ij,l-1)+mass(ij,l))
      END DO
   END DO
   ! We need a barrier here because we compute W_etadot above and do a vertical difference below
   !$OMP BARRIER
   IF (ll_begin==1) THEN
      !DIR$ SIMD
      DO ij=ij_begin, ij_end
         dW(ij,1) = dW(ij,1) - W_etadot(ij,1)
      END DO
   END IF
   DO l = ll_beginp1, ll_end
      !DIR$ SIMD
      DO ij=ij_begin, ij_end
         dW(ij,l) = dW(ij,l) + W_etadot(ij,l-1) - W_etadot(ij,l)
      END DO
   END DO
   IF(ll_endp1==llm+1) THEN
      !DIR$ SIMD
      DO ij=ij_begin, ij_end
         dW(ij,llm+1) = dW(ij,llm+1) + W_etadot(ij,llm+1 -1)
      END DO
   END IF
   !---------------------------- caldyn_vert_NH ----------------------------------
   !--------------------------------------------------------------------------