Changeset 657 for codes/icosagcm/devel/src/kernels_hex/caldyn_solver.k90

Timestamp:

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

Author:

dubos

Message:

devel/hex : updated kernels

File:

• codes/icosagcm/devel/src/kernels_hex/caldyn_solver.k90 (modified) (4 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