Changeset 5716

Timestamp:

Jun 18, 2025, 11:46:38 AM (5 weeks ago)

Author:

aborella

Message:

Minor bug fixes (temperature and water conservation)

Location:

LMDZ6/branches/contrails/libf/phylmd

Files:

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

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

@@ -968 +968 @@
                         ENDIF
 
+                        IF ( ok_ice_sedim ) THEN
+                          qice_sedim = flsed(i) / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
+                          ! Add the ice that was sedimented, as it is not included in zqn
+                          qlbef(i) = qlbef(i) + qice_sedim
+                        ENDIF
+
                         num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
                               +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
@@ -1086 +1092 @@
             ENDIF
 
-            ! Remove the ice that was sedimentated. As it is not included in zqn,
-            ! we only remove it from the total water
             IF ( ok_ice_sedim ) THEN
-              zq(i) = zq(i) - flsed(i) / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
+              qice_sedim = flsed(i) / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
+              ! Remove the ice that was sedimented. As it is not included in zqn,
+              ! we only remove it from the total water
+              zq(i) = zq(i) - qice_sedim
+              ! Temperature update due to phase change (sedimented ice was condensed)
+              zt(i) = zt(i) + qice_sedim &
+                    * RLSTT / RCPD / ( 1. + RVTMP2 * ( zq(i) + zmqc(i) + zcond(i) ) )
             ENDIF
-            
-            
+
             ! temperature update due to phase change
             zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & 

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

@@ -133 +133 @@
 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, chi_mixing
 USE lmdz_lscp_ini,   ONLY: aspect_ratio_cirrus, cooling_rate_ice_thresh
-USE lmdz_lscp_ini,   ONLY: ok_ice_sedim, fallice_sedim
+USE lmdz_lscp_ini,   ONLY: ok_ice_sedim, fallice_sedim, cice_velo, dice_velo
+USE lmdz_lscp_ini,   ONLY: chi_sedim
 
 USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, aspect_ratio_lincontrails
@@ -289 +290 @@
 ! for sedimentation
 REAL, DIMENSION(klon) :: cfsed_abv, qised_abv
-REAL :: qice_sedim
+REAL :: iwc, icefall_velo, qice_sedim
 REAL :: sedfra_abv, sedfra1, sedfra2, sedfra3
 !
@@ -732 +733 @@
         !--New PDF
         rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
-        rhl_clr = MAX(0., MIN(150., rhl_clr))
+        !--rhl_clr is limited to 80%, because the fit is not valid above.
+        !--This corresponds to rhi_clr ~= 140% at 210K, and ~= 105% at 245K
+        rhl_clr = MAX(0., MIN(80., rhl_clr))
 
         !--Calculation of the properties of the PDF
@@ -740 +743 @@
         !--Coefficient for standard deviation:
         !--  tuning coef * (clear sky area**0.25) * (function of temperature)
-        pdf_e1 = beta_pdf_lscp * ( clrfra(i) * cell_area(i) )**0.25 &
-                * MAX( MAX(205., MIN(250., temp(i))) - temp_thresh_pdf_lscp, 0. )
+        !--MAX on the sqrt of the area because the parameterization was derived for
+        !-- L > 50 km
+        pdf_e1 = beta_pdf_lscp * SQRT(MAX(5e4, SQRT( clrfra(i) * cell_area(i) ) )) &
+                * MAX( MAX(210., MIN(245., temp(i))) - temp_thresh_pdf_lscp, 0. )
         IF ( rhl_clr .GT. 50. ) THEN
           pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp
@@ -748 +753 @@
         ENDIF
         pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * &
-                    MAX( temp_nowater - MAX(205., MIN(250., temp(i))), 0. )
+                    MAX( temp_nowater - MAX(210., MIN(245., temp(i))), 0. )
         pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3
         pdf_alpha(i) = MIN(10., pdf_alpha(i)) !--Avoid overflows
@@ -972 +977 @@
           !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing
           !--decreases the size of the clouds
-          sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim )
+          !sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim )
+          !--The chance that the air surrounding cirrus clouds is moist is higher then
+          !--the actual distribution of humidity in the gridbox. 
+          !--This is quantified with chi_mixing. For chi_mixing=1, the two are equal.
+          !--If chi_mixing > 1, the chance that moist air is around cirrus clouds is increased
+          sigma_mix = MAX(1. - pdf_fra_below_lim / clrfra_mix, &
+                      MIN(pdf_fra_above_lim / clrfra_mix, &
+                  chi_mixing * pdf_fra_above_lim &
+                  / ( pdf_fra_below_lim + chi_mixing * pdf_fra_above_lim )))
 
           IF ( pdf_fra_above_lim .GT. eps ) THEN
-            dcf_mix(i)  = clrfra_mix * sigma_mix
-            dqvc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
+            IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
+              dcf_mix(i)  = clrfra_mix * sigma_mix
+              dqvc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
+            ELSE
+              dcf_mix(i)  = clrfra_mix * sigma_mix
+              dqvc_mix(i) = clrfra_mix * sigma_mix * qsat(i)
+              dqi_mix(i)  = clrfra_mix * sigma_mix * (pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
+            ENDIF
           ENDIF
 
@@ -1130 +1149 @@
           !--For aviation, we increase the chance that the air surrounding contrails
           !--is moist. This is quantified with chi_mixing_lincontrails
-          sigma_mix = chi_mixing_lincontrails * pdf_fra_above_lim &
-                    / ( pdf_fra_below_lim + chi_mixing_lincontrails * pdf_fra_above_lim )
+          sigma_mix = MAX(1. - pdf_fra_below_lim / clrfra_mix, &
+                      MIN(pdf_fra_above_lim / clrfra_mix, &
+                  chi_mixing_lincontrails * pdf_fra_above_lim &
+                  / ( pdf_fra_below_lim + chi_mixing_lincontrails * pdf_fra_above_lim )))
 
           IF ( pdf_fra_above_lim .GT. eps ) THEN
@@ -1290 +1311 @@
           !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing
           !--decreases the size of the clouds
-          sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim )
+          sigma_mix = MAX(1. - pdf_fra_below_lim / clrfra_mix, &
+                      MIN(pdf_fra_above_lim / clrfra_mix, &
+                  chi_mixing * pdf_fra_above_lim &
+                  / ( pdf_fra_below_lim + chi_mixing * pdf_fra_above_lim )))
 
           IF ( pdf_fra_above_lim .GT. eps ) THEN
@@ -1362 +1386 @@
         !--in-cloud ice water content is constant.
         !
-        qice_sedim = 0.
         IF ( cldfra(i) .GT. eps ) THEN
-          dzsed(i) = MIN(fallice_sedim * dtime, dz)
+          iwc = rho * ( qcld(i) - qvc(i) ) / cldfra(i)
+          icefall_velo = fallice_sedim * cice_velo * iwc**dice_velo
+          dzsed(i) = icefall_velo * dtime * MAX(0., MIN(1., (2. - icefall_velo * dtime / dz)))
           cfsed(i) = cldfra(i)
           qice_sedim = ( qcld(i) - qvc(i) ) * dzsed(i) / dz
@@ -1433 +1458 @@
                                 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
             ENDIF
+            
+            sedfra1 = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, &
+                    sedfra1 * chi_sedim * pdf_fra_above_lim &
+                    / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim )))
 
             IF ( pdf_fra_above_lim .GT. eps ) THEN
-              sedfra1 = sedfra1 * pdf_fra_above_lim / clrfra(i)
-              dcf_sed(i)  = dcf_sed(i)  + sedfra1
-              dqvc_sed(i) = dqvc_sed(i) + sedfra1 * pdf_q_above_lim / pdf_fra_above_lim
-              dqi_sed(i)  = dqi_sed(i)  + sedfra1 * qiceincld
+              IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
+                dcf_sed(i)  = dcf_sed(i)  + sedfra1
+                dqvc_sed(i) = dqvc_sed(i) + sedfra1 * pdf_q_above_lim / pdf_fra_above_lim
+                dqi_sed(i)  = dqi_sed(i)  + sedfra1 * qiceincld
+              ELSE
+                dcf_sed(i)  = dcf_sed(i)  + sedfra1
+                dqvc_sed(i) = dqvc_sed(i) + sedfra1 * qsat(i)
+                dqi_sed(i)  = dqi_sed(i)  + sedfra1 &
+                                 * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
+              ENDIF
             ENDIF
           ENDIF ! ( sedfra1 .GT. eps .AND. ( clrfra(i) .GT. eps )
@@ -1446 +1481 @@
 
           !--For region (3)
-          IF ( sedfra3 .GT. eps ) THEN
+          IF ( cldfra(i) .GT. eps ) THEN
             dcf_sed(i)  = dcf_sed(i)  + sedfra3
             dqvc_sed(i) = dqvc_sed(i) + sedfra3 * qvc(i) / cldfra(i)
@@ -1452 +1487 @@
           ENDIF
         ENDIF ! qised_abv(i) .GT. 0.
-        !PRINT *, 'A', cldfra(i), dcf_sed(i), clrfra(i)
-        !PRINT *, 'B', qcld(i) - qvc(i), qvc(i), dqvc_sed(i), dqi_sed(i)
 
         !--Add tendencies

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

@@ -221 +221 @@
   REAL, SAVE, PROTECTED :: coef_shear_lscp=0.72              ! [-] additional coefficient for the shearing process (subprocess of the mixing process)
   !$OMP THREADPRIVATE(coef_shear_lscp)
+  
+  REAL, SAVE, PROTECTED :: chi_mixing=1.5                    ! [-] factor for increasing the chance that moist air is surrounding  cirrus clouds
+  !$OMP THREADPRIVATE(chi_mixing)
   !--End of the parameters for condensation and ice supersaturation
 
@@ -373 +376 @@
   REAL, SAVE, PROTECTED :: fallice_sedim=1.                 ! Tuning factor for ice fallspeed velocity for sedimentation [-]
   !$OMP THREADPRIVATE(fallice_sedim)
+  
+  REAL, SAVE, PROTECTED :: chi_sedim=1E5                    ! [-] factor for increasing the chance that sedimented ice falls into moist air
+  !$OMP THREADPRIVATE(chi_sedim)
   !--End of the parameters for poprecip
 
@@ -494 +500 @@
     CALL getin_p('ok_ice_sedim',ok_ice_sedim)
     CALL getin_p('fallice_sedim',fallice_sedim)
+    CALL getin_p('chi_sedim',chi_sedim)
     ! for condensation and ice supersaturation
     CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds)
@@ -517 +524 @@
     CALL getin_p('coef_mixing_lscp',coef_mixing_lscp)
     CALL getin_p('coef_shear_lscp',coef_shear_lscp)
+    CALL getin_p('chi_mixing',chi_mixing)
     ! for aviation
     CALL getin_p('ok_precip_contrails',ok_precip_contrails)
@@ -599 +607 @@
     WRITE(lunout,*) 'lscp_ini, ok_ice_sedim:', ok_ice_sedim
     WRITE(lunout,*) 'lscp_ini, fallice_sedim:', fallice_sedim
+    WRITE(lunout,*) 'lscp_ini, chi_sedim:', chi_sedim
     ! for condensation and ice supersaturation
     WRITE(lunout,*) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat
@@ -622 +631 @@
     WRITE(lunout,*) 'lscp_ini, coef_mixing_lscp:', coef_mixing_lscp
     WRITE(lunout,*) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp
+    WRITE(lunout,*) 'lscp_ini, chi_mixing:', chi_mixing
     ! for aviation
     WRITE(lunout,*) 'lscp_ini, ok_precip_contrails:', ok_precip_contrails