Changeset 5495

Timestamp:

Jan 20, 2025, 4:01:44 PM (7 hours ago)

Author:

evignon

Message:

correction du diagnostique de l'humidite relative en ciel clair.
n'affecte que le coating des aerosols.
Audran et Etienne

File:

• LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90 (modified) (2 diffs)

LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90

@@ -841 +841 @@
             IF (mpc_bl_points(i,k) .GT. 0) THEN
                 zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
-                ! following line is very strange and probably wrong
-                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
+                !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
+                !--This is slighly approximated, the actual formula is
+                !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
+                !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
+                !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
+                rhcl(i,k) = ( zq(i) - qincloud_mpc(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
                 ! water vapor update and partition function if thermals
                 zq(i) = zq(i) - zcond(i)        
@@ -872 +876 @@
                       zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
                     ENDIF
-                    ! following line is very strange and probably wrong:
-                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
+                    !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
+                    !--This is slighly approximated, the actual formula is
+                    !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
+                    !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
+                    !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
+                    rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
                     ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param) 
                     IF (iflag_icefrac .GE. 1) THEN