Changeset 875 for trunk/LMDZ.GENERIC

Timestamp:

Feb 11, 2013, 11:06:56 AM (12 years ago)

Author:

jleconte

Message:

11/02/2013 == JL

• Updated moist convection scheme to handle situations with a large water vapor content
• Added a keyword to enable ocean runoff in callphys.def (activerunoff)

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/phystd/hydrol.F90 (modified) (3 diffs)
• libf/phystd/largescale.F90 (modified) (3 diffs)
• libf/phystd/moistadj.F90 (modified) (6 diffs)
• libf/phystd/su_watercycle.F90 (modified) (2 diffs)
• libf/phystd/watercommon_h.F90 (modified) (4 diffs)

trunk/LMDZ.GENERIC/README

@@ -896 +896 @@
 - Cosmetic changes in optcv (mostly spaces and line breaks) 
     ... so that comparisons with optci are easy e.g. through vimdiff
+    
+== 11/02/2013 == JL
+- Updated moist convection scheme to handle situations with a large water vapor content 
+- Added a keyword to enable ocean runoff in callphys.def (activerunoff)
+
+ 
+

trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90

@@ -3 +3 @@
      albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
 
+  use ioipsl_getincom 
   use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
   USE surfdat_h
@@ -55 +56 @@
       parameter (oceantime=10*24*3600)
 
-      logical oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
-      logical activerunoff ! enable simple runoff scheme?
-      logical oceanalbvary ! ocean albedo varies with the diurnal cycle?
+      logical,save :: oceanbulkavg ! relax ocean temperatures to a GLOBAL mean value?
+      logical,save :: activerunoff ! enable simple runoff scheme?
+      logical,save :: oceanalbvary ! ocean albedo varies with the diurnal cycle?
 
 !     Arguments
@@ -97 +98 @@
       data firstcall /.true./
 
-      oceanbulkavg=.false.
-      activerunoff=.false.
-      oceanalbvary=.false.
 
       if(firstcall)then
 
+         oceanbulkavg=.false.
+         oceanalbvary=.false.
+         write(*,*)"Activate runnoff into oceans?"
+         activerunoff=.false.
+         call getin("activerunoff",activerunoff)
+         write(*,*)" activerunoff = ",activerunoff
+         
+         
+         
          ALLOCATE(runoff(ngrid))
 

trunk/LMDZ.GENERIC/libf/phystd/largescale.F90

@@ -50 +50 @@
       EXTERNAL CBRT
       INTEGER i, k , nn
-      INTEGER,PARAMETER :: nitermax=1000
-      REAL,PARAMETER :: alpha=.5,qthreshold=1.e-6
+      INTEGER,PARAMETER :: nitermax=2000
+      REAL,PARAMETER :: alpha=.5,qthreshold=1.e-6 
+      ! JL13: if "careful, T<Tmin in psat water" appears often, you may want to stabilise the model by
+      !                   decreasing alpha and increasing nitermax accordingly
       REAL zt(ngrid), zq(ngrid)
       REAL zcond(ngrid),zcond_iter
       REAL zdelq(ngrid)
       REAL zqs(ngrid), zdqs(ngrid)
-      REAL psat_tmp
+      REAL psat_tmp,dlnpsat_tmp
       
 ! evaporation calculations
@@ -115 +117 @@
             Do nn=1,nitermax  
                call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
-               call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
+               call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)
                zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i))    
                   !zcond can be negative here
@@ -131 +133 @@
             zcond(i) = 0.0
             Do nn=1,nitermax  
-               call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
+               call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)
                zcond_iter = MAX(0.0,alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i)))           
                   !zcond always postive! cannot evaporate clouds!

trunk/LMDZ.GENERIC/libf/phystd/moistadj.F90

@@ -1 +1 @@
 subroutine moistadj(ngrid, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
 
-  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, Psat_water, Lcpdqsat_water
+  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, RCPV, Psat_water, Lcpdqsat_water
   USE tracer_h
 
@@ -59 +59 @@
       REAL cp_delta_t(nlayermx)
       REAL v_cptj(nlayermx), v_cptjk1, v_ssig
-      REAL v_cptt(ngrid,nlayermx), v_p, v_t, zpsat
-      REAL zqs(ngrid,nlayermx), zdqs(ngrid,nlayermx)
+      REAL v_cptt(ngrid,nlayermx), v_p, v_t, v_zqs,v_cptt2,v_pratio,v_dlnpsat
+      REAL zqs(ngrid,nlayermx), zdqs(ngrid,nlayermx),zpsat(ngrid,nlayermx),zdlnpsat(ngrid,nlayermx)
       REAL zq1(ngrid), zq2(ngrid)
       REAL gamcpdz(ngrid,2:nlayermx)
@@ -122 +122 @@
             v_p = pplay(i,k)
 
-            call Psat_water(v_t,v_p,zpsat,zqs(i,k))
-            call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
+            call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
+            call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
          ENDDO
       ENDDO
@@ -132 +132 @@
             zdp = pplev(i,k)-pplev(i,k+1)
             zdpm = pplev(i,k-1)-pplev(i,k)
-            gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
-                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
-                * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
-                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
+!            gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
+!                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
+!               * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
+!                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
+! general case where water is not a trace gas (JL13)
+            v_zqs     = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp)
+            v_cptt2   = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp)
+            v_pratio  = ((1./(1.+zpsat(i,k-1)/pplay(i,k-1)))*zdpm + (1./(1.+zpsat(i,k)/pplay(i,k)))*zdp)/(zdpm+zdp)
+            v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
+            gamcpdz(i,k) = v_pratio*( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
+                / ((1.- v_zqs) + v_zqs * RCPV/RCPD + v_zqs * v_pratio * v_dlnpsat)                    
          ENDDO
       ENDDO
@@ -200 +207 @@
          v_p = pplay(i,k)
          
-         call Psat_water(v_t,v_p,zpsat,zqs(i,k))
-         call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
+         call Psat_water(v_t,v_p,zpsat(i,k),zqs(i,k))
+         call Lcpdqsat_water(v_t,v_p,zpsat(i,k),zqs(i,k),zdqs(i,k),zdlnpsat(i,k))
 
 
@@ -226 +233 @@
          zdpm = pplev(i,k-1) - pplev(i,k)
          zdp = pplev(i,k) - pplev(i,k+1)
-         gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)             &
-                +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
-                * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
-                / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
+!         gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)          &
+!             +  RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) )                           &
+!             * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) )                                          &
+!             / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) )                    
+! general case where water is not a trace gas (JL13)
+         v_zqs     = (zqs(i,k-1)*zdpm + zqs(i,k)*zdp)/(zdpm+zdp)
+         v_cptt2   = (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp)/(zdpm+zdp)
+         v_pratio  = ((1./(1.+zpsat(i,k-1)/pplay(i,k-1)))*zdpm + (1./(1.+zpsat(i,k)/pplay(i,k)))*zdp)/(zdpm+zdp)
+         v_dlnpsat = (zdlnpsat(i,k-1)*zdpm + zdlnpsat(i,k)*zdp)/(zdpm+zdp)
+         gamcpdz(i,k) = v_pratio*( (R/RCPD*v_cptt2*(1.- v_zqs) + RLVTT*v_zqs) * (pplay(i,k-1)-pplay(i,k))/pplev(i,k) )  &
+                / ((1.- v_zqs) + v_zqs * RCPV/RCPD + v_zqs * v_pratio * v_dlnpsat)                    
       ENDDO
 

trunk/LMDZ.GENERIC/libf/phystd/su_watercycle.F90

@@ -1 +1 @@
       subroutine su_watercycle
 
-! to use  'getin'
-      use ioipsl_getincom 
       use watercommon_h
       implicit none
@@ -28 +26 @@
       RV = 1000.*8.314/mH2O ! caution! R is R/mugaz already!
 
-      RCPV   = 4. *RV    ! could be more precise...
-      !RCPV   = 5/2 *8.31*1000/mH2O ! BC modif here: doesn't appear to be correct
+      RCPV   = 1.88e3 ! specific heat capacity of water vapor at 350K
 
       RVTMP2 = RCPV/RCPD-1. ! not currently used...

trunk/LMDZ.GENERIC/libf/phystd/watercommon_h.F90

@@ -77 +77 @@
 
 !==================================================================
-      subroutine Lcpdqsat_water(T,p,psat,qsat,dqsat)
+      subroutine Lcpdqsat_water(T,p,psat,qsat,dqsat,dlnpsat)
 
          implicit none
@@ -84 +84 @@
 !     Purpose
 !     -------
-!     Compute L/cp*d (q_sat)/d T
+!     Compute dqsat=L/cp*d (q_sat)/d T and dlnpsat=L/cp d(ln Psat)/d T
 !     for a given temperature (K)! 
 !     Based on the Tetens formula from L.Li physical parametrization manual
@@ -98 +98 @@
   
 !        output
-         real dqsat
+         real dqsat,dlnpsat
 
 ! JL12 variables for tetens formula
@@ -125 +125 @@
 
          dqsat=RLVTT/RCPD*qsat*(p/(p-(1.-epsi)*psat))*dummy
+         dlnpsat=RLVTT/RCPD*dummy
          return
       end subroutine Lcpdqsat_water