Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90	(revision 5407)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90	(revision 5408)
@@ -108,12 +108,12 @@
 USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
 USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
-USE lmdz_lscp_poprecip, ONLY : poprecip_precld, poprecip_postcld
+USE lmdz_lscp_precip, ONLY : histprecip_precld, poprecip_precld, poprecip_postcld
 
 ! Use du module lmdz_lscp_ini contenant les constantes
 USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
-USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
+USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, DDT0, rain_int_min
 USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
 USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min 
-USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
+USE lmdz_lscp_ini, ONLY : cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
 USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
 USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
@@ -267,5 +267,4 @@
   ! LOCAL VARIABLES:
   !----------------
-  REAL,DIMENSION(klon) :: qsl, qsi                                ! saturation threshold at current vertical level
   REAL,DIMENSION(klon) :: qice_ini
   REAL :: zct, zcl,zexpo
@@ -284,8 +283,6 @@
   REAL, DIMENSION(klon) :: zfice_th
   REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
-  REAL, DIMENSION(klon) :: zrfl, zrfln
-  REAL :: zqev, zqevt 
-  REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
-  REAL :: zqev0,zqevi, zqevti
+  REAL, DIMENSION(klon) :: zrfl
+  REAL, DIMENSION(klon) :: zifl, ziflprev
   REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
   REAL, DIMENSION(klon) :: zoliql, zoliqi
@@ -294,11 +291,10 @@
   REAL, DIMENSION(klon) :: zdz,iwc
   REAL :: zchau,zfroi
-  REAL, DIMENSION(klon) :: zfice,zneb,znebprecip
-  REAL :: zmelt,zrain,zsnow,zprecip
+  REAL, DIMENSION(klon) :: zfice,zneb
+  REAL :: zrain,zsnow,zprecip
   REAL, DIMENSION(klon) :: dzfice
   REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
   REAL :: zsolid
   REAL, DIMENSION(klon) :: qtot, qzero
-  REAL, DIMENSION(klon) :: dqsl,dqsi
   REAL :: smallestreal
   !  Variables for Bergeron process
@@ -310,6 +306,4 @@
   REAL :: zfrac_lessi
   REAL, DIMENSION(klon) :: zprec_cond
-  REAL :: zmair
-  REAL :: zcpair, zcpeau
   REAL, DIMENSION(klon) :: zlh_solid
   REAL, DIMENSION(klon) :: ztupnew
@@ -318,5 +312,5 @@
   REAL, DIMENSION(klon) :: zrflclr, zrflcld
   REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
-  REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr 
+  REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
   REAL, DIMENSION(klon) :: ziflclr, ziflcld
   REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
@@ -369,5 +363,4 @@
 ! Few initialisations
 
-znebprecip(:)=0.0
 ctot_vol(1:klon,1:klev)=0.0
 rneblsvol(1:klon,1:klev)=0.0
@@ -488,4 +481,6 @@
 DO k = klev, 1, -1
 
+    qice_ini = qi_seri(:,k)
+
     IF (k.LE.klev-1) THEN
        iftop=.false. 
@@ -511,257 +506,34 @@
         !c_iso init of iso
     ENDDO
+
+    ! --------------------------------------------------------------------
+    ! P1> Precipitation processes, before cloud formation:
+    !     Thermalization of precipitation falling from the overlying layer AND
+    !     Precipitation evaporation/sublimation/melting
+    !---------------------------------------------------------------------
    
     !================================================================
     ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
-    IF (ok_poprecip) THEN
-
-            CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
-                              zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
-                              zqvapclr, zqupnew, &
-                              zrfl, zrflclr, zrflcld, &
-                              zifl, ziflclr, ziflcld, &
-                              dqreva(:,k),dqssub(:,k) &
-                             )
+    IF ( ok_poprecip .OR. ( iflag_evap_prec .EQ. 6 ) ) THEN
+
+      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
+                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
+                        zqvapclr, zqupnew, &
+                        zrfl, zrflclr, zrflcld, &
+                        zifl, ziflclr, ziflcld, &
+                        dqreva(:,k), dqssub(:,k) &
+                        )
 
     !================================================================
     ELSE
 
-    ! --------------------------------------------------------------------
-    ! P1> Thermalization of precipitation falling from the overlying layer
-    ! --------------------------------------------------------------------
-    ! Computes air temperature variation due to enthalpy transported by
-    ! precipitation. Precipitation is then thermalized with the air in the
-    ! layer. 
-    ! The precipitation should remain thermalized throughout the different
-    ! thermodynamical transformations. 
-    ! The corresponding water mass should
-    ! be added when calculating the layer's enthalpy change with 
-    ! temperature
-    ! See lmdzpedia page todoan
-    ! todoan: check consistency with ice phase
-    ! todoan: understand why several steps
-    ! ---------------------------------------------------------------------
-    
-    IF (iftop) THEN
-
-        DO i = 1, klon
-            zmqc(i) = 0.
-        ENDDO
-
-    ELSE
-
-        DO i = 1, klon
- 
-            zmair=(paprs(i,k)-paprs(i,k+1))/RG
-            ! no condensed water so cp=cp(vapor+dry air)
-            ! RVTMP2=rcpv/rcpd-1
-            zcpair=RCPD*(1.0+RVTMP2*zq(i)) 
-            zcpeau=RCPD*RVTMP2
-
-            ! zmqc: precipitation mass that has to be thermalized with 
-            ! layer's air so that precipitation at the ground has the
-            ! same temperature as the lowermost layer
-            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
-            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
-            zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) &
-             / (zcpair + zmqc(i)*zcpeau)
-
-        ENDDO
- 
-    ENDIF
-
-    qice_ini = qi_seri(:,k)
-    ! --------------------------------------------------------------------
-    ! P2> Precipitation evaporation/sublimation/melting
-    ! --------------------------------------------------------------------
-    ! A part of the precipitation coming from above is evaporated/sublimated/melted.
-    ! --------------------------------------------------------------------
-    
-    IF (iflag_evap_prec.GE.1) THEN 
-
-        ! Calculation of saturation specific humidity
-        ! depending on temperature:
-        CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,0,.false.,zqs,zdqs)
-        ! wrt liquid water
-        CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsl,dqsl)
-        ! wrt ice
-        CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsi,dqsi)
-
-        DO i = 1, klon
-
-            ! if precipitation
-            IF (zrfl(i)+zifl(i).GT.0.) THEN
-
-            ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
-            ! c_iso: likely important to distinguish cs from neb isotope precipitation
-
-                IF (iflag_evap_prec.GE.4) THEN
-                    zrfl(i) = zrflclr(i)
-                    zifl(i) = ziflclr(i)
-                ENDIF
-
-                IF (iflag_evap_prec.EQ.1) THEN
-                    znebprecip(i)=zneb(i)
-                ELSE
-                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
-                ENDIF
-
-                IF (iflag_evap_prec.GT.4) THEN
-                ! Max evaporation not to saturate the clear sky precip fraction
-                ! i.e. the fraction where evaporation occurs
-                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
-                ELSEIF (iflag_evap_prec .EQ. 4) THEN
-                ! Max evaporation not to saturate the whole mesh
-                ! Pay attention -> lead to unrealistic and excessive evaporation
-                    zqev0 = MAX(0.0, zqs(i)-zq(i))
-                ELSE
-                ! Max evap not to saturate the fraction below the cloud
-                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
-                ENDIF
-
-                ! Evaporation of liquid precipitation coming from above
-                ! dP/dz=beta*(1-q/qsat)*sqrt(P)
-                ! formula from Sundquist 1988, Klemp & Wilhemson 1978
-                ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
-
-                IF (iflag_evap_prec.EQ.3) THEN
-                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
-                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
-                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-                ELSE IF (iflag_evap_prec.GE.4) THEN
-                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
-                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
-                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-                ELSE
-                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
-                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-                ENDIF
-
-                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
-                    *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 
-
-                ! sublimation of the solid precipitation coming from above
-                IF (iflag_evap_prec.EQ.3) THEN
-                    zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
-                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
-                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-                ELSE IF (iflag_evap_prec.GE.4) THEN
-                     zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
-                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
-                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-                ELSE
-                    zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
-                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-                ENDIF
-
-                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
-                *RG*dtime/(paprs(i,k)-paprs(i,k+1))
-
-                ! A. JAM
-                ! Evaporation limit: we ensure that the layer's fraction below
-                ! the cloud or the whole mesh (depending on iflag_evap_prec)
-                ! does not reach saturation. In this case, we 
-                ! redistribute zqev0 conserving the ratio liquid/ice 
-
-                IF (zqevt+zqevti.GT.zqev0) THEN
-                    zqev=zqev0*zqevt/(zqevt+zqevti)
-                    zqevi=zqev0*zqevti/(zqevt+zqevti)
-                ELSE
-                    zqev=zqevt
-                    zqevi=zqevti
-                ENDIF
-
-                !--Diagnostics
-                dqreva(i,k) = - zqev / dtime
-                dqssub(i,k) = - zqevti / dtime
-
-                ! New solid and liquid precipitation fluxes after evap and sublimation
-                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
-                         /RG/dtime)
-                zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
-                         /RG/dtime)
-                
-
-                ! vapor, temperature, precip fluxes update
-                ! vapor is updated after evaporation/sublimation (it is increased)
-                zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
-                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-                ! zmqc is the total condensed water in the precip flux (it is decreased)
-                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
-                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-                ! air and precip temperature (i.e., gridbox temperature)
-                ! is updated due to latent heat cooling
-                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
-                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
-                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
-                + (zifln(i)-zifl(i))                      &
-                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
-                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
-
-                ! New values of liquid and solid precipitation 
-                zrfl(i) = zrfln(i)
-                zifl(i) = zifln(i)
-
-                ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
-                !                         due to evap + sublim                                     
-                                           
-
-                IF (iflag_evap_prec.GE.4) THEN
-                    zrflclr(i) = zrfl(i)
-                    ziflclr(i) = zifl(i)
-                    IF(zrflclr(i) + ziflclr(i).LE.0) THEN
-                        znebprecipclr(i) = 0.0
-                    ENDIF
-                    zrfl(i) = zrflclr(i) + zrflcld(i)
-                    zifl(i) = ziflclr(i) + ziflcld(i)
-                ENDIF
-
-                ! c_iso duplicate for isotopes or loop on isotopes
-
-                ! Melting:
-                zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
-                ! precip fraction that is melted
-                zmelt = MIN(MAX(zmelt,0.),1.)
-
-                ! update of rainfall and snowfall due to melting
-                IF (iflag_evap_prec.GE.4) THEN
-                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
-                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
-                    zrfl(i)=zrflclr(i)+zrflcld(i)
-                ELSE
-                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
-                ENDIF
-
-
-                ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
-
-                ! Latent heat of melting because of precipitation melting
-                ! NB: the air + precip temperature is simultaneously updated
-                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
-                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
-               
-                IF (iflag_evap_prec.GE.4) THEN
-                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
-                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
-                    zifl(i)=ziflclr(i)+ziflcld(i)
-                ELSE
-                    zifl(i)=zifl(i)*(1.-zmelt)
-                ENDIF
-
-            ELSE
-                ! if no precip, we reinitialize the cloud fraction used for the precip to 0 
-                znebprecip(i)=0.
-
-            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
-
-        ENDDO ! loop on klon
-
-    ENDIF ! (iflag_evap_prec>=1)
+      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
+                        zt, ztupnew, zq, zmqc, zneb, znebprecipclr, &
+                        zrfl, zrflclr, zrflcld, &
+                        zifl, ziflclr, ziflcld, &
+                        dqreva(:,k), dqssub(:,k) &
+                        )
 
     ENDIF ! (ok_poprecip)
-
-    ! --------------------------------------------------------------------
-    ! End precip evaporation
-    ! --------------------------------------------------------------------
     
     ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
@@ -782,5 +554,5 @@
     
     ! --------------------------------------------------------------------
-    ! P3> Cloud formation
+    ! P2> Cloud formation
     !---------------------------------------------------------------------
     !
@@ -793,5 +565,5 @@
     ! -------------------------------------------------------------------
 
-        ! P3.1> With the PDFs (log-normal, bigaussian) 
+        ! P2.1> With the PDFs (log-normal, bigaussian) 
         ! cloud properties calculation with the initial values of t and q
         ! ----------------------------------------------------------------
@@ -1038,5 +810,5 @@
         ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
 
-        ! P3.2> Final quantities calculation 
+        ! P2.2> Final quantities calculation 
         ! Calculated variables: 
         !   rneb : cloud fraction
@@ -1136,8 +908,10 @@
 
     ! ----------------------------------------------------------------
-    ! P4> Precipitation formation 
+    ! P3> Precipitation processes, after cloud formation 
+    !     - precipitation formation, melting/freezing
     ! ----------------------------------------------------------------
 
     !================================================================
+    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
     IF (ok_poprecip) THEN
 
@@ -1468,10 +1242,10 @@
     ENDIF ! ok_poprecip
 
-    ! End of precipitation formation
-    ! --------------------------------
+    ! End of precipitation processes after cloud formation
+    ! ----------------------------------------------------
 
 
     !----------------------------------------------------------------------
-    ! P5> Calculation of cloud condensates amount seen by radiative scheme
+    ! P4> Calculation of cloud condensates amount seen by radiative scheme
     !----------------------------------------------------------------------
 
@@ -1506,5 +1280,5 @@
 
     ! ----------------------------------------------------------------
-    ! P6> Wet scavenging
+    ! P5> Wet scavenging
     ! ----------------------------------------------------------------
 
@@ -1563,5 +1337,5 @@
     
     !------------------------------------------------------------
-    ! P7 > write diagnostics and outputs
+    ! P6 > write diagnostics and outputs
     !------------------------------------------------------------
    
Index: DZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.f90	(revision 5407)
+++ 	(revision )
@@ -1,980 +1,0 @@
-MODULE lmdz_lscp_poprecip
-!----------------------------------------------------------------
-! Module for the process-oriented treament of precipitation
-! that are called in LSCP
-! Authors: Atelier Nuage (G. Riviere, L. Raillard, M. Wimmer,
-! N. Dutrievoz, E. Vignon, A. Borella, et al.)
-! Jan. 2024
-
-
-IMPLICIT NONE
-
-CONTAINS
-
-!----------------------------------------------------------------
-! Computes the processes-oriented precipitation formulations for
-! evaporation and sublimation
-! 
-SUBROUTINE poprecip_precld( &
-           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
-           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
-           rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
-           )
-
-USE lmdz_lscp_ini, ONLY : prt_level, lunout
-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_tools, ONLY : calc_qsat_ecmwf
-
-IMPLICIT NONE
-
-
-INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
-REAL,    INTENT(IN)                     :: dtime          !--time step [s]
-LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
-
-
-REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
-REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
-REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
-
-REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
-REAL,    INTENT(INOUT), DIMENSION(klon) :: tempupnew      !--updated temperature of the overlying layer [K]
-
-REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
-REAL,    INTENT(INOUT), DIMENSION(klon) :: qprecip        !--specific humidity in the precipitation falling from the upper layer [kg/kg]
-
-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 [-]
-
-REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
-REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [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) :: 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(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
-
-
-!--Integer for interating over klon
-INTEGER :: i
-!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
-REAL, DIMENSION(klon) :: dhum_to_dflux
-!--
-REAL, DIMENSION(klon) :: rho, dz
-
-!--Saturation values
-REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
-!--Vapor in the clear sky
-REAL :: qvapclr
-!--Fluxes tendencies because of evaporation and sublimation
-REAL :: dprecip_evasub_max, draineva, dsnowsub, dprecip_evasub_tot
-!--Specific humidity tendencies because of evaporation and sublimation
-REAL :: dqrevap, dqssubl
-!--Specific heat constant
-REAL :: cpair, cpw
-
-!--Initialisation
-qzero(:)  = 0.
-dqreva(:) = 0.
-dqssub(:) = 0.
-dqrevap   = 0.
-dqssubl   = 0.
-
-!-- 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
-!--depending on temperature:
-CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:))
-!--wrt liquid water
-CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
-!--wrt ice
-CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
-
-
-
-!--First step consists in "thermalizing" the layer:
-!--as the flux of precip from layer above "advects" some heat (as the precip is at the temperature 
-!--of the overlying layer) we recalculate a mean temperature that both the air and the precip in the 
-!--layer have.
-
-IF (iftop) THEN
-
-  DO i = 1, klon
-    qprecip(i) = 0.
-  ENDDO
-
-ELSE
-
-  DO i = 1, klon
-    !--No condensed water so cp=cp(vapor+dry air)
-    !-- RVTMP2=rcpv/rcpd-1
-    cpair = RCPD * ( 1. + RVTMP2 * qvap(i) )
-    cpw = RCPD * RVTMP2
-    !--qprecip has to be thermalized with 
-    !--layer's air so that precipitation at the ground has the
-    !--same temperature as the lowermost layer
-    !--we convert the flux into a specific quantity qprecip
-    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) ) &
-            / ( cpair + qprecip(i) * cpw )
-  ENDDO
-
-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.
-
-
-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
-      !--the clear and the cloudy sky as in the layer above. This
-      !--extends the assumption that the cloud fraction is the same
-      !--as the layer above. This is assumed only for the evap / subl
-      !--process
-      !--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) )
-    ELSE
-      !--Legacy version from Ludo - we use the total specific humidity
-      !--for the evap / subl process
-      qvapclr = qvap(i)
-    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
-
-
-    !--New solid and liquid precipitation fluxes after evap and sublimation
-    dqrevap = draineva / dhum_to_dflux(i)
-    dqssubl = dsnowsub / dhum_to_dflux(i)
-
-
-    !--Vapor is updated after evaporation/sublimation (it is increased)
-    qvap(i) = qvap(i) - dqrevap - dqssubl
-    !--qprecip is the total condensed water in the precip flux (it is decreased)
-    qprecip(i) = qprecip(i) + dqrevap + dqssubl
-    !--Air and precip temperature (i.e., gridbox temperature)
-    !--is updated due to latent heat cooling
-    temp(i) = temp(i) &
-            + dqrevap * RLVTT / RCPD &
-            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
-            + dqssubl * RLSTT / RCPD &
-            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
-
-    !--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
-    IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
-
-    !--Calculation of the total fluxes
-    rain(i) = rainclr(i) + raincld(i)
-    snow(i) = snowclr(i) + snowcld(i)
-
-  ELSEIF ( ( rain(i) + snow(i) ) .LE. 0. ) THEN
-    !--If no precip, we reinitialize the cloud fraction used for the precip to 0 
-    precipfraccld(i) = 0.
-    precipfracclr(i) = 0.
-
-  ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
-
-  !--Diagnostic tendencies
-  dqssub(i) = dqssubl / dtime
-  dqreva(i) = dqrevap / dtime
-
-ENDDO ! loop on klon
-
-
-END SUBROUTINE poprecip_precld
-
-!----------------------------------------------------------------
-! Computes the processes-oriented precipitation formulations for
-! - autoconversion (auto) via a deposition process
-! - aggregation (agg)
-! - riming (rim)
-! - collection (col)
-! - melting (melt)
-! - freezing (freez)
-! 
-SUBROUTINE poprecip_postcld( &
-           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
-           temp, qvap, qliq, qice, icefrac, cldfra, &
-           precipfracclr, precipfraccld, &
-           rain, rainclr, raincld, snow, snowclr, snowcld, &
-           qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
-           dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
-
-USE lmdz_lscp_ini, ONLY : prt_level, lunout
-USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
-USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
-
-USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
-                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & 
-                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
-                          rho_rain, r_rain, r_snow, rho_ice,                   &
-                          tau_auto_snow_min, tau_auto_snow_max,                &
-                          thresh_precip_frac, eps,                             &
-                          gamma_melt, alpha_freez, beta_freez, temp_nowater,   &
-                          iflag_cloudth_vert, iflag_rain_incloud_vol,          &
-                          cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez,       &
-                          rain_fallspeed_clr, rain_fallspeed_cld,              &
-                          snow_fallspeed_clr, snow_fallspeed_cld
-
-
-IMPLICIT NONE
-
-INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
-REAL,    INTENT(IN)                     :: dtime          !--time step [s]
-
-REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
-REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
-REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
-
-REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
-LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
-
-REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
-REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity [kg/kg]
-REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--current liquid water specific humidity [kg/kg]
-REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--current ice water specific humidity [kg/kg]
-REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
-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) :: 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(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(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) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsrim         !--snow tendency due to riming [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
-REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
-
-
-
-!--Local variables
-
-INTEGER :: i
-!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
-REAL, DIMENSION(klon) :: dhum_to_dflux
-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
-REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
-REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
-REAL :: rho_snow
-REAL :: dqlcol           !--loss of liquid cloud content due to collection by rain [kg/kg/s]
-REAL :: dqiagg           !--loss of ice cloud content due to collection by aggregation [kg/kg/s]
-REAL :: dqlrim           !--loss of liquid cloud content due to riming on snow [kg/kg/s]
-
-!--Autoconversion
-REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
-REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
-REAL :: dqlauto          !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
-REAL :: dqiauto          !--loss of ice cloud content due to autoconversion to snow [kg/kg/s]
-
-!--Melting
-REAL :: dqsmelt_max, air_thermal_conduct
-REAL :: nb_snowflake_clr, nb_snowflake_cld
-REAL :: capa_snowflake, temp_wetbulb
-REAL :: rho, r_ice
-REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
-REAL, DIMENSION(klon) :: qzero, qsat, dqsat
-
-!--Freezing
-REAL :: dqrfreez_max
-REAL :: tau_freez
-REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
-REAL :: coef_freez
-REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
-REAL :: Eff_rain_ice
-
-
-!--Initialisation of variables
-
-
-qzero(:)    = 0.
-
-dqrcol(:)   = 0.
-dqsagg(:)   = 0.
-dqsauto(:)  = 0.
-dqrauto(:)  = 0.
-dqsrim(:)   = 0.
-dqrmelt(:)  = 0.
-dqsmelt(:)  = 0.
-dqrfreez(:) = 0.
-dqsfreez(:) = 0.
-
-
-DO i = 1, klon
-
-  !--Variables initialisation
-  dqlcol  = 0.
-  dqiagg  = 0.
-  dqiauto = 0.
-  dqlauto = 0.
-  dqlrim  = 0.
-
-  !-- 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) ) / 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
-  !--the "true" volume fraction instead of a modified 
-  !--surface fraction (which is larger and artificially
-  !--reduces the in-cloud water).
-  IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
-    eff_cldfra = ctot_vol(i)
-  ELSE
-    eff_cldfra = cldfra(i)
-  ENDIF 
-
-
-  !--Start precipitation growth processes
-
-  !--If the cloud is big enough, the precipitation processes activate
-  ! TODO met on seuil_neb ici ?
-  IF ( cldfra(i) .GE. seuil_neb ) THEN
-
-    !---------------------------------------------------------
-    !--           COLLECTION AND AGGREGATION
-    !---------------------------------------------------------
-    !--Collection: processus through which rain collects small liquid droplets
-    !--in suspension, and add it to the rain flux
-    !--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
-    !--Those processes are treated before autoconversion because we do not
-    !--want to collect/aggregate the newly formed fluxes, which already
-    !--"saw" the cloud as they come from it
-    !--The formulas come from Muench and Lohmann 2020
-
-    !--gamma_col: tuning coefficient [-]
-    !--rho_rain: volumic mass of rain [kg/m3]
-    !--r_rain: size of the rain droplets [m]
-    !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
-    !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
-    !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
-    !--get in-cloud mean quantities. The two divisions by eff_cldfra are
-    !--then simplified.
-
-    !--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. ) THEN
-      !--Exact explicit version, which does not need a barrier because of
-      !--the exponential decrease
-      dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
-
-      !--Add tendencies
-      qliq(i) = qliq(i) + dqlcol
-      raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i)
-
-      !--Diagnostic tendencies
-      dqrcol(i)  = - dqlcol  / dtime
-    ENDIF
-
-    !--Same as for aggregation
-    !--Eff_snow_liq formula:
-    !--it s a product of a collection efficiency and a sticking efficiency
-    ! Milbrandt and Yau formula that gives very low values:
-    ! Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
-    ! Lin 1983's formula 
-    Eff_snow_ice = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) ) 
-    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
-    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
-    coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
-    IF ( snowcld(i) .GT. 0. ) THEN
-      !--Exact explicit version, which does not need a barrier because of
-      !--the exponential decrease
-      dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
-
-      !--Add tendencies
-      qice(i) = qice(i) + dqiagg
-      snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
-
-      !--Diagnostic tendencies
-      dqsagg(i)  = - dqiagg  / dtime
-    ENDIF
-
-
-    !---------------------------------------------------------
-    !--                  AUTOCONVERSION
-    !---------------------------------------------------------
-    !--Autoconversion converts liquid droplets/ice crystals into
-    !--rain drops/snowflakes. It relies on the formulations by
-    !--Sundqvist 1978.
-
-    !--If we are in a convective point, we have different parameters
-    !--for the autoconversion
-    IF ( ptconv(i) ) THEN
-        qthresh_auto_rain = cld_lc_con
-        qthresh_auto_snow = cld_lc_con_snow
-
-        tau_auto_rain = cld_tau_con
-        !--tau for snow depends on the ice fraction in mixed-phase clouds     
-        tau_auto_snow = tau_auto_snow_max &
-                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
-
-        expo_auto_rain = cld_expo_con
-        expo_auto_snow = cld_expo_con
-    ELSE
-        qthresh_auto_rain = cld_lc_lsc
-        qthresh_auto_snow = cld_lc_lsc_snow
-
-        tau_auto_rain = cld_tau_lsc
-        !--tau for snow depends on the ice fraction in mixed-phase clouds     
-        tau_auto_snow = tau_auto_snow_max &
-                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
-
-        expo_auto_rain = cld_expo_lsc
-        expo_auto_snow = cld_expo_lsc
-    ENDIF
-
-
-    ! Liquid water quantity to remove according to (Sundqvist, 1978)
-    ! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
-    !
-    !--And same formula for ice
-    !
-    !--We first treat the second term (with exponential) in an explicit way
-    !--and then treat the first term (-q/tau) in an exact way
-
-    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
-               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
-
-    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
-               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
-
-
-    !--Barriers so that we don't create more rain/snow
-    !--than there is liquid/ice
-    dqlauto = MAX( - qliq(i), dqlauto )
-    dqiauto = MAX( - qice(i), dqiauto )
-
-    !--Add tendencies
-    qliq(i) = qliq(i) + dqlauto
-    qice(i) = qice(i) + dqiauto
-    raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
-    snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
-
-    !--Diagnostic tendencies
-    dqsauto(i) = - dqiauto / dtime
-    dqrauto(i) = - dqlauto / dtime
-
-
-    !---------------------------------------------------------
-    !--                  RIMING
-    !---------------------------------------------------------
-    !--Process which converts liquid droplets in suspension into
-    !--snow because of the collision between
-    !--those and falling snowflakes.
-    !--The formula comes from Muench and Lohmann 2020
-    !--NB.: this process needs a temperature adjustment
-
-    !--Eff_snow_liq formula: following Ferrier 1994,
-    !--assuming 1
-    Eff_snow_liq = 1.0
-    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
-    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
-    coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
-    IF ( snowcld(i) .GT. 0. ) THEN
-      !--Exact version, which does not need a barrier because of
-      !--the exponential decrease
-      dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. )
-
-      !--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) = snowcld(i) - dqlrim * dhum_to_dflux(i)
-
-      !--Temperature adjustment with the release of latent
-      !--heat because of solid condensation
-      temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
-                        / ( 1. + RVTMP2 * qtot(i) )
-
-      !--Diagnostic tendencies
-      dqsrim(i)  = - dqlrim  / dtime
-    ENDIF
-
-  ENDIF ! cldfra .GE. seuil_neb
-
-ENDDO ! loop on klon
-
-
-!--Re-calculation of saturation specific humidity
-!--because riming changed temperature
-CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
-
-DO i = 1, klon
-
-  !---------------------------------------------------------
-  !--                  MELTING
-  !---------------------------------------------------------
-  !--Process through which snow melts into rain. 
-  !--The formula is homemade.
-  !--NB.: this process needs a temperature adjustment
-
-  !--dqsmelt_max         : maximum snow melting so that temperature
-  !--                      stays higher than 273 K [kg/kg]
-  !--capa_snowflake      : capacitance of a snowflake, equal to
-  !--                      the radius if the snowflake is a sphere [m]
-  !--temp_wetbulb        : wet-bulb temperature [K]
-  !--snow_fallspeed      : snow fall velocity (in clear/cloudy sky) [m/s]
-  !--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
-  !--gamma_melt          : tuning parameter for melting [-]
-  !--nb_snowflake        : number of snowflakes (in clear/cloudy air) [-]
-
-  IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
-    !--Computed according to
-    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
-    dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
-                        * ( 1. + RVTMP2 * qtot(i) ))
-   
-    !--Initialisation 
-    dqsclrmelt = 0.
-    dqscldmelt = 0.
-
-    !--We assume that the snowflakes are spherical
-    capa_snowflake = r_snow
-    !--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
-    air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184    
-    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
-    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
-
-    !--In clear air
-    IF ( ( snowclr(i) .GT. 0. ) .AND.  ( precipfracclr(i) .GT. 0. ) ) THEN
-      !--Formula for the wet-bulb temperature from ECMWF (IFS)
-      !--The vapor used is the vapor in the clear sky
-      temp_wetbulb = temp(i) &
-                   - ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
-                   * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
-                   - 40.637 * ( temp(i) - 275. ) )
-      !--Calculated according to
-      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
-      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
-                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
-      dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
-                   * capa_snowflake / RLMLT * gamma_melt &
-                   * MAX(0., temp_wetbulb - RTT) * dtime
-      
-      !--Barrier to limit the melting flux to the clr snow flux in the mesh
-      dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
-    ENDIF
-
-
-    !--In cloudy air
-    IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
-      !--As the air is saturated, the wet-bulb temperature is equal to the
-      !--temperature
-      temp_wetbulb = temp(i)
-      !--Calculated according to
-      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
-      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
-                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
-      dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
-                   * capa_snowflake / RLMLT * gamma_melt &
-                   * MAX(0., temp_wetbulb - RTT) * dtime
-
-      !--Barrier to limit the melting flux to the cld snow flux in the mesh
-      dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
-   ENDIF
-    
-
-    !--Barrier on temperature. If the total melting flux leads to a
-    !--positive temperature, it is limited to keep temperature above 0 degC.
-    !--It is activated if the total is LOWER than the max
-    !--because everything is negative
-    dqstotmelt = dqsclrmelt + dqscldmelt
-    IF ( dqstotmelt .LT. dqsmelt_max ) THEN
-      !--We redistribute the max melted snow keeping
-      !--the clear/cloud partition of the melted snow
-      dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
-      dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
-      dqstotmelt = dqsmelt_max
-
-    ENDIF
-
-    !--Add tendencies
-    !--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))
-
-    !--Temperature adjustment with the release of latent
-    !--heat because of melting
-    temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
-                      / ( 1. + RVTMP2 * qtot(i) )
-
-    !--Diagnostic tendencies
-    dqrmelt(i) = - dqstotmelt / dtime
-    dqsmelt(i) = dqstotmelt / dtime
-
-  ENDIF
-
-
-  !---------------------------------------------------------
-  !--                  FREEZING
-  !---------------------------------------------------------
-  !--Process through which rain freezes into snow. 
-  !-- We parameterize it as a 2 step process:
-  !--first: freezing following collision with ice crystals
-  !--second: immersion freezing following (inspired by Bigg 1953)
-  !--the latter is parameterized as an exponential decrease of the rain
-  !--water content with a homemade formulya
-  !--This is based on a caracteritic time of freezing, which
-  !--exponentially depends on temperature so that it is
-  !--equal to 1 for temp_nowater (see below) and is close to
-  !--0 for RTT (=273.15 K).
-  !--NB.: this process needs a temperature adjustment
-  !--dqrfreez_max : maximum rain freezing so that temperature
-  !--              stays lower than 273 K [kg/kg]
-  !--tau_freez    : caracteristic time of freezing [s]
-  !--gamma_freez  : tuning parameter [s-1]
-  !--alpha_freez  : tuning parameter for the shape of the exponential curve [-]
-  !--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
-
-  IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
-
- 
-    !--1st step: freezing following collision with ice crystals
-    !--Sub-process of freezing which quantifies the collision between
-    !--ice crystals in suspension and falling rain droplets.
-    !--The rain droplets freeze, becoming graupel, and carrying
-    !--the ice crystal (which acted as an ice nucleating particle).
-    !--The formula is adapted from the riming formula.
-    !--it works only in the cloudy part
-
-    dqifreez = 0.
-    dqrtotfreez_step1 = 0.
-
-    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
-         ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
-      dqrclrfreez = 0.
-      dqrcldfreez = 0.
-
-      !--Computed according to
-      !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
-      dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
-                         * ( 1. + RVTMP2 * qtot(i) ))
- 
-
-      !--The collision efficiency is assumed unity 
-      Eff_rain_ice = 1.
-      coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
-      !--Exact version, which does not need a barrier because of
-      !--the exponential decrease.
-      dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
-
-      !--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)
-      !--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)
-      !--Dry density [kg/m3]
-      rho = pplay(i) / temp(i) / RD
-      !--r_ice formula from Sun and Rikus (1999)
-      r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
-            + 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
-      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 * 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) )
-
-    ENDIF
-    
-    !-- Second step immersion freezing of rain drops
-    !-- with a homemade timeconstant depending on temperature
-
-    dqrclrfreez = 0.
-    dqrcldfreez = 0.
-    dqrtotfreez_step2 = 0.
-    !--Computed according to
-    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
-
-    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
-                         * ( 1. + RVTMP2 * qtot(i) ))
-
-    
-    tau_freez = 1. / ( beta_freez &
-              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
-
-
-    !--In clear air
-    IF ( rainclr(i) .GT. 0. ) THEN
-      !--Exact solution of dqrain/dt = -qrain/tau_freez
-      dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
-    ENDIF
-
-    !--In cloudy air
-    IF ( raincld(i) .GT. 0. ) THEN
-      !--Exact solution of dqrain/dt = -qrain/tau_freez
-      dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
-    ENDIF
-
-    !--temperature barrier step 2
-    !--It is activated if the total is LOWER than the max
-    !--because everything is negative
-    dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
-  
-    IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
-      !--We redistribute the max freezed rain keeping
-      !--the clear/cloud partition of the freezing rain
-      dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
-      dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
-      dqrtotfreez_step2 = dqrfreez_max
-    ENDIF
-
-
-    !--Add tendencies
-    !--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))
-
-
-    !--Temperature adjustment with the uptake of latent
-    !--heat because of freezing
-    temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
-                      / ( 1. + RVTMP2 * qtot(i) )
-
-    !--Diagnostic tendencies
-    dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2          
-    dqrfreez(i) = dqrtotfreez / dtime
-    dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
-
-  ENDIF
-
-
-
-  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
-  !--we reduce the precipiration fraction in the clear/cloudy air so that the new
-  !--local flux of rain+snow is equal to rain_int_min.
-  !--Here, rain+snow is the gridbox-mean flux of precip.
-  !--Therefore, (rain+snow)/precipfrac is the local flux of precip.
-  !--If the local flux of precip is lower than rain_int_min, i.e.,
-  !-- (rain+snow)/precipfrac < rain_int_min , i.e., 
-  !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
-  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
-  !--Note that this is physically different than what is proposed in LTP thesis.
-  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
-  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
-
-  !--Calculate outputs
-  rain(i) = rainclr(i) + raincld(i)
-  snow(i) = snowclr(i) + snowcld(i)
-
-  !--Diagnostics
-  !--BEWARE this is indeed a diagnostic: this is an estimation from
-  !--the value of the flux at the bottom interface of the mesh and
-  !--and assuming an upstream numerical calculation
-  !--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
-  !--used for computing the total ice water content in the mesh
-  !--for radiation only
-  qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
-                 + raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
-  qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
-                 + snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
-
-
-ENDDO ! loop on klon
-
-
-END SUBROUTINE poprecip_postcld
-
-END MODULE lmdz_lscp_poprecip
Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90	(revision 5408)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90	(revision 5408)
@@ -0,0 +1,1274 @@
+MODULE lmdz_lscp_precip
+!----------------------------------------------------------------
+! Module for the process-oriented treament of precipitation
+! that are called in LSCP
+! Authors: Atelier Nuage (G. Riviere, L. Raillard, M. Wimmer,
+! N. Dutrievoz, E. Vignon, A. Borella, et al.)
+! Jan. 2024
+
+
+IMPLICIT NONE
+
+CONTAINS
+
+!----------------------------------------------------------------
+! historical (till CMIP6) version of the pre-cloud formation
+! precipitation scheme containing precip evaporation and melting
+
+SUBROUTINE histprecip_precld( &
+           klon, dtime, iftop, paprsdn, paprsup, pplay, zt, ztupnew, zq, &
+           zmqc, zneb, znebprecipclr, &
+           zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub &
+           )
+
+USE lmdz_lscp_ini, ONLY : iflag_evap_prec
+USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, ztfondue
+USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
+USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
+
+IMPLICIT NONE
+
+
+INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
+REAL,    INTENT(IN)                     :: dtime          !--time step [s]
+LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
+
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: zt             !--current temperature [K]
+REAL,    INTENT(IN),    DIMENSION(klon) :: ztupnew        !--updated temperature of the overlying layer [K]
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: zq             !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: znebprecipclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: zrfl           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflcld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: zifl           !--flux of snow gridbox-mean coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--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]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
+
+
+REAL :: zmair
+REAL :: zcpair, zcpeau
+REAL, DIMENSION(klon) :: znebprecip
+REAL, DIMENSION(klon) :: qzero
+REAL, DIMENSION(klon) :: zqs, zdqs
+REAL, DIMENSION(klon) :: qsl, qsi
+REAL, DIMENSION(klon) :: dqsl, dqsi
+
+REAL :: zqev0
+REAL :: zqev, zqevt
+REAL :: zqevi, zqevti
+
+REAL, DIMENSION(klon) :: zrfln, zifln
+
+REAL :: zmelt
+INTEGER :: i
+
+qzero(:) = 0.
+
+
+! --------------------------------------------------------------------
+! P1> Thermalization of precipitation falling from the overlying layer
+! --------------------------------------------------------------------
+! Computes air temperature variation due to enthalpy transported by
+! precipitation. Precipitation is then thermalized with the air in the
+! layer. 
+! The precipitation should remain thermalized throughout the different
+! thermodynamical transformations. 
+! The corresponding water mass should
+! be added when calculating the layer's enthalpy change with 
+! temperature
+! ---------------------------------------------------------------------
+
+IF (iftop) THEN
+
+  DO i = 1, klon
+    zmqc(i) = 0.
+  ENDDO
+
+ELSE
+
+  DO i = 1, klon
+
+    zmair=(paprsdn(i)-paprsup(i))/RG
+    ! no condensed water so cp=cp(vapor+dry air)
+    ! RVTMP2=rcpv/rcpd-1
+    zcpair=RCPD*(1.0+RVTMP2*zq(i)) 
+    zcpeau=RCPD*RVTMP2
+
+    ! zmqc: precipitation mass that has to be thermalized with 
+    ! layer's air so that precipitation at the ground has the
+    ! same temperature as the lowermost layer
+    zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
+    ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
+    zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) &
+          / (zcpair + zmqc(i)*zcpeau)
+
+  ENDDO
+
+ENDIF
+
+! --------------------------------------------------------------------
+! P2> Precipitation evaporation/sublimation/melting
+! --------------------------------------------------------------------
+! A part of the precipitation coming from above is evaporated/sublimated/melted.
+! --------------------------------------------------------------------
+
+IF (iflag_evap_prec.GE.1) THEN 
+
+  ! Calculation of saturation specific humidity
+  ! depending on temperature:
+  CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,0,.false.,zqs,zdqs)
+  ! wrt liquid water
+  CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,1,.false.,qsl,dqsl)
+  ! wrt ice
+  CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,2,.false.,qsi,dqsi)
+
+  DO i = 1, klon
+
+    ! if precipitation
+    IF (zrfl(i)+zifl(i).GT.0.) THEN
+
+    ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
+    ! c_iso: likely important to distinguish cs from neb isotope precipitation
+
+      IF (iflag_evap_prec.GE.4) THEN
+        zrfl(i) = zrflclr(i)
+        zifl(i) = ziflclr(i)
+      ENDIF
+
+      IF (iflag_evap_prec.EQ.1) THEN
+        znebprecip(i)=zneb(i)
+      ELSE
+        znebprecip(i)=MAX(zneb(i),znebprecip(i))
+      ENDIF
+
+      IF (iflag_evap_prec.GT.4) THEN
+      ! Max evaporation not to saturate the clear sky precip fraction
+      ! i.e. the fraction where evaporation occurs
+        zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
+      ELSEIF (iflag_evap_prec .EQ. 4) THEN
+      ! Max evaporation not to saturate the whole mesh
+      ! Pay attention -> lead to unrealistic and excessive evaporation
+        zqev0 = MAX(0.0, zqs(i)-zq(i))
+      ELSE
+      ! Max evap not to saturate the fraction below the cloud
+        zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
+      ENDIF
+
+      ! Evaporation of liquid precipitation coming from above
+      ! dP/dz=beta*(1-q/qsat)*sqrt(P)
+      ! formula from Sundquist 1988, Klemp & Wilhemson 1978
+      ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
+
+      IF (iflag_evap_prec.EQ.3) THEN
+        zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
+          *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
+          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
+      ELSE IF (iflag_evap_prec.GE.4) THEN
+        zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
+          *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
+          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
+      ELSE
+        zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
+          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
+      ENDIF
+
+      zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * RG*dtime/(paprsdn(i)-paprsup(i)) 
+
+      ! sublimation of the solid precipitation coming from above
+      IF (iflag_evap_prec.EQ.3) THEN
+        zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
+          *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
+          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
+      ELSE IF (iflag_evap_prec.GE.4) THEN
+        zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
+          *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
+          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
+      ELSE
+        zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
+          *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG
+      ENDIF
+
+      zqevti = MAX(0.0,MIN(zqevti,zifl(i))) * RG*dtime/(paprsdn(i)-paprsup(i))
+
+      ! A. JAM
+      ! Evaporation limit: we ensure that the layer's fraction below
+      ! the cloud or the whole mesh (depending on iflag_evap_prec)
+      ! does not reach saturation. In this case, we 
+      ! redistribute zqev0 conserving the ratio liquid/ice 
+
+      IF (zqevt+zqevti.GT.zqev0) THEN
+        zqev=zqev0*zqevt/(zqevt+zqevti)
+        zqevi=zqev0*zqevti/(zqevt+zqevti)
+      ELSE
+        zqev=zqevt
+        zqevi=zqevti
+      ENDIF
+
+      !--Diagnostics
+      dqreva(i) = - zqev / dtime
+      dqssub(i) = - zqevti / dtime
+
+      ! New solid and liquid precipitation fluxes after evap and sublimation
+      zrfln(i) = Max(0.,zrfl(i) - zqev*(paprsdn(i)-paprsup(i))/RG/dtime)
+      zifln(i) = Max(0.,zifl(i) - zqevi*(paprsdn(i)-paprsup(i))/RG/dtime)
+      
+
+      ! vapor, temperature, precip fluxes update
+      ! vapor is updated after evaporation/sublimation (it is increased)
+      zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
+        * (RG/(paprsdn(i)-paprsup(i)))*dtime
+      ! zmqc is the total condensed water in the precip flux (it is decreased)
+      zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
+        * (RG/(paprsdn(i)-paprsup(i)))*dtime
+      ! air and precip temperature (i.e., gridbox temperature)
+      ! is updated due to latent heat cooling
+      zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
+        * (RG/(paprsdn(i)-paprsup(i)))*dtime    &
+        * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
+        + (zifln(i)-zifl(i))                      &
+        * (RG/(paprsdn(i)-paprsup(i)))*dtime    &
+        * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+
+      ! New values of liquid and solid precipitation 
+      zrfl(i) = zrfln(i)
+      zifl(i) = zifln(i)
+
+      ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
+      !                         due to evap + sublim                                     
+                                 
+
+      IF (iflag_evap_prec.GE.4) THEN
+        zrflclr(i) = zrfl(i)
+        ziflclr(i) = zifl(i)
+        IF(zrflclr(i) + ziflclr(i).LE.0) THEN
+          znebprecipclr(i) = 0.0
+        ENDIF
+        zrfl(i) = zrflclr(i) + zrflcld(i)
+        zifl(i) = ziflclr(i) + ziflcld(i)
+      ENDIF
+
+      ! c_iso duplicate for isotopes or loop on isotopes
+
+      ! Melting:
+      zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
+      ! precip fraction that is melted
+      zmelt = MIN(MAX(zmelt,0.),1.)
+
+      ! update of rainfall and snowfall due to melting
+      IF (iflag_evap_prec.GE.4) THEN
+        zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
+        zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
+        zrfl(i)=zrflclr(i)+zrflcld(i)
+      ELSE
+        zrfl(i)=zrfl(i)+zmelt*zifl(i)
+      ENDIF
+
+
+      ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
+
+      ! Latent heat of melting because of precipitation melting
+      ! NB: the air + precip temperature is simultaneously updated
+      zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprsdn(i)-paprsup(i)) &
+        *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
+      
+      IF (iflag_evap_prec.GE.4) THEN
+        ziflclr(i)=ziflclr(i)*(1.-zmelt)
+        ziflcld(i)=ziflcld(i)*(1.-zmelt)
+        zifl(i)=ziflclr(i)+ziflcld(i)
+      ELSE
+        zifl(i)=zifl(i)*(1.-zmelt)
+      ENDIF
+
+    ELSE
+      ! if no precip, we reinitialize the cloud fraction used for the precip to 0
+      ! AB note that this assignment is useless, as znebprecip is not re-used
+      znebprecip(i)=0.
+
+    ENDIF ! (zrfl(i)+zifl(i).GT.0.)
+
+  ENDDO ! loop on klon
+
+ENDIF ! (iflag_evap_prec>=1)
+
+END SUBROUTINE histprecip_precld
+
+
+!----------------------------------------------------------------
+! Computes the processes-oriented precipitation formulations for
+! evaporation and sublimation
+! 
+SUBROUTINE poprecip_precld( &
+           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
+           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
+           rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
+           )
+
+USE lmdz_lscp_ini, ONLY : prt_level, lunout
+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_tools, ONLY : calc_qsat_ecmwf
+
+IMPLICIT NONE
+
+
+INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
+REAL,    INTENT(IN)                     :: dtime          !--time step [s]
+LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
+
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
+REAL,    INTENT(IN),    DIMENSION(klon) :: tempupnew      !--updated temperature of the overlying layer [K]
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qprecip        !--specific humidity in the precipitation falling from the upper layer [kg/kg]
+
+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 [-]
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
+REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [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) :: 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(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
+
+
+!--Integer for interating over klon
+INTEGER :: i
+!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
+REAL, DIMENSION(klon) :: dhum_to_dflux
+!--
+REAL, DIMENSION(klon) :: rho, dz
+
+!--Saturation values
+REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
+!--Vapor in the clear sky
+REAL :: qvapclr
+!--Fluxes tendencies because of evaporation and sublimation
+REAL :: dprecip_evasub_max, draineva, dsnowsub, dprecip_evasub_tot
+!--Specific humidity tendencies because of evaporation and sublimation
+REAL :: dqrevap, dqssubl
+!--Specific heat constant
+REAL :: cpair, cpw
+
+!--Initialisation
+qzero(:)  = 0.
+dqreva(:) = 0.
+dqssub(:) = 0.
+dqrevap   = 0.
+dqssubl   = 0.
+
+!-- 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
+!--depending on temperature:
+CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:))
+!--wrt liquid water
+CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
+!--wrt ice
+CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
+
+
+
+!--First step consists in "thermalizing" the layer:
+!--as the flux of precip from layer above "advects" some heat (as the precip is at the temperature 
+!--of the overlying layer) we recalculate a mean temperature that both the air and the precip in the 
+!--layer have.
+
+IF (iftop) THEN
+
+  DO i = 1, klon
+    qprecip(i) = 0.
+  ENDDO
+
+ELSE
+
+  DO i = 1, klon
+    !--No condensed water so cp=cp(vapor+dry air)
+    !-- RVTMP2=rcpv/rcpd-1
+    cpair = RCPD * ( 1. + RVTMP2 * qvap(i) )
+    cpw = RCPD * RVTMP2
+    !--qprecip has to be thermalized with 
+    !--layer's air so that precipitation at the ground has the
+    !--same temperature as the lowermost layer
+    !--we convert the flux into a specific quantity qprecip
+    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) ) &
+            / ( cpair + qprecip(i) * cpw )
+  ENDDO
+
+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.
+
+
+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
+      !--the clear and the cloudy sky as in the layer above. This
+      !--extends the assumption that the cloud fraction is the same
+      !--as the layer above. This is assumed only for the evap / subl
+      !--process
+      !--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) )
+    ELSE
+      !--Legacy version from Ludo - we use the total specific humidity
+      !--for the evap / subl process
+      qvapclr = qvap(i)
+    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
+
+
+    !--New solid and liquid precipitation fluxes after evap and sublimation
+    dqrevap = draineva / dhum_to_dflux(i)
+    dqssubl = dsnowsub / dhum_to_dflux(i)
+
+
+    !--Vapor is updated after evaporation/sublimation (it is increased)
+    qvap(i) = qvap(i) - dqrevap - dqssubl
+    !--qprecip is the total condensed water in the precip flux (it is decreased)
+    qprecip(i) = qprecip(i) + dqrevap + dqssubl
+    !--Air and precip temperature (i.e., gridbox temperature)
+    !--is updated due to latent heat cooling
+    temp(i) = temp(i) &
+            + dqrevap * RLVTT / RCPD &
+            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
+            + dqssubl * RLSTT / RCPD &
+            / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
+
+    !--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
+    IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
+
+    !--Calculation of the total fluxes
+    rain(i) = rainclr(i) + raincld(i)
+    snow(i) = snowclr(i) + snowcld(i)
+
+  ELSEIF ( ( rain(i) + snow(i) ) .LE. 0. ) THEN
+    !--If no precip, we reinitialize the cloud fraction used for the precip to 0 
+    precipfraccld(i) = 0.
+    precipfracclr(i) = 0.
+
+  ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
+
+  !--Diagnostic tendencies
+  dqssub(i) = dqssubl / dtime
+  dqreva(i) = dqrevap / dtime
+
+ENDDO ! loop on klon
+
+END SUBROUTINE poprecip_precld
+
+
+!----------------------------------------------------------------
+! Computes the processes-oriented precipitation formulations for
+! - autoconversion (auto) via a deposition process
+! - aggregation (agg)
+! - riming (rim)
+! - collection (col)
+! - melting (melt)
+! - freezing (freez)
+! 
+SUBROUTINE poprecip_postcld( &
+           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
+           temp, qvap, qliq, qice, icefrac, cldfra, &
+           precipfracclr, precipfraccld, &
+           rain, rainclr, raincld, snow, snowclr, snowcld, &
+           qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
+           dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
+
+USE lmdz_lscp_ini, ONLY : prt_level, lunout
+USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI
+USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
+
+USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
+                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & 
+                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
+                          rho_rain, r_rain, r_snow, rho_ice,                   &
+                          tau_auto_snow_min, tau_auto_snow_max,                &
+                          thresh_precip_frac, eps,                             &
+                          gamma_melt, alpha_freez, beta_freez, temp_nowater,   &
+                          iflag_cloudth_vert, iflag_rain_incloud_vol,          &
+                          cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez,       &
+                          rain_fallspeed_clr, rain_fallspeed_cld,              &
+                          snow_fallspeed_clr, snow_fallspeed_cld
+
+
+IMPLICIT NONE
+
+INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
+REAL,    INTENT(IN)                     :: dtime          !--time step [s]
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
+LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity [kg/kg]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--current liquid water specific humidity [kg/kg]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--current ice water specific humidity [kg/kg]
+REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
+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) :: 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(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(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) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsrim         !--snow tendency due to riming [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
+
+
+
+!--Local variables
+
+INTEGER :: i
+!--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation
+REAL, DIMENSION(klon) :: dhum_to_dflux
+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
+REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp
+REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq
+REAL :: rho_snow
+REAL :: dqlcol           !--loss of liquid cloud content due to collection by rain [kg/kg/s]
+REAL :: dqiagg           !--loss of ice cloud content due to collection by aggregation [kg/kg/s]
+REAL :: dqlrim           !--loss of liquid cloud content due to riming on snow [kg/kg/s]
+
+!--Autoconversion
+REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
+REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
+REAL :: dqlauto          !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
+REAL :: dqiauto          !--loss of ice cloud content due to autoconversion to snow [kg/kg/s]
+
+!--Melting
+REAL :: dqsmelt_max, air_thermal_conduct
+REAL :: nb_snowflake_clr, nb_snowflake_cld
+REAL :: capa_snowflake, temp_wetbulb
+REAL :: rho, r_ice
+REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
+REAL, DIMENSION(klon) :: qzero, qsat, dqsat
+
+!--Freezing
+REAL :: dqrfreez_max
+REAL :: tau_freez
+REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
+REAL :: coef_freez
+REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
+REAL :: Eff_rain_ice
+
+
+!--Initialisation of variables
+
+
+qzero(:)    = 0.
+
+dqrcol(:)   = 0.
+dqsagg(:)   = 0.
+dqsauto(:)  = 0.
+dqrauto(:)  = 0.
+dqsrim(:)   = 0.
+dqrmelt(:)  = 0.
+dqsmelt(:)  = 0.
+dqrfreez(:) = 0.
+dqsfreez(:) = 0.
+
+
+DO i = 1, klon
+
+  !--Variables initialisation
+  dqlcol  = 0.
+  dqiagg  = 0.
+  dqiauto = 0.
+  dqlauto = 0.
+  dqlrim  = 0.
+
+  !-- 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) ) / 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
+  !--the "true" volume fraction instead of a modified 
+  !--surface fraction (which is larger and artificially
+  !--reduces the in-cloud water).
+  IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
+    eff_cldfra = ctot_vol(i)
+  ELSE
+    eff_cldfra = cldfra(i)
+  ENDIF 
+
+
+  !--Start precipitation growth processes
+
+  !--If the cloud is big enough, the precipitation processes activate
+  ! TODO met on seuil_neb ici ?
+  IF ( cldfra(i) .GE. seuil_neb ) THEN
+
+    !---------------------------------------------------------
+    !--           COLLECTION AND AGGREGATION
+    !---------------------------------------------------------
+    !--Collection: processus through which rain collects small liquid droplets
+    !--in suspension, and add it to the rain flux
+    !--Aggregation: same for snow (precip flux) and ice crystals (in suspension)
+    !--Those processes are treated before autoconversion because we do not
+    !--want to collect/aggregate the newly formed fluxes, which already
+    !--"saw" the cloud as they come from it
+    !--The formulas come from Muench and Lohmann 2020
+
+    !--gamma_col: tuning coefficient [-]
+    !--rho_rain: volumic mass of rain [kg/m3]
+    !--r_rain: size of the rain droplets [m]
+    !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1)
+    !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are
+    !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to
+    !--get in-cloud mean quantities. The two divisions by eff_cldfra are
+    !--then simplified.
+
+    !--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. ) THEN
+      !--Exact explicit version, which does not need a barrier because of
+      !--the exponential decrease
+      dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. )
+
+      !--Add tendencies
+      qliq(i) = qliq(i) + dqlcol
+      raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i)
+
+      !--Diagnostic tendencies
+      dqrcol(i)  = - dqlcol  / dtime
+    ENDIF
+
+    !--Same as for aggregation
+    !--Eff_snow_liq formula:
+    !--it s a product of a collection efficiency and a sticking efficiency
+    ! Milbrandt and Yau formula that gives very low values:
+    ! Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) )
+    ! Lin 1983's formula 
+    Eff_snow_ice = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) ) 
+    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
+    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
+    coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice
+    IF ( snowcld(i) .GT. 0. ) THEN
+      !--Exact explicit version, which does not need a barrier because of
+      !--the exponential decrease
+      dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. )
+
+      !--Add tendencies
+      qice(i) = qice(i) + dqiagg
+      snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i)
+
+      !--Diagnostic tendencies
+      dqsagg(i)  = - dqiagg  / dtime
+    ENDIF
+
+
+    !---------------------------------------------------------
+    !--                  AUTOCONVERSION
+    !---------------------------------------------------------
+    !--Autoconversion converts liquid droplets/ice crystals into
+    !--rain drops/snowflakes. It relies on the formulations by
+    !--Sundqvist 1978.
+
+    !--If we are in a convective point, we have different parameters
+    !--for the autoconversion
+    IF ( ptconv(i) ) THEN
+        qthresh_auto_rain = cld_lc_con
+        qthresh_auto_snow = cld_lc_con_snow
+
+        tau_auto_rain = cld_tau_con
+        !--tau for snow depends on the ice fraction in mixed-phase clouds     
+        tau_auto_snow = tau_auto_snow_max &
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+
+        expo_auto_rain = cld_expo_con
+        expo_auto_snow = cld_expo_con
+    ELSE
+        qthresh_auto_rain = cld_lc_lsc
+        qthresh_auto_snow = cld_lc_lsc_snow
+
+        tau_auto_rain = cld_tau_lsc
+        !--tau for snow depends on the ice fraction in mixed-phase clouds     
+        tau_auto_snow = tau_auto_snow_max &
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+
+        expo_auto_rain = cld_expo_lsc
+        expo_auto_snow = cld_expo_lsc
+    ENDIF
+
+
+    ! Liquid water quantity to remove according to (Sundqvist, 1978)
+    ! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) )
+    !
+    !--And same formula for ice
+    !
+    !--We first treat the second term (with exponential) in an explicit way
+    !--and then treat the first term (-q/tau) in an exact way
+
+    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
+               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
+
+    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
+               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
+
+
+    !--Barriers so that we don't create more rain/snow
+    !--than there is liquid/ice
+    dqlauto = MAX( - qliq(i), dqlauto )
+    dqiauto = MAX( - qice(i), dqiauto )
+
+    !--Add tendencies
+    qliq(i) = qliq(i) + dqlauto
+    qice(i) = qice(i) + dqiauto
+    raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i)
+    snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i)
+
+    !--Diagnostic tendencies
+    dqsauto(i) = - dqiauto / dtime
+    dqrauto(i) = - dqlauto / dtime
+
+
+    !---------------------------------------------------------
+    !--                  RIMING
+    !---------------------------------------------------------
+    !--Process which converts liquid droplets in suspension into
+    !--snow because of the collision between
+    !--those and falling snowflakes.
+    !--The formula comes from Muench and Lohmann 2020
+    !--NB.: this process needs a temperature adjustment
+
+    !--Eff_snow_liq formula: following Ferrier 1994,
+    !--assuming 1
+    Eff_snow_liq = 1.0
+    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
+    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
+    coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq
+    IF ( snowcld(i) .GT. 0. ) THEN
+      !--Exact version, which does not need a barrier because of
+      !--the exponential decrease
+      dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. )
+
+      !--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) = snowcld(i) - dqlrim * dhum_to_dflux(i)
+
+      !--Temperature adjustment with the release of latent
+      !--heat because of solid condensation
+      temp(i) = temp(i) - dqlrim * RLMLT / RCPD &
+                        / ( 1. + RVTMP2 * qtot(i) )
+
+      !--Diagnostic tendencies
+      dqsrim(i)  = - dqlrim  / dtime
+    ENDIF
+
+  ENDIF ! cldfra .GE. seuil_neb
+
+ENDDO ! loop on klon
+
+
+!--Re-calculation of saturation specific humidity
+!--because riming changed temperature
+CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat)
+
+DO i = 1, klon
+
+  !---------------------------------------------------------
+  !--                  MELTING
+  !---------------------------------------------------------
+  !--Process through which snow melts into rain. 
+  !--The formula is homemade.
+  !--NB.: this process needs a temperature adjustment
+
+  !--dqsmelt_max         : maximum snow melting so that temperature
+  !--                      stays higher than 273 K [kg/kg]
+  !--capa_snowflake      : capacitance of a snowflake, equal to
+  !--                      the radius if the snowflake is a sphere [m]
+  !--temp_wetbulb        : wet-bulb temperature [K]
+  !--snow_fallspeed      : snow fall velocity (in clear/cloudy sky) [m/s]
+  !--air_thermal_conduct : thermal conductivity of the air [J/m/K/s]
+  !--gamma_melt          : tuning parameter for melting [-]
+  !--nb_snowflake        : number of snowflakes (in clear/cloudy air) [-]
+
+  IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
+    !--Computed according to
+    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
+    dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD &
+                        * ( 1. + RVTMP2 * qtot(i) ))
+   
+    !--Initialisation 
+    dqsclrmelt = 0.
+    dqscldmelt = 0.
+
+    !--We assume that the snowflakes are spherical
+    capa_snowflake = r_snow
+    !--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
+    air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184    
+    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
+    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
+
+    !--In clear air
+    IF ( ( snowclr(i) .GT. 0. ) .AND.  ( precipfracclr(i) .GT. 0. ) ) THEN
+      !--Formula for the wet-bulb temperature from ECMWF (IFS)
+      !--The vapor used is the vapor in the clear sky
+      temp_wetbulb = temp(i) &
+                   - ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
+                   * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
+                   - 40.637 * ( temp(i) - 275. ) )
+      !--Calculated according to
+      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
+      nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr &
+                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
+      dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct &
+                   * capa_snowflake / RLMLT * gamma_melt &
+                   * MAX(0., temp_wetbulb - RTT) * dtime
+      
+      !--Barrier to limit the melting flux to the clr snow flux in the mesh
+      dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i))
+    ENDIF
+
+
+    !--In cloudy air
+    IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
+      !--As the air is saturated, the wet-bulb temperature is equal to the
+      !--temperature
+      temp_wetbulb = temp(i)
+      !--Calculated according to
+      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
+      nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld &
+                       / ( 4. / 3. * RPI * r_snow**3. * rho_snow )
+      dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct &
+                   * capa_snowflake / RLMLT * gamma_melt &
+                   * MAX(0., temp_wetbulb - RTT) * dtime
+
+      !--Barrier to limit the melting flux to the cld snow flux in the mesh
+      dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i))
+   ENDIF
+    
+
+    !--Barrier on temperature. If the total melting flux leads to a
+    !--positive temperature, it is limited to keep temperature above 0 degC.
+    !--It is activated if the total is LOWER than the max
+    !--because everything is negative
+    dqstotmelt = dqsclrmelt + dqscldmelt
+    IF ( dqstotmelt .LT. dqsmelt_max ) THEN
+      !--We redistribute the max melted snow keeping
+      !--the clear/cloud partition of the melted snow
+      dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt
+      dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
+      dqstotmelt = dqsmelt_max
+
+    ENDIF
+
+    !--Add tendencies
+    !--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))
+
+    !--Temperature adjustment with the release of latent
+    !--heat because of melting
+    temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD &
+                      / ( 1. + RVTMP2 * qtot(i) )
+
+    !--Diagnostic tendencies
+    dqrmelt(i) = - dqstotmelt / dtime
+    dqsmelt(i) = dqstotmelt / dtime
+
+  ENDIF
+
+
+  !---------------------------------------------------------
+  !--                  FREEZING
+  !---------------------------------------------------------
+  !--Process through which rain freezes into snow. 
+  !-- We parameterize it as a 2 step process:
+  !--first: freezing following collision with ice crystals
+  !--second: immersion freezing following (inspired by Bigg 1953)
+  !--the latter is parameterized as an exponential decrease of the rain
+  !--water content with a homemade formulya
+  !--This is based on a caracteritic time of freezing, which
+  !--exponentially depends on temperature so that it is
+  !--equal to 1 for temp_nowater (see below) and is close to
+  !--0 for RTT (=273.15 K).
+  !--NB.: this process needs a temperature adjustment
+  !--dqrfreez_max : maximum rain freezing so that temperature
+  !--              stays lower than 273 K [kg/kg]
+  !--tau_freez    : caracteristic time of freezing [s]
+  !--gamma_freez  : tuning parameter [s-1]
+  !--alpha_freez  : tuning parameter for the shape of the exponential curve [-]
+  !--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC)
+
+  IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
+
+ 
+    !--1st step: freezing following collision with ice crystals
+    !--Sub-process of freezing which quantifies the collision between
+    !--ice crystals in suspension and falling rain droplets.
+    !--The rain droplets freeze, becoming graupel, and carrying
+    !--the ice crystal (which acted as an ice nucleating particle).
+    !--The formula is adapted from the riming formula.
+    !--it works only in the cloudy part
+
+    dqifreez = 0.
+    dqrtotfreez_step1 = 0.
+
+    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
+         ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
+      dqrclrfreez = 0.
+      dqrcldfreez = 0.
+
+      !--Computed according to
+      !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
+      dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
+                         * ( 1. + RVTMP2 * qtot(i) ))
+ 
+
+      !--The collision efficiency is assumed unity 
+      Eff_rain_ice = 1.
+      coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
+      !--Exact version, which does not need a barrier because of
+      !--the exponential decrease.
+      dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
+
+      !--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)
+      !--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)
+      !--Dry density [kg/m3]
+      rho = pplay(i) / temp(i) / RD
+      !--r_ice formula from Sun and Rikus (1999)
+      r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
+            + 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
+      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 * 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) )
+
+    ENDIF
+    
+    !-- Second step immersion freezing of rain drops
+    !-- with a homemade timeconstant depending on temperature
+
+    dqrclrfreez = 0.
+    dqrcldfreez = 0.
+    dqrtotfreez_step2 = 0.
+    !--Computed according to
+    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
+
+    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
+                         * ( 1. + RVTMP2 * qtot(i) ))
+
+    
+    tau_freez = 1. / ( beta_freez &
+              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
+
+
+    !--In clear air
+    IF ( rainclr(i) .GT. 0. ) THEN
+      !--Exact solution of dqrain/dt = -qrain/tau_freez
+      dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
+    ENDIF
+
+    !--In cloudy air
+    IF ( raincld(i) .GT. 0. ) THEN
+      !--Exact solution of dqrain/dt = -qrain/tau_freez
+      dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. )
+    ENDIF
+
+    !--temperature barrier step 2
+    !--It is activated if the total is LOWER than the max
+    !--because everything is negative
+    dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
+  
+    IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
+      !--We redistribute the max freezed rain keeping
+      !--the clear/cloud partition of the freezing rain
+      dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
+      dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
+      dqrtotfreez_step2 = dqrfreez_max
+    ENDIF
+
+
+    !--Add tendencies
+    !--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))
+
+
+    !--Temperature adjustment with the uptake of latent
+    !--heat because of freezing
+    temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
+                      / ( 1. + RVTMP2 * qtot(i) )
+
+    !--Diagnostic tendencies
+    dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2          
+    dqrfreez(i) = dqrtotfreez / dtime
+    dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
+
+  ENDIF
+
+
+
+  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
+  !--we reduce the precipiration fraction in the clear/cloudy air so that the new
+  !--local flux of rain+snow is equal to rain_int_min.
+  !--Here, rain+snow is the gridbox-mean flux of precip.
+  !--Therefore, (rain+snow)/precipfrac is the local flux of precip.
+  !--If the local flux of precip is lower than rain_int_min, i.e.,
+  !-- (rain+snow)/precipfrac < rain_int_min , i.e., 
+  !-- (rain+snow)/rain_int_min < precipfrac , then we want to reduce
+  !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/rain_int_min.
+  !--Note that this is physically different than what is proposed in LTP thesis.
+  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
+  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
+
+  !--Calculate outputs
+  rain(i) = rainclr(i) + raincld(i)
+  snow(i) = snowclr(i) + snowcld(i)
+
+  !--Diagnostics
+  !--BEWARE this is indeed a diagnostic: this is an estimation from
+  !--the value of the flux at the bottom interface of the mesh and
+  !--and assuming an upstream numerical calculation
+  !--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is
+  !--used for computing the total ice water content in the mesh
+  !--for radiation only
+  qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr &
+                 + raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld )
+  qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr &
+                 + snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld )
+
+
+ENDDO ! loop on klon
+
+END SUBROUTINE poprecip_postcld
+
+END MODULE lmdz_lscp_precip
