Changeset 5615 for LMDZ6/branches

Timestamp:

Apr 14, 2025, 11:29:13 PM (3 months ago)

Author:

aborella

Message:

Comments + small bugfixes

Location:

LMDZ6/branches/contrails/libf/phylmd

Files:

• lmdz_aviation.f90 (modified) (2 diffs)
• lmdz_lscp.f90 (modified) (2 diffs)
• lmdz_lscp_condensation.f90 (modified) (18 diffs)
• lmdz_lscp_ini.f90 (modified) (6 diffs)

LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90

@@ -179 +179 @@
         pdf_q_above_qcritcont = 0.
       ELSEIF ( pdf_y .LT. -10. ) THEN
-        pdf_fra_above_qcritcont = 1.
-        pdf_q_above_qcritcont = qclr(i) / clrfra(i)
+        pdf_fra_above_qcritcont = clrfra(i)
+        pdf_q_above_qcritcont = qclr(i)
       ELSE
         pdf_y = EXP( pdf_y )
@@ -196 +196 @@
         pdf_q_above_qsat = 0.
       ELSEIF ( pdf_y .LT. -10. ) THEN
-        pdf_fra_above_qsat = 1.
-        pdf_q_above_qsat = qclr(i) / clrfra(i)
+        pdf_fra_above_qsat = clrfra(i)
+        pdf_q_above_qsat = qclr(i)
       ELSE
         pdf_y = EXP( pdf_y )

LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90

@@ -334 +334 @@
   REAL, DIMENSION(klon) :: contfra, perscontfra, qcont
   LOGICAL, DIMENSION(klon) :: pt_pron_clds
-  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrail)
-  ! Constants used for calculating ratios that are advected (using a parent-child
-  ! formalism). This is not done in the dynamical core because at this moment,
-  ! only isotopes can use this parent-child formalism. Note that the two constants
-  ! are the same as the one use in the dynamical core, being also defined in
-  ! dyn3d_common/infotrac.F90
-  REAL :: min_qParent, min_ratio
   !--for Lamquin et al 2012 diagnostics
   REAL, DIMENSION(klon) :: issrfra100to150UP, issrfra150to200UP, issrfra200to250UP
@@ -446 +439 @@
 qvc(:)          = 0.
 shear(:)        = 0.
-min_qParent     = 1.e-30
-min_ratio       = 1.e-16
+pt_pron_clds(:) = .FALSE.
 
 !--for Lamquin et al (2012) diagnostics

LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90

@@ -124 +124 @@
 
 USE lmdz_lscp_ini,   ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
-USE lmdz_lscp_ini,   ONLY: mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp
+USE lmdz_lscp_ini,   ONLY: N_ice_volume, corr_incld_depsub, nu_iwc_pdf_lscp
+USE lmdz_lscp_ini,   ONLY: beta_pdf_lscp, temp_thresh_pdf_lscp
 USE lmdz_lscp_ini,   ONLY: std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp
-USE lmdz_lscp_ini,   ONLY: coef_mixing_lscp, coef_shear_lscp
+USE lmdz_lscp_ini,   ONLY: coef_mixing_lscp, coef_shear_lscp, incld_ice_thresh
 
 USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, aspect_ratio_contrails
@@ -241 +242 @@
 REAL :: pres_sat, kappa_depsub, tauinv_depsub
 REAL :: air_thermal_conduct, water_vapor_diff
-REAL :: N_ice
 !
 ! for dissipation
@@ -402 +402 @@
       !--HYPOTHESIS 1: the ice crystals are spherical, therefore
       !-- mi = 4/3 * pi * rm_ice**3 * rho_ice
-      !--HYPOTHESIS 2: the ice crystals are monodisperse with the
-      !--initial radius rm_ice_0.
-      !--NB. this is notably different than the assumption
-      !--of a distributed qice in the cloud made in the sublimation process;
-      !--should it be consistent?
+      !--HYPOTHESIS 2: the ice crystals concentration is constant in the cloud
+      !
+      !--The equation in terms of q_ice is valide locally, and the local ice water content
+      !--follows a Gamma distribution with a factor nu_iwc_pdf_lscp. Therefore, by
+      !--integrating the local equation over the PDF (entire cloud), a correcting factor
+      !--must be included, equal to
+      !-- corr_incld_depsub = GAMMA(nu + 1/3) / GAMMA(nu) / nu**(1/3)
+      !--NB. this is equal to about 0.9, hence the correction is not big
+      !--NB. to lighten the calculated, corr_incld_depsub is calculated in lmdz_lscp_ini
       !
       !--As the deposition process does not create new ice crystals,
@@ -412 +416 @@
       !--therefore the sublimation process does not destroy ice crystals
       !--(or, in a limit case, it destroys all ice crystals), then
-      !--Nice is a constant during the sublimation/deposition process.
-      !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
+      !--Nice is a constant during the sublimation/deposition process
+      !--hence dmi = dqi
       !
-      !--The deposition equation then reads:
-      !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
-      !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
+      !--The deposition equation then reads for qi the in-cloud ice water content:
+      !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat * corr_incld_depsub &
+      !--            / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
+      !-- dqi/dt = alpha*4pi*capa_cond_cirrus*Nice*corr_incld_depsub &
+      !--             / ( 4pi/3 N_ice rho_ice )**(1/3) &
       !--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
-      !--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
-      !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
-      !--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
+      !--             qi**(1/3) * (qvc - qsat) / qsat
       !--and we have
-      !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
-      !-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
-      !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
+      !-- dqvc/dt = - alpha * kappa(T) * qi**(1/3) * (qvc - qsat)
+      !-- dqi/dt  =   alpha * kappa(T) * qi**(1/3) * (qvc - qsat)
       !
       !--This system of equations can be resolved with an exact
       !--explicit numerical integration, having one variable resolved
-      !--explicitly, the other exactly. The exactly resolved variable is
-      !--the one decreasing, so it is qvc if the process is
-      !--condensation, qi if it is sublimation.
+      !--explicitly, the other exactly. qvc is always the variable solved exactly.
       !
       !--kappa is computed as an initialisation constant, as it depends only
@@ -440 +441 @@
       !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
       water_vapor_diff = 0.211 * ( temp(i) / RTT )**1.94 * ( 101325. / pplay(i) ) * 1.e-4
-      !--Median ice crystal concentration [#/m3] in cirrus clouds from Kramer et al (2020)
-      N_ice = 0.003 * 1e6
-      !--NB. the greater kappa_depsub, the lower the efficiency of the
+      !--NB. the greater kappa_depsub, the more efficient is the
       !--deposition/sublimation process
-      kappa_depsub = ( 4. / 3. * RPI * N_ice / rho * rho_ice )**(1./3.) &
-                   * qsat(i) / ( 4. * RPI * capa_cond_cirrus * N_ice / rho ) &
-                   * ( RV * temp(i) / water_vapor_diff / pres_sat &
+      kappa_depsub = 4. * RPI * capa_cond_cirrus * N_ice_volume / rho * corr_incld_depsub &
+                   / qsat(i) / ( 4. / 3. * RPI * N_ice_volume / rho * rho_ice )**(1./3.) &
+                   / ( RV * temp(i) / water_vapor_diff / pres_sat &
                      + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
 
@@ -488 +487 @@
         
         !--If the ice water content is too low, the cloud is purely sublimated
-        IF ( qiceincld .LT. ( 0.005 * qsat(i) ) ) THEN
+        IF ( qiceincld .LT. eps ) THEN
           dcfa_sub(i) = - contfra(i)
           dqia_sub(i) = - qiceincld * contfra(i)
@@ -530 +529 @@
         ELSE
 
-          !IF ( qiceincld .LT. ( 0.005 * qsat(i) ) ) THEN
-          !  dcf_sub(i) = ( qiceincld / ( 0.005 * qsat(i) ) - 1. ) * cldfra(i)
-          !  dqvc_sub(i) = qvapincld * dcf_sub(i)
-
-          !  cldfra(i) = cldfra(i) + dcf_sub(i)
-          !  qcld(i) = qcld(i) + dqvc_sub(i)
-          !  qvc(i) = qvc(i) + dqvc_sub(i)
-          !  clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
-          !  qclr(i) = qclr(i) - dqvc_sub(i)
-          !ENDIF
+          !--Cirrus clouds cannot have an in-cloud ice water content lower than
+          !--incld_ice_thresh times the saturation
+          IF ( qiceincld .LT. ( incld_ice_thresh * qsat(i) ) ) THEN
+            dcf_sub(i) = ( qiceincld / ( incld_ice_thresh * qsat(i) ) - 1. ) * cldfra(i)
+            dqvc_sub(i) = qvapincld * dcf_sub(i)
+
+            cldfra(i) = cldfra(i) + dcf_sub(i)
+            qcld(i) = qcld(i) + dqvc_sub(i)
+            qvc(i) = qvc(i) + dqvc_sub(i)
+            clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
+            qclr(i) = qclr(i) - dqvc_sub(i)
+          ENDIF
 
           IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
@@ -571 +572 @@
 
           !--If the dissipation process must be activated
-          !IF ( cldfra(i) .GT. eps ) THEN
           IF ( qvapincld_new .GT. qvapincld ) THEN
             !--Gamma distribution starting at qvapincld
-            pdf_shape = ( mu_subl_pdf_lscp + 1. ) / qiceincld
+            pdf_shape = nu_iwc_pdf_lscp / qiceincld
             pdf_y = pdf_shape * ( qvapincld_new - qvapincld )
-            pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
-            pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
+            pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp      , pdf_y )
+            pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y )
 
             !--Tendencies and diagnostics
@@ -649 +649 @@
         !--Coefficient for standard deviation:
         !--  tuning coef * (clear sky area**0.25) * (function of temperature)
-        !pdf_e1 = beta_pdf_lscp &
-        !       * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) &
-        !       * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
-        !--  tuning coef * (cell area**0.25) * (function of temperature)
         pdf_e1 = beta_pdf_lscp * ( clrfra(i) * cell_area(i) )**0.25 &
                 * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
@@ -873 +869 @@
             pdf_q_above_lim = 0.
           ELSEIF ( pdf_y .LT. -10. ) THEN
-            pdf_fra_above_lim = 1.
-            pdf_q_above_lim = qclr(i) / clrfra(i)
+            pdf_fra_above_lim = clrfra(i)
+            pdf_q_above_lim = qclr(i)
           ELSE
             pdf_y = EXP( pdf_y )
@@ -1030 +1026 @@
             pdf_q_above_lim = 0.
           ELSEIF ( pdf_y .LT. -10. ) THEN
-            pdf_fra_above_lim = 1.
-            pdf_q_above_lim = qclr(i) / clrfra(i)
+            pdf_fra_above_lim = clrfra(i)
+            pdf_q_above_lim = qclr(i)
           ELSE
             pdf_y = EXP( pdf_y )
@@ -1103 +1099 @@
       !
       !--If ice supersaturation is activated, the cloud properties are prognostic.
-      !--The falling ice is then considered a new cloud in the gridbox.
+      !--The falling ice is then partly considered a new cloud in the gridbox.
       !--BEWARE with this parameterisation, we can create a new cloud with a much
       !--different ice water content and water vapor content than the existing cloud
@@ -1134 +1130 @@
             pdf_q_above_lim = 0.
           ELSEIF ( pdf_y .LT. -10. ) THEN
-            pdf_fra_above_lim = 1.
-            pdf_q_above_lim = qclr(i) / clrfra(i)
+            pdf_fra_above_lim = clrfra(i)
+            pdf_q_above_lim = qclr(i)
           ELSE
             pdf_y = EXP( pdf_y )
@@ -1159 +1155 @@
           dqi_sed(i) = qice_sedim(i)
 
-        ENDIF ! clrfra(i) .GT. eps
+        ENDIF ! ( clrfra_sed .GT. eps .AND. ( clrfra(i) .GT. eps )
 
         !--Add tendencies
@@ -1169 +1165 @@
         qclr(i)    = qclr(i) - dqvc_sed(i) + qice_sedim(i) - dqi_sed(i)
 
-      ENDIF ! qice_sedim(i) .GT. 0.
+      ENDIF ! icesed_flux(i) .GT. 0.
 
 
@@ -1188 +1184 @@
           qissr(i) = 0.
         ELSEIF ( pdf_y .LT. -10. ) THEN
-          issrfra(i) = 1.
-          qissr(i) = qclr(i) / clrfra(i)
+          issrfra(i) = clrfra(i)
+          qissr(i) = qclr(i)
         ELSE
           pdf_y = EXP( pdf_y )
@@ -1232 +1228 @@
       dqvc_sed(i) = dqvc_sed(i) / dtime
 
-    ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
+    ENDIF ! pt_pron_clds(i)
 
   ENDIF   ! end keepgoing
@@ -1320 +1316 @@
 
     ENDIF ! keepgoing
-  ENDDO
+  ENDDO ! loop on klon
 ENDIF ! ok_plane_contrail
 

LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90

@@ -165 +165 @@
   !$OMP THREADPRIVATE(ffallv_issr)
 
+  REAL, SAVE, PROTECTED :: incld_ice_thresh=1.E-3            ! [-] minimum in-cloud ice water content, relative to saturation
+  !$OMP THREADPRIVATE(incld_ice_thresh)
+
   REAL, SAVE, PROTECTED :: depo_coef_cirrus=.7               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
   !$OMP THREADPRIVATE(depo_coef_cirrus)
@@ -171 +174 @@
   !$OMP THREADPRIVATE(capa_cond_cirrus)
 
-  REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp=1./3.            ! [-] factor for the ice distribution inside cirrus clouds
-  !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
-  
-  REAL, SAVE, PROTECTED :: beta_pdf_lscp=8.75E-4             ! [SI] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
+  REAL, SAVE, PROTECTED :: N_ice_volume=3.E4                 ! [#/m3] ice crystal concentration in cirrus clouds (default value from  Kramer et al, 2020)
+  !$OMP THREADPRIVATE(N_ice_volume)
+
+  REAL, SAVE, PROTECTED :: nu_iwc_pdf_lscp=4./3.             ! [-] shape factor for the ice distribution inside cirrus clouds
+  !$OMP THREADPRIVATE(nu_iwc_pdf_lscp)
+
+  REAL, SAVE, PROTECTED :: corr_incld_depsub                 ! [-] correction factor for in-cloud IWC rather than local IWC in dep / sub process BEWARE NO GETIN (calculated automatically)
+  !$OMP THREADPRIVATE(corr_incld_depsub)
+  
+  REAL, SAVE, PROTECTED :: beta_pdf_lscp=8.46E-4             ! [SI] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
   !$OMP THREADPRIVATE(beta_pdf_lscp)
   
@@ -183 +192 @@
   !$OMP THREADPRIVATE(k0_pdf_lscp)
   
-  REAL, SAVE, PROTECTED :: kappa_pdf_lscp=0.0192             ! [] factor for the PDF fit of water vapor in UTLS
+  REAL, SAVE, PROTECTED :: kappa_pdf_lscp=0.0192             ! [K-1] factor for the PDF fit of water vapor in UTLS
   !$OMP THREADPRIVATE(kappa_pdf_lscp)
   
@@ -475 +484 @@
     ffallv_issr=ffallv_lsc
     CALL getin_p('ffallv_issr',ffallv_issr)
+    CALL getin_p('incld_ice_thresh',incld_ice_thresh)
     CALL getin_p('depo_coef_cirrus',depo_coef_cirrus)
     CALL getin_p('capa_cond_cirrus',capa_cond_cirrus)
-    CALL getin_p('mu_subl_pdf_lscp',mu_subl_pdf_lscp)
+    CALL getin_p('nu_iwc_pdf_lscp',nu_iwc_pdf_lscp)
     CALL getin_p('beta_pdf_lscp',beta_pdf_lscp)
     CALL getin_p('temp_thresh_pdf_lscp',temp_thresh_pdf_lscp)
@@ -573 +583 @@
     WRITE(lunout,*) 'lscp_ini, ok_weibull_warm_clouds:', ok_weibull_warm_clouds
     WRITE(lunout,*) 'lscp_ini, ffallv_issr', ffallv_issr
+    WRITE(lunout,*) 'lscp_ini, incld_ice_thresh', incld_ice_thresh
     WRITE(lunout,*) 'lscp_ini, depo_coef_cirrus:', depo_coef_cirrus
     WRITE(lunout,*) 'lscp_ini, capa_cond_cirrus:', capa_cond_cirrus
-    WRITE(lunout,*) 'lscp_ini, mu_subl_pdf_lscp:', mu_subl_pdf_lscp
+    WRITE(lunout,*) 'lscp_ini, nu_iwc_pdf_lscp:', nu_iwc_pdf_lscp
     WRITE(lunout,*) 'lscp_ini, beta_pdf_lscp:', beta_pdf_lscp
     WRITE(lunout,*) 'lscp_ini, temp_thresh_pdf_lscp:', temp_thresh_pdf_lscp
@@ -629 +640 @@
     ENDIF
 
+    !--Calculated here to lighten calculations
+    corr_incld_depsub = GAMMA(nu_iwc_pdf_lscp + 1./3.) / GAMMA(nu_iwc_pdf_lscp) &
+                      / nu_iwc_pdf_lscp**(1./3.)
+
 
     !AA Temporary initialisation