source: codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90 @ 935

   !--------------------------------------------------------------------------
   !---------------------------- caldyn_fast ----------------------------------
   !
   SELECT CASE(caldyn_thermo)
   CASE(thermo_boussinesq)
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         !DIR$ SIMD
         DO l = 1, llm
            berni(l,ij) = pk(l,ij)
            ! from now on pk contains the vertically-averaged geopotential
            pk(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij))
         END DO
      END DO
      !$OMP END DO
   CASE(thermo_theta)
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         !DIR$ SIMD
         DO l = 1, llm
            Phi_ik = .5*(geopot(l,ij)+geopot(l+1,ij))
            berni(l,ij) = Phi_ik
         END DO
      END DO
      !$OMP END DO
   CASE(thermo_entropy)
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         !DIR$ SIMD
         DO l = 1, llm
            Phi_ik = .5*(geopot(l,ij)+geopot(l+1,ij))
            berni(l,ij) = Phi_ik + pk(l,ij)*(cpp-theta(l,ij,1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
         END DO
      END DO
      !$OMP END DO
   CASE(thermo_variable_Cp)
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         !DIR$ SIMD
         DO l = 1, llm
            Phi_ik = .5*(geopot(l,ij)+geopot(l+1,ij))
            ! thermodynamics with variable Cp
            ! Cp(T) = Cp0 * (T/T0)^nu
            ! => h = Cp(T).T/(nu+1)
            cp_ik = cpp*(pk(ij,l)/Treff)**nu
            berni(l,ij) = Phi_ik + pk(l,ij)*(cp_ik/(nu+1.)-theta(l,ij,1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
         END DO
      END DO
      !$OMP END DO
   CASE DEFAULT
      PRINT *, 'Unsupported value of caldyn_thermo : ',caldyn_thermo ! FIXME
      STOP
   END SELECT
   !
   !$OMP DO SCHEDULE(STATIC)
   DO edge = 1, edge_num
      ij_left = left(edge)
      ij_right = right(edge)
      !DIR$ SIMD
      DO l = 1, llm
         due = .5*(theta(l,ij_left,1)+theta(l,ij_right,1))*(pk(l,ij_right)-pk(l,ij_left)) + berni(l,ij_right)-berni(l,ij_left)
         du(l,edge) = du(l,edge) - 1.*due
         u(l,edge) = u(l,edge) + tau*du(l,edge)
      END DO
   END DO
   !$OMP END DO
   !
   !---------------------------- caldyn_fast ----------------------------------
   !--------------------------------------------------------------------------