Changeset 5634

Timestamp:

Apr 29, 2025, 12:07:04 AM (6 weeks ago)

Author:

aborella

Message:

Small fixes condensation

Location:

LMDZ6/branches/contrails/libf/phylmd

Files:

• lmdz_aviation.f90 (modified) (2 diffs)
• lmdz_lscp_condensation.f90 (modified) (10 diffs)

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

@@ -174 +174 @@
       pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
       pdf_x = qcritcont(i) / qsatl(i) * 100.
-      pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+      pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
       IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
         pdf_fra_above_qcritcont = 0.
@@ -191 +191 @@
     
       pdf_x = qsat(i) / qsatl(i) * 100.
-      pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+      pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
       IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
         pdf_fra_above_qsat = 0.

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

@@ -556 +556 @@
               !--the old one, it is comprised between 0 and 1
               tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
-              qice_ratio = tauinv_depsub * dtime / 1.5 / qiceinmix * ( qsat(i) - qvapincld )
+              qice_ratio = tauinv_depsub * dtime / 1.5 / qiceincld * ( qsat(i) - qvapincld )
               !--The new vapor in the cloud is increased with the
               !--sublimated ice
@@ -683 +683 @@
         pdf_gamma(i) = GAMMA(1. + 1. / pdf_alpha(i))
         !--Barrier to avoid overflows
-        pdf_scale(i) = MAX(0.01, pdf_std / SQRT( &
-                                    GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 ))
+        pdf_scale(i) = MAX(eps, MIN(rhl_clr / pdf_gamma(i), pdf_std / SQRT( &
+                                    GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 )))
         pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
 
@@ -692 +692 @@
 
         pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
-        pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+        pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
         IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
           pdf_fra_above_nuc = 0.
@@ -747 +747 @@
         pdf_alpha(i) = 1.
         pdf_gamma(i) = 1.
-        pdf_scale(i) = 0.01
+        pdf_scale(i) = eps
       ENDIF ! clrfra(i) .GT. eps
 
@@ -786 +786 @@
         !--are spherical.
         !-- b / a = bovera = MAX(0.1, cldfra)
-        bovera = MAX(0.1, cldfra(i))
+        !bovera = MAX(0.1, cldfra(i))
+        bovera = 0.111111111
         !--P / a is a function of b / a only, that we can calculate
         !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
@@ -877 +878 @@
 
           rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+          pdf_x = qvapinclr_lim / qsatl(i) * 100.
           pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
-          pdf_x = qvapinclr_lim / qsatl(i) * 100.
-          pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+          pdf_x = qsat(i) / qsatl(i) * 100.
+          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
           IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
             pdf_fra_above_lim = 0.
@@ -1034 +1036 @@
 
           rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+          pdf_x = qvapinclr_lim / qsatl(i) * 100.
           pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
-          pdf_x = qvapinclr_lim / qsatl(i) * 100.
-          pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+          pdf_x = qsat(i) / qsatl(i) * 100.
+          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
           IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
             pdf_fra_above_lim = 0.
@@ -1138 +1141 @@
 
           rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
+          pdf_x = qvapinclr_lim / qsatl(i) * 100.
           pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
-          pdf_x = qvapinclr_lim / qsatl(i) * 100.
-          pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+          pdf_x = qsat(i) / qsatl(i) * 100.
+          pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
           IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
             pdf_fra_above_lim = 0.
@@ -1193 +1197 @@
         pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
         pdf_x = qsat(i) / qsatl(i) * 100.
-        pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i)
+        pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
         IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
           issrfra(i) = 0.
@@ -1206 +1210 @@
           issrfra(i) = EXP( - pdf_y ) * clrfra(i)
           qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
+        ENDIF
+        IF ( issrfra(i) .LT. eps ) THEN
+          issrfra(i) = 0.
+          qissr(i) = 0.
         ENDIF
       ELSE