Changeset 5430 for LMDZ6/trunk/libf

Timestamp:

Dec 19, 2024, 3:41:51 PM (16 hours ago)

Author:

evignon

Message:

ajout de la deposition de vapeur sur les precipitations
glacees (croissance) dans poprecip, A borella

Location:

LMDZ6/trunk/libf/phylmd

Files:

• lmdz_lscp.f90 (modified) (1 diff)
• lmdz_lscp_precip.f90 (modified) (12 diffs)

LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90

@@ -494 +494 @@
                         zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
                         zqvapclr, zqupnew, &
+                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
                         zrfl, zrflclr, zrflcld, &
                         zifl, ziflclr, ziflcld, &

LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90

@@ -714 +714 @@
            klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
            qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
+           cldfra, rvc_seri, qliq, qice, &
            rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
            )
@@ -720 +721 @@
 USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
-USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub
+USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds
+USE lmdz_lscp_ini, ONLY : eps, temp_nowater
 USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
 
@@ -747 +749 @@
 REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
 
+REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
+REAL,    INTENT(IN),    DIMENSION(klon) :: rvc_seri       !--cloud water vapor at the beginning of lscp (ratio wrt total water) - used only if the cloud properties are advected [kg/kg]
+REAL,    INTENT(IN),    DIMENSION(klon) :: qliq           !--liquid water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
+REAL,    INTENT(IN),    DIMENSION(klon) :: qice           !--ice water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
+
+
 REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
 REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
-REAL,    INTENT(IN),    DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
 REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
 REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
-REAL,    INTENT(IN),    DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
 
 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
@@ -762 +770 @@
 !--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
 REAL, DIMENSION(klon) :: dhum_to_dflux
-!--
+!--Air density [kg/m3] and layer thickness [m]
 REAL, DIMENSION(klon) :: rho, dz
+!--Temporary precip fractions and fluxes for the evaporation
+REAL, DIMENSION(klon) :: precipfracclr_tmp, precipfraccld_tmp
+REAL, DIMENSION(klon) :: rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp
 
 !--Saturation values
 REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
-!--Vapor in the clear sky
-REAL :: qvapclr
+!--Vapor in the clear sky and cloud
+REAL :: qvapclr, qvapcld
 !--Fluxes tendencies because of evaporation and sublimation
-REAL :: dprecip_evasub_max, draineva, dsnowsub, dprecip_evasub_tot
+REAL :: dprecip_evasub_max, dprecip_evasub_tot
+REAL :: drainclreva, draincldeva, dsnowclrsub, dsnowcldsub
 !--Specific humidity tendencies because of evaporation and sublimation
 REAL :: dqrevap, dqssubl
@@ -828 +840 @@
 ENDIF
 
-! TODO Probleme : on utilise qvap total dans la maille pour l'evap / sub
-! alors qu'on n'evap / sub que dans le ciel clair
-! deux options pour cette routine :
-! - soit on diagnostique le nuage AVANT l'evap / sub et on estime donc
-! la fraction precipitante ciel clair dans la maille, ce qui permet de travailler
-! avec des fractions, des fluxs et surtout un qvap dans le ciel clair
-! - soit on pousse la param de Ludo au bout, et on prend un qvap de k+1
-! dans le ciel clair, avec un truc comme :
-!   qvapclr(k) = qvapclr(k+1)/qtot(k+1) * qtot(k)
-! UPDATE : on code la seconde version. A voir si on veut mettre la premiere version.
-
-
+
+!--Initialise the precipitation fractions and fluxes
 DO i = 1, klon
+  precipfracclr_tmp(i) = precipfracclr(i)
+  precipfraccld_tmp(i) = precipfraccld(i)
+  rainclr_tmp(i) = rainclr(i)
+  snowclr_tmp(i) = snowclr(i)
+  raincld_tmp(i) = raincld(i)
+  snowcld_tmp(i) = snowcld(i)
+ENDDO
+
+!--If we relax the LTP assumption that the cloud is the same than in the
+!--gridbox above, we can change the tmp quantities
+IF ( ok_ice_supersat ) THEN
+  !--Update the precipitation fraction using advected cloud fraction
+  CALL poprecip_fracupdate( &
+             klon, cldfra, precipfracclr_tmp, precipfraccld_tmp, &
+             rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp)
+
+! here we could put an ELSE statement and do one iteration of the condensation scheme
+! (but it would only diagnose cloud from lognormal - link with thermals?)
+
+  DO i = 1, klon
+    !--We cannot have a total fraction greater than the one in the gridbox above
+    !--because new cloud formation has not yet occured.
+    !--If this happens, we reduce the precip frac in the cloud, because it can
+    !--only mean that the cloud fraction at this level is greater than the total
+    !--precipitation fraction of the level above, and the precip frac in clear sky
+    !--is zero.
+    !--NB. this only works because we assume a max-random cloud and precip overlap
+    precipfraccld_tmp(i) = MIN( precipfraccld_tmp(i), precipfracclr(i) + precipfraccld(i) )
+  ENDDO
+
+ENDIF
+!--At this stage, we guarantee that
+!-- precipfracclr + precipfraccld = precipfracclr_tmp + precipfraccld_tmp
+
+
+DO i = 1, klon
 
   !--If there is precipitation from the layer above
-  ! NOTE TODO here we could replace the condition on precipfracclr(i) by a condition
-  ! such as eps or thresh_precip_frac, to remove the senseless barrier in the formulas
-  ! of evap / sublim
-  IF ( ( ( rain(i) + snow(i) ) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN
-
-    IF ( ok_corr_vap_evasub ) THEN
-      !--Corrected version - we use the same water ratio between
+  IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
+
+    !--Init
+    drainclreva = 0.
+    draincldeva = 0.
+    dsnowclrsub = 0.
+    dsnowcldsub = 0.
+
+    !--Init the properties of the air where reevaporation occurs
+    IF ( ok_ice_supersat ) THEN
+      !--We can diagnose the water vapor in the cloud using the advected
+      !--water vapor in the cloud and cloud fraction
+      IF ( ok_unadjusted_clouds .AND. ( cldfra(i) .GT. eps ) ) THEN
+        qvapcld = rvc_seri(i) * qvap(i) / cldfra(i)
+      ENDIF
+      !--We can diagnose completely the water vapor in clear sky, because all
+      !--the needed variables (ice, liq, vapor in cld, cloud fraction) are
+      !--advected
+      !--Note that qvap(i) is the total water in the gridbox
+      IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
+        qvapclr = ( qvap(i) - qice(i) - qliq(i) - rvc_seri(i) * qvap(i) ) / ( 1. - cldfra(i) )
+      ENDIF
+    ELSEIF ( ok_corr_vap_evasub ) THEN
+      !--Corrected version from Ludo - we use the same water ratio between
       !--the clear and the cloudy sky as in the layer above. This
       !--extends the assumption that the cloud fraction is the same
@@ -856 +911 @@
       !--Note that qvap(i) is the total water in the gridbox, and
       !--precipfraccld(i) is the cloud fraction in the layer above
-      qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
+      IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN
+        qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
+      ENDIF
     ELSE
       !--Legacy version from Ludo - we use the total specific humidity
@@ -863 +920 @@
     ENDIF
 
-    !--Evaporation of liquid precipitation coming from above
-    !--in the clear sky only
-    !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
-    !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
-    !--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. - qvapclr / 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)
-    !--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. - qvapclr / qsati(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., dsnowsub)
-
-    !--Evaporation limit: we ensure that the layer's fraction below
-    !--the clear sky does not reach saturation. In this case, we 
-    !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice 
-    !--Max evaporation is computed not to saturate the clear sky precip fraction
-    !--(i.e., the fraction where evaporation occurs)
-    !--It is expressed as a max flux dprecip_evasub_max
-    
-    dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr(i)) &
-                     * dhum_to_dflux(i)
-    dprecip_evasub_tot = draineva + dsnowsub
-
-    !--Barriers
-    !--If activates if the total is LOWER than the max because
-    !--everything is negative
-    IF ( dprecip_evasub_tot .LT. dprecip_evasub_max ) THEN
-      draineva = dprecip_evasub_max * draineva / dprecip_evasub_tot
-      dsnowsub = dprecip_evasub_max * dsnowsub / dprecip_evasub_tot
-    ENDIF
+    !---------------------------
+    !-- EVAP / SUBL IN CLEAR SKY
+    !---------------------------
+
+    IF ( precipfracclr_tmp(i) .GT. eps ) THEN
+
+      !--Evaporation of liquid precipitation coming from above
+      !--in the clear sky only
+      !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
+      !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
+      !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
+      !--which does not need a barrier on rainclr, because included in the formula
+      drainclreva = precipfracclr_tmp(i) * MAX(0., &
+                  - coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) &
+                  + ( rainclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_eva ) &
+                  )**( 1. / ( 1. - expo_eva ) ) - rainclr_tmp(i)
+               
+      !--Evaporation is limited by 0
+      !--NB. with ok_ice_supersat activated, this barrier should be useless
+      drainclreva = MIN(0., drainclreva)
+
+
+      !--Sublimation of the solid precipitation coming from above
+      !--(same formula as for liquid precip)
+      !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly)
+      !--which does not need a barrier on snowclr, because included in the formula
+      dsnowclrsub = precipfracclr_tmp(i) * MAX(0., &
+                  - coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) &
+                  + ( snowclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_sub ) &
+                  )**( 1. / ( 1. - expo_sub ) ) - snowclr_tmp(i)
+
+      !--If ice supersaturation is activated, we allow vapor to condense on falling snow
+      !--i.e., having a positive dsnowclrsub
+      IF ( .NOT. ok_ice_supersat ) THEN
+        !--Sublimation is limited by 0
+        dsnowclrsub = MIN(0., dsnowclrsub)
+      ENDIF
+
+
+      !--Evaporation limit: we ensure that the layer's fraction below
+      !--the clear sky does not reach saturation. In this case, we 
+      !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice 
+      !--Max evaporation is computed not to saturate the clear sky precip fraction
+      !--(i.e., the fraction where evaporation occurs)
+      !--It is expressed as a max flux dprecip_evasub_max
+      
+      IF ( .NOT. ok_ice_supersat ) THEN
+        dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr_tmp(i)) &
+                         * dhum_to_dflux(i)
+      ELSE
+        dprecip_evasub_max = ( qvapclr - qsat(i) ) * precipfracclr_tmp(i) &
+                         * dhum_to_dflux(i)
+      ENDIF
+      dprecip_evasub_tot = drainclreva + dsnowclrsub
+
+      !--Barriers
+      !--If activates if the total is LOWER than the max because
+      !--everything is negative
+      IF ( (( dprecip_evasub_max .LT. 0. ) .AND. &
+            ( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. &
+           (( dprecip_evasub_max .GT. 0. ) .AND. &
+            ( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN
+        drainclreva = dprecip_evasub_max * drainclreva / dprecip_evasub_tot
+        dsnowclrsub = dprecip_evasub_max * dsnowclrsub / dprecip_evasub_tot
+      ENDIF
+
+    ENDIF ! precipfracclr_tmp .GT. eps
+
+
+    !---------------------------
+    !-- EVAP / SUBL IN THE CLOUD
+    !---------------------------
+
+    IF ( ok_unadjusted_clouds .AND. ( temp(i) .LE. temp_nowater ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
+      !--Evaporation of liquid precipitation coming from above
+      !--in the cloud only
+      !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
+      !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
+      !--Exact explicit formulation (raincld is resolved exactly, qvap explicitly)
+      !--which does not need a barrier on raincld, because included in the formula
+      !draincldeva = precipfraccld_tmp(i) * MAX(0., &
+      !            - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
+      !            + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
+      !            )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
+               
+      !--Evaporation is limited by 0
+      !--NB. with ok_ice_supersat activated, this barrier should be useless
+      !draincldeva = MIN(0., draincldeva)
+      draincldeva = 0.
+
+
+      !--Sublimation of the solid precipitation coming from above
+      !--(same formula as for liquid precip)
+      !--Exact explicit formulation (snowcld is resolved exactly, qvap explicitly)
+      !--which does not need a barrier on snowcld, because included in the formula
+      dsnowcldsub = precipfraccld_tmp(i) * MAX(0., &
+                  - coef_sub * ( 1. - expo_sub ) * (1. - qvapcld / qsati(i)) * dz(i) &
+                  + ( snowcld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_sub ) &
+                  )**( 1. / ( 1. - expo_sub ) ) - snowcld_tmp(i)
+
+      !--There is no barrier because we want deposition to occur if it is possible
+
+
+      !--Evaporation limit: we ensure that the layer's fraction below
+      !--the clear sky does not reach saturation. In this case, we 
+      !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice 
+      !--Max evaporation is computed not to saturate the clear sky precip fraction
+      !--(i.e., the fraction where evaporation occurs)
+      !--It is expressed as a max flux dprecip_evasub_max
+      
+      dprecip_evasub_max = ( qvapcld - qsat(i) ) * precipfraccld_tmp(i) * dhum_to_dflux(i)
+      dprecip_evasub_tot = draincldeva + dsnowcldsub
+
+      !--Barriers
+      !--If activates if the total is LOWER than the max because
+      !--everything is negative
+      IF ( (( dprecip_evasub_max .LT. 0. ) .AND. &
+            ( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. &
+           (( dprecip_evasub_max .GT. 0. ) .AND. &
+            ( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN
+        draincldeva = dprecip_evasub_max * draincldeva / dprecip_evasub_tot
+        dsnowcldsub = dprecip_evasub_max * dsnowcldsub / dprecip_evasub_tot
+      ENDIF
+
+    ENDIF ! ok_unadjusted_clouds
 
 
     !--New solid and liquid precipitation fluxes after evap and sublimation
-    dqrevap = draineva / dhum_to_dflux(i)
-    dqssubl = dsnowsub / dhum_to_dflux(i)
-
+    dqrevap = ( drainclreva + draincldeva ) / dhum_to_dflux(i)
+    dqssubl = ( dsnowclrsub + dsnowcldsub ) / dhum_to_dflux(i)
 
     !--Vapor is updated after evaporation/sublimation (it is increased)
@@ -929 +1064 @@
 
     !--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)
-
-    !--If there is no more precip fluxes, the precipitation fraction in clear
-    !--sky is set to 0
+    !--The MAX is needed because in some cases, the flux can be slightly
+    !--negative (numerical precision)
+    rainclr_tmp(i) = MAX(0., rainclr_tmp(i) + drainclreva)
+    snowclr_tmp(i) = MAX(0., snowclr_tmp(i) + dsnowclrsub)
+    raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva)
+    snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub)
+
+
+    !--We reattribute fluxes according to initial fractions, using proportionnality
+    !--Note that trough this process, we conserve total precip fluxes
+    !-- rainclr_tmp + raincld_tmp = rainclr + raincld (same for snow)
+    !--If the cloud has increased, a part of the precip in clear sky falls in clear sky,
+    !--the rest falls in the cloud
+    IF ( ( precipfraccld(i) .GT. eps ) .AND. ( precipfraccld_tmp(i) .GT. precipfraccld(i) ) ) THEN
+      !--All the precip from cloud falls in the cloud
+      raincld(i) = MAX(0., raincld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i))
+      rainclr(i) = MAX(0., rainclr_tmp(i) + raincld_tmp(i) - raincld(i))
+      snowcld(i) = MAX(0., snowcld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i))
+      snowclr(i) = MAX(0., snowclr_tmp(i) + snowcld_tmp(i) - snowcld(i))
+
+    !--If the cloud has narrowed, a part of the precip in cloud falls in clear sky,
+    !--the rest falls in the cloud
+    ELSEIF ( ( precipfracclr(i) .GT. eps ) .AND. ( precipfracclr_tmp(i) .GT. precipfracclr(i) ) ) THEN
+      !--All the precip from clear sky falls in clear sky
+      rainclr(i) = MAX(0., rainclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i))
+      raincld(i) = MAX(0., raincld_tmp(i) + rainclr_tmp(i) - rainclr(i))
+      snowclr(i) = MAX(0., snowclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i))
+      snowcld(i) = MAX(0., snowcld_tmp(i) + snowclr_tmp(i) - snowclr(i))
+
+    ELSE
+      !--If the cloud stays the same or if there is no cloud above and
+      !--in the current layer, there is no reattribution
+      rainclr(i) = rainclr_tmp(i)
+      raincld(i) = raincld_tmp(i)
+      snowclr(i) = snowclr_tmp(i)
+      snowcld(i) = snowcld_tmp(i)                                     
+    ENDIF
+
+    !--If there is no more precip fluxes, the precipitation fraction is set to 0
     IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
+    IF ( ( raincld(i) + snowcld(i) ) .LE. 0. ) precipfraccld(i) = 0.
 
     !--Calculation of the total fluxes
@@ -955 +1124 @@
 
 END SUBROUTINE poprecip_precld
+
+
+!----------------------------------------------------------------
+! Computes the precipitation fraction update
+! The goal of this routine is to reattribute precipitation fractions
+! and fluxes to clear or cloudy air, depending on the variation of
+! the cloud fraction on the vertical dimension. We assume a
+! maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
+! and LTP thesis, 2021)
+! NB. in fact, we assume a maximum-random overlap of the total precip. frac
+! 
+SUBROUTINE poprecip_fracupdate( &
+           klon, cldfra, precipfracclr, precipfraccld, &
+           rainclr, raincld, snowclr, snowcld)
+
+USE lmdz_lscp_ini, ONLY : eps
+
+IMPLICIT NONE
+
+INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction [-]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
+                                                          !--NB. at the end of the routine, becomes the fraction of precip
+                                                          !--in the current layer
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
+
+!--Local variables
+INTEGER :: i
+REAL :: dcldfra
+REAL :: precipfractot
+REAL :: dprecipfracclr, dprecipfraccld
+REAL :: drainclr, dsnowclr
+REAL :: draincld, dsnowcld
+
+
+DO i = 1, klon
+
+  !--Initialisation
+  precipfractot = precipfracclr(i) + precipfraccld(i)
+
+  !--Instead of using the cloud cover which was use in LTP thesis, we use the
+  !--total precip. fraction to compute the maximum-random overlap. This is
+  !--because all the information of the cloud cover is embedded into
+  !--precipfractot, and this allows for taking into account the potential
+  !--reduction of the precipitation fraction because either the flux is too
+  !--small (see barrier at the end of poprecip_postcld) or the flux is completely
+  !--evaporated (see barrier at the end of poprecip_precld)
+  !--NB. precipfraccld(i) is here the cloud fraction of the layer above
+  !precipfractot = 1. - ( 1. - precipfractot ) * &
+  !               ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
+  !             / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
+
+
+  IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
+    precipfractot = 1.
+  ELSE
+    precipfractot = 1. - ( 1. - precipfractot ) * &
+                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
+               / ( 1. - precipfraccld(i) )
+  ENDIF
+
+  !--precipfraccld(i) is here the cloud fraction of the layer above
+  dcldfra = cldfra(i) - precipfraccld(i)
+  !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
+  !--calculation of the current CS precip. frac.
+  !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
+  !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
+  !--if precipfractot < cldfra)
+  dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i)
+  !--Tendency of the cloudy precipitation fraction. We add a MAX on the
+  !--calculation of the current CS precip. frac.
+  !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) ) 
+  !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
+  !--if cldfra < 0)
+  dprecipfraccld = dcldfra
+
+
+  !--If the cloud extends
+  IF ( dprecipfraccld .GT. 0. ) THEN
+    !--If there is no CS precip, nothing happens.
+    !--If there is, we reattribute some of the CS precip flux
+    !--to the cloud precip flux, proportionnally to the
+    !--decrease of the CS precip fraction
+    IF ( precipfracclr(i) .LT. eps ) THEN
+      drainclr = 0.
+      dsnowclr = 0.
+    ELSE
+      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i) 
+      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i) 
+    ENDIF
+  !--If the cloud narrows
+  ELSEIF ( dprecipfraccld .LT. 0. ) THEN
+    !--We reattribute some of the cloudy precip flux
+    !--to the CS precip flux, proportionnally to the
+    !--decrease of the cloud precip fraction
+    IF ( precipfraccld(i) .LT. eps ) THEN
+      draincld = 0.
+      dsnowcld = 0.
+    ELSE
+      draincld = dprecipfraccld / precipfraccld(i) * raincld(i) 
+      dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i) 
+    ENDIF
+    drainclr = - draincld
+    dsnowclr = - dsnowcld
+  !--If the cloud stays the same or if there is no cloud above and
+  !--in the current layer, nothing happens
+  ELSE
+    drainclr = 0.
+    dsnowclr = 0.
+  ENDIF
+
+  !--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)
+
+ENDDO
+
+END SUBROUTINE poprecip_fracupdate
 
 
@@ -1043 +1341 @@
 REAL, DIMENSION(klon) :: qtot                             !--includes vap, liq, ice and precip
 
-!--Partition of the fluxes
-REAL :: dcldfra
-REAL :: precipfractot
-REAL :: dprecipfracclr, dprecipfraccld
-REAL :: drainclr, dsnowclr
-REAL :: draincld, dsnowcld
-
 !--Collection, aggregation and riming
 REAL :: eff_cldfra
@@ -1098 +1389 @@
 
 
+!--Update the precipitation fraction following cloud formation
+CALL poprecip_fracupdate( &
+           klon, cldfra, precipfracclr, precipfraccld, &
+           rainclr, raincld, snowclr, snowcld)
+
+
 DO i = 1, klon
 
@@ -1111 +1408 @@
   qtot(i) = qvap(i) + qliq(i) + qice(i) &
           + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i)
-
-  !------------------------------------------------------------
-  !--     PRECIPITATION FRACTIONS UPDATE
-  !------------------------------------------------------------
-  !--The goal of this routine is to reattribute precipitation fractions
-  !--and fluxes to clear or cloudy air, depending on the variation of
-  !--the cloud fraction on the vertical dimension. We assume a
-  !--maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
-  !--and LTP thesis, 2021)
-  !--NB. in fact, we assume a maximum-random overlap of the total precip. frac
-
-  !--Initialisation
-  precipfractot = precipfracclr(i) + precipfraccld(i)
-
-  !--Instead of using the cloud cover which was use in LTP thesis, we use the
-  !--total precip. fraction to compute the maximum-random overlap. This is
-  !--because all the information of the cloud cover is embedded into
-  !--precipfractot, and this allows for taking into account the potential
-  !--reduction of the precipitation fraction because either the flux is too
-  !--small (see barrier at the end of poprecip_postcld) or the flux is completely
-  !--evaporated (see barrier at the end of poprecip_precld)
-  !--NB. precipfraccld(i) is here the cloud fraction of the layer above
-  !precipfractot = 1. - ( 1. - precipfractot ) * &
-  !               ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
-  !             / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
-
-
-  IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
-    precipfractot = 1.
-  ELSE
-    precipfractot = 1. - ( 1. - precipfractot ) * &
-                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
-               / ( 1. - precipfraccld(i) )
-  ENDIF
-
-  !--precipfraccld(i) is here the cloud fraction of the layer above
-  dcldfra = cldfra(i) - precipfraccld(i)
-  !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
-  !--calculation of the current CS precip. frac.
-  !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
-  !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
-  !--if precipfractot < cldfra)
-  dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i)
-  !--Tendency of the cloudy precipitation fraction. We add a MAX on the
-  !--calculation of the current CS precip. frac.
-  !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) ) 
-  !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
-  !--if cldfra < 0)
-  dprecipfraccld = dcldfra
-
-
-  !--If the cloud extends
-  IF ( dprecipfraccld .GT. 0. ) THEN
-    !--If there is no CS precip, nothing happens.
-    !--If there is, we reattribute some of the CS precip flux
-    !--to the cloud precip flux, proportionnally to the
-    !--decrease of the CS precip fraction
-    IF ( precipfracclr(i) .LE. 0. ) THEN
-      drainclr = 0.
-      dsnowclr = 0.
-    ELSE
-      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i) 
-      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i) 
-    ENDIF
-  !--If the cloud narrows
-  ELSEIF ( dprecipfraccld .LT. 0. ) THEN
-    !--We reattribute some of the cloudy precip flux
-    !--to the CS precip flux, proportionnally to the
-    !--decrease of the cloud precip fraction
-    draincld = dprecipfraccld / precipfraccld(i) * raincld(i) 
-    dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i) 
-    drainclr = - draincld
-    dsnowclr = - dsnowcld
-  !--If the cloud stays the same or if there is no cloud above and
-  !--in the current layer, nothing happens
-  ELSE
-    drainclr = 0.
-    dsnowclr = 0.
-  ENDIF
-
-  !--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)
   
   !--If vertical heterogeneity is taken into account, we use