Changeset 659

Timestamp:

12/30/17 02:01:51 (6 years ago)

Author:

dubos

Message:

devel/codegen : update NH kernels

File:

• codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin (modified) (6 diffs)

codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin

@@ -14 +14 @@
   g2=g*g
   gm2 = 1./g2
+  vreff = Treff*cpp/preff*kappa
   gamma = 1./(1.-kappa)
     
@@ -37 +38 @@
     CASE(thermo_entropy)
       {% call() compute_p_and_Aik() %}
-        X_ij = Treff*exp(theta(CELL)/cpp) ! theta = Tref.exp(s/Cp)
-        X_ij = (cpp/preff)*kappa*X_ij*rho_ij
-        p_ik(CELL) = preff*(X_ij**gamma)
+        X_ij = log(vreff*rho_ij) + theta(CELL)/cpp
+        p_ik(CELL) = preff*exp(X_ij*gamma)
         c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
       {% endcall %}
@@ -125 +125 @@
 END_BLOCK
 
-KERNEL(caldyn_solver)
+KERNEL(caldyn_mil)
   FORALL_CELLS_EXT('1', 'llm+1')
     ON_PRIMAL
@@ -133 +133 @@
     END_BLOCK
   END_BLOCK
-  !
-  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
+END_BLOCK
+
+KERNEL(caldyn_solver)
   ! 
   ! Compute pressure (pres) and Exner function (pk)
@@ -142 +141 @@
   ! 1-kappa = Cv/Cp
   ! Cp/Cv = 1/(1-kappa)
-  gamma = 1./(1.-kappa)
+  gamma    = 1./(1.-kappa)
   vreff    = Rd*Treff/preff ! reference specific volume
-  Cvd   = cpp-Rd
+  Cvd      = 1./(cpp-Rd)
+  Rd_preff = kappa*cpp/preff
   FORALL_CELLS_EXT()
     SELECT CASE(caldyn_thermo)
     CASE(thermo_theta)
       ON_PRIMAL
-        rho_ij = g*rhodz(CELL)/(geopot(UP(CELL))-geopot(CELL))
-        X_ij = (cpp/preff)*kappa*theta(CELL,1)*rho_ij
+        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
+        rho_ij = rho_ij*g*rhodz(CELL)
+        X_ij = Rd_preff*theta(CELL,1)*rho_ij
         ! kappa.theta.rho = p/exner
         ! => X = (p/p0)/(exner/Cp)
@@ -161 +162 @@
     CASE(thermo_entropy)
       ON_PRIMAL
-        rho_ij = g*rhodz(CELL)/(geopot(UP(CELL))-geopot(CELL))
-        T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))/Cvd )
+        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
+        rho_ij = rho_ij*g*rhodz(CELL)
+        T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))*Cvd )
         pres(CELL) = rho_ij*Rd*T_ij
         pk(CELL)   = T_ij