Changeset 786 for trunk/LMDZ.GENERIC/libf

Timestamp:

Sep 19, 2012, 12:27:11 PM (12 years ago)

Author:

jleconte

Message:

19/09/2012 == JL

• Correction in largescale to improve robustness when large water vapor amount
• Correction in soil_setting to allow change of the number of subsurface layers

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• ave_stelspec.F90 (modified) (1 diff)
• dimphys.h (modified) (1 diff)
• largescale.F90 (modified) (2 diffs)
• mass_redistribution.F90 (modified) (1 diff)
• soil_settings.F (modified) (2 diffs)
• watercommon_h.F90 (modified) (4 diffs)

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

@@ -152 +152 @@
          band=1
          Do while(lam(1).lt. real(10000.0/BWNV(band+1)))
-            if (band+1.ge.L_NSPECTV+1) exit
+            if (band.gt.L_NSPECTV-1) exit
             band=band+1
          enddo

trunk/LMDZ.GENERIC/libf/phystd/dimphys.h

@@ -9 +9 @@
 ! nsoilmx : number of subterranean layers
       !integer, parameter :: nsoilmx = 4 ! for a test 
-      integer, parameter :: nsoilmx = 18 ! for z1=0.02 cm, depth = 104.8 m 
+      integer, parameter :: nsoilmx = 18 ! for z1=0.0002 cm, depth = 18 m => mars case 
+      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 cm, depth = 104.8 m => earth case
 !-----------------------------------------------------------------------

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

@@ -51 +51 @@
       EXTERNAL CBRT
       INTEGER i, k , nn
-      INTEGER,PARAMETER :: niter=4
+      INTEGER,PARAMETER :: nitermax=1000
+      REAL,PARAMETER :: alpha=.5,qthreshold=1.e-6
       REAL zt(ngridmx), zq(ngridmx)
       REAL zcond(ngridmx),zcond_iter
@@ -92 +93 @@
 
          if(zt(i).le.15.) then
-            print*,'in lsc',i,zt(i)
-            zt(i)=15.   ! check too low temperatures
+            print*,'in lsc',i,k,zt(i)
+!           zt(i)=15.   ! check too low temperatures
          endif
-!         call watersat(zt(i),pplay(i,k),zqs(i))
          call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
  
-         zdelq(i) = ratqs * zq(i)
-         if(zq(i)+zdelq(i).lt.0.999) then
-            rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i))
+         zdelq(i) = MAX(MIN(ratqs * zq(i),1.-zq(i)),1.e-12)
+         rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i))
+!        print*,zq(i),zdelq(i),zqs(i),rneb(i,k)
+         if (rneb(i,k).lt.0.) then  !no clouds
+
+            rneb(i,k)=0.
+            zcond(i)=0.
+
+         else if (rneb(i,k).gt.1.) then    !complete cloud cover, we start without evaporating
+
+            rneb(i,k)=1.
+            zt(i)=pt(i,k)+pdt(i,k)*ptimestep
+            zx_q(i) = pq(i,k,igcm_h2o_vap)+pdq(i,k,igcm_h2o_vap)*ptimestep
+            dqevap(i,k)=0.
+!           iterative process to stabilize the scheme when large water amounts JL12
+            zcond(i) = 0.0
+            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))
+               zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i))    
+                  !zcond can be negative here
+               zx_q(i) = zx_q(i) - zcond_iter
+               zcond(i) = zcond(i) + zcond_iter
+               zt(i) = zt(i) + zcond_iter*RLVTT/RCPD
+               if (ABS(zcond_iter/alpha).lt.qthreshold) exit
+            End do ! niter
+            zcond(i)=MAX(zcond(i),-(pq(i,k,igcm_h2o_ice)+pdq(i,k,igcm_h2o_ice)*ptimestep))
+
+         else   !standard case      
+
             zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0 !water vapor in cloudy sky
-            if (rneb(i,k) .LE. 0.0) zx_q(i) = 0.0
-            if (rneb(i,k) .GE. 1.0) zx_q(i) = zq(i)
-            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
-         else
-            if(zq(i).gt.zqs(i)) then
-               rneb(i,k)=1.
-               zx_q(i)=zq(i)
-            else
-               rneb(i,k)=0.
-               zx_q(i)=zqs(i)  !no condensation needed
-            Endif
+!           iterative process to stabilize the scheme when large water amounts JL12
+            zcond(i) = 0.0
+            Do nn=1,nitermax  
+               call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
+               zcond_iter = MAX(0.0,alpha*(zx_q(i)-zqs(i))/(1.+zdqs(i)))           
+                  !zcond always postive! cannot evaporate clouds!
+                  !this is why we must reevaporate before largescale
+               zx_q(i) = zx_q(i) - zcond_iter
+               zcond(i) = zcond(i) + zcond_iter
+               if (ABS(zcond_iter/alpha).lt.qthreshold) exit
+               zt(i) = zt(i) + zcond_iter*RLVTT/RCPD
+               call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
+            End do ! niter
+
          Endif
 
-!        iterative process to stabilize the scheme when large water amounts JL12
-         zcond(i) = 0.0
-         Do nn=1,niter  
-!            call watersat_grad(zt(i),zqs(i),zdqs(i))
-            call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
-            zcond_iter = MAX(0.0,(zx_q(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i)))          
-               !zcond always postive! cannot evaporate clouds!
-               !this is why we must reevaporate before largescale
-            zx_q(i) = zx_q(i) - zcond_iter
-            zcond(i) = zcond(i) + zcond_iter
-            zt(i) = zt(i) + zcond_iter*RLVTT/RCPD
-            call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
-         End do ! niter
+         zcond(i) = zcond(i)*rneb(i,k)/ptimestep ! JL12
 
-
-         zcond(i) = zcond(i)/ptimestep ! added by RDW
       ENDDO
 

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

@@ -146 +146 @@
                   zmassboil(ig)=pqs(ig,igcm_h2o_vap)/ptimestep
                endif
-               zmassboil(ig)=zmassboil(ig)*0.3 !momentary, should be 1. JL12
+               zmassboil(ig)=zmassboil(ig)*0.0 !momentary, should be 1. JL12
                pdqsmr(ig,igcm_h2o_vap)=-zmassboil(ig)
                pdtsrfmr(ig)=-zmassboil(ig)*RLVTT/pcapcal(ig)

trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F

@@ -133 +133 @@
         ! read <depth> coordinate
         if (interpol) then !put values in oldmlayer
+          if (.not.allocated(oldmlayer)) then
+             allocate(oldmlayer(dimlen),stat=ierr)
+          endif
 #ifdef NC_DOUBLE
            ierr = NF_GET_VAR_DOUBLE(nid,nvarid,oldmlayer)
@@ -159 +162 @@
       ! default mlayer distribution, following a power law:
       !  mlayer(k)=lay1*alpha**(k-1/2)
-        lay1=2.e-4
+        lay1=2.e-4  !mars case (nsoilmx=18)
+!        lay1=3.e-2  !earth case (nsoilmx=13)
         alpha=2
         do iloop=0,nsoil-1

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

@@ -51 +51 @@
          real,parameter :: Trefvaporization=35.86
          real,parameter :: Trefsublimation=7.66
+         real,parameter :: Tmin=8.
          real,parameter :: r3vaporization=17.269
          real,parameter :: r3sublimation=21.875
@@ -56 +57 @@
 ! checked vs. old watersat data 14/05/2012 by JL.
 
-         if (T.lt.T_h2O_ice_liq) then ! solid / vapour
-            psat = Pref_solid_liquid*Exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation))
-         else                 ! liquid / vapour
-            psat = Pref_solid_liquid*Exp(r3vaporization*(T-T_h2O_ice_liq)/(T-Trefvaporization))
+         if (T.gt.T_h2O_ice_liq) then 
+            psat = Pref_solid_liquid*Exp(r3vaporization*(T-T_h2O_ice_liq)/(T-Trefvaporization)) ! liquid / vapour
+         else if (T.lt.Tmin) then
+            print*, "careful, T<Tmin in psat water"
+            psat = Pref_solid_liquid*Exp(r3sublimation*(Tmin-T_h2O_ice_liq)/(Tmin-Trefsublimation)) ! min psat  
+         else                 
+            psat = Pref_solid_liquid*Exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation)) ! solid / vapour
          endif
          if(psat.gt.p) then
@@ -99 +103 @@
          real,parameter :: Pref_solid_liquid=611.14
          real,parameter :: Trefvaporization=35.86
+         real,parameter :: Tmin=8.
          real,parameter :: Trefsublimation=7.66
          real,parameter :: r3vaporization=17.269
@@ -110 +115 @@
          endif
 
-         if (T.lt.T_h2O_ice_liq) then ! solid / vapour
-            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2
-         else                 ! liquid / vapour
-            dummy = r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2
+         if (T.gt.T_h2O_ice_liq) then
+            dummy = r3vaporization*(T_h2O_ice_liq-Trefvaporization)/(T-Trefvaporization)**2  ! liquid / vapour
+         else if (T.lt.Tmin) then
+            print*, "careful, T<Tmin in Lcp psat water"
+            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(Tmin-Trefsublimation)**2  ! solid / vapour
+         else               
+            dummy = r3sublimation*(T_h2O_ice_liq-Trefsublimation)/(T-Trefsublimation)**2  ! solid / vapour
          endif
 
-         dqsat=RLVTT/RCPD*qsat**2*(p/(epsi*psat))*dummy
+         dqsat=RLVTT/RCPD*qsat*(p/(p-(1.-epsi)*psat))*dummy
          return
       end subroutine Lcpdqsat_water