Changeset 657 for codes/icosagcm/devel/src/kernels_hex

Timestamp:

12/30/17 01:56:49 (6 years ago)

Author:

dubos

Message:

devel/hex : updated kernels

Location:

codes/icosagcm/devel/src/kernels_hex

Files:

• caldyn_mil.k90 (added)
• caldyn_solver.k90 (modified) (4 diffs)
• compute_NH_geopot.k90 (modified) (2 diffs)

codes/icosagcm/devel/src/kernels_hex/caldyn_solver.k90

@@ -1 +1 @@
    !--------------------------------------------------------------------------
    !---------------------------- 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)
@@ -30 +8 @@
    gamma = 1./(1.-kappa)
    vreff = Rd*Treff/preff ! reference specific volume
-   Cvd = cpp-Rd
+   Cvd = 1./(cpp-Rd)
+   Rd_preff = kappa*cpp/preff
    DO l = ll_begin, ll_end
       !DIR$ SIMD
@@ -36 +15 @@
          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
+            rho_ij = 1./(geopot(ij,l+1)-geopot(ij,l))
+            rho_ij = rho_ij*g*rhodz(ij,l)
+            X_ij = Rd_preff*theta(ij,l,1)*rho_ij
             ! kappa.theta.rho = p/exner
             ! => X = (p/p0)/(exner/Cp)
@@ -46 +26 @@
             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 )
+            rho_ij = 1./(geopot(ij,l+1)-geopot(ij,l))
+            rho_ij = rho_ij*g*rhodz(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

codes/icosagcm/devel/src/kernels_hex/compute_NH_geopot.k90

@@ -4 +4 @@
    g2=g*g
    gm2 = 1./g2
+   vreff = Treff*cpp/preff*kappa
    gamma = 1./(1.-kappa)
    !$OMP BARRIER
@@ -30 +31 @@
             DO ij=ij_omp_begin_ext,ij_omp_end_ext
                rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l))
-               X_ij = Treff*exp(theta(ij,l)/cpp) ! theta = Tref.exp(s/Cp)
-               X_ij = (cpp/preff)*kappa*X_ij*rho_ij
-               p_ik(ij,l) = preff*(X_ij**gamma)
+               X_ij = log(vreff*rho_ij) + theta(ij,l)/cpp
+               p_ik(ij,l) = preff*exp(X_ij*gamma)
                c2_mik = gamma*p_ik(ij,l)/(rho_ij*m_ik(ij,l)) ! c^2 = gamma*R*T = gamma*p/rho
                A_ik(ij,l) = c2_mik*(tau/g*rho_ij)**2