Changeset 401 for codes/icosagcm/trunk/src/caldyn_kernels_base.f90

Timestamp:

06/08/16 01:51:21 (8 years ago)

Author:

dubos

Message:

Introduced entropy as prognostic variable - tested with JW06

File:

• codes/icosagcm/trunk/src/caldyn_kernels_base.f90 (modified) (4 diffs)

codes/icosagcm/trunk/src/caldyn_kernels_base.f90

@@ -37 +37 @@
 
     INTEGER :: i,j,ij,l
-    REAL(rstd) :: p_ik, exner_ik
+    REAL(rstd) :: Rd, p_ik, exner_ik, temp_ik
     INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext
-
 
     CALL trace_start("compute_geopot")
@@ -47 +46 @@
     ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
 
+    Rd = kappa*cpp
+
     IF(caldyn_eta==eta_mass .AND. .NOT. DEC) THEN
 !!! Compute exner function and geopotential
+       STOP
        DO l = 1,llm
           !DIR$ SIMD
@@ -89 +91 @@
           END DO
           ! now pk contains the Lagrange multiplier (pressure)
-       ELSE ! non-Boussinesq, compute geopotential and Exner pressure
+       ELSE ! non-Boussinesq, compute pressure, Exner pressure or temperature, then geopotential
           ! uppermost layer
 
@@ -109 +111 @@
              END DO
           END IF
-          ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+
           DO l = 1,llm
-             !DIR$ SIMD
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext         
-                p_ik = pk(ij,l)
-                exner_ik = cpp * (p_ik/preff) ** kappa
-                pk(ij,l) = exner_ik
-                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
-             ENDDO
+             SELECT CASE(caldyn_thermo)
+             CASE(thermo_theta)
+                !DIR$ SIMD
+                DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                   p_ik = pk(ij,l)
+                   exner_ik = cpp * (p_ik/preff) ** kappa
+                   pk(ij,l) = exner_ik
+                   ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+                   geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
+                ENDDO
+             CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff)
+                !DIR$ SIMD
+                DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                   p_ik = pk(ij,l)
+                   temp_ik = Treff*exp((theta(ij,l) + Rd*log(p_ik/preff))/cpp)
+                   pk(ij,l) = temp_ik
+                   ! specific volume v = Rd*T/p = dphi/g/rhodz
+                   geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik
+                ENDDO
+             CASE DEFAULT
+                STOP
+             END SELECT
           ENDDO
        END IF