Changeset 5614 for LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90

Timestamp:

Apr 14, 2025, 9:21:07 PM (8 weeks ago)

Author:

evignon

Message:

Commission liée à un update majeur de la routine de condensation grande echelle suite au travail
de Lea, Audran et Etienne
Elle inclue une restructuration des routines pour clarifier le role "moniteur" de la routine lscp_main,
une mise à jour de la parametrisation de partitionnement de phase de Lea pour inclure les nuages de couche limite,
ainsi que des corrections des routines de precipitations "poprecip".
Convergence numerique verifiee en prod et debug pour les physiques NPv6.3 et 7.0.1c

File:

• LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90 (modified) (24 diffs)

LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90

@@ -330 +330 @@
 LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
 
+
 REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
 REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
@@ -540 +541 @@
           IF (zoliqi(i) .GT. 0.) THEN
             zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & 
-              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
+              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)/zneb(i)*zrho(i)*ffallv)**(-1./dice_velo))
           ELSE
             zfroi=0.
@@ -618 +619 @@
     ! Computation of DT if all the liquid precip freezes
     DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
+    
+
     ! T should not exceed the freezing point
     ! that is Delta > RTT-zt(i)
@@ -722 +725 @@
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
 USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds
-USE lmdz_lscp_ini, ONLY : eps, temp_nowater
+USE lmdz_lscp_ini, ONLY : eps, temp_nowater, ok_growth_precip_deposition
 USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
 
@@ -792 +795 @@
 dqreva(:) = 0.
 dqssub(:) = 0.
-dqrevap   = 0.
-dqssubl   = 0.
 
 !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
@@ -880 +881 @@
 DO i = 1, klon
 
+  dqrevap   = 0.
+  dqssubl   = 0.
   !--If there is precipitation from the layer above
   IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
@@ -913 +916 @@
       IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN
         qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
+      ELSE
+        qvapclr = 0.
+      ENDIF
+      IF (  precipfraccld(i)  .GT. eps ) THEN
+        qvapcld = MAX(qtotupnew(i)-qvapclrup(i) , 0.) / qtotupnew(i) * qvap(i) /  precipfraccld(i) 
+      ELSE
+        qvapcld = 0.
       ENDIF
     ELSE
@@ -918 +928 @@
       !--for the evap / subl process
       qvapclr = qvap(i)
+      qvapcld = qvap(i)
     ENDIF
 
@@ -940 +951 @@
       !--NB. with ok_ice_supersat activated, this barrier should be useless
       drainclreva = MIN(0., drainclreva)
-
+     
+      ! we set it to 0 as not sufficiently tested
+      drainclreva = 0.
 
       !--Sublimation of the solid precipitation coming from above
@@ -986 +999 @@
       ENDIF
 
+    ELSE
+            
+    !--All the precipitation is sublimated if the fraction is zero
+       drainclreva = - rainclr_tmp(i)
+       dsnowclrsub = - snowclr_tmp(i)
+
     ENDIF ! precipfracclr_tmp .GT. eps
 
@@ -993 +1012 @@
     !---------------------------
 
-    IF ( ok_unadjusted_clouds .AND. ( temp(i) .LE. temp_nowater ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
+    IF ( ok_growth_precip_deposition .AND. ( temp(i) .LE. RTT ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
       !--Evaporation of liquid precipitation coming from above
       !--in the cloud only
@@ -1000 +1019 @@
       !--Exact explicit formulation (raincld is resolved exactly, qvap explicitly)
       !--which does not need a barrier on raincld, because included in the formula
-      !draincldeva = precipfraccld_tmp(i) * MAX(0., &
-      !            - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
-      !            + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
-      !            )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
-               
+      
+      draincldeva = precipfraccld_tmp(i) * MAX(0., &
+                  - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
+                  + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
+                  )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
+
       !--Evaporation is limited by 0
       !--NB. with ok_ice_supersat activated, this barrier should be useless
-      !draincldeva = MIN(0., draincldeva)
-      draincldeva = 0.
+      draincldeva = MIN(0., draincldeva)
 
 
@@ -1277 +1296 @@
 
 USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
-                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & 
+                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc,               & 
                           thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
                           rho_rain, r_rain, r_snow, rho_ice,                   &
+                          expo_tau_auto_snow,                                  &
                           tau_auto_snow_min, tau_auto_snow_max,                &
-                          thresh_precip_frac, eps,                             &
+                          thresh_precip_frac, eps, rain_int_min,               &
                           gamma_melt, alpha_freez, beta_freez, temp_nowater,   &
                           iflag_cloudth_vert, iflag_rain_incloud_vol,          &
@@ -1340 +1360 @@
 REAL, DIMENSION(klon) :: dhum_to_dflux
 REAL, DIMENSION(klon) :: qtot                             !--includes vap, liq, ice and precip
+REAL                  :: min_precip                       !--minimum precip flux below which precip fraction decreases
 
 !--Collection, aggregation and riming
@@ -1502 +1523 @@
         !--tau for snow depends on the ice fraction in mixed-phase clouds     
         tau_auto_snow = tau_auto_snow_max &
-                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
 
         expo_auto_rain = cld_expo_con
@@ -1513 +1534 @@
         !--tau for snow depends on the ice fraction in mixed-phase clouds     
         tau_auto_snow = tau_auto_snow_max &
-                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
 
         expo_auto_rain = cld_expo_lsc
@@ -1534 +1555 @@
                - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
 
-
     !--Barriers so that we don't create more rain/snow
     !--than there is liquid/ice
@@ -1543 +1563 @@
     qliq(i) = qliq(i) + dqlauto
     qice(i) = qice(i) + dqiauto
+
     raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
     snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
@@ -1710 +1731 @@
   !--second: immersion freezing following (inspired by Bigg 1953)
   !--the latter is parameterized as an exponential decrease of the rain
-  !--water content with a homemade formulya
+  !--water content with a homemade formula
   !--This is based on a caracteritic time of freezing, which
   !--exponentially depends on temperature so that it is
@@ -1717 +1738 @@
   !--NB.: this process needs a temperature adjustment
   !--dqrfreez_max : maximum rain freezing so that temperature
-  !--              stays lower than 273 K [kg/kg]
+  !--               stays lower than 273 K [kg/kg]
   !--tau_freez    : caracteristic time of freezing [s]
   !--gamma_freez  : tuning parameter [s-1]
@@ -1798 +1819 @@
               * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
 
-
     !--In clear air
     IF ( rainclr(i) .GT. 0. ) THEN
@@ -1827 +1847 @@
     !--Add tendencies
     !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
+    
     rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
     raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
@@ -1833 +1854 @@
 
 
+
     !--Temperature adjustment with the uptake of latent
     !--heat because of freezing
+
     temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
                       / ( 1. + RVTMP2 * qtot(i) )
-
     !--Diagnostic tendencies
     dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2          
@@ -1847 +1869 @@
 
 
-  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
-  !--we reduce the precipiration fraction in the clear/cloudy air so that the new
-  !--local flux of rain+snow is equal to rain_int_min.
+  !--If the local flux of rain+snow in clear air is lower than min_precip,
+  !--we reduce the precipiration fraction in the clear air so that the new
+  !--local flux of rain+snow is equal to min_precip.
+  !--we apply the minimum only on the clear-sky fraction because the cloudy precip fraction
+  !--already decreases out of clouds
   !--Here, rain+snow is the gridbox-mean flux of precip.
   !--Therefore, (rain+snow)/precipfrac is the local flux of precip.
-  !--If the local flux of precip is lower than rain_int_min, i.e.,
-  !-- (rain+snow)/precipfrac < rain_int_min , i.e., 
-  !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
-  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
+  !--If the local flux of precip is lower than min_precip, i.e.,
+  !-- (rain+snow)/precipfrac < min_precip , i.e., 
+  !-- (rain+snow)/min_precip < precipfrac , then we want to reduce
+  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/min_precip.
   !--Note that this is physically different than what is proposed in LTP thesis.
-  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
-  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
+  !--min_precip is either equal to rain_int_min or calculated as a very small fraction
+  !--of the minimum precip flux estimated as the flux associated with the 
+  !--autoconversion threshold mass content
+  !min_precip=1.e-6*(pplay(i)/RD/temp(i))*MIN(rain_fallspeed_clr*cld_lc_lsc,snow_fallspeed_clr*cld_lc_lsc_snow)
+  min_precip=rain_int_min
+  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / min_precip )
 
   !--Calculate outputs