Changeset 836 for codes/icosagcm/devel

Timestamp:

05/03/19 13:25:21 (5 years ago)

Author:

dubos

Message:

devel : Cp(T) thermodynamics (TBC)

Location:

codes/icosagcm/devel

Files:

• Python/src/kernels_caldyn_hydro.jin (modified) (1 diff)
• src/base/earth_const.f90 (modified) (1 diff)
• src/dynamics/caldyn_kernels_base.F90 (modified) (2 diffs)
• src/kernels_hex/compute_geopot.k90 (modified) (1 diff)
• src/kernels_unst/compute_geopot.k90 (modified) (1 diff)
• src/unstructured/caldyn_unstructured.F90 (modified) (2 diffs)
• src/unstructured/data_unstructured.F90 (modified) (4 diffs)

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

@@ -65 +65 @@
       gv = (g*Rd)*temp_ik/p_ik   ! specific volume v = Rd*T/p
     END_GEOPOT
+  CASE(thermo_variable_Cp)
+    ! thermodynamics with variable Cp
+    !      Cp.dT = dh = Tds + vdp
+    !              pv = RT
+    ! =>           ds = (dh+v.dp)/T = Cp.dT/T - R dp/p
+    !           Cp(T) = Cp0 * (T/T0)^nu
+    ! =>       s(p,T) = Cp(T)/nu - R log(p/preff) 
+    !               h = Cp(T).T/(nu+1)
+    BALANCE( rhodz(THECELL) )
+    GEOPOT('temp_ik')
+      Cp_ik = nu*( theta(CELL,1) + Rd*log(p_ik/preff) )
+      temp_ik = Treff* (Cp_ik/cpp)**(1./nu)
+      gv = (g*Rd)*temp_ik/p_ik   ! specific volume v = Rd*T/p
+    END_GEOPOT
   CASE(thermo_moist)
     BALANCE( rhodz(THECELL)*(1.+theta(THECELL,2)) )

codes/icosagcm/devel/src/base/earth_const.f90

@@ -7 +7 @@
 ! hence the variables set here are not THREADPRIVATE
 
-  REAL(rstd) :: g=9.80616
-  REAL(rstd) :: radius=6.37122E6
+! Some variable are BIND(C) to be accessible from Python
+
+  REAL(rstd), BIND(C) :: g=9.80616    ! acceleration of gravity assumed constant
+  REAL(rstd):: radius=6.37122E6
   REAL(rstd) :: omega=7.292E-5
   REAL(rstd) :: scale_factor=1.
   REAL(rstd),PARAMETER :: daysec=86400
 
-  REAL(rstd) :: kappa=0.2857143
-  REAL(rstd) :: cpp=1004.70885
-  REAL(rstd) :: nu=0.35 ! exponent in variable-Cp law Cp=cpp*(T/Treff)^nu
-  REAL(rstd) :: cppv=1860.
-  REAL(rstd) :: Rv=461.5
-  REAL(rstd) :: Treff=273.
-  REAL(rstd) :: preff=101325.
-  REAL(rstd) :: pa=50000. ! default value set to preff/2 by disvert_std
+  REAL(rstd), BIND(C) ::   Treff=273.,     &   ! reference temperature, used in definition of entropy
+                           preff=101325.,  &   ! reference pressure, used in definition of entropy and theta
+                           cpp=1004.70885, &   ! specific heat at constant pressure at t=Treff for dry air
+                           cppv=1860.,     &   ! specific heat at constant pressure at t=Treff for water vapor
+                           Rd,             &   ! specific perfect gas constant for dry air
+                           Rv=461.5            ! specific perfect gas constant for water vapor
+
+  REAL(rstd) :: nu=0.35     ! exponent in variable-Cp law Cp=cpp*(T/Treff)^nu
+  REAL(rstd) :: pa=50000.   ! default value set to preff/2 by disvert_std
   REAL(rstd) :: scale_height=8000. ! atmospheric scale height (m)
   REAL(rstd) :: gas_constant = 8.3144621 
-  REAL(rstd) :: Rd, mu             ! specific perfect gas constant and molar mass (?)
+  REAL(rstd) :: mu             ! molar mass (?)
 
-  INTEGER, PUBLIC, PARAMETER :: thermo_none=-99, thermo_theta=1, thermo_entropy=2, thermo_variable_Cp=3, &
+  REAL(rstd) :: kappa=0.2857143
+
+  INTEGER, PARAMETER :: thermo_none=-99, thermo_theta=1, thermo_entropy=2, thermo_variable_Cp=3, &
        thermo_moist=4, thermo_boussinesq=5, thermo_dry=10, thermo_fake_moist=11, thermo_moist_debug=100
-  LOGICAL, PUBLIC :: boussinesq
-  INTEGER, PUBLIC :: caldyn_thermo, physics_thermo
+  LOGICAL :: boussinesq
+  INTEGER, BIND(C) :: caldyn_thermo, physics_thermo
 
 CONTAINS

codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90

@@ -24 +24 @@
 
     INTEGER :: i,j,ij,l
-    REAL(rstd) :: Rd, p_ik, exner_ik, temp_ik, qv, chi, Rmix, gv
+    REAL(rstd) :: p_ik, exner_ik, Cp_ik, temp_ik, qv, chi, Rmix, gv
     INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext
 
@@ -32 +32 @@
 
     CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
-
-    Rd = kappa*cpp
 
     IF(dysl_geopot) THEN

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

@@ -82 +82 @@
          END DO
       END DO
+   CASE(thermo_variable_Cp)
+      ! thermodynamics with variable Cp
+      ! Cp.dT = dh = Tds + vdp
+      ! pv = RT
+      ! => ds = (dh+v.dp)/T = Cp.dT/T - R dp/p
+      ! Cp(T) = Cp0 * (T/T0)^nu
+      ! => s(p,T) = Cp(T)/nu - R log(p/preff)
+      ! h = Cp(T).T/(nu+1)
+      !DIR$ SIMD
+      DO ij=ij_omp_begin_ext,ij_omp_end_ext
+         pk(ij,llm) = ptop + .5*g* rhodz(ij,llm)
+      END DO
+      DO l = llm-1,1,-1
+         !DIR$ SIMD
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            pk(ij,l) = pk(ij,l+1) + (.5*g)*( rhodz(ij,l) + rhodz(ij,l+1) )
+         END DO
+      END DO
+      IF(caldyn_eta == eta_lag) THEN
+         !DIR$ SIMD
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            ps(ij) = pk(ij,1) + .5*g* rhodz(ij,1)
+         END DO
+      END IF
+      DO l = 1,llm
+         !DIR$ SIMD
+         DO ij=ij_omp_begin_ext,ij_omp_end_ext
+            p_ik = pk(ij,l)
+            Cp_ik = nu*( theta(ij,l,1) + Rd*log(p_ik/preff) )
+            temp_ik = Treff* (Cp_ik/cpp)**(1./nu)
+            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
+            pk(ij,l) = temp_ik
+            geopot(ij,l+1) = geopot(ij,l) + gv*rhodz(ij,l)
+         END DO
+      END DO
    CASE(thermo_moist)
       !DIR$ SIMD

codes/icosagcm/devel/src/kernels_unst/compute_geopot.k90

@@ -70 +70 @@
       END DO
       !$OMP END DO
+   CASE(thermo_variable_Cp)
+      ! thermodynamics with variable Cp
+      ! Cp.dT = dh = Tds + vdp
+      ! pv = RT
+      ! => ds = (dh+v.dp)/T = Cp.dT/T - R dp/p
+      ! Cp(T) = Cp0 * (T/T0)^nu
+      ! => s(p,T) = Cp(T)/nu - R log(p/preff)
+      ! h = Cp(T).T/(nu+1)
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)
+         DO l = llm-1,1,-1
+            pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
+         END DO
+         IF(caldyn_eta == eta_lag) THEN
+            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)
+         END IF
+      END DO
+      !$OMP END DO
+      !$OMP DO SCHEDULE(STATIC)
+      DO ij=1,primal_num
+         DO l = 1,llm
+            p_ik = pk(l,ij)
+            Cp_ik = nu*( theta(l,ij,1) + Rd*log(p_ik/preff) )
+            temp_ik = Treff* (Cp_ik/cpp)**(1./nu)
+            gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
+            pk(l,ij) = temp_ik
+            geopot(l+1,ij) = geopot(l,ij) + gv*rhodz(l,ij)
+         END DO
+      END DO
+      !$OMP END DO
    CASE(thermo_moist)
       !$OMP DO SCHEDULE(STATIC)

codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90

@@ -2 +2 @@
   USE ISO_C_BINDING
   USE OMP_LIB
+  USE earth_const
   USE data_unstructured_mod
   IMPLICIT NONE
@@ -119 +120 @@
   FIELD_PS     :: ps        ! OUT
   DECLARE_INDICES
-  NUM :: gdz, ke, uu, chi, gv, exner_ik, temp_ik, p_ik, qv, Rmix
+  NUM :: gdz, ke, uu, chi, gv, exner_ik, temp_ik, p_ik, qv, Rmix, Cp_ik
   START_TRACE(id_geopot, 3,0,3)
 #include "../kernels_unst/compute_geopot.k90"

codes/icosagcm/devel/src/unstructured/data_unstructured.F90

@@ -1 +1 @@
 MODULE data_unstructured_mod
   USE ISO_C_BINDING
+  USE earth_const, ONLY : thermo_theta
   USE mpipara, ONLY : is_mpi_master
   USE grid_param, ONLY : llm, nqdyn
@@ -13 +14 @@
 
   INTEGER, PARAMETER :: eta_mass=1, eta_lag=2, &
-       thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4, &
+!       thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4, &
        caldyn_vert_cons=1, max_nb_stage=5
-  INDEX,  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, &
+  INDEX,  BIND(C) :: caldyn_eta=eta_lag, &
        caldyn_vert_variant=caldyn_vert_cons, nb_threads=0, nb_stage=0
+!  INDEX,  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, &
+!       caldyn_vert_variant=caldyn_vert_cons, nb_threads=0, nb_stage=0
   LOGICAL(C_BOOL), BIND(C) :: hydrostatic=.TRUE., debug_on=.FALSE.
   LOGICAL(C_BOOL), BIND(C, NAME='debug_hevi_solver') :: debug_hevi_solver_=.TRUE.
@@ -40 +43 @@
   TIME, PARAMETER :: print_trace_interval = 1.
   TIME, BIND(C) :: elapsed
-  NUM, BIND(C) :: g, ptop, cpp, cppv, Rd, Rv, preff, Treff, pbot, Phi_bot, rho_bot
-  NUM :: kappa
+  NUM, BIND(C) :: ptop, pbot, Phi_bot, rho_bot
+!  NUM, BIND(C) :: g, ptop, cpp, cppv, Rd, Rv, preff, Treff, pbot, Phi_bot, rho_bot
+!  NUM :: kappa
   NUM1(max_nb_stage), BIND(C)              :: tauj       ! diagonal of fast Butcher tableau
   NUM2(max_nb_stage,max_nb_stage), BIND(C) :: cslj, cflj ! slow and fast modified Butcher tableaus
@@ -252 +256 @@
   !
   SUBROUTINE init_params() BINDC(init_params)
+    USE earth_const
     kappa = Rd/cpp
     IF(is_mpi_master) THEN