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, znebprecip, 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) :: znebprecip     !--fraction of precipitation 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) :: 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.
         END DO

      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)

         END DO

      END IF

! --------------------------------------------------------------------
! 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)
               END IF

               IF (iflag_evap_prec .EQ. 1) THEN
                  znebprecip(i) = zneb(i)
               ELSE
                  znebprecip(i) = MAX(zneb(i), znebprecip(i))
               END IF

               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))
               END IF

               ! 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
               END IF

               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
               END IF

               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
               END IF

               !--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
                  END IF
                  zrfl(i) = zrflclr(i) + zrflcld(i)
                  zifl(i) = ziflclr(i) + ziflcld(i)
               END IF

               ! 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)
               END IF

               ! 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)
               END IF

            ELSE
               ! if no precip, we reinitialize the cloud fraction used for the precip to 0
               znebprecip(i) = 0.

            END IF ! (zrfl(i)+zifl(i).GT.0.)

         END DO ! loop on klon

      END IF ! (iflag_evap_prec>=1)

   END SUBROUTINE histprecip_precld

   SUBROUTINE histprecip_postcld( &
      klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, zdqsdT_raw, &
      zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
      rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
      zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
      zradocond, zradoice, dqrauto, dqsauto &
      )

      USE lmdz_lscp_ini, ONLY: RD, RG, RCPD, RVTMP2, RLSTT, RLMLT, RTT
      USE lmdz_lscp_ini, ONLY: cld_lc_con, cld_tau_con, cld_expo_con, ffallv_con, &
                               cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, ffallv_lsc, &
                               seuil_neb, rain_int_min, cice_velo, dice_velo
      USE lmdz_lscp_ini, ONLY: iflag_evap_prec, iflag_cloudth_vert, iflag_rain_incloud_vol, &
                               iflag_autoconversion, ok_radocond_snow, ok_bug_phase_lscp, &
                               niter_lscp
      USE lmdz_lscp_tools, ONLY: fallice_velocity

      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(IN), DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
      LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv         !--true if we are in a convective point
      REAL, INTENT(IN), DIMENSION(klon) :: zdqsdT_raw     !--derivative of qsat w.r.t. temperature times L/Cp [SI]

      REAL, INTENT(INOUT), DIMENSION(klon) :: zt             !--current temperature [K]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zq             !--current water vapor specific humidity [kg/kg]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zoliq          !--current liquid+ice water specific humidity [kg/kg]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zoliql         !--current liquid water specific humidity [kg/kg]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zoliqi         !--current ice water specific humidity [kg/kg]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zcond          !--liquid+ice water specific humidity AFTER CONDENSATION (no precip process) [kg/kg]
      REAL, INTENT(IN), DIMENSION(klon) :: zfice          !--ice fraction AFTER CONDENSATION [-]
      REAL, INTENT(IN), DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]

      REAL, INTENT(IN), DIMENSION(klon) :: rneb           !--cloud fraction [-]
      REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
      REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipcld  !--fraction of precipitation in the cloud IN THE LAYER ABOVE [-]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
      REAL, INTENT(INOUT), DIMENSION(klon) :: tot_zneb       !--total cloud cover above the current mesh [-]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zrho_up        !--air density IN THE LAYER ABOVE [kg/m3]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zvelo_up       !--ice fallspeed velocity IN THE LAYER ABOVE [m/s]

      REAL, INTENT(INOUT), DIMENSION(klon) :: zrfl           !--flux of rain gridbox-mean [kg/s/m2]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zrflclr        !--flux of rain gridbox-mean in clear sky [kg/s/m2]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zrflcld        !--flux of rain gridbox-mean in cloudy air [kg/s/m2]
      REAL, INTENT(INOUT), DIMENSION(klon) :: zifl           !--flux of snow gridbox-mean [kg/s/m2]
      REAL, INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky [kg/s/m2]
      REAL, INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air [kg/s/m2]

      REAL, INTENT(OUT), DIMENSION(klon) :: zradocond      !--condensed water used in the radiation scheme [kg/kg]
      REAL, INTENT(OUT), DIMENSION(klon) :: zradoice       !--condensed ice water used in the radiation scheme [kg/kg]
      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]

! Local variables for precip fraction update
      REAL :: smallestreal
      REAL, DIMENSION(klon) :: tot_znebn, d_tot_zneb
      REAL, DIMENSION(klon) :: d_znebprecip_cld_clr, d_znebprecip_clr_cld
      REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
      REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld

! Local variables for autoconversion
      REAL :: zct, zcl, zexpo, ffallv
      REAL :: zchau, zfroi
      REAL :: zrain, zsnow, zprecip
      REAL :: effective_zneb
      REAL, DIMENSION(klon) :: zrho, zvelo
      REAL, DIMENSION(klon) :: zdz, iwc

! Local variables for Bergeron process
      REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
      REAL, DIMENSION(klon) :: zqpreci, zqprecl

! Local variables for calculation of quantity of condensates seen by the radiative scheme
      REAL, DIMENSION(klon) :: zradoliq
      REAL, DIMENSION(klon) :: ziflprev

! Misc
      INTEGER :: i, n

! Initialisation
      smallestreal = 1.e-9
      zqprecl(:) = 0.
      zqpreci(:) = 0.
      ziflprev(:) = 0.

      IF (iflag_evap_prec .GE. 4) THEN

         !Partitionning between precipitation coming from clouds and that coming from CS

         !0) Calculate tot_zneb, total cloud fraction above the cloud
         !assuming a maximum-random overlap (voir Jakob and Klein, 2000)

         DO i = 1, klon
            tot_znebn(i) = 1.-(1.-tot_zneb(i))*(1 - max(rneb(i), zneb(i))) &
                           /(1.-min(zneb(i), 1.-smallestreal))
            d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
            tot_zneb(i) = tot_znebn(i)

            !1) Cloudy to clear air
            d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i), znebprecipcld(i))
            IF (znebprecipcld(i) .GT. 0.) THEN
               d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
               d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
            ELSE
               d_zrfl_cld_clr(i) = 0.
               d_zifl_cld_clr(i) = 0.
            END IF

            !2) Clear to cloudy air
            d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i) - d_tot_zneb(i) - zneb(i)))
            IF (znebprecipclr(i) .GT. 0) THEN
               d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
               d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
            ELSE
               d_zrfl_clr_cld(i) = 0.
               d_zifl_clr_cld(i) = 0.
            END IF

            !Update variables
            znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
            znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
            zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
            ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
            zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
            ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)

            ! c_iso  do the same thing for isotopes precip
         END DO
      END IF

! Autoconversion
! -------------------------------------------------------------------------------

! Initialisation
      DO i = 1, klon
         zrho(i) = pplay(i)/zt(i)/RD
         zdz(i) = (paprsdn(i) - paprsup(i))/(zrho(i)*RG)
         iwc(i) = 0.
         zneb(i) = MAX(rneb(i), seuil_neb)

         ! variables for calculation of quantity of condensates seen by the radiative scheme
         ! NB. only zradocond and zradoice are outputed, but zradoliq is used if ok_radocond_snow
         ! is TRUE
         zradocond(i) = zoliq(i)/REAL(niter_lscp + 1)
         zradoliq(i) = zoliq(i)*(1.-zfice(i))/REAL(niter_lscp + 1)
         zradoice(i) = zoliq(i)*zfice(i)/REAL(niter_lscp + 1)
      END DO

      DO n = 1, niter_lscp

         DO i = 1, klon
            IF (rneb(i) .GT. 0.0) THEN
               iwc(i) = zrho(i)*zoliqi(i)/zneb(i) ! in-cloud ice condensate content
            END IF
         END DO

         CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, zvelo)

         DO i = 1, klon

            IF (rneb(i) .GT. 0.0) THEN

               ! Initialization of zrain, zsnow and zprecip:
               zrain = 0.
               zsnow = 0.
               zprecip = 0.
               ! c_iso same init for isotopes. Externalisation?

               IF (zneb(i) .EQ. seuil_neb) THEN
                  zprecip = 0.0
                  zsnow = 0.0
                  zrain = 0.0
               ELSE

                  IF (ptconv(i)) THEN ! if convective point
                     zcl = cld_lc_con
                     zct = 1./cld_tau_con
                     zexpo = cld_expo_con
                     ffallv = ffallv_con
                  ELSE
                     zcl = cld_lc_lsc
                     zct = 1./cld_tau_lsc
                     zexpo = cld_expo_lsc
                     ffallv = ffallv_lsc
                  END IF

                  ! 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
                     effective_zneb = ctot_vol(i)
                  ELSE
                     effective_zneb = zneb(i)
                  END IF

                  ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
                  ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
                  !.........................................................
                  IF (iflag_autoconversion .EQ. 2) THEN
                     ! two-steps resolution with niter_lscp=1 sufficient
                     ! we first treat the second term (with exponential) in an explicit way
                     ! and then treat the first term (-q/tau) in an exact way
                     zchau = zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct &
                                               *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo))))
                  ELSE
                     ! old explicit resolution with subtimesteps
                     zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) &
                             *(1.0 - EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo))
                  END IF

                  ! Ice water quantity to remove (Zender & Kiehl, 1997)
                  ! dqice/dt=1/rho*d(rho*wice*qice)/dz
                  !....................................
                  IF (iflag_autoconversion .EQ. 2) THEN
                     ! exact resolution, niter_lscp=1 is sufficient but works only
                     ! with iflag_vice=0
                     IF (zoliqi(i) .GT. 0.) THEN
                        zfroi = (zoliqi(i) - ((zoliqi(i)**(-dice_velo)) &
                                       + dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)/zneb(i)*zrho(i)*ffallv)**(-1./dice_velo))
                     ELSE
                        zfroi = 0.
                     END IF
                  ELSE
                     ! old explicit resolution with subtimesteps
                     zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*zvelo(i)
                  END IF

                  zrain = MIN(MAX(zchau, 0.0), zoliql(i))
                  zsnow = MIN(MAX(zfroi, 0.0), zoliqi(i))
                  zprecip = MAX(zrain + zsnow, 0.0)

               END IF

               IF (iflag_autoconversion .GE. 1) THEN
                  ! debugged version with phase conservation through the autoconversion process
                  zoliql(i) = MAX(zoliql(i) - 1.*zrain, 0.0)
                  zoliqi(i) = MAX(zoliqi(i) - 1.*zsnow, 0.0)
                  zoliq(i) = MAX(zoliq(i) - zprecip, 0.0)
               ELSE
                  ! bugged version with phase resetting
                  zoliql(i) = MAX(zoliq(i)*(1.-zfice(i)) - 1.*zrain, 0.0)
                  zoliqi(i) = MAX(zoliq(i)*zfice(i) - 1.*zsnow, 0.0)
                  zoliq(i) = MAX(zoliq(i) - zprecip, 0.0)
               END IF

               ! c_iso: call isotope_conversion (for readibility) or duplicate

               ! variables for calculation of quantity of condensates seen by the radiative scheme
               zradocond(i) = zradocond(i) + zoliq(i)/REAL(niter_lscp + 1)
               zradoliq(i) = zradoliq(i) + zoliql(i)/REAL(niter_lscp + 1)
               zradoice(i) = zradoice(i) + zoliqi(i)/REAL(niter_lscp + 1)

               !--Diagnostics
               dqrauto(i) = dqrauto(i) + zrain/dtime
               dqsauto(i) = dqsauto(i) + zsnow/dtime

            END IF ! rneb >0

         END DO  ! i = 1,klon

      END DO      ! n = 1,niter

! Precipitation flux calculation

      DO i = 1, klon

         IF (iflag_evap_prec .GE. 4) THEN
            ziflprev(i) = ziflcld(i)
         ELSE
            ziflprev(i) = zifl(i)*zneb(i)
         END IF

         IF (rneb(i) .GT. 0.0) THEN

            ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
            ! If T<0C, liquid precip are converted into ice, which leads to an increase in
            ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
            ! taken into account through a linearization of the equations and by approximating
            ! the liquid precipitation process with a "threshold" process. We assume that
            ! condensates are not modified during this operation. Liquid precipitation is
            ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
            ! Water vapor increases as well
            ! Note that compared to fisrtilp, we always assume iflag_bergeron=2

            IF (ok_bug_phase_lscp) THEN
               zqpreci(i) = (zcond(i) - zoliq(i))*zfice(i)
               zqprecl(i) = (zcond(i) - zoliq(i))*(1.-zfice(i))
            ELSE
               zqpreci(i) = zcond(i)*zfice(i) - zoliqi(i)
               zqprecl(i) = zcond(i)*(1.-zfice(i)) - zoliql(i)
            END IF
            zcp = RCPD*(1.0 + RVTMP2*(zq(i) + zmqc(i) + zcond(i)))
            coef1 = rneb(i)*RLSTT/zcp*zdqsdT_raw(i)
            ! Computation of DT if all the liquid precip freezes
            DeltaT = RLMLT*zqprecl(i)/(zcp*(1.+coef1))

            ! T should not exceed the freezing point
            ! that is Delta > RTT-zt(i)
            DeltaT = max(min(RTT - zt(i), DeltaT), 0.)
            zt(i) = zt(i) + DeltaT
            ! water vaporization due to temp. increase
            Deltaq = rneb(i)*zdqsdT_raw(i)*DeltaT
            ! we add this vaporized water to the total vapor and we remove it from the precip
            zq(i) = zq(i) + Deltaq
            ! The three "max" lines herebelow protect from rounding errors
            zcond(i) = max(zcond(i) - Deltaq, 0.)
            ! liquid precipitation converted to ice precip
            Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
            zqprecl(i) = max(zqprecl(i) + Deltaqprecl, 0.)
            ! iced water budget
            zqpreci(i) = max(zqpreci(i) - Deltaqprecl - Deltaq, 0.)

            ! c_iso : mv here condensation of isotopes + redispatchage en precip

            IF (iflag_evap_prec .GE. 4) THEN
               zrflcld(i) = zrflcld(i) + zqprecl(i) &
                            *(paprsdn(i) - paprsup(i))/(RG*dtime)
               ziflcld(i) = ziflcld(i) + zqpreci(i) &
                            *(paprsdn(i) - paprsup(i))/(RG*dtime)
               znebprecipcld(i) = rneb(i)
               zrfl(i) = zrflcld(i) + zrflclr(i)
               zifl(i) = ziflcld(i) + ziflclr(i)
            ELSE
               zrfl(i) = zrfl(i) + zqprecl(i) &
                         *(paprsdn(i) - paprsup(i))/(RG*dtime)
               zifl(i) = zifl(i) + zqpreci(i) &
                         *(paprsdn(i) - paprsup(i))/(RG*dtime)
            END IF
            ! c_iso : same for isotopes

         END IF ! rneb>0

      END DO

! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
! if iflag_evap_prec>=4
      IF (iflag_evap_prec .GE. 4) THEN

         DO i = 1, klon

            IF ((zrflclr(i) + ziflclr(i)) .GT. 0.) THEN
               znebprecipclr(i) = min(znebprecipclr(i), max( &
                                      zrflclr(i)/(MAX(znebprecipclr(i), seuil_neb)*rain_int_min), &
                                      ziflclr(i)/(MAX(znebprecipclr(i), seuil_neb)*rain_int_min)))
            ELSE
               znebprecipclr(i) = 0.0
            END IF

            IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
               znebprecipcld(i) = min(znebprecipcld(i), max( &
                                      zrflcld(i)/(MAX(znebprecipcld(i), seuil_neb)*rain_int_min), &
                                      ziflcld(i)/(MAX(znebprecipcld(i), seuil_neb)*rain_int_min)))
            ELSE
               znebprecipcld(i) = 0.0
            END IF
         END DO

      END IF

! we recalculate zradoice to account for contributions from the precipitation flux
! if ok_radocond_snow is true
      IF (ok_radocond_snow) THEN
         IF (.NOT. iftop) THEN
            ! for the solid phase (crystals + snowflakes)
            ! we recalculate zradoice by summing
            ! the ice content calculated in the mesh
            ! + the contribution of the non-evaporated snowfall
            ! from the overlying layer
            DO i = 1, klon
               IF (ziflprev(i) .NE. 0.0) THEN
                  zradoice(i) = zoliqi(i) + zqpreci(i) + ziflprev(i)/zrho_up(i)/zvelo_up(i)
               ELSE
                  zradoice(i) = zoliqi(i) + zqpreci(i)
               END IF
               zradocond(i) = zradoliq(i) + zradoice(i)
            END DO
         END IF
         ! keep in memory air density and ice fallspeed velocity
         zrho_up(:) = zrho(:)
         zvelo_up(:) = zvelo(:)
      END IF

   END SUBROUTINE histprecip_postcld

!----------------------------------------------------------------
! 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, &
      cldfra, rvc_seri, qliq, qice, &
      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, ok_ice_supersat, ok_unadjusted_clouds
      USE lmdz_lscp_ini, ONLY: eps, temp_nowater, ok_growth_precip_deposition
      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(IN), DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
      REAL, INTENT(IN), DIMENSION(klon) :: rvc_seri       !--cloud water vapor at the beginning of lscp (ratio wrt total water) - used only if the cloud properties are advected [kg/kg]
      REAL, INTENT(IN), DIMENSION(klon) :: qliq           !--liquid water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
      REAL, INTENT(IN), DIMENSION(klon) :: qice           !--ice water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]

      REAL, INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
      REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
      REAL, INTENT(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) :: 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
!--Air density [kg/m3] and layer thickness [m]
      REAL, DIMENSION(klon) :: rho, dz
!--Temporary precip fractions and fluxes for the evaporation
      REAL, DIMENSION(klon) :: precipfracclr_tmp, precipfraccld_tmp
      REAL, DIMENSION(klon) :: rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp

!--Saturation values
      REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
!--Vapor in the clear sky and cloud
      REAL :: qvapclr, qvapcld
!--Fluxes tendencies because of evaporation and sublimation
      REAL :: dprecip_evasub_max, dprecip_evasub_tot
      REAL :: drainclreva, draincldeva, dsnowclrsub, dsnowcldsub
!--Specific humidity tendencies because of evaporation and sublimation
      REAL :: dqrevap, dqssubl
!--Specific heat constant
      REAL :: cpair, cpw

!--Initialisation
      qzero(:) = 0.
      dqreva(:) = 0.
      dqssub(:) = 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.
         END DO

      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)
         END DO

      END IF

!--Initialise the precipitation fractions and fluxes
      DO i = 1, klon
         precipfracclr_tmp(i) = precipfracclr(i)
         precipfraccld_tmp(i) = precipfraccld(i)
         rainclr_tmp(i) = rainclr(i)
         snowclr_tmp(i) = snowclr(i)
         raincld_tmp(i) = raincld(i)
         snowcld_tmp(i) = snowcld(i)
      END DO

!--If we relax the LTP assumption that the cloud is the same than in the
!--gridbox above, we can change the tmp quantities
      IF (ok_ice_supersat) THEN
         !--Update the precipitation fraction using advected cloud fraction
         CALL poprecip_fracupdate( &
            klon, cldfra, precipfracclr_tmp, precipfraccld_tmp, &
            rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp)

! here we could put an ELSE statement and do one iteration of the condensation scheme
! (but it would only diagnose cloud from lognormal - link with thermals?)

         DO i = 1, klon
            !--We cannot have a total fraction greater than the one in the gridbox above
            !--because new cloud formation has not yet occured.
            !--If this happens, we reduce the precip frac in the cloud, because it can
            !--only mean that the cloud fraction at this level is greater than the total
            !--precipitation fraction of the level above, and the precip frac in clear sky
            !--is zero.
            !--NB. this only works because we assume a max-random cloud and precip overlap
            precipfraccld_tmp(i) = MIN(precipfraccld_tmp(i), precipfracclr(i) + precipfraccld(i))
         END DO

      END IF
!--At this stage, we guarantee that
!-- precipfracclr + precipfraccld = precipfracclr_tmp + precipfraccld_tmp

      DO i = 1, klon

         dqrevap = 0.
         dqssubl = 0.
         !--If there is precipitation from the layer above
         IF ((rain(i) + snow(i)) .GT. 0.) THEN

            !--Init
            drainclreva = 0.
            draincldeva = 0.
            dsnowclrsub = 0.
            dsnowcldsub = 0.

            !--Init the properties of the air where reevaporation occurs
            IF (ok_ice_supersat) THEN
               !--We can diagnose the water vapor in the cloud using the advected
               !--water vapor in the cloud and cloud fraction
               IF (ok_unadjusted_clouds .AND. (cldfra(i) .GT. eps)) THEN
                  qvapcld = rvc_seri(i)*qvap(i)/cldfra(i)
               END IF
               !--We can diagnose completely the water vapor in clear sky, because all
               !--the needed variables (ice, liq, vapor in cld, cloud fraction) are
               !--advected
               !--Note that qvap(i) is the total water in the gridbox
               IF ((1.-cldfra(i)) .GT. eps) THEN
                  qvapclr = (qvap(i) - qice(i) - qliq(i) - rvc_seri(i)*qvap(i))/(1.-cldfra(i))
               END IF
            ELSEIF (ok_corr_vap_evasub) THEN
               !--Corrected version from Ludo - we use the same water ratio between
               !--the clear and the cloudy sky as in the layer above. This
               !--extends the assumption that the cloud fraction is the same
               !--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
               IF ((1.-precipfraccld(i)) .GT. eps) THEN
                  qvapclr = qvapclrup(i)/qtotupnew(i)*qvap(i)/(1.-precipfraccld(i))
               ELSE
                  qvapclr = 0.
               END IF
               IF (precipfraccld(i) .GT. eps) THEN
                  qvapcld = MAX(qtotupnew(i) - qvapclrup(i), 0.)/qtotupnew(i)*qvap(i)/precipfraccld(i)
               ELSE
                  qvapcld = 0.
               END IF
            ELSE
               !--Legacy version from Ludo - we use the total specific humidity
               !--for the evap / subl process
               qvapclr = qvap(i)
               qvapcld = qvap(i)
            END IF

            !---------------------------
            !-- EVAP / SUBL IN CLEAR SKY
            !---------------------------

            IF (precipfracclr_tmp(i) .GT. eps) THEN

               !--Evaporation of liquid precipitation coming from above
               !--in the clear sky only
               !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
               !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
               !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
               !--which does not need a barrier on rainclr, because included in the formula
               drainclreva = precipfracclr_tmp(i)*MAX(0., &
                                                      -coef_eva*(1.-expo_eva)*(1.-qvapclr/qsatl(i))*dz(i) &
                                                      + (rainclr_tmp(i)/precipfracclr_tmp(i))**(1.-expo_eva) &
                                                      )**(1./(1.-expo_eva)) - rainclr_tmp(i)

               !--Evaporation is limited by 0
               !--NB. with ok_ice_supersat activated, this barrier should be useless
               drainclreva = MIN(0., drainclreva)

               !--Sublimation of the solid precipitation coming from above
               !--(same formula as for liquid precip)
               !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly)
               !--which does not need a barrier on snowclr, because included in the formula
               dsnowclrsub = precipfracclr_tmp(i)*MAX(0., &
                                                      -coef_sub*(1.-expo_sub)*(1.-qvapclr/qsati(i))*dz(i) &
                                                      + (snowclr_tmp(i)/precipfracclr_tmp(i))**(1.-expo_sub) &
                                                      )**(1./(1.-expo_sub)) - snowclr_tmp(i)

               !--If ice supersaturation is activated, we allow vapor to condense on falling snow
               !--i.e., having a positive dsnowclrsub
               IF (.NOT. ok_growth_precip_deposition) THEN
                  !--Sublimation is limited by 0
                  dsnowclrsub = MIN(0., dsnowclrsub)
               END IF

               !--Evaporation limit: we ensure that the layer's fraction below
               !--the clear sky does not reach saturation. In this case, we
               !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
               !--Max evaporation is computed not to saturate the clear sky precip fraction
               !--(i.e., the fraction where evaporation occurs)
               !--It is expressed as a max flux dprecip_evasub_max

               IF (.NOT. ok_ice_supersat) THEN
                  dprecip_evasub_max = MIN(0., (qvapclr - qsat(i))*precipfracclr_tmp(i)) &
                                       *dhum_to_dflux(i)
               ELSE
                  dprecip_evasub_max = (qvapclr - qsat(i))*precipfracclr_tmp(i) &
                                       *dhum_to_dflux(i)
               END IF
               dprecip_evasub_tot = drainclreva + dsnowclrsub

               !--Barriers
               !--If activates if the total is LOWER than the max because
               !--everything is negative
               IF (((dprecip_evasub_max .LT. 0.) .AND. &
                    (dprecip_evasub_tot .LT. dprecip_evasub_max)) .OR. &
                   ((dprecip_evasub_max .GT. 0.) .AND. &
                    (dprecip_evasub_tot .GT. dprecip_evasub_max))) THEN
                  drainclreva = dprecip_evasub_max*drainclreva/dprecip_evasub_tot
                  dsnowclrsub = dprecip_evasub_max*dsnowclrsub/dprecip_evasub_tot
               END IF

            ELSE

               !--All the precipitation is sublimated if the fraction is zero
               drainclreva = -rainclr_tmp(i)
               dsnowclrsub = -snowclr_tmp(i)

            END IF ! precipfracclr_tmp .GT. eps

            !---------------------------
            !-- EVAP / SUBL IN THE CLOUD
            !---------------------------

            IF (ok_growth_precip_deposition .AND. (temp(i) .LE. RTT) .AND. (precipfraccld_tmp(i) .GT. eps)) THEN
               !--Evaporation of liquid precipitation coming from above
               !--in the cloud only
               !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
               !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
               !--Exact explicit formulation (raincld is resolved exactly, qvap explicitly)
               !--which does not need a barrier on raincld, because included in the formula

               draincldeva = precipfraccld_tmp(i)*MAX(0., &
                                                      -coef_eva*(1.-expo_eva)*(1.-qvapcld/qsatl(i))*dz(i) &
                                                      + (raincld_tmp(i)/precipfraccld_tmp(i))**(1.-expo_eva) &
                                                      )**(1./(1.-expo_eva)) - raincld_tmp(i)

               !--Evaporation is limited by 0
               !--NB. with ok_ice_supersat activated, this barrier should be useless
               draincldeva = MIN(0., draincldeva)

               !--Sublimation of the solid precipitation coming from above
               !--(same formula as for liquid precip)
               !--Exact explicit formulation (snowcld is resolved exactly, qvap explicitly)
               !--which does not need a barrier on snowcld, because included in the formula
               dsnowcldsub = precipfraccld_tmp(i)*MAX(0., &
                                                      -coef_sub*(1.-expo_sub)*(1.-qvapcld/qsati(i))*dz(i) &
                                                      + (snowcld_tmp(i)/precipfraccld_tmp(i))**(1.-expo_sub) &
                                                      )**(1./(1.-expo_sub)) - snowcld_tmp(i)

               !--There is no barrier because we want deposition to occur if it is possible

               !--Evaporation limit: we ensure that the layer's fraction below
               !--the clear sky does not reach saturation. In this case, we
               !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice
               !--Max evaporation is computed not to saturate the clear sky precip fraction
               !--(i.e., the fraction where evaporation occurs)
               !--It is expressed as a max flux dprecip_evasub_max

               dprecip_evasub_max = (qvapcld - qsat(i))*precipfraccld_tmp(i)*dhum_to_dflux(i)
               dprecip_evasub_tot = draincldeva + dsnowcldsub

               !--Barriers
               !--If activates if the total is LOWER than the max because
               !--everything is negative
               IF (((dprecip_evasub_max .LT. 0.) .AND. &
                    (dprecip_evasub_tot .LT. dprecip_evasub_max)) .OR. &
                   ((dprecip_evasub_max .GT. 0.) .AND. &
                    (dprecip_evasub_tot .GT. dprecip_evasub_max))) THEN
                  draincldeva = dprecip_evasub_max*draincldeva/dprecip_evasub_tot
                  dsnowcldsub = dprecip_evasub_max*dsnowcldsub/dprecip_evasub_tot
               END IF

            END IF ! ok_unadjusted_clouds

            !--New solid and liquid precipitation fluxes after evap and sublimation
            dqrevap = (drainclreva + draincldeva)/dhum_to_dflux(i)
            dqssubl = (dsnowclrsub + dsnowcldsub)/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_tmp(i) = MAX(0., rainclr_tmp(i) + drainclreva)
            snowclr_tmp(i) = MAX(0., snowclr_tmp(i) + dsnowclrsub)
            raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva)
            snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub)

            !--We reattribute fluxes according to initial fractions, using proportionnality
            !--Note that trough this process, we conserve total precip fluxes
            !-- rainclr_tmp + raincld_tmp = rainclr + raincld (same for snow)
            !--If the cloud has increased, a part of the precip in clear sky falls in clear sky,
            !--the rest falls in the cloud
            IF ((precipfraccld(i) .GT. eps) .AND. (precipfraccld_tmp(i) .GT. precipfraccld(i))) THEN
               !--All the precip from cloud falls in the cloud
               raincld(i) = MAX(0., raincld_tmp(i)*precipfraccld(i)/precipfraccld_tmp(i))
               rainclr(i) = MAX(0., rainclr_tmp(i) + raincld_tmp(i) - raincld(i))
               snowcld(i) = MAX(0., snowcld_tmp(i)*precipfraccld(i)/precipfraccld_tmp(i))
               snowclr(i) = MAX(0., snowclr_tmp(i) + snowcld_tmp(i) - snowcld(i))

               !--If the cloud has narrowed, a part of the precip in cloud falls in clear sky,
               !--the rest falls in the cloud
            ELSEIF ((precipfracclr(i) .GT. eps) .AND. (precipfracclr_tmp(i) .GT. precipfracclr(i))) THEN
               !--All the precip from clear sky falls in clear sky
               rainclr(i) = MAX(0., rainclr_tmp(i)*precipfracclr(i)/precipfracclr_tmp(i))
               raincld(i) = MAX(0., raincld_tmp(i) + rainclr_tmp(i) - rainclr(i))
               snowclr(i) = MAX(0., snowclr_tmp(i)*precipfracclr(i)/precipfracclr_tmp(i))
               snowcld(i) = MAX(0., snowcld_tmp(i) + snowclr_tmp(i) - snowclr(i))

            ELSE
               !--If the cloud stays the same or if there is no cloud above and
               !--in the current layer, there is no reattribution
               rainclr(i) = rainclr_tmp(i)
               raincld(i) = raincld_tmp(i)
               snowclr(i) = snowclr_tmp(i)
               snowcld(i) = snowcld_tmp(i)
            END IF

            !--If there is no more precip fluxes, the precipitation fraction is set to 0
            IF ((rainclr(i) + snowclr(i)) .LE. 0.) precipfracclr(i) = 0.
            IF ((raincld(i) + snowcld(i)) .LE. 0.) precipfraccld(i) = 0.

            !--Calculation of the total fluxes
            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.

         END IF ! ( ( rain(i) + snow(i) ) .GT. 0. )

         !--Diagnostic tendencies
         dqssub(i) = dqssubl/dtime
         dqreva(i) = dqrevap/dtime

      END DO ! loop on klon

   END SUBROUTINE poprecip_precld

!----------------------------------------------------------------
! Computes the precipitation fraction update
! The goal of this routine is to reattribute precipitation fractions
! and fluxes to clear or cloudy air, depending on the variation of
! the cloud fraction on the vertical dimension. We assume a
! maximum-random overlap of the cloud cover (see Jakob and Klein, 2000,
! and LTP thesis, 2021)
! NB. in fact, we assume a maximum-random overlap of the total precip. frac
!
   SUBROUTINE poprecip_fracupdate( &
      klon, cldfra, precipfracclr, precipfraccld, &
      rainclr, raincld, snowclr, snowcld)

      USE lmdz_lscp_ini, ONLY: eps

      IMPLICIT NONE

      INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]

      REAL, INTENT(IN), DIMENSION(klon) :: cldfra         !--cloud fraction [-]
      REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
      REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
      !--NB. at the end of the routine, becomes the fraction of precip
      !--in the current layer

      REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2]
      REAL, INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
      REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2]
      REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]

!--Local variables
      INTEGER :: i
      REAL :: dcldfra
      REAL :: precipfractot
      REAL :: dprecipfracclr, dprecipfraccld
      REAL :: drainclr, dsnowclr
      REAL :: draincld, dsnowcld

      DO i = 1, klon

         !--Initialisation
         precipfractot = precipfracclr(i) + precipfraccld(i)

         !--Instead of using the cloud cover which was use in LTP thesis, we use the
         !--total precip. fraction to compute the maximum-random overlap. This is
         !--because all the information of the cloud cover is embedded into
         !--precipfractot, and this allows for taking into account the potential
         !--reduction of the precipitation fraction because either the flux is too
         !--small (see barrier at the end of poprecip_postcld) or the flux is completely
         !--evaporated (see barrier at the end of poprecip_precld)
         !--NB. precipfraccld(i) is here the cloud fraction of the layer above
         !precipfractot = 1. - ( 1. - precipfractot ) * &
         !               ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
         !             / ( 1. - MIN( precipfraccld(i), 1. - eps ) )

         IF (precipfraccld(i) .GT. (1.-eps)) THEN
            precipfractot = 1.
         ELSE
            precipfractot = 1.-(1.-precipfractot)* &
                            (1.-MAX(cldfra(i), precipfraccld(i))) &
                            /(1.-precipfraccld(i))
         END IF

         !--precipfraccld(i) is here the cloud fraction of the layer above
         dcldfra = cldfra(i) - precipfraccld(i)
         !--Tendency of the clear-sky precipitation fraction. We add a MAX on the
         !--calculation of the current CS precip. frac.
         !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
         !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated
         !--if precipfractot < cldfra)
         dprecipfracclr = (precipfractot - cldfra(i)) - precipfracclr(i)
         !--Tendency of the cloudy precipitation fraction. We add a MAX on the
         !--calculation of the current CS precip. frac.
         !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) )
         !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated
         !--if cldfra < 0)
         dprecipfraccld = dcldfra

         !--If the cloud extends
         IF (dprecipfraccld .GT. 0.) THEN
            !--If there is no CS precip, nothing happens.
            !--If there is, we reattribute some of the CS precip flux
            !--to the cloud precip flux, proportionnally to the
            !--decrease of the CS precip fraction
            IF (precipfracclr(i) .LT. eps) THEN
               drainclr = 0.
               dsnowclr = 0.
            ELSE
               drainclr = dprecipfracclr/precipfracclr(i)*rainclr(i)
               dsnowclr = dprecipfracclr/precipfracclr(i)*snowclr(i)
            END IF
            !--If the cloud narrows
         ELSEIF (dprecipfraccld .LT. 0.) THEN
            !--We reattribute some of the cloudy precip flux
            !--to the CS precip flux, proportionnally to the
            !--decrease of the cloud precip fraction
            IF (precipfraccld(i) .LT. eps) THEN
               draincld = 0.
               dsnowcld = 0.
            ELSE
               draincld = dprecipfraccld/precipfraccld(i)*raincld(i)
               dsnowcld = dprecipfraccld/precipfraccld(i)*snowcld(i)
            END IF
            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.
         END IF

         !--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)

      END DO

   END SUBROUTINE poprecip_fracupdate

!----------------------------------------------------------------
! Computes the processes-oriented precipitation formulations for
! - autoconversion (auto) via a deposition (ice) and coalescence (liquid) 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, &
                               thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
                               rho_rain, r_rain, r_snow, rho_ice, &
                               expo_tau_auto_snow, &
                               tau_auto_snow_min, tau_auto_snow_max, &
                               eps, rain_int_min, &
                               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
      REAL                  :: min_precip                       !--minimum precip flux below which precip fraction decreases

!--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.

!--Update the precipitation fraction following cloud formation
      CALL poprecip_fracupdate( &
         klon, cldfra, precipfracclr, precipfraccld, &
         rainclr, raincld, snowclr, snowcld)

      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)

         !--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)
         END IF

         !--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
            END IF

            !--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
            END IF

            !---------------------------------------------------------
            !--                  AUTOCONVERSION
            !---------------------------------------------------------
            !--Autoconversion converts liquid droplets/ice crystals into
            !--rain drops/snowflakes following a Kessler type scheme.
            !-- 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_tau_auto_snow)

               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_tau_auto_snow)

               expo_auto_rain = cld_expo_lsc
               expo_auto_snow = cld_expo_lsc
            END IF

            ! 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
            END IF

         END IF ! cldfra .GE. seuil_neb

      END DO ! 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))
            END IF

            !--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))
            END IF

            !--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

            END IF

            !--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

         END IF

         !---------------------------------------------------------
         !--                  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 formula
         !--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))

            END IF

            !-- 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.)
            END IF

            !--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.)
            END IF

            !--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
            END IF

            !--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

         END IF

         !--If the local flux of rain+snow in clear air is lower than min_precip,
         !--we reduce the precipitation fraction in the clear air so that the new
         !--local flux of rain+snow is equal to min_precip.
         !--we apply the minimum only on the clear-sky fraction because the cloudy precip fraction
         !--already decreases out of clouds
         !--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 min_precip, i.e.,
         !-- (rain+snow)/precipfrac < min_precip , i.e.,
         !-- (rain+snow)/min_precip < precipfrac , then we want to reduce
         !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/min_precip.
         !--Note that this is physically different than what is proposed in LTP thesis.
         !--min_precip is either equal to rain_int_min or calculated as a very small fraction
         !--of the minimum precip flux estimated as the flux associated with the
         !--autoconversion threshold mass content
         !min_precip=1.e-6*(pplay(i)/RD/temp(i))*MIN(rain_fallspeed_clr*cld_lc_lsc,snow_fallspeed_clr*cld_lc_lsc_snow)
         min_precip = rain_int_min
         precipfracclr(i) = MIN(precipfracclr(i), (rainclr(i) + snowclr(i))/min_precip)

         !--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)

      END DO ! loop on klon

   END SUBROUTINE poprecip_postcld

END MODULE lmdz_lscp_precip
