!--------------------------------------------------------------------------
   !---------------------------- caldyn_solver ----------------------------------
   IF (ll_begin==1) THEN
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         m_il(ij,1) = .5*rhodz(ij,1)
      END DO
   END IF
   DO l = ll_beginp1, ll_end
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         m_il(ij,l) = .5*(rhodz(ij,l)+rhodz(ij,l-1))
      END DO
   END DO
   IF(ll_endp1==llm+1) THEN
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         m_il(ij,llm+1) = .5*rhodz(ij,llm+1 -1)
      END DO
   END IF
   !
   IF(tau>0) THEN ! solve implicit problem for geopotential
      CALL compute_NH_geopot(tau,PHI_BOT_VAR, rhodz, m_il, theta, W, geopot)
   END IF
   !
   ! Compute pressure (pres) and Exner function (pk)
   ! kappa = R/Cp
   ! 1-kappa = Cv/Cp
   ! Cp/Cv = 1/(1-kappa)
   gamma = 1./(1.-kappa)
   vreff = Rd*Treff/preff ! reference specific volume
   Cvd = cpp-Rd
   DO l = ll_begin, ll_end
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         SELECT CASE(caldyn_thermo)
         CASE(thermo_theta)
            rho_ij = g*rhodz(ij,l)/(geopot(ij,l+1)-geopot(ij,l))
            X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij
            ! kappa.theta.rho = p/exner
            ! => X = (p/p0)/(exner/Cp)
            ! = (p/p0)^(1-kappa)
            pres(ij,l) = preff*(X_ij**gamma) ! pressure
            ! Compute Exner function (needed by compute_caldyn_fast)
            ! other formulae possible if exponentiation is slow
            pk(ij,l) = cpp*((pres(ij,l)/preff)**kappa) ! Exner
         CASE(thermo_entropy)
            rho_ij = g*rhodz(ij,l)/(geopot(ij,l+1)-geopot(ij,l))
            T_ij = Treff*exp( (theta(ij,l,1)+Rd*log(vreff*rho_ij))/Cvd )
            pres(ij,l) = rho_ij*Rd*T_ij
            pk(ij,l) = T_ij
         CASE DEFAULT
            STOP
         END SELECT
      END DO
   END DO
   ! We need a barrier here because we compute pres above and do a vertical difference below
   !$OMP BARRIER
   IF (ll_begin==1) THEN
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         ! Lower BC
         dW(ij,1) = (1./g)*(pbot-rho_bot*(geopot(ij,1)-PHI_BOT(ij))-pres(ij,1)) - m_il(ij,1)
         W(ij,1) = W(ij,1)+tau*dW(ij,1) ! update W
         dPhi(ij,1) = g*g*W(ij,1)/m_il(ij,1)
      END DO
   END IF
   DO l = ll_beginp1, ll_end
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         dW(ij,l) = (1./g)*(pres(ij,l-1)-pres(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)
      END DO
   END DO
   IF(ll_endp1==llm+1) THEN
      !DIR$ SIMD
      DO ij=ij_begin_ext, ij_end_ext
         ! Top BC
         dW(ij,llm+1) = (1./g)*(pres(ij,llm+1 -1)-ptop) - m_il(ij,llm+1)
         W(ij,llm+1) = W(ij,llm+1)+tau*dW(ij,llm+1) ! update W
         dPhi(ij,llm+1) = g*g*W(ij,llm+1)/m_il(ij,llm+1)
      END DO
   END IF
   ! We need a barrier here because we update W 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 du = -0.5*g^2.grad(w^2)
         berni(ij,l) = (-.25*g*g)*((W(ij,l)/m_il(ij,l))**2 + (W(ij,l+1)/m_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
         du(ij+u_right,l) = ne_right*(berni(ij,l)-berni(ij+t_right,l))
         du(ij+u_lup,l) = ne_lup*(berni(ij,l)-berni(ij+t_lup,l))
         du(ij+u_ldown,l) = ne_ldown*(berni(ij,l)-berni(ij+t_ldown,l))
      END DO
   END DO
   !---------------------------- caldyn_solver ----------------------------------
   !--------------------------------------------------------------------------