!--------------------------------------------------------------------------
   !---------------------------- caldyn_solver ----------------------------------
   !$OMP DO SCHEDULE(STATIC)
   DO ij = 1, primal_num
      l=1
      m_il(l,ij) = .5*rhodz(l,ij)
      DO l = 2, llm
         m_il(l,ij) = .5*(rhodz(l,ij)+rhodz(l-1,ij))
      END DO
      l=llm+1
      m_il(l,ij) = .5*rhodz(l-1,ij)
   END DO
   !$OMP END DO
   !
   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
   SELECT CASE(caldyn_thermo)
   CASE(thermo_theta)
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         DO l = 1, llm
            rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij))
            X_ij = (cpp/preff)*kappa*theta(l,ij,1)*rho_ij
            ! kappa.theta.rho = p/exner
            ! => X = (p/p0)/(exner/Cp)
            ! = (p/p0)^(1-kappa)
            pres(l,ij) = preff*(X_ij**gamma) ! pressure
            ! Compute Exner function (needed by compute_caldyn_fast)
            ! other formulae possible if exponentiation is slow
            pk(l,ij) = cpp*((pres(l,ij)/preff)**kappa) ! Exner
         END DO
      END DO
      !$OMP END DO
   CASE(thermo_entropy)
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         DO l = 1, llm
            rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij))
            T_ij = Treff*exp( (theta(l,ij,1)+Rd*log(vreff*rho_ij))/Cvd )
            pres(l,ij) = rho_ij*Rd*T_ij
            pk(l,ij) = T_ij
         END DO
      END DO
      !$OMP END DO
   CASE DEFAULT
      STOP
   END SELECT
   ! We need a barrier here because we compute pres above and do a vertical difference below
   !$OMP BARRIER
   !$OMP DO SCHEDULE(STATIC)
   DO ij = 1, primal_num
      l=1
      ! Lower BC
      dW(l,ij) = (1./g)*(pbot-rho_bot*(geopot(l,ij)-PHI_BOT(ij))-pres(l,ij)) - m_il(l,ij)
      W(l,ij) = W(l,ij)+tau*dW(l,ij) ! update W
      dPhi(l,ij) = g*g*W(l,ij)/m_il(l,ij)
      DO l = 2, llm
         dW(l,ij) = (1./g)*(pres(l-1,ij)-pres(l,ij)) - m_il(l,ij)
         W(l,ij) = W(l,ij)+tau*dW(l,ij) ! update W
         dPhi(l,ij) = g*g*W(l,ij)/m_il(l,ij)
      END DO
      l=llm+1
      ! Top BC
      dW(l,ij) = (1./g)*(pres(l-1,ij)-ptop) - m_il(l,ij)
      W(l,ij) = W(l,ij)+tau*dW(l,ij) ! update W
      dPhi(l,ij) = g*g*W(l,ij)/m_il(l,ij)
   END DO
   !$OMP END DO
   ! We need a barrier here because we update W above and do a vertical average below
   !$OMP BARRIER
   !$OMP DO SCHEDULE(STATIC)
   DO ij = 1, primal_num
      DO l = 1, llm
         ! compute du = -0.5*g^2.grad(w^2)
         berni(l,ij) = (-.25*g*g)*((W(l,ij)/m_il(l,ij))**2 + (W(l+1,ij)/m_il(l+1,ij))**2 )
      END DO
   END DO
   !$OMP END DO
   !$OMP DO SCHEDULE(STATIC)
   DO edge = 1, edge_num
      ij_left = left(edge)
      ij_right = right(edge)
      DO l = 1, llm
         du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right))
      END DO
   END DO
   !$OMP END DO
   !---------------------------- caldyn_solver ----------------------------------
   !--------------------------------------------------------------------------