Changeset 2330 for trunk/ICOSA_LMDZ

Timestamp:

May 25, 2020, 4:21:38 PM (4 years ago)

Author:

slebonnois

Message:

SL: modifications to be able to turn the physics with Cp(T)

File:

• trunk/ICOSA_LMDZ/src/phyvenus/interface_icosa_lmdz.f90 (modified) (11 diffs)

trunk/ICOSA_LMDZ/src/phyvenus/interface_icosa_lmdz.f90

@@ -326 +326 @@
     ! value in physics daysec=10087066.76s, important for solar radiation, 
     ! for the rest ptime/timeofday is what matters. To well estimate them,
-    ! day_length must be multiple of dt*itau_physics. Error order:1e-4.
+    ! day_length must be multiple of dt*itau_physics. Error order:1e-5.
     day_length=10087000   
     CALL getin('day_length',day_length)
@@ -368 +368 @@
           ! we reset date to midnight at lon=0 in the physics
           ptime=0.0
+          day_ini=0
         else
           write(*,*)"Error: physics time step in startphy.nc is different"
@@ -407 +408 @@
 
     ! Initialize some physical constants
-    CALL suphec
+    CALL suphec(cpp)
 
     ! Initialize cpdet_phy module
-    nu_venus=0.35
-    t0_venus=460.
+! FLAG CP CONSTANT OR NOT... (in LMDZ.COMMON: cpofT in gcm.def)
+! here, it allows to run the physics with a variable Cp together with dynamics
+! with constant Cp. 
+! Conversion Theta/T takes this into account (see below)
+! IF CP(T)
+!    nu_venus=0.35
+!    t0_venus=460.
+! IF CP CONSTANT
+    nu_venus=0.
+    t0_venus=0.
     CALL init_cpdet_phy(cpp,nu_venus,t0_venus)
   
@@ -733 +742 @@
     USE wind_mod
     USE omp_para
+    USE cpdet_phy_mod, only: tpot2t
     IMPLICIT NONE
     
     REAL(rstd) :: uc(3)
     INTEGER :: i,j,ij,l
+    REAL    :: tmpT(1),tmptheta(1),tmppk(1)
     
 
@@ -776 +787 @@
 
 !   compute theta, temperature and pression at layer
+! Conversion Theta/T takes Cp(T) into account
     DO    l    = ll_begin, ll_end
       DO j=jj_begin,jj_end
@@ -781 +793 @@
           ij=(j-1)*iim+i
           theta(ij,l) = theta_rhodz(ij,l,1) / ((p(ij,l)-p(ij,l+1))/g)
-          Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp
+          tmptheta(1) = theta(ij,l)
+          tmppk(1)    = pk(ij,l)
+          call tpot2t(1,tmptheta,tmpT,tmppk)
+          Temp(ij,l)  = tmpT(1)
           p_layer(ij,l)=preff*(pk(ij,l)/cpp)**(1./kappa) 
         ENDDO
@@ -789 +804 @@
 
 !!! Compute geopotential
+!!  This computation (with teta = cp T / pk !) is identical to 
+!!     delta phi = R/RMD T/p delta p         (r=R/RMD=cpp*kappa)
        
   ! for first layer
@@ -795 +812 @@
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i
-        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
+        phi( ij,1 ) = phis( ij ) + ( cpp*Temp(ij,1)/pk(ij,1) ) &
+                                 * ( pks(ij) - pk(ij,1) )
       ENDDO
     ENDDO
@@ -806 +824 @@
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i
-        phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
-                         * (  pk(ij,l-1) -  pk(ij,l)    )
+        phi(ij,l) =  0.5 *cpp*( Temp(ij,l)/pk(ij,l)+Temp(ij,l-1)/pk(ij,l-1) ) & 
+                         * ( pk(ij,l-1) - pk(ij,l) )
       ENDDO
     ENDDO
@@ -855 +873 @@
     USE theta2theta_rhodz_mod
     USE omp_para
+    USE cpdet_phy_mod, only: t2tpot
     IMPLICIT NONE
       INTEGER :: i,j,ij,l,iq
+      REAL    :: tmptheta(1),tmpT(1),tmppk(1)
           
       DO l=ll_begin,ll_end
@@ -948 +968 @@
 
 !   compute theta, temperature and pression at layer
+! Conversion Theta/T takes Cp(T) into account
     DO    l    = ll_begin, ll_end
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
           ij=(j-1)*iim+i
-          theta_rhodz(ij,l,1) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
+          tmpT(1)  = Temp(ij,l)
+          tmppk(1) = pk(ij,l)
+          call t2tpot(1,tmpT,tmptheta,tmppk)
+          theta_rhodz(ij,l,1) = tmptheta(1) * ((p(ij,l)-p(ij,l+1))/g)
         ENDDO
       ENDDO