Changeset 401

Timestamp:

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

Author:

dubos

Message:

Introduced entropy as prognostic variable - tested with JW06

Location:

codes/icosagcm/trunk/src

Files:

• caldyn_gcm.f90 (modified) (1 diff)
• caldyn_kernels_base.f90 (modified) (4 diffs)
• caldyn_kernels_hevi.f90 (modified) (1 diff)
• earth_const.f90 (modified) (3 diffs)
• etat0.f90 (modified) (1 diff)

codes/icosagcm/trunk/src/caldyn_gcm.f90

@@ -35 +35 @@
     IF (is_master) PRINT *, 'caldyn_conserv=',def
 
+    def='theta'
+    CALL getin('caldyn_thermo',def)
+    SELECT CASE(TRIM(def))
+    CASE('theta')
+       caldyn_thermo=thermo_theta
+    CASE('entropy')
+       caldyn_thermo=thermo_entropy
+    CASE DEFAULT
+       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_thermo : <', &
+            TRIM(def),'> options are <theta>, <entropy>'
+       STOP
+    END SELECT
+    
     CALL allocate_caldyn
 

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

codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90

@@ -401 +401 @@
              ! from now on pk contains the vertically-averaged geopotential
              pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-          ENDDO
-       ENDDO
-
+          END DO
+       END DO
     ELSE ! compressible
 
        DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-          ENDDO
-       ENDDO
+          SELECT CASE(caldyn_thermo)
+          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+             END DO
+          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s) 
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end         
+                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                     + pk(ij,l)*(cpp-theta(ij,l)) ! pk=temperature, theta=entropy
+             END DO
+          END SELECT
+       END DO
 
     END IF ! Boussinesq/compressible

codes/icosagcm/trunk/src/earth_const.f90

@@ -9 +9 @@
   REAL(rstd),SAVE :: kappa=0.2857143
   REAL(rstd),SAVE :: cpp=1004.70885
+  REAL(rstd),SAVE :: cppv=1004.70885 ! FIXME
+  REAL(rstd),SAVE :: Rv=461.5
+  REAL(rstd),SAVE :: Treff=273.
   REAL(rstd),SAVE :: preff=101325.
   REAL(rstd),SAVE :: pa=50000.
@@ -15 +18 @@
   REAL(rstd),SAVE :: gas_constant = 8.3144621 
   REAL(rstd),SAVE :: mu                 ! molar mass of the atmosphere
+
+  INTEGER, PARAMETER,PUBLIC :: thermo_theta=1, thermo_entropy=2, thermo_moist=3
+  INTEGER, PUBLIC :: caldyn_thermo
+  !$OMP THREADPRIVATE(caldyn_thermo) 
 
   LOGICAL, SAVE :: boussinesq, hydrostatic
@@ -33 +40 @@
     CALL getin("cpp",cpp)  
     CALL getin("preff",preff)  
+    CALL getin("Treff",Treff)  
     CALL getin("scale_height",scale_height)
     

codes/icosagcm/trunk/src/etat0.f90

@@ -172 +172 @@
          CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q)
       ENDIF
+
+      IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1)
+    
     ENDDO
     
-    IF( TRIM(etat0_type)/='williamson91.6' ) CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
-    
     CALL deallocate_field(f_temp)
     
   END SUBROUTINE etat0_collocated
+
+  SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset)
+    USE icosa
+    USE pression_mod
+    USE exner_mod
+    USE omp_para
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
+    REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: q(iim*jjm,llm,nqtot)
+    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
+    INTEGER,INTENT(IN) :: offset
+
+    REAL(rstd) :: p(iim*jjm,llm+1)
+    REAL(rstd) :: cppd,Rd, mass, p_ij, q_ij,r_ij, chi,nu, entropy, theta
+    INTEGER :: i,j,ij,l
+
+    cppd=cpp
+    Rd=kappa*cppd
+
+    CALL compute_pression(ps,p,offset)
+    ! flush p
+    !$OMP BARRIER
+    DO    l    = ll_begin, ll_end
+       DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+             ij=(j-1)*iim+i
+             mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass
+             p_ij = .5*(p(ij,l)+p(ij,l+1))  ! pressure at full level
+             SELECT CASE(caldyn_thermo)
+             CASE(thermo_theta)
+                theta = temp(ij,l)*(p_ij/preff)**(-kappa) 
+                theta_rhodz(ij,l) = mass * theta
+             CASE(thermo_entropy)
+                nu = log(p_ij/preff)
+                chi = log(temp(ij,l)/Treff)
+                entropy = cppd*chi-Rd*nu
+                theta_rhodz(ij,l) = mass * entropy
+!             CASE(thermo_moist)
+!                q_ij=q(ij,l,1)
+!                r_ij=1.-q_ij
+!                mass=mass*(1-q_ij) ! dry mass
+!                nu = log(p_ij/preff)
+!                chi = log(temp(ij,l)/Treff)
+!                entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu)
+!                theta_rhodz(ij,l) = mass * entropy                
+                CASE DEFAULT
+                   STOP
+             END SELECT
+          ENDDO
+       ENDDO
+    ENDDO
+    !$OMP BARRIER  
+  END SUBROUTINE compute_temperature2entropy
 
   SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q)