Changeset 4833

Timestamp:

Feb 27, 2024, 5:48:56 PM (10 months ago)

Author:

evignon

Message:

commission poprecip liee a l'atelier nuage du jour:
traitement exact de l'evaporation et de la sublimation
validation de la nouvelle formulation du freezing de la pluie

File:

• LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90 (modified) (26 diffs)

LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90

@@ -63 +63 @@
 !--Integer for interating over klon
 INTEGER :: i
-!--hum_to_flux: coef to convert a specific quantity to a flux 
-REAL, DIMENSION(klon) :: hum_to_flux
+!--dhum_to_dflux: coef to convert a specific quantity to a flux 
+REAL, DIMENSION(klon) :: dhum_to_dflux
+!--
+REAL, DIMENSION(klon) :: rho, dz
 
 !--Saturation values
@@ -82 +84 @@
 dqssubl   = 0.
 
-!-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt
-hum_to_flux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime
+!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
+dhum_to_dflux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime
+rho(:) = pplay(:) / temp(:) / RD
+dz(:) = ( paprsdn(:) - paprsup(:) ) / RG / rho(:)
 
 !--Calculation of saturation specific humidity
@@ -117 +121 @@
     !--same temperature as the lowermost layer
     !--we convert the flux into a specific quantity qprecip
-    qprecip(i) = ( rain(i) + snow(i) ) / hum_to_flux(i)
+    qprecip(i) = ( rain(i) + snow(i) ) / dhum_to_dflux(i)
     !-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer
     temp(i) = ( tempupnew(i) * qprecip(i) * cpw + cpair * temp(i) ) &
@@ -133 +137 @@
     !--Evaporation of liquid precipitation coming from above
     !--in the clear sky only
-    !--dP/dz=beta*(1-q/qsat)*(P**expo_eva) (lines 1-2)
-    !--multiplying by dz = - dP / g / rho (line 3-4)
+    !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
     !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
-    draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) &
-             * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva &
-             * temp(i) * RD / pplay(i) &
-             * ( paprsdn(i) - paprsup(i) ) / RG
-
-    !--Evaporation is limited by 0 and by the total water amount in
-    !--the precipitation
-    draineva = MIN(0., MAX(draineva, -rainclr(i)))
+    !--Explicit formulation
+    !draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) * dz(i) &
+    !         * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva &
+    !draineva = MAX( - rainclr(i), draineva)
+    !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
+    !--which does not need a barrier on rainclr, because included in the formula
+    draineva = precipfracclr(i) * ( MAX(0., &
+             - coef_eva * ( 1. - expo_eva ) * (1. - qvap(i) / qsatl(i)) * dz(i) &
+             + ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_eva ) &
+               ) )**( 1. / ( 1. - expo_eva ) ) - rainclr(i)
+             
+    !--Evaporation is limited by 0
+    draineva = MIN(0., draineva)
 
 
     !--Sublimation of the solid precipitation coming from above
     !--(same formula as for liquid precip)
-    dsnowsub = - precipfracclr(i) * coef_sub * (1. - qvap(i) / qsati(i)) &
-             * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_sub &
-             * temp(i) * RD / pplay(i) &
-             * ( paprsdn(i) - paprsup(i) ) / RG
-
-    !--Sublimation is limited by 0 and by the total water amount in
-    !--the precipitation
+    !--Explicit formulation
+    !dsnowsub = - precipfracclr(i) * coef_sub * (1. - qvap(i) / qsatl(i)) * dz(i) &
+    !         * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_sub &
+    !dsnowsub = MAX( - snowclr(i), dsnowsub)
+    !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly)
+    !--which does not need a barrier on snowclr, because included in the formula
+    dsnowsub = precipfracclr(i) * ( MAX(0., &
+             - coef_sub * ( 1. - expo_sub ) * (1. - qvap(i) / qsatl(i)) * dz(i) &
+             + ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_sub ) &
+             ) )**( 1. / ( 1. - expo_sub ) ) - snowclr(i)
+
+    !--Sublimation is limited by 0
     ! TODO: change max when we will allow for vapor deposition in supersaturated regions
-    dsnowsub = MIN(0., MAX(dsnowsub, -snowclr(i)))
+    dsnowsub = MIN(0., dsnowsub)
 
     !--Evaporation limit: we ensure that the layer's fraction below
@@ -166 +179 @@
     
     dprecip_evasub_max = MIN(0., ( qvap(i) - qsat(i) ) * precipfracclr(i)) &
-                     * hum_to_flux(i)
+                     * dhum_to_dflux(i)
     dprecip_evasub_tot = draineva + dsnowsub
 
@@ -179 +192 @@
 
     !--New solid and liquid precipitation fluxes after evap and sublimation
-    dqrevap = draineva / hum_to_flux(i)
-    dqssubl = dsnowsub / hum_to_flux(i)
+    dqrevap = draineva / dhum_to_dflux(i)
+    dqssubl = dsnowsub / dhum_to_dflux(i)
 
 
@@ -196 +209 @@
 
     !--Add tendencies
-
+    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
     rainclr(i) = MAX(0., rainclr(i) + draineva)
     snowclr(i) = MAX(0., snowclr(i) + dsnowsub)
@@ -306 +319 @@
 
 INTEGER :: i
-REAL, DIMENSION(klon) :: hum_to_flux
+REAL, DIMENSION(klon) :: dhum_to_dflux
 REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip
 
@@ -371 +384 @@
   dqifreez= 0.
 
-  !--hum_to_flux: coef to convert a delta in specific quantity to a flux 
-  !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt
-  hum_to_flux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
+  !--dhum_to_dflux: coef to convert a delta in specific quantity to a flux 
+  !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
+  dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
   qtot(i) = qvap(i) + qliq(i) + qice(i) &
-          + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / hum_to_flux(i)
+          + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i)
 
   !------------------------------------------------------------
@@ -457 +470 @@
 
   !--We add the tendencies
+  !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
   precipfraccld(i) = precipfraccld(i) + dprecipfraccld
   precipfracclr(i) = precipfracclr(i) + dprecipfracclr
-
-  rainclr(i) = MAX( 0. , rainclr(i) + drainclr )
-  snowclr(i) = MAX( 0. , snowclr(i) + dsnowclr )
-  raincld(i) = MAX( 0. , raincld(i) - drainclr )
-  snowcld(i) = MAX( 0. , snowcld(i) - dsnowclr )
+  rainclr(i) = MAX(0., rainclr(i) + drainclr)
+  snowclr(i) = MAX(0., snowclr(i) + dsnowclr)
+  raincld(i) = MAX(0., raincld(i) - drainclr)
+  snowcld(i) = MAX(0., snowcld(i) - dsnowclr)
   
   !--If vertical heterogeneity is taken into account, we use
@@ -501 +514 @@
     !--then simplified.
 
-    !--The sticking efficacy is perfect.
+    !--The collection efficiency is perfect.
     Eff_rain_liq = 1.
     coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq
-    IF ((raincld(i) .GT. 0.) .AND. (coef_col .GT. 0.)) THEN
+    IF ( raincld(i) .GT. 0. ) THEN
       !-- ATTENTION Double implicit version
       !-- BEWARE the formule below is FALSE (because raincld is a flux, not a delta flux)
-      !qrain_tmp = raincld(i) / hum_to_flux(i)
+      !qrain_tmp = raincld(i) / dhum_to_dflux(i)
       !coef_tmp = coef_col * dtime * ( qrain_tmp / precipfraccld(i) + qliq(i) / eff_cldfra )
       !dqlcol = qliq(i) * ( 1. / ( 1. + 0.5 * ( coef_tmp - 1. + SQRT( &
@@ -521 +534 @@
       !--Add tendencies
       qliq(i) = qliq(i) + dqlcol
-      raincld(i) =  MAX( 0. , raincld(i) - dqlcol * hum_to_flux(i) )
+      raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i)
 
       !--Diagnostic tendencies
@@ -532 +545 @@
     Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
     coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
-    IF ((snowcld(i) .GT. 0.) .AND. (coef_agg .GT. 0.)) THEN
+    IF ( snowcld(i) .GT. 0. ) THEN
       !-- ATTENTION Double implicit version?
       !--Barriers so that the processes do not consume more liquid/ice than
@@ -543 +556 @@
       !--Add tendencies
       qice(i) = qice(i) + dqiagg
-      snowcld(i) = MAX( 0. , snowcld(i) - dqiagg * hum_to_flux(i) )
+      snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
 
       !--Diagnostic tendencies
@@ -607 +620 @@
     qliq(i) = qliq(i) + dqlauto
     qice(i) = qice(i) + dqiauto
-    raincld(i) = MAX( 0. , raincld(i) - dqlauto * hum_to_flux(i) )
-    snowcld(i) = MAX( 0. , snowcld(i) - dqiauto * hum_to_flux(i) )
+    raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
+    snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
 
     !--Diagnostic tendencies
@@ -628 +641 @@
     Eff_snow_liq = 0.2
     coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
-    IF ((snowcld(i) .GT. 0.) .AND. (coef_rim .GT. 0.)) THEN
+    IF ( snowcld(i) .GT. 0. ) THEN
       !-- ATTENTION Double implicit version?
       !--Barriers so that the processes do not consume more liquid than
@@ -638 +651 @@
 
       !--Add tendencies
+      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
       qliq(i) = qliq(i) + dqlrim
-      snowcld(i) = MAX( 0. , snowcld(i) - dqlrim * hum_to_flux(i) )
+      snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i)
 
       !--Temperature adjustment with the release of latent
@@ -706 +720 @@
 
       !--Barrier to limit the melting flux to the clr snow flux in the mesh
-      dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / hum_to_flux(i))
+      dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
     ENDIF
 
@@ -720 +734 @@
 
       !--Barrier to limit the melting flux to the cld snow flux in the mesh
-      dqscldmelt = MAX( dqscldmelt , -snowcld(i) / hum_to_flux(i))
+      dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
     ENDIF
     
@@ -738 +752 @@
 
     !--Add tendencies
-
-    rainclr(i) = MAX( 0. , rainclr(i) - dqsclrmelt * hum_to_flux(i) )
-    raincld(i) = MAX( 0. , raincld(i) - dqscldmelt * hum_to_flux(i) )
-    snowclr(i) = MAX( 0. , snowclr(i) + dqsclrmelt * hum_to_flux(i) )
-    snowcld(i) = MAX( 0. , snowcld(i) + dqscldmelt * hum_to_flux(i) )
+    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
+    rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i))
+    raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i))
+    snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i))
+    snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i))
 
 
@@ -792 +806 @@
     dqrtotfreez_step1 = 0.
 
-    IF ((qice(i) .GT. 0.) .AND. (raincld(i) .GT. 0.)) THEN
+    IF ( ( qice(i) .GT. 0. ) .AND. ( raincld(i) .GT. 0. ) ) THEN
       dqrclrfreez = 0.
       dqrcldfreez = 0.
@@ -815 +829 @@
 
       !--We add the part of rain water that freezes, limited by a temperature barrier
-      !-- this quantity is calculated assuming that the number of drop that freeze correspond to the number
-      !-- of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
-      dqrcldfreez = MAX(dqifreez*rho_rain*r_rain/(rho_ice*r_ice),-raincld(i)/hum_to_flux(i))
+      !--this quantity is calculated assuming that the number of drop that freeze correspond to the number
+      !--of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
+      !--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3
+      !--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3
+      !--The assumption above corresponds to dNi = dNr, i.e.,
+      !-- dqr = dqi * (4/3 PI rho_rain * r_rain**3) / (4/3 PI rho_ice * r_ice**3)
+      dqrcldfreez = dqifreez * rho_rain * r_rain**3. / ( rho_ice * r_ice**3. )
+      dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i))
       dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max) 
       dqrtotfreez_step1 = dqrcldfreez
 
       !--Add tendencies
+      !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
       qice(i) = qice(i) + dqifreez
-      raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i) )
-      snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) - dqifreez * hum_to_flux(i) )
+      raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
+      snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i))
       temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
-                      / ( 1. + RVTMP2 * qtot(i) )
+                        / ( 1. + RVTMP2 * qtot(i) )
 
     ENDIF
@@ -850 +870 @@
     IF ( rainclr(i) .GT. 0. ) THEN
       !--Exact solution of dqrain/dt = -qrain/tau_freez
-      dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
+      dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
     ENDIF
 
@@ -856 +876 @@
     IF ( raincld(i) .GT. 0. ) THEN
       !--Exact solution of dqrain/dt = -qrain/tau_freez
-      dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
+      dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
     ENDIF
 
@@ -874 +894 @@
 
     !--Add tendencies
-    rainclr(i) = MAX( 0. , rainclr(i) + dqrclrfreez * hum_to_flux(i) )
-    raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i) )
-    snowclr(i) = MAX( 0. , snowclr(i) - dqrclrfreez * hum_to_flux(i) )
-    snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) )
+    !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
+    rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i))
+    raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i))
+    snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i))
+    snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i))
 
 
@@ -885 +906 @@
                       / ( 1. + RVTMP2 * qtot(i) )
 
-    dqrtotfreez=dqrtotfreez_step1+dqrtotfreez_step2          
     !--Diagnostic tendencies
+    dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2          
     dqrfreez(i) = dqrtotfreez / dtime
     dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime