Changeset 4830 for LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90

Timestamp:

Feb 22, 2024, 5:29:02 PM (3 months ago)

Author:

evignon

Message:

changements suite à l'atelier nuage d'aujourd'hui.

File:

• LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90 (modified) (25 diffs)

LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90

@@ -23 +23 @@
 
 USE lmdz_lscp_ini, ONLY : prt_level, lunout
-USE lmdz_lscp_ini, ONLY : coef_eva, coef_eva_i, expo_eva, expo_eva_i, thresh_precip_frac
+USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
 USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
@@ -149 +149 @@
     !--Sublimation of the solid precipitation coming from above
     !--(same formula as for liquid precip)
-    dsnowsub = - precipfracclr(i) * coef_eva_i * (1. - qvap(i) / qsati(i)) &
-             * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva_i &
+    dsnowsub = - precipfracclr(i) * coef_sub * (1. - qvap(i) / qsati(i)) &
+             * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_sub &
              * temp(i) * RD / pplay(i) &
              * ( paprsdn(i) - paprsup(i) ) / RG
@@ -229 +229 @@
 ! - aggregation (agg)
 ! - riming (rim)
-! - collection (coll)
+! - collection (col)
 ! - melting (melt)
-! - freezing (free)
+! - freezing (freez)
 ! 
 SUBROUTINE poprecip_postcld( &
@@ -238 +238 @@
            precipfracclr, precipfraccld, &
            rain, rainclr, raincld, snow, snowclr, snowcld, &
-           qrain, qsnow, dqrauto, dqrcol, dqrmelt, dqrfreez, &
+           qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
            dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
 
@@ -251 +251 @@
                           tau_auto_snow_min, tau_auto_snow_max,                &
                           thresh_precip_frac, eps, air_thermal_conduct,        &
-                          coef_ventil, alpha_freez, gamma_freez, temp_nowater, &
-                          iflag_cloudth_vert, iflag_rain_incloud_vol
+                          coef_ventil, alpha_freez, beta_freez, temp_nowater,  &
+                          iflag_cloudth_vert, iflag_rain_incloud_vol,          &
+                          cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez,       &
+                          rain_fallspeed_clr, rain_fallspeed_cld,              &
+                          snow_fallspeed_clr, snow_fallspeed_cld
 
 
@@ -286 +289 @@
 REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
 
-REAL,    INTENT(OUT),   DIMENSION(klon) :: qrain          !--specific rain content [kg/kg]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnow          !--specific snow content [kg/kg]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
@@ -306 +309 @@
 REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip
 
-! partition of the fluxes
+!--Partition of the fluxes
 REAL :: dcldfra
 REAL :: precipfractot
@@ -313 +316 @@
 REAL :: draincld, dsnowcld
 
-! collection, aggregation and riming
+!--Collection, aggregation and riming
 REAL :: eff_cldfra
 REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
 REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
-REAL :: dqlcol           ! loss of liquid cloud content due to collection by rain [kg/kg/s]
-REAL :: dqiagg           ! loss of ice cloud content due to collection by aggregation [kg/kg/s]
-REAL :: dqlrim           ! loss of liquid cloud content due to riming on snow[kg/kg/s]
-
-! autoconversion
+REAL :: dqlcol           !--loss of liquid cloud content due to collection by rain [kg/kg/s]
+REAL :: dqiagg           !--loss of ice cloud content due to collection by aggregation [kg/kg/s]
+REAL :: dqlrim           !--loss of liquid cloud content due to riming on snow [kg/kg/s]
+
+!--Autoconversion
 REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
 REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
-REAL :: dqlauto          ! loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
-REAL :: dqiauto          ! loss of ice cloud content due to autoconversion to snow [kg/kg/s]
-
-! melting
+REAL :: dqlauto          !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
+REAL :: dqiauto          !--loss of ice cloud content due to autoconversion to snow [kg/kg/s]
+
+!--Melting
 REAL :: dqsmelt_max
-REAL :: fallice_clr, fallice_cld
 REAL :: nb_snowflake_clr, nb_snowflake_cld
 REAL :: capa_snowflake, temp_wetbulb
@@ -335 +337 @@
 REAL, DIMENSION(klon) :: qzero, qsat, dqsat
 
-! freezing
+!--Freezing
 REAL :: dqrfreez_max
 REAL :: tau_freez
 REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez
+REAL :: coef_freez
+REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
+REAL :: Eff_rain_ice
 
 
@@ -358 +363 @@
 DO i = 1, klon
 
-  ! variables initialisation
-  dqlrim  = 0.
+  !--Variables initialisation
   dqlcol  = 0.
   dqiagg  = 0.
   dqiauto = 0.
   dqlauto = 0.
+  dqlrim  = 0.
+  dqifreez= 0.
 
   !--hum_to_flux: coef to convert a specific quantity to a flux 
@@ -401 +407 @@
   !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
   !--calculation of the current CS precip. frac.
-  dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
+  !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
+  !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
+  !--if precipfractot < cldfra)
+  dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i)
   !--Tendency of the cloudy precipitation fraction. We add a MAX on the
   !--calculation of the current CS precip. frac.
   !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) ) 
-  !--We remove it, because cldfra is guaranteed to be > O (the MAX is activated
+  !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
   !--if cldfra < 0)
   dprecipfraccld = dcldfra
@@ -448 +457 @@
   
 
-  ! if vertical heterogeneity is taken into account, we use
-  ! the "true" volume fraction instead of a modified 
-  ! surface fraction (which is larger and artificially
-  ! reduces the in-cloud water).
+  !--If vertical heterogeneity is taken into account, we use
+  !--the "true" volume fraction instead of a modified 
+  !--surface fraction (which is larger and artificially
+  !--reduces the in-cloud water).
   IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
     eff_cldfra = ctot_vol(i)
@@ -459 +468 @@
 
 
-  ! Start precipitation growth processes
+  !--Start precipitation growth processes
 
   !--If the cloud is big enough, the precipitation processes activate
@@ -483 +492 @@
     !--get in-cloud mean quantities. The two divisions by eff_cldfra are
     !--then simplified.
+
+    !--The sticking efficacy is perfect.
     Eff_rain_liq = 1.
     coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
     IF ((raincld(i) .GT. 0.) .AND. (coef_col .GT. 0.)) THEN
       !-- ATTENTION Double implicit version
+      !-- BEWARE the formule below is FALSE (because raincld is a flux, not a delta flux)
       !qrain_tmp = raincld(i) / hum_to_flux(i)
       !coef_tmp = coef_col * dtime * ( qrain_tmp / precipfraccld(i) + qliq(i) / eff_cldfra )
@@ -495 +507 @@
       !--available.
       !dqlcol = MAX( - qliq(i), dqlcol )
-      !--Exact version, which does not need a barrier because of
+      !--Exact explicit version, which does not need a barrier because of
       !--the exponential decrease
       dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
@@ -517 +529 @@
       !--available.
       !dqiagg = MAX( - qice(i), dqiagg )
-      !--Exact version, which does not need a barrier because of
+      !--Exact explicit version, which does not need a barrier because of
       !--the exponential decrease
       dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
@@ -536 +548 @@
     !--rain drops/snowflakes. It relies on the formulations by
     !--Sundqvist 1978.
-    ! TODO what else?
 
     !--If we are in a convective point, we have different parameters
     !--for the autoconversion
     IF ( ptconv(i) ) THEN
-        ! ATTENTION voir les constantes ici qui ne devraient pas etre identiques
         qthresh_auto_rain = cld_lc_con
-        qthresh_auto_snow = cld_lc_con
+        qthresh_auto_snow = cld_lc_con_snow
 
         tau_auto_rain = cld_tau_con
@@ -553 +563 @@
         expo_auto_snow = cld_expo_con
     ELSE
-        ! ATTENTION voir les constantes ici qui ne devraient pas etre identiques
         qthresh_auto_rain = cld_lc_lsc
-        qthresh_auto_snow = cld_lc_lsc
+        qthresh_auto_snow = cld_lc_lsc_snow
 
         tau_auto_rain = cld_tau_lsc
@@ -568 +577 @@
 
     ! Liquid water quantity to remove according to (Sundqvist, 1978)
-    ! dqliq/dt = -qliq/tau * ( 1-exp(-(qcin/clw)**2) )
-    !.........................................................
-    ! we first treat the second term (with exponential) in an explicit way
-    ! and then treat the first term (-q/tau) in an exact way
+    ! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
+    !
+    !--And same formula for ice
+    !
+    !--We first treat the second term (with exponential) in an explicit way
+    !--and then treat the first term (-q/tau) in an exact way
 
     dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
@@ -602 +613 @@
     !--snow (graupel in fact) because of the collision between
     !--those and falling snowflakes.
-    !--The formula come from Muench and Lohmann 2020
+    !--The formula comes from Muench and Lohmann 2020
     !--NB.: this process needs a temperature adjustment
 
@@ -631 +642 @@
     ENDIF
 
-  ! ATTENTION veut on faire un processus similaire et symetrique qui
-  ! convertirait la pluie en surfusion qui tombe sur des cristaux de glace
-  ! en flux de neige ? (si la temperature est negative)
-
   ENDIF ! cldfra .GE. seuil_neb
 
@@ -658 +665 @@
   !--                the radius if the snowflake is a sphere [m]
   !--temp_wetbulb: wet-bulb temperature [K]
-  !--fallice: snow fall velocity (in clear/cloudy sky) [m/s]
+  !--snow_fallspeed: snow fall velocity (in clear/cloudy sky) [m/s]
   !--air_thermal_conduct: thermal conductivity of the air [J/m/K/s]
   !--coef_ventil: ventilation coefficient [-]
@@ -682 +689 @@
     !--In clear air
     IF ( snowclr(i) .GT. 0. ) THEN
-      ! ATTENTION fallice a definir
-      fallice_clr = 1.
       !--Calculated according to
       !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
-      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / fallice_clr &
+      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
                        / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
       dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
                    * capa_snowflake / RLMLT * coef_ventil &
                    * MAX(0., temp_wetbulb - RTT) * dtime
+
+      !--Barrier to limit the melting flux to the clr snow flux in the mesh
       dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / hum_to_flux(i))
-      ELSE
-      dqsclrmelt = 0.
    ENDIF
 
     !--In cloudy air
     IF ( snowcld(i) .GT. 0. ) THEN
-      ! ATTENTION fallice a definir
-      fallice_cld = 1.
       !--Calculated according to
       !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
-      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / fallice_cld &
+      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
                        / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
       dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
                    * capa_snowflake / RLMLT * coef_ventil &
                    * MAX(0., temp_wetbulb - RTT) * dtime
+
+      !--Barrier to limit the melting flux to the cld snow flux in the mesh
       dqscldmelt = MAX( dqscldmelt , -snowcld(i) / hum_to_flux(i))
-    ELSE
-      dqscldmelt = 0.
     ENDIF
     
-    !--Barriers
+    !--Barrier on temperature. If the total melting flux leads to a
+    !--positive temperature, it is limited to keep temperature above 0 degC.
     !--It is activated if the total is LOWER than the max
     !--because everything is negative
@@ -769 +773 @@
                          * ( 1. + RVTMP2 * qtot(i) ))
  
-    tau_freez = 1. / ( gamma_freez &
+    tau_freez = 1. / ( beta_freez &
               * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
+
+    !--Initialisation 
+    dqrclrfreez = 0.
+    dqrcldfreez = 0.
 
     !--In clear air
     IF ( rainclr(i) .GT. 0. ) THEN
-       !--Exact solution of dqrain/dt = -qrain/tau_freez
-       dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
-       dqrclrfreez = MAX(dqrclrfreez, -rainclr(i) / hum_to_flux(i))
-    ELSE
-       dqrclrfreez = 0.
+      !--Exact solution of dqrain/dt = -qrain/tau_freez
+      dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
     ENDIF
 
     !--In cloudy air
     IF ( raincld(i) .GT. 0. ) THEN
-       !--Exact solution of dqrain/dt = -qrain/tau_freez
-       dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
-       dqrcldfreez = MAX(dqrcldfreez, -raincld(i) / hum_to_flux(i))
-    ELSE
-       dqrcldfreez = 0.
+      !--Exact solution of dqrain/dt = -qrain/tau_freez
+      dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
+
+      !---------------------------------------------------------
+      !--    FREEZING OF RAIN WITH CONTACT OF ICE CRYSTALS
+      !---------------------------------------------------------
+      !--Sub-process of freezing which quantifies the collision between
+      !--ice crystals in suspension and falling rain droplets.
+      !--The rain droplets freeze, becoming graupel, and carrying
+      !--the ice crystal (which acted as an ice nucleating particle).
+      !--The formula is adapted from the riming formula.
+
+      !--The sticking efficacy is perfect.
+      Eff_rain_ice = 1.
+      coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
+      !-- ATTENTION Double implicit version? + barriers if needed
+      !--Exact version, which does not need a barrier because of
+      !--the exponential decrease
+      dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
+
+      !--We add the two freezing processes in cloud
+      dqrcldfreez = dqrcldfreez - dqifreez
     ENDIF
 
@@ -838 +860 @@
 
   !--Diagnostics
-  ! ATTENTION A REPRENDRE
-  qrain(i) = rain(i) / hum_to_flux(i)
-  qsnow(i) = snow(i) / hum_to_flux(i)
+  !--BEWARE this is indeed a diagnostic: this is an estimation from
+  !--the value of the flux at the bottom interface of the mesh and
+  !--and assuming an upstream numerical calculation
+  !--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
+  !--used for computing the total ice water content in the mesh
+  !--for radiation only
+  qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
+                 + raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
+  qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
+                 + snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
 
 ENDDO ! loop on klon