Changeset 1993 for trunk

Timestamp:

Aug 29, 2018, 4:41:34 PM (6 years ago)

Author:

jleconte

Message:

29/08/2018 == JL

-watersat was used only in vdifc and thus it was not consistent with other routines (turbdiff, rain, largescale...) which used Psat_water from watercommon.
This is now harmonized. ALl routines use Psat_water. Watersat.F has been removed, but the routine is now in watercommon for archival purpose. It is not used anymore.
-also changed the number of chars for tname in the dyn3D/infotrac.F90 to be able to run rcm1d.

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/dyn3d/infotrac.F90 (modified) (1 diff)
• libf/phystd/callcorrk.F90 (modified) (3 diffs)
• libf/phystd/largescale.F90 (modified) (7 diffs)
• libf/phystd/rain.F90 (modified) (2 diffs)
• libf/phystd/turbdiff.F90 (modified) (3 diffs)
• libf/phystd/vdifc.F (modified) (4 diffs)
• libf/phystd/watercommon_h.F90 (modified) (1 diff)

trunk/LMDZ.GENERIC/README

@@ -1389 +1389 @@
 correct bug on rain initialization at each timestep in physiqu_mod so that mass_redist can work without rain (if precipitations are taken care with sedimentation for example)
 change a table that was used as a float to a float in gfluxv for speedup. Does not change results bit for bit
+
+== 29/08/2018 == JL
+-watersat was used only in vdifc and thus it was not consistent with other routines (turbdiff, rain, largescale...) which used Psat_water from watercommon.
+This is now harmonized. ALl routines use Psat_water. Watersat.F has been removed, but the routine is now in watercommon for archival purpose. It is not used anymore.
+-also changed the number of chars for tname in the dyn3D/infotrac.F90 to be able to run rcm1d.

trunk/LMDZ.GENERIC/libf/dyn3d/infotrac.F90

@@ -12 +12 @@
 
 ! Name variables
-  CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics
-  CHARACTER(len=23), ALLOCATABLE, DIMENSION(:), SAVE :: ttext ! tracer long name for diagnostics
+  CHARACTER(len=30), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics
+  CHARACTER(len=33), ALLOCATABLE, DIMENSION(:), SAVE :: ttext ! tracer long name for diagnostics
 
 ! iadv  : index of trasport schema for each tracer

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

@@ -167 +167 @@
       real*8 pq_temp(nlayer)
 ! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
-      real ptemp, Ttemp, qsat
+      real psat,qsat
 
       logical OLRz
@@ -533 +533 @@
             if(RH.lt.0.0) RH=0.0
             
-            ptemp=pplay(ig,l)
-            Ttemp=pt(ig,l)
-            call watersat(Ttemp,ptemp,qsat)
+            call Psat_water(pt(ig,l),pplay(ig,l),psat,qsat)
 
             !pq_temp(l) = qsat      ! fully saturated everywhere
@@ -551 +549 @@
          RH = satval * (1 - 0.02) / 0.98
          if(RH.lt.0.0) RH=0.0
-
-!         ptemp = pplev(ig,1)
-!         Ttemp = tsurf(ig)
-!         call watersat(Ttemp,ptemp,qsat)
 
          qvar(2*nlayer+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)

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

@@ -5 +5 @@
       use ioipsl_getin_p_mod, only: getin_p
       use watercommon_h, only : RLVTT, RCPD, RVTMP2,  &
-          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_waterDP,Lcpdqsat_waterDP
+          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_water,Lcpdqsat_water
       USE tracer_h
       IMPLICIT none
@@ -50 +50 @@
       ! JL13: if "careful, T<Tmin in psat water" appears often, you may want to stabilise the model by
       !                   decreasing alpha and increasing nitermax accordingly
-      DOUBLE PRECISION zt(ngrid), zq(ngrid)
+      DOUBLE PRECISION zq(ngrid)
       DOUBLE PRECISION zcond(ngrid),zcond_iter
       DOUBLE PRECISION zdelq(ngrid)
-      DOUBLE PRECISION zqs(ngrid), zdqs(ngrid)
-      DOUBLE PRECISION local_p,psat_tmp,dlnpsat_tmp,Lcp
+      DOUBLE PRECISION zqs(ngrid)
+      real zt(ngrid),local_p,psat_tmp,dlnpsat_tmp,Lcp,zqs_temp,zdqs
       
 ! evaporation calculations
@@ -106 +106 @@
 !           zt(i)=15.   ! check too low temperatures
          endif
-         call Psat_waterDP(zt(i),local_p,psat_tmp,zqs(i))
+         call Psat_water(zt(i),local_p,psat_tmp,zqs_temp)
+         zqs(i)=zqs_temp
  
          zdelq(i) = MAX(MIN(ratqs * zq(i),1.-zq(i)),1.d-12)
@@ -127 +128 @@
             zcond(i) = 0.0d0
             Do nn=1,nitermax  
-               call Psat_waterDP(zt(i),local_p,psat_tmp,zqs(i))
-               call Lcpdqsat_waterDP(zt(i),local_p,psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)
-               zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs(i))          
+               call Psat_water(zt(i),local_p,psat_tmp,zqs_temp)
+               zqs(i)=zqs_temp
+               call Lcpdqsat_water(zt(i),local_p,psat_tmp,zqs_temp,zdqs,dlnpsat_tmp)
+               zcond_iter = alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs)     
                   !zcond can be negative here
                zx_q(i) = zx_q(i) - zcond_iter
@@ -150 +152 @@
             zcond(i) = 0.0d0
             Do nn=1,nitermax  
-               call Lcpdqsat_waterDP(zt(i),local_p,psat_tmp,zqs(i),zdqs(i),dlnpsat_tmp)
-               zcond_iter = MAX(0.0d0,alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs(i)))       
+               call Lcpdqsat_water(zt(i),local_p,psat_tmp,zqs(i),zdqs,dlnpsat_tmp)
+               zcond_iter = MAX(0.0d0,alpha*(zx_q(i)-zqs(i))/(1.d0+zdqs))          
                   !zcond always postive! cannot evaporate clouds!
                   !this is why we must reevaporate before largescale
@@ -159 +161 @@
 !              if (ABS(zcond_iter/alpha).lt.qthreshold) exit
                zt(i) = zt(i) + zcond_iter*Lcp*rneb(i,k)
-               call Psat_waterDP(zt(i),local_p,psat_tmp,zqs(i))
+               call Psat_water(zt(i),local_p,psat_tmp,zqs_temp)
+               zqs(i)=zqs_temp
                if (nn.eq.nitermax) print*,'itermax in largescale'
             End do ! niter
@@ -172 +175 @@
          pdqvaplsc(1:ngrid,k)  = dqevap(1:ngrid,k) - zcond(1:ngrid)
          pdqliqlsc(1:ngrid,k) = - pdqvaplsc(1:ngrid,k)
-         pdtlsc(1:ngrid,k)  = pdqliqlsc(1:ngrid,k)*real(Lcp)
+         pdtlsc(1:ngrid,k)  = pdqliqlsc(1:ngrid,k)*Lcp
 
    Enddo ! k= nlayer, 1, -1

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

@@ -82 +82 @@
       real reffh2oliq(ngrid,nlayer),reffh2oice(ngrid,nlayer)
   
-      real ttemp, ptemp, psat_tmp
+      real ptemp, psat_tmp
       real tnext(ngrid,nlayer)
 
@@ -183 +183 @@
       DO k = 1, nlayer
          DO i = 1, ngrid
-            ttemp = zt(i,k) 
             ptemp = pplay(i,k)
-!            call watersat(ttemp,ptemp,zqs(i,k))
-            call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
+            call Psat_water(zt(i,k) ,ptemp,psat_tmp,zqs(i,k))
             call Tsat_water(ptemp,Tsat(i,k))
          ENDDO

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

@@ -7 +7 @@
           pdqdif,pdqevap,pdqsdif,flux_u,flux_v,lastcall)
 
-      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
+      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water, Lcpdqsat_water
       use radcommon_h, only : sigma, glat
       use surfdat_h, only: dryness
@@ -116 +116 @@
       REAL zdmassevap(ngrid)
       REAL rho(ngrid)         ! near-surface air density
-      REAL qsat(ngrid),psat(ngrid)
       REAL kmixmin
 
 !     Variables added for implicit latent heat inclusion
 !     --------------------------------------------------
-      real dqsat(ngrid), qsat_temp1, qsat_temp2,psat_temp
+      real dqsat(ngrid),psat_temp,qsat(ngrid),psat(ngrid)
 
       integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
@@ -518 +517 @@
                 ! Calculate the value of qsat at the surface (water)
                 call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
-                call Psat_water(ptsrf(ig)-0.0001,pplev(ig,1),psat_temp,qsat_temp1)
-                call Psat_water(ptsrf(ig)+0.0001,pplev(ig,1),psat_temp,qsat_temp2)
-                dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
-                ! calculate dQsat / dT by finite differences
-                ! we cannot use the updated temperature value yet...
+                call Lcpdqsat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig),dqsat(ig),psat_temp)
+                dqsat(ig)=dqsat(ig)*RCPD/RLVTT
                enddo
 

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

@@ -7 +7 @@
      &     pdqdif,pdqsdif,lastcall)
 
-      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
+      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol 
+     &    ,Psat_water, Lcpdqsat_water
       use radcommon_h, only : sigma
       USE surfdat_h
@@ -112 +113 @@
       REAL zq1temp(ngrid)
       REAL rho(ngrid)         ! near-surface air density
-      REAL qsat(ngrid)
       DATA firstcall/.true./
       REAL kmixmin
@@ -118 +118 @@
 !     Variables added for implicit latent heat inclusion
 !     --------------------------------------------------
-      real latconst, dqsat(ngrid), qsat_temp1, qsat_temp2
-      real z1_Tdry(ngrid), z2_Tdry(ngrid)
-      real z1_T(ngrid), z2_T(ngrid)
-      real zb_T(ngrid)
-      real zc_T(ngrid,nlay)
-      real zd_T(ngrid,nlay)
-      real lat1(ngrid), lat2(ngrid)
-      real surfh2otot
-      logical surffluxdiag
-      integer isub ! sub-loop for precision
+      real dqsat(ngrid),psat_temp,qsat(ngrid),psat(ngrid)
 
       integer ivap, iice ! also make liq for clarity on surface...
@@ -520 +511 @@
 
                do ig=1,ngrid
-
                 ! Calculate the value of qsat at the surface (water)
-                call watersat(ptsrf(ig),pplev(ig,1),qsat(ig))
-                call watersat(ptsrf(ig)-0.0001,pplev(ig,1),qsat_temp1)
-                call watersat(ptsrf(ig)+0.0001,pplev(ig,1),qsat_temp2)
-                dqsat(ig)=(qsat_temp2-qsat_temp1)/0.0002
-                ! calculate dQsat / dT by finite differences
-                ! we cannot use the updated temperature value yet...
- 
+                call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
+                call Lcpdqsat_water(ptsrf(ig),pplev(ig,1),psat(ig),
+     &             qsat(ig),dqsat(ig),psat_temp)
+                dqsat(ig)=dqsat(ig)*RCPD/RLVTT
                enddo
 

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

@@ -208 +208 @@
          return
       end subroutine Tsat_water
-
-!==================================================================
-      subroutine Psat_waterDP(T,p,psat,qsat)
-
-         implicit none
-
-!==================================================================
-!     Purpose
-!     -------
-!     Compute the saturation vapor pressure and mass mixing ratio at saturation (kg/kg)
+!==================================================================
+
+
+!==================================================================
+subroutine watersat(T,p,qsat)
+  implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the water mass mixing ratio at saturation (kg/kg)
 !     for a given pressure (Pa) and temperature (K)
-!     DOUBLE PRECISION
-!     Based on the Tetens formula from L.Li physical parametrization manual
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-
-!        input
-         double precision, intent(in) :: T, p
-  
-!        output
-         double precision psat,qsat
-
-! JL12 variables for tetens formula
-         double precision,parameter :: Pref_solid_liquid=611.14d0
-         double precision,parameter :: Trefvaporization=35.86D0
-         double precision,parameter :: Trefsublimation=7.66d0
-         double precision,parameter :: Tmin=8.d0
-         double precision,parameter :: r3vaporization=17.269d0
-         double precision,parameter :: r3sublimation=21.875d0
-
-! checked vs. old watersat data 14/05/2012 by JL.
-
-         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  
-         ! Ehouarn: gfortran says: Error: Result of EXP underflows its kind,
-         !          so set psat to the smallest possible value instead
-            psat=tiny(psat)
-         else                 
-            psat = Pref_solid_liquid*Exp(r3sublimation*(T-T_h2O_ice_liq)/(T-Trefsublimation)) ! solid / vapour
-         endif
-         if(psat.gt.p) then
-            qsat=1.d0
-         else
-            qsat=epsi*psat/(p-(1.d0-epsi)*psat)
-         endif
-         return
-      end subroutine Psat_waterDP
-
-
-
-
-!==================================================================
-      subroutine Lcpdqsat_waterDP(T,p,psat,qsat,dqsat,dlnpsat)
-
-         implicit none
-
-!==================================================================
-!     Purpose
-!     -------
-!     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
-!
-!     Authors
-!     -------
-!     Jeremy Leconte (2012)
-!
-!==================================================================
-
-!        input
-         double precision T, p, psat, qsat
-  
-!        output
-         double precision dqsat,dlnpsat
-
-! JL12 variables for tetens formula
-         double precision,parameter :: Pref_solid_liquid=611.14d0
-         double precision,parameter :: Trefvaporization=35.86d0
-         double precision,parameter :: Tmin=8.d0
-         double precision,parameter :: Trefsublimation=7.66d0
-         double precision,parameter :: r3vaporization=17.269d0
-         double precision,parameter :: r3sublimation=21.875d0
-
-         double precision :: dummy
-
-         if (psat.gt.p) then
-            dqsat=0.d0
-            return
-         endif
-
-         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*(p/(p-(1.d0-epsi)*psat))*dummy
-         dlnpsat=RLVTT/RCPD*dummy
-         return
-      end subroutine Lcpdqsat_waterDP
+!     A replacement for the old watersat.F in the Martian GCM.
+!     Based on FCTTRE.h in the LMDTERRE model.
+!
+!     JL18 watersat was used only in vdifc and thus it was not consistent with other routines (turbdiff, rain, largescale...)
+!        which used Psat_water. This is now harmonized
+!        we put it here for archival purpose, but it is not used anymore.
+!
+!     Authors
+!     -------
+!     Robin Wordsworth (2010)
+!
+!==================================================================
+
+!   input
+  real T, p
+  
+!   output
+  real qsat
+
+! checked vs. NIST data 22/06/2010 by RW.
+! / by p gives partial pressure
+! x by epsi converts to mass mixing ratio
+
+  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
+     qsat = 100.0 * 10**(2.07023 - 0.00320991             &
+          * T - 2484.896 / T + 3.56654 * alog10(T))
+  else                 ! liquid / vapour
+     qsat = 100.0 * 10**(23.8319 - 2948.964 / T - 5.028  &
+          * alog10(T) - 29810.16 * exp( -0.0699382 * T)  &
+          + 25.21935 * exp(-2999.924/T))
+  endif
+!  qsat=epsi*qsat/p
+  if(qsat.gt.p) then
+     qsat=1.
+  else
+     qsat=epsi*qsat/(p-(1.-epsi)*qsat)
+  endif
+  return
+end subroutine watersat
+
+subroutine watersat_grad(T,qsat,dqsat)
+
+  implicit none
+
+!==================================================================
+!     Purpose
+!     -------
+!     Compute the L/cp*d (q_sat)/d T
+!     for a given temperature (K)
+!
+!     JL18: see watersat
+!
+!     Authors
+!     -------
+!     Robin Wordsworth (2010)
+!
+!==================================================================
+
+!   input
+  real T,qsat
+  
+!   output
+  real dqsat
+
+!  if (T.lt.T_coup) then ! solid / vapour !why use T_coup?????????? JL12
+  if (T.lt.T_h2O_ice_liq) then ! solid / vapour
+     dqsat = RLVTT/RCPD*qsat*(3.56654/T             &
+          +2484.896*LOG(10.)/T**2                   &
+          -0.00320991*LOG(10.))
+  else                 ! liquid / vapour
+     dqsat = RLVTT/RCPD*qsat*LOG(10.)*              &
+          (2948.964/T**2-5.028/LOG(10.)/T           &
+          +25.21935*2999.924/T**2*EXP(-2999.924/T)  &
+          +29810.16*0.0699382*EXP(-0.0699382*T))
+  end if
+
+  return
+end subroutine watersat_grad
+