Changeset 5625

Timestamp:

Apr 22, 2025, 6:03:13 PM (2 months ago)

Author:

aborella

Message:

Correction to coupling with deep convection + correction to cloud dissipation

Location:

LMDZ6/branches/contrails/libf/phylmd

Files:

• lmdz_lscp.f90 (modified) (2 diffs)
• lmdz_lscp_condensation.f90 (modified) (8 diffs)
• lmdz_lscp_ini.f90 (modified) (3 diffs)

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

@@ -1204 +1204 @@
       pcf_seri(:,k) = perscontfra(:)
       qva_seri(:,k) = zqs(:) * contfra(:)
-      qia_seri(:,k) = qcont(:) - zqs(:) * contfra(:)
       DO i = 1, klon
+        IF ( zoliqi(i) .GT. 0. ) THEN
+          !--The quantity of ice in linear contrails has been reduced by precipitation
+          !--with the same factor as the rest of the clouds
+          qia_seri(i,k) = ( qcont(i) - zqs(i) * contfra(i) ) / zoliqi(i) * radocond(i,k)
+        ELSE
+          qia_seri(i,k) = 0.
+        ENDIF
         IF ( ( rneb(i,k) - cfa_seri(i,k) ) .GT. eps ) THEN
+          !--The in-cloud quantity of ice in persistent contrail cirrus is the same as
+          !--the one from all natural cirrus clouds
           qice_perscont(i,k) = ( radocond(i,k) - qia_seri(i,k) ) &
               * perscontfra(i) / ( rneb(i,k) - cfa_seri(i,k) )
@@ -1219 +1227 @@
       DO i = 1, klon
 
+        !--We save the cloud properties that will be advected
         cf_seri(i,k) = rneb(i,k)
-
+        qvc_seri(i,k) = qvc(i)
+
+        !--Deep convection clouds properties are removed from radiative properties
+        !--outputed from lscp (NB. rneb and radocond are only used for the radiative
+        !--properties and are NOT prognostics)
+        !--We must have iflag_coupl == 5 for this coupling to work
+        IF ( ptconv(i,k) .AND. pt_pron_clds(i) ) THEN
+          rneb(i,k) = rneb(i,k) - cfcon(i,k)
+          radocond(i,k) = radocond(i,k) - qccon(i,k) * cfcon(i,k)
+        ENDIF
+
+        !--If everything was precipitated, the remaining empty cloud is dissipated
+        !--and everything is transfered to the subsaturated clear sky region
+        !--NB. we do not change rneb, as it is a diagnostic only
         IF ( zoliq(i) .LE. 0. ) THEN
-          !--If everything was precipitated, the remaining empty cloud is dissipated
-          !--and everything is transfered to the subsaturated clear sky region
-          !--NB. we do not change rneb, as it is a diagnostic only
           cf_seri(i,k) = 0.
+          qvc_seri(i,k) = 0.
           qvc(i) = 0.
         ENDIF
-
-        qvc_seri(i,k) = qvc(i)
 
         !--Diagnostics

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

@@ -127 +127 @@
 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, incld_ice_thresh
+USE lmdz_lscp_ini,   ONLY: coef_mixing_lscp, coef_shear_lscp, cooling_rate_ice_thresh
 
 USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, aspect_ratio_contrails
@@ -244 +244 @@
 !
 ! for dissipation
-REAL :: pdf_shape
-REAL :: dcf_sub_tmp, dqvc_sub_tmp
+REAL, DIMENSION(klon) :: temp_diss, qsati_diss
+REAL :: pdf_shape, qiceincld_min
 !
 ! for condensation
-REAL, DIMENSION(klon) :: qsatl, dqsatl
+REAL, DIMENSION(klon) :: qsatl, dqsat_tmp
 REAL, DIMENSION(klon) :: pdf_alpha, pdf_scale, pdf_gamma
 REAL :: rhl_clr, pdf_loc
@@ -278 +278 @@
 
 !--Calculation of qsat w.r.t. liquid
-CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsatl)
+CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 1, .FALSE., qsatl, dqsat_tmp)
+!--Calculation of qsat max for dissipation
+temp_diss(:) = temp(:) + cooling_rate_ice_thresh * dtime
+CALL calc_qsat_ecmwf(klon, temp_diss, qzero, pplay, RTT, 2, .FALSE., qsati_diss, dqsat_tmp)
 
 !
@@ -530 +533 @@
         ELSE
 
-          !--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_tmp = ( qiceincld / ( incld_ice_thresh * qsat(i) ) - 1. ) * cldfra(i)
-            dqvc_sub_tmp = qvapincld * dcf_sub_tmp
-
-            cldfra(i) = cldfra(i) + dcf_sub_tmp
-            qcld(i) = qcld(i) + dqvc_sub_tmp
-            qvc(i) = qvc(i) + dqvc_sub_tmp
-            clrfra(i) = MIN(1., clrfra(i) - dcf_sub_tmp)
-            qclr(i) = qclr(i) - dqvc_sub_tmp
-          ELSE
-            dcf_sub_tmp = 0.
-            dqvc_sub_tmp = 0.
-          ENDIF
-
           IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
             IF ( qvapincld .GE. qsat(i) ) THEN
@@ -574 +561 @@
           !--    DISSIPATION OF THE CLOUD    --
           !------------------------------------
+          !--Additionally to a minimum in cloud water vapor, we impose a minimum
+          !--in in-cloud ice water content. It is calculated following
+          !--Marti and Mauersberger (1993), see also Schiller et al. (2008)
+          qiceincld_min = qsati_diss(i) - qsat(i)
 
           !--If the dissipation process must be activated
-          IF ( qvapincld_new .GT. qvapincld ) THEN
+          IF ( ( qvapincld_new + qiceincld_min ) .GT. qvapincld ) THEN
             !--Gamma distribution starting at qvapincld
             pdf_shape = nu_iwc_pdf_lscp / qiceincld
-            pdf_y = pdf_shape * ( qvapincld_new - qvapincld )
+            pdf_y = pdf_shape * ( qvapincld_new + qiceincld_min - qvapincld )
             pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp      , pdf_y )
             pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y )
@@ -610 +601 @@
             qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
           ENDIF ! qvapincld_new .GT. qvapincld
-
-          !--We readd previously calculated dcf_sub_tmp and dqvc_sub_tmp to the diagnostics
-          dcf_sub(i) = dcf_sub(i) + dcf_sub_tmp
-          dqvc_sub(i) = dqvc_sub(i) + dqvc_sub_tmp
 
 
@@ -849 +836 @@
           qiceinmix = ( qcld(i) - qvc(i) ) / cldfra(i) / ( 1. + clrfra_mix / cldfra_mix )
           tauinv_depsub = qiceinmix**(1./3.) * kappa_depsub
-          qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime / 2. * tauinv_depsub ) )
+          qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) )
           qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) &
                         - qvc(i) / cldfra(i) * cldfra_mix / clrfra_mix
@@ -1126 +1113 @@
           IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
             tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
-            qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime / 2. * tauinv_depsub ) )
+            qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) )
           ELSE
             qvapinclr_lim = qsat(i) - qiceincld

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 :: cooling_rate_ice_thresh=2.e-5     ! [K/s] minimum virtual cooling rate before ice is sublimated
+  !$OMP THREADPRIVATE(cooling_rate_ice_thresh)
 
   REAL, SAVE, PROTECTED :: depo_coef_cirrus=.7               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
@@ -484 +484 @@
     ffallv_issr=ffallv_lsc
     CALL getin_p('ffallv_issr',ffallv_issr)
-    CALL getin_p('incld_ice_thresh',incld_ice_thresh)
+    CALL getin_p('cooling_rate_ice_thresh',cooling_rate_ice_thresh)
     CALL getin_p('depo_coef_cirrus',depo_coef_cirrus)
     CALL getin_p('capa_cond_cirrus',capa_cond_cirrus)
@@ -583 +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, cooling_rate_ice_thresh', cooling_rate_ice_thresh
     WRITE(lunout,*) 'lscp_ini, depo_coef_cirrus:', depo_coef_cirrus
     WRITE(lunout,*) 'lscp_ini, capa_cond_cirrus:', capa_cond_cirrus