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.
  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
      znebprecip(i)=0.

    ENDIF ! (zrfl(i)+zifl(i).GT.0.)

  ENDDO ! loop on klon

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

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

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


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

   
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
    ENDIF
  ENDDO

  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
        ENDIF


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


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


        ! 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.
          ENDIF
        ELSE
        ! old explicit resolution with subtimesteps
          zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*zvelo(i) 
        ENDIF

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

      ENDIF


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

      ! 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

    ENDIF ! rneb >0

  ENDDO  ! i = 1,klon

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

  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)
    ENDIF
    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)
    ENDIF
    ! c_iso : same for isotopes

  ENDIF ! rneb>0

ENDDO

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

    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
    ENDIF
  ENDDO

ENDIF

! 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)
      ENDIF
      zradocond(i)=zradoliq(i)+zradoice(i)
    ENDDO
  ENDIF
  ! keep in memory air density and ice fallspeed velocity
  zrho_up(:) = zrho(:)
  zvelo_up(:) = zvelo(:)
ENDIF

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


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

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

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

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

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


DO i = 1, klon

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

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

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

      !--Evaporation of liquid precipitation coming from above
      !--in the clear sky only
      !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva)
      !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
      !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly)
      !--which does not need a barrier on rainclr, because included in the formula
      drainclreva = precipfracclr_tmp(i) * MAX(0., &
                  - coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) &
                  + ( rainclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_eva ) &
                  )**( 1. / ( 1. - expo_eva ) ) - rainclr_tmp(i)
               
      !--Evaporation is limited by 0
      !--NB. with ok_ice_supersat activated, this barrier should be useless
      drainclreva = MIN(0., drainclreva)
     
      ! we set it to 0 as not sufficiently tested
      drainclreva = 0.

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

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


      !--Evaporation limit: we ensure that the layer's fraction below
      !--the clear sky does not reach saturation. In this case, we 
      !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice 
      !--Max evaporation is computed not to saturate the clear sky precip fraction
      !--(i.e., the fraction where evaporation occurs)
      !--It is expressed as a max flux dprecip_evasub_max
      
      IF ( .NOT. ok_ice_supersat ) THEN
        dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr_tmp(i)) &
                         * dhum_to_dflux(i)
      ELSE
        dprecip_evasub_max = ( qvapclr - qsat(i) ) * precipfracclr_tmp(i) &
                         * dhum_to_dflux(i)
      ENDIF
      dprecip_evasub_tot = drainclreva + dsnowclrsub

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

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

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

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

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

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

USE lmdz_lscp_ini, ONLY : eps

IMPLICIT NONE

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

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

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

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


DO i = 1, klon

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

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


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

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


  !--If the cloud extends
  IF ( dprecipfraccld .GT. 0. ) THEN
    !--If there is no CS precip, nothing happens.
    !--If there is, we reattribute some of the CS precip flux
    !--to the cloud precip flux, proportionnally to the
    !--decrease of the CS precip fraction
    IF ( precipfracclr(i) .LT. eps ) THEN
      drainclr = 0.
      dsnowclr = 0.
    ELSE
      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i) 
      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i) 
    ENDIF
  !--If the cloud narrows
  ELSEIF ( dprecipfraccld .LT. 0. ) THEN
    !--We reattribute some of the cloudy precip flux
    !--to the CS precip flux, proportionnally to the
    !--decrease of the cloud precip fraction
    IF ( precipfraccld(i) .LT. eps ) THEN
      draincld = 0.
      dsnowcld = 0.
    ELSE
      draincld = dprecipfraccld / precipfraccld(i) * raincld(i) 
      dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i) 
    ENDIF
    drainclr = - draincld
    dsnowclr = - dsnowcld
  !--If the cloud stays the same or if there is no cloud above and
  !--in the current layer, nothing happens
  ELSE
    drainclr = 0.
    dsnowclr = 0.
  ENDIF

  !--We add the tendencies
  !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision)
  precipfraccld(i) = precipfraccld(i) + dprecipfraccld
  precipfracclr(i) = precipfracclr(i) + dprecipfracclr
  rainclr(i) = MAX(0., rainclr(i) + drainclr)
  snowclr(i) = MAX(0., snowclr(i) + dsnowclr)
  raincld(i) = MAX(0., raincld(i) - drainclr)
  snowcld(i) = MAX(0., snowcld(i) - dsnowclr)

ENDDO

END SUBROUTINE poprecip_fracupdate


!----------------------------------------------------------------
! 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,               & 
                          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,                &
                          thresh_precip_frac, 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)
  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_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
    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 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) )

    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 air is lower than min_precip,
  !--we reduce the precipiration 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 )


ENDDO ! loop on klon

END SUBROUTINE poprecip_postcld

END MODULE lmdz_lscp_precip
