Changeset 2871 for trunk/LMDZ.GENERIC/libf

Timestamp:

Jan 16, 2023, 4:47:11 PM (2 years ago)

Author:

jleconte

Message:

Improved rain evaporation for les + correction in callkeys

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• callkeys_mod.F90 (modified) (2 diffs)
• largescale.F90 (modified) (3 diffs)
• physiq_mod.F90 (modified) (8 diffs)
• rain.F90 (modified) (11 diffs)

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

@@ -42 +42 @@
       logical,save :: generic_rain
 !$OMP THREADPRIVATE(varactive,varfixed,radfixed,sedimentation,generic_condensation,generic_rain)
-      logical,save :: water,watercond,waterrain
-!$OMP THREADPRIVATE(water,watercond,waterrain)
-      logical,save :: aeroco2,aeroh2o,aeroh2so4,aeroback2lay
-!$OMP THREADPRIVATE(aeroco2,aeroh2o,aeroh2so4,aeroback2lay)
+      logical,save :: water ,watercond, waterrain, moistadjustment
+!$OMP THREADPRIVATE(water, watercond, waterrain, moistadjustment)
+      logical,save :: aeroco2, aeroh2o, aeroh2so4, aeroback2lay
+!$OMP THREADPRIVATE(aeroco2, aeroh2o, aeroh2so4, aeroback2lay)
       logical,save :: aeronh3, aeronlay, aeroaurora
 !$OMP THREADPRIVATE(aeronh3,aeronlay,aeroaurora)
@@ -52 +52 @@
 !$OMP THREADPRIVATE(aerovenus)
       ! detailed sub-options when with "Venus-like" aerosol additions
-      logical,save :: aerovenus1,aerovenus2,aerovenus2p,aerovenus3,aerovenusUV
-!$OMP THREADPRIVATE(aerovenus1,aerovenus2,aerovenus2p,aerovenus3,aerovenusUV)
+      logical,save :: aerovenus1, aerovenus2, aerovenus2p, aerovenus3, aerovenusUV
+!$OMP THREADPRIVATE(aerovenus1, aerovenus2, aerovenus2p, aerovenus3, aerovenusUV)
 
-      logical,save :: aerofixco2,aerofixh2o
-!$OMP THREADPRIVATE(aerofixco2,aerofixh2o)
+      logical,save :: aerofixco2, aerofixh2o
+!$OMP THREADPRIVATE(aerofixco2, aerofixh2o)
       integer,save :: aerogeneric ! number of aerosols of "generic" kind
 !$OMP THREADPRIVATE(aerogeneric)

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

@@ -40 +40 @@
 !     Options du programme
       REAL, SAVE :: ratqs   ! determine largeur de la distribution de vapeur
-!$OMP THREADPRIVATE(ratqs)
+      REAL, SAVE :: qvap_deep   ! deep mixing ratio of water vapor when simulating bottom less planets
+!$OMP THREADPRIVATE(ratqs, qvap_deep)
 
 !     Variables locales
@@ -72 +73 @@
          call getin_p("ratqs",ratqs)
          write(*,*) " ratqs = ",ratqs
+
+         write(*,*) "Deep water vapor mixing ratio ? (no effect if negative) "
+         qvap_deep=-1. ! default value
+         call getin_p("qvap_deep",qvap_deep)
+         write(*,*) " qvap_deep = ",qvap_deep
 
          firstcall = .false.
@@ -181 +187 @@
 
    Enddo ! k= nlayer, 1, -1
- 
+
+   if (qvap_deep >= 0.) then
+      !brings lower vapor ratio to a fixed value.
+      ! tau=3600. seems too fast
+      pdqvaplsc(1:ngrid,1) = (qvap_deep - pq(1:ngrid,1,igcm_h2o_vap))/14400. - pdq(1:ngrid,1,igcm_h2o_vap)
+   endif
+
 
       end

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

@@ -61 +61 @@
                               startphy_file, testradtimes, tlocked, &
                               tracer, UseTurbDiff, water, watercond, &
-                              waterrain, generic_rain, global1d, szangle
+                              waterrain, generic_rain, global1d, moistadjustment, szangle
       use generic_tracer_index_mod, only: generic_tracer_index
       use nonoro_gwd_ran_mod, only: nonoro_gwd_ran
@@ -85 +85 @@
          comm_FLUXTOP_LW,comm_FLUXSURF_SW,          &
          comm_FLUXSURF_LW,comm_FLXGRD,              &
-         comm_DTRAIN,comm_DTLSC,comm_H2OICE_REFF
+         comm_DTRAIN,comm_DTLSC,comm_H2OICE_REFF,   &
+         comm_LATENT_HF
+
 #endif
 
@@ -1425 +1427 @@
             if(watercond.and.(RLVTT.gt.1.e-8))then
                
-               if (.not.calltherm) then
+               if ((.not.calltherm).and.moistadjustment) then
                   dqmoist(1:ngrid,1:nlayer,1:nq)=0.
                   dtmoist(1:ngrid,1:nlayer)=0.
@@ -1506 +1508 @@
                zdqssnow(1:ngrid)    = 0.0
 
-               call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
+               call rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pphi,pt,pdt,pq,pdq,            &
                          zdtrain,zdqrain,zdqsrain,zdqssnow,reevap_precip,cloudfrac)
 
@@ -2448 +2450 @@
 
             call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
-            call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
+            if ((.not.calltherm).and.moistadjustment) then
+               call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
+            endif
             call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
              
@@ -2571 +2575 @@
    comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer)
    comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer)
+   comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
+   if (.not.calldifv) comm_LATENT_HF(:)=0.0
    if ((tracer).and.(water)) then
      comm_CLOUDFRAC(1:ngrid,1:nlayer)=cloudfrac(1:ngrid,1:nlayer)
@@ -2576 +2582 @@
      comm_SURFRAIN(1:ngrid)=zdqsrain(1:ngrid)
      comm_DQVAP(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_h2o_vap)
-     comm_ALBEQ(1:ngrid)=albedo_equivalent(1:ngrid)
+     comm_DQICE(1:ngrid,1:nlayer)=pdq(1:ngrid,1:nlayer,igcm_h2o_ice)
      comm_H2OICE_REFF(1:ngrid,1:nlayer)=reffrad(1:ngrid,1:nlayer,iaero_h2o)
      comm_REEVAP(1:ngrid)=reevap_precip(1:ngrid)
@@ -2582 +2588 @@
      comm_DTLSC(1:ngrid,1:nlayer)=dtlscale(1:ngrid,1:nlayer)
      comm_RH(1:ngrid,1:nlayer)=RH(1:ngrid,1:nlayer)
+   else
+      comm_CLOUDFRAC(1:ngrid,1:nlayer)=0.
+      comm_TOTCLOUDFRAC(1:ngrid)=0.
+      comm_SURFRAIN(1:ngrid)=0.
+      comm_DQVAP(1:ngrid,1:nlayer)=0.
+      comm_DQICE(1:ngrid,1:nlayer)=0.
+      comm_H2OICE_REFF(1:ngrid,1:nlayer)=0.
+      comm_REEVAP(1:ngrid)=0.
+      comm_DTRAIN(1:ngrid,1:nlayer)=0.
+      comm_DTLSC(1:ngrid,1:nlayer)=0.
+      comm_RH(1:ngrid,1:nlayer)=0.
    endif
-   !comm_DQICE(1:ngrid,1:nlayer)=zdqdyn(1:ngrid,1:nlayer)
    comm_FLUXTOP_DN(1:ngrid)=fluxtop_dn(1:ngrid)
    comm_FLUXABS_SW(1:ngrid)=fluxabs_sw(1:ngrid)

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

@@ -1 @@
-subroutine rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,reevap_precip,rneb)
+subroutine rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,pphi,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,reevap_precip,rneb)
 
 
@@ -30 +30 @@
       real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
       real,intent(in) :: pplay(ngrid,nlayer)   ! mid-layer pressure (Pa)
+      real,intent(in) :: pphi(ngrid,nlayer)   ! mid-layer geopotential
       real,intent(in) :: t(ngrid,nlayer) ! input temperature (K)
       real,intent(in) :: pdt(ngrid,nlayer) ! input tendency on temperature (K/s)      
@@ -74 +75 @@
       INTEGER i, k, n
       REAL zqs(ngrid,nlayer),Tsat(ngrid,nlayer), zdelta, zcor
-      REAL zrfl(ngrid), zrfln(ngrid), zqev, zqevt
+      REAL precip_rate(ngrid), precip_rate_tmp(ngrid), zqev, zqevt
 
       REAL zoliq(ngrid)
@@ -85 +86 @@
       real tnext(ngrid,nlayer)
 
-      real l2c(ngrid,nlayer)
+      real dmass(ngrid,nlayer)
       real dWtot
 
@@ -177 +178 @@
       d_q(1:ngrid,1:nlayer) = 0.0
       d_ql(1:ngrid,1:nlayer) = 0.0
-      zrfl(1:ngrid) = 0.0
-      zrfln(1:ngrid) = 0.0
+      precip_rate(1:ngrid) = 0.0
+      precip_rate_tmp(1:ngrid) = 0.0
 
       ! calculate saturation mixing ratio
@@ -192 +193 @@
       DO k = 1, nlayer
          DO i = 1, ngrid
-            l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
+            dmass(i,k)=(pplev(i,k)-pplev(i,k+1))/g
          ENDDO
       ENDDO
@@ -203 +204 @@
          IF (evap_prec) THEN ! note no rneb dependence!
             DO i = 1, ngrid
-               IF (zrfl(i) .GT.0.) THEN
+               IF (precip_rate(i) .GT.0.) THEN
 
                   if(zt(i,k).gt.Tsat(i,k))then
-!!                   treat the case where all liquid water should boil
-                     zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*l2c(i,k)/RLVTT/ptimestep,zrfl(i))
-                     zrfl(i)=MAX(zrfl(i)-zqev,0.)
-                     d_q(i,k)=zqev/l2c(i,k)*ptimestep
-                     d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
-                  else
-                     zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here
+!!                   treat the case where all liquid water should boil
+                     zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*dmass(i,k)/RLVTT/ptimestep,precip_rate(i))
+                     precip_rate(i)=MAX(precip_rate(i)-zqev,0.)
+                     d_q(i,k)=zqev/dmass(i,k)*ptimestep
+                     d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
+                  else
+                     zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmass(i,k)/ptimestep !there was a bug here
                      zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k))    & !default was 2.e-5
-                        *sqrt(zrfl(i))*l2c(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here
+                        *sqrt(precip_rate(i))*dmass(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here
                      zqevt = MAX (zqevt, 0.0)
                      zqev  = MIN (zqev, zqevt)
                      zqev  = MAX (zqev, 0.0)
-                     zrfln(i)= zrfl(i) - zqev
-                     zrfln(i)= max(zrfln(i),0.0)
-
-                     d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
+                     precip_rate_tmp(i)= precip_rate(i) - zqev
+                     precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0)
+
+                     d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmass(i,k)*ptimestep
                      !d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
                      d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
-                     zrfl(i)  = zrfln(i)
-                  end if
-                     
-
-               ENDIF ! of IF (zrfl(i) .GT.0.)
+                     precip_rate(i)  = precip_rate_tmp(i)
+                  end if
+#ifdef MESOSCALE 
+                  d_t(i,k) = d_t(i,k)+(pphi(i,k+1)-pphi(i,k))*precip_rate(i)*ptimestep/(RCPD*dmass(i,k))
+                  ! JL22. Accounts for gravitational energy of falling precipitations (probably not to be used in the GCM
+                  !   where the counterpart is not included in the dynamics.)
+#endif
+
+               ENDIF ! of IF (precip_rate(i) .GT.0.)
             ENDDO
          ENDIF ! of IF (evap_prec)
@@ -253 +258 @@
                   IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
                      d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
-                     zrfl(i)   = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
+                     precip_rate(i)   = precip_rate(i) - d_ql(i,k)*dmass(i,k)/ptimestep
                   ENDIF
                ENDIF
@@ -350 +355 @@
                IF (rneb(i,k).GT.0.0) THEN
                   d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
-                  zrfl(i)   = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
+                  precip_rate(i)   = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmass(i,k)/ptimestep
                ENDIF
             ENDDO
@@ -361 +366 @@
 !     Rain or snow on the ground
       DO i = 1, ngrid
-         if(zrfl(i).lt.0.0)then
+         if(precip_rate(i).lt.0.0)then
             print*,'Droplets of negative rain are falling...'
             call abort
          endif
          IF (t(i,1) .LT. T_h2O_ice_liq) THEN
-            dqssnow(i) = zrfl(i)
+            dqssnow(i) = precip_rate(i)
             dqsrain(i) = 0.0
          ELSE
             dqssnow(i) = 0.0
-            dqsrain(i) = zrfl(i) ! liquid water = ice for now
+            dqsrain(i) = precip_rate(i) ! liquid water = ice for now
          ENDIF
       ENDDO
@@ -381 +386 @@
            reevap_precip(i)=0.
            do k=1,nlayer
-              reevap_precip(i)=reevap_precip(i)+dqrain(i,k,i_vap)*l2c(i,k)
+              reevap_precip(i)=reevap_precip(i)+dqrain(i,k,i_vap)*dmass(i,k)
            enddo
         enddo