Changeset 4898

Timestamp:

Apr 11, 2024, 5:50:38 PM (7 months ago)

Author:

evignon

Message:

nettoyage poprecip + changements de valeur par defaut du seuil
pour l'autoconversion solide dans poprecip
Lea, Audran, Etienne

Location:

LMDZ6/trunk/libf/phylmd

Files:

• lmdz_lscp_ini.F90 (modified) (4 diffs)
• lmdz_lscp_poprecip.F90 (modified) (24 diffs)

LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90

@@ -140 +140 @@
   !$OMP THREADPRIVATE(ok_poprecip)
 
-  REAL, SAVE, PROTECTED :: cld_lc_lsc_snow               ! snow autoconversion coefficient, stratiform rain
+  REAL, SAVE, PROTECTED :: cld_lc_lsc_snow=2.e-5            ! snow autoconversion coefficient, stratiform. default from  Chaboureau and PInty 2006
   !$OMP THREADPRIVATE(cld_lc_lsc_snow)
 
-  REAL, SAVE, PROTECTED :: cld_lc_con_snow                ! snow autoconversion coefficient, convective rain
+  REAL, SAVE, PROTECTED :: cld_lc_con_snow=2.e-5            ! snow autoconversion coefficient, convective 
   !$OMP THREADPRIVATE(cld_lc_con_snow)
 
@@ -166 +166 @@
   REAL, SAVE, PROTECTED :: rho_ice=920.                    ! A COMMENTER TODO [kg/m3]
   !$OMP THREADPRIVATE(rho_ice)
-
-  REAL, SAVE, PROTECTED :: rho_snow                        ! A COMMENTER TODO [kg/m3]
-  !$OMP THREADPRIVATE(rho_snow)
 
   REAL, SAVE, PROTECTED :: r_rain=500.E-6                   ! A COMMENTER TODO [m]
@@ -279 +276 @@
     CALL getin_p('cld_lc_lsc',cld_lc_lsc)
     CALL getin_p('cld_lc_con',cld_lc_con)
-    cld_lc_lsc_snow=cld_lc_lsc
-    cld_lc_con_snow=cld_lc_con
     CALL getin_p('cld_lc_lsc_snow',cld_lc_lsc_snow)
     CALL getin_p('cld_lc_con_snow',cld_lc_con_snow)
@@ -393 +388 @@
 
 
-    !--Initialisations of constants depending on other constants
-    !--rho_snow formula from r_snow (Brandes et al. 2007 - JAMC)
-    rho_snow = 1.e3 * 0.178 * (r_snow * 2. * 1000.)**(-0.922)
-
-
     !AA Temporary initialisation
     a_tr_sca(1) = -0.5

LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90

@@ -63 +63 @@
 !--Integer for interating over klon
 INTEGER :: i
-!--dhum_to_dflux: coef to convert a specific quantity to a flux 
+!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
 REAL, DIMENSION(klon) :: dhum_to_dflux
 !--
@@ -139 +139 @@
     !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
     !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
-    !--Explicit formulation
-    !draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) * dz(i) &
-    !         * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva &
-    !draineva = MAX( - rainclr(i), draineva)
     !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
     !--which does not need a barrier on rainclr, because included in the formula
@@ -149 +145 @@
              + ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_eva ) &
                ) )**( 1. / ( 1. - expo_eva ) ) - rainclr(i)
+    ! TODO veut on vraiment mettre qvap ici, qui est en fait qtot dans la maille ?
              
     !--Evaporation is limited by 0
@@ -156 +153 @@
     !--Sublimation of the solid precipitation coming from above
     !--(same formula as for liquid precip)
-    !--Explicit formulation
-    !dsnowsub = - precipfracclr(i) * coef_sub * (1. - qvap(i) / qsatl(i)) * dz(i) &
-    !         * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_sub &
-    !dsnowsub = MAX( - snowclr(i), dsnowsub)
     !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly)
     !--which does not need a barrier on snowclr, because included in the formula
@@ -166 +159 @@
              + ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_sub ) &
              ) )**( 1. / ( 1. - expo_sub ) ) - snowclr(i)
+    ! TODO veut on vraiment mettre qvap ici, qui est en fait qtot dans la maille ?
 
     !--Sublimation is limited by 0
@@ -178 +172 @@
     !--It is expressed as a max flux dprecip_evasub_max
     
+    ! TODO veut on vraiment mettre qvap ici, qui est en fait qtot dans la maille ?
     dprecip_evasub_max = MIN(0., ( qvap(i) - qsat(i) ) * precipfracclr(i)) &
                      * dhum_to_dflux(i)
@@ -261 +256 @@
                           cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & 
                           thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
-                          rho_rain, rho_snow, r_rain, r_snow, rho_ice, r_ice,  &
+                          rho_rain, r_rain, r_snow, rho_ice, r_ice,            &
                           tau_auto_snow_min, tau_auto_snow_max,                &
                           thresh_precip_frac, eps,                             &
@@ -319 +314 @@
 
 INTEGER :: i
+!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
 REAL, DIMENSION(klon) :: dhum_to_dflux
-REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip
+REAL, DIMENSION(klon) :: qtot                             !--includes vap, liq, ice and precip
 
 !--Partition of the fluxes
@@ -333 +329 @@
 REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
 REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
+REAL :: rho_snow
 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]
@@ -384 +381 @@
   dqlrim  = 0.
 
-  !--dhum_to_dflux: coef to convert a delta in specific quantity to a flux 
   !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
   dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
@@ -492 +488 @@
 
   !--If the cloud is big enough, the precipitation processes activate
+  ! TODO met on seuil_neb ici ?
   IF ( cldfra(i) .GE. seuil_neb ) THEN
 
@@ -518 +515 @@
     coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
     IF ( raincld(i) .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) / dhum_to_dflux(i)
-      !coef_tmp = coef_col * dtime * ( qrain_tmp / precipfraccld(i) + qliq(i) / eff_cldfra )
-      !dqlcol = qliq(i) * ( 1. / ( 1. + 0.5 * ( coef_tmp - 1. + SQRT( &
-      !    ( 1. - coef_tmp )**2. + 4. * coef_col * dtime * qrain_tmp / precipfraccld(i) ) &
-      !                   ) ) - 1. )
-      !--Barriers so that the processes do not consume more liquid/ice than
-      !--available.
-      !dqlcol = MAX( - qliq(i), dqlcol )
       !--Exact explicit version, which does not need a barrier because of
       !--the exponential decrease
@@ -544 +531 @@
     !--it s a product of a collection efficiency and a sticking efficiency
     Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
+    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
+    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
     coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
     IF ( snowcld(i) .GT. 0. ) THEN
-      !-- ATTENTION Double implicit version?
-      !--Barriers so that the processes do not consume more liquid/ice than
-      !--available.
-      !dqiagg = MAX( - qice(i), dqiagg )
       !--Exact explicit version, which does not need a barrier because of
       !--the exponential decrease
@@ -612 +597 @@
 
 
-    !--Barriers so that we don t create more rain/snow
+    !--Barriers so that we don't create more rain/snow
     !--than there is liquid/ice
     dqlauto = MAX( - qliq(i), dqlauto )
@@ -632 +617 @@
     !---------------------------------------------------------
     !--Process which converts liquid droplets in suspension into
-    !--snow (graupel in fact) because of the collision between
+    !--snow because of the collision between
     !--those and falling snowflakes.
     !--The formula comes from Muench and Lohmann 2020
@@ -640 +625 @@
     !--assuming a cloud droplet diameter of 20 microns.
     Eff_snow_liq = 0.2
+    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
+    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
     coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
     IF ( snowcld(i) .GT. 0. ) THEN
-      !-- ATTENTION Double implicit version?
-      !--Barriers so that the processes do not consume more liquid than
-      !--available.
-      !dqlrim = MAX( - qliq(i), dqlrim )
       !--Exact version, which does not need a barrier because of
       !--the exponential decrease
@@ -682 +665 @@
   !--NB.: this process needs a temperature adjustment
 
-  !--dqsmelt_max: maximum snow melting so that temperature
-  !--             stays higher than 273 K [kg/kg]
-  !--capa_snowflake: capacitance of a snowflake, equal to
-  !--                the radius if the snowflake is a sphere [m]
-  !--temp_wetbulb: wet-bulb temperature [K]
-  !--snow_fallspeed: snow fall velocity (in clear/cloudy sky) [m/s]
-  !--air_thermal_conduct: thermal conductivity of the air [J/m/K/s]
-  !--gamma_melt: melting tuning parameter [-]
-  !--nb_snowflake: number of snowflakes (in clear/cloudy air) [-]
+  !--dqsmelt_max         : maximum snow melting so that temperature
+  !--                      stays higher than 273 K [kg/kg]
+  !--capa_snowflake      : capacitance of a snowflake, equal to
+  !--                      the radius if the snowflake is a sphere [m]
+  !--temp_wetbulb        : wet-bulb temperature [K]
+  !--snow_fallspeed      : snow fall velocity (in clear/cloudy sky) [m/s]
+  !--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
+  !--gamma_melt          : tuning parameter for melting [-]
+  !--nb_snowflake        : number of snowflakes (in clear/cloudy air) [-]
 
   IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
@@ -705 +688 @@
     ! ATTENTION POUR L'INSTANT C'EST UN WBULB SELON FORMULE ECMWF
     ! ATTENTION EST-CE QVAP ?????
+    ! TODO veut on vraiment mettre qvap ici, on pourrait vouloir qvap dans le nuage / clr sky ?
+    ! (a mettre dans les deux IF)
     temp_wetbulb = temp(i) - ( qsat(i) - qvap(i) ) &
                  * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
                  - 40.637 * ( temp(i) - 275. ) )
-    air_thermal_conduct=(5.69+0.017*(temp(i)-RTT))*1.e-3*4.184                ! thermal conductivity of the air, SI
-
+    
+    !--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
+    air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184    
+
+    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
+    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
 
     !--In clear air
@@ -779 +768 @@
   !--Process through which rain freezes into snow. 
   !-- We parameterize it as a 2 step process:
-  !-- first: freezing following collision with ice crystals
-  !-- 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
+  !--first: freezing following collision with ice crystals
+  !--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
   !--This is based on a caracteritic time of freezing, which
   !--exponentially depends on temperature so that it is
@@ -788 +777 @@
   !--0 for RTT (=273.15 K).
   !--NB.: this process needs a temperature adjustment
-  !--dqrfreez_max: maximum rain freezing so that temperature
+  !--dqrfreez_max : maximum rain freezing so that temperature
   !--              stays lower than 273 K [kg/kg]
-  !--tau_freez: caracteristic time of freezing [s]
-  !--gamma_freez: tuning parameter [s-1]
-  !--alpha_freez: tuning parameter for the shape of the exponential curve [-]
-  !--temp_nowater: temperature below which no liquid water exists [K] (about -40 degC)
+  !--tau_freez    : caracteristic time of freezing [s]
+  !--gamma_freez  : tuning parameter [s-1]
+  !--alpha_freez  : tuning parameter for the shape of the exponential curve [-]
+  !--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
 
   IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
 
-      
+ 
     !--1st step: freezing following collision with ice crystals
     !--Sub-process of freezing which quantifies the collision between
@@ -804 +793 @@
     !--the ice crystal (which acted as an ice nucleating particle).
     !--The formula is adapted from the riming formula.
-    !-- it works only in the cloudy part
+    !--it works only in the cloudy part
 
     dqifreez = 0.
@@ -818 +807 @@
                          * ( 1. + RVTMP2 * qtot(i) ))
  
-      !--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 collision efficiency is assumed unity 
@@ -832 +816 @@
 
       !--We add the part of rain water that freezes, limited by a temperature barrier
-      !--this quantity is calculated assuming that the number of drop that freeze correspond to the number
+      !--This quantity is calculated assuming that the number of drop that freeze correspond to the number
       !--of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
       !--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3
@@ -882 +866 @@
     ENDIF
 
-    !--temperature barrie step 2
+    !--temperature barrier step 2
     !--It is activated if the total is LOWER than the max
     !--because everything is negative