MODULE lmdz_lscp

IMPLICIT NONE

CONTAINS

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
     paprs,pplay,temp,qt,qice_save,ptconv,ratqs,        &
     d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri,  &
     pfraclr, pfracld,                                  &
     cldfraliq, sigma2_icefracturb,mean_icefracturb,    &
     radocond, radicefrac, rain, snow,                  &
     frac_impa, frac_nucl, beta,                        &
     prfl, psfl, rhcl, qta, fraca,                      &
     tv, pspsk, tla, thl, iflag_cld_th,                 &
     iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
     distcltop, temp_cltop, tke, tke_dissip,            &
     qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
     Tcontr, qcontr, qcontr2, fcontrN, fcontrP,         &
     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
     dqsmelt, dqsfreez)

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
!--------------------------------------------------------------------------------
! Date: 01/2021
!--------------------------------------------------------------------------------
! Aim: Large Scale Clouds and Precipitation (LSCP)

! This code is a new version of the fisrtilp.F90 routine, which is itself a 
! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
! and 'ilp' (il pleut, L. Li)

! Compared to the original fisrtilp code, lscp
! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
! -> consider always precipitation thermalisation (fl_cor_ebil>0)
! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T) 
! -> option "oldbug" by JYG has been removed
! -> iflag_t_glace >0 always
! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
! -> rectangular distribution from L. Li no longer available
! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt) 
!--------------------------------------------------------------------------------
! References: 

! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
! - Touzze-Peifert Ludo, PhD thesis, p117-124
! -------------------------------------------------------------------------------
! Code structure:

! P0>     Thermalization of the precipitation coming from the overlying layer
! P1>     Evaporation of the precipitation (falling from the k+1 level)
! P2>     Cloud formation (at the k level) 
! P2.A.1> With the PDFs, calculation of cloud properties using the inital
!         values of T and Q
! P2.A.2> Coupling between condensed water and temperature
! P2.A.3> Calculation of final quantities associated with cloud formation
! P2.B>   Release of Latent heat after cloud formation
! P3>     Autoconversion to precipitation (k-level)
! P4>     Wet scavenging
!------------------------------------------------------------------------------
! Some preliminary comments (JBM) :

! The cloud water that the radiation scheme sees is not the same that the cloud
! water used in the physics and the dynamics 

! During the autoconversion to precipitation (P3 step), radocond (cloud water used
! by the radiation scheme) is calculated as an average of the water that remains
! in the cloud during the precipitation and not the water remaining at the end
! of the time step. The latter is used in the rest of the physics and advected
! by the dynamics. 

! In summary:

! Radiation:
! xflwc(newmicro)+xfiwc(newmicro) =
!   radocond=lwcon(nc)+iwcon(nc)

! Notetheless, be aware of:

! radocond .NE. ocond(nc)
! i.e.:
! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
! but oliq+(ocond-oliq) .EQ. ocond
! (which is not trivial)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

! USE de modules contenant des fonctions.
USE lmdz_cloudth, ONLY: cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat
USE lmdz_lscp_tools, ONLY: icefrac_lscp, icefrac_lscp_turb
USE lmdz_lscp_tools, ONLY: fallice_velocity, distance_to_cloud_top
USE ice_sursat_mod, ONLY: ice_sursat
USE lmdz_lscp_poprecip, ONLY: poprecip_precld, poprecip_postcld

! Use du module lmdz_lscp_ini contenant les constantes
USE lmdz_lscp_ini, ONLY: prt_level, lunout
USE lmdz_lscp_ini, ONLY: seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
USE lmdz_lscp_ini, ONLY: ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
USE lmdz_lscp_ini, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
USE lmdz_lscp_ini, ONLY: coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
USE lmdz_lscp_ini, ONLY: iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
USE lmdz_lscp_ini, ONLY: iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
USE lmdz_lscp_ini, ONLY: ok_poprecip
USE lmdz_lscp_ini, ONLY: iflag_icefrac
IMPLICIT NONE

!===============================================================================
! VARIABLES DECLARATION
!===============================================================================

! INPUT VARIABLES:
!-----------------

  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output

  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg] 
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qice_save       ! ice specific from previous time step [kg/kg]
  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
  LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated

  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active

  !Inputs associated with thermal plumes

  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp) 
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]

  ! INPUT/OUTPUT variables
  !------------------------
  
  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl              ! liquid potential temperature [K]
  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs            ! function of pressure that sets the large-scale


  ! Input sursaturation en glace
  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rneb_seri           ! fraction nuageuse en memoire
  
  ! OUTPUT variables
  !-----------------

  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t                 ! temperature increment [K]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q                 ! specific humidity increment [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql                ! liquid water increment [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi                ! cloud ice mass increment [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb                ! cloud fraction [-]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol           ! cloud fraction per unit volume [-]  
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr             ! precip fraction clear-sky part [-]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld             ! precip fraction cloudy part [-]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq           ! liquid fraction of cloud [-]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond            ! condensed water used in the radiation scheme [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac          ! ice fraction of condensed water for radiation scheme
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl                ! clear-sky relative humidity [-]
  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain                ! surface large-scale rainfall [kg/s/m2] 
  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow                ! surface large-scale snowfall [kg/s/m2]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl               ! saturation specific humidity wrt liquid [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsats               ! saturation specific humidity wrt ice [kg/kg]  
  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl                ! large-scale rainfall flux in the column [kg/s/m2]
  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl                ! large-scale snowfall flux in the column [kg/s/m2]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop           ! distance to cloud top [m]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop          ! temperature of cloud top [K]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta                ! conversion rate of condensed water

  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
  
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
  
  ! for supersaturation and contrails parameterisation
  
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qclr                ! specific total water content in clear sky region [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld                ! specific total water content in cloudy region [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qss                 ! specific total water content in supersat region [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qvc                 ! specific vapor content in clouds [kg/kg] 
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebclr             ! mesh fraction of clear sky [-]   
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebss              ! mesh fraction of ISSR [-]  
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_ss            ! coefficient governing the ice nucleation RHi threshold [-]      
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr              ! threshold temperature for contrail formation [K]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr              ! threshold humidity for contrail formation [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2             ! // (2nd expression more consistent with LMDZ expression of q)          
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN             ! fraction of grid favourable to non-persistent contrails
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP             ! fraction of grid favourable to persistent contrails
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth         ! mean saturation deficit in thermals
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv        ! mean saturation deficit in environment
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath     ! std of saturation deficit in thermals
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv    ! std of saturation deficit in environment

  ! for POPRECIP

  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]


  ! LOCAL VARIABLES:
  !----------------
  REAL,DIMENSION(klon) :: qsl, qsi                                ! saturation threshold at current vertical level
  REAL :: zct, zcl,zexpo
  REAL, DIMENSION(klon,klev) :: ctot
  REAL, DIMENSION(klon,klev) :: ctot_vol
  REAL, DIMENSION(klon) :: zqs, zdqs
  REAL :: zdelta, zcor, zcvm5
  REAL, DIMENSION(klon) :: zdqsdT_raw
  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
  REAL, DIMENSION(klon) :: Tbef,qlbef,DT                          ! temperature, humidity and temp. variation during lognormal iteration
  REAL :: num,denom
  REAL :: cste
  REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta             ! lognormal parameters
  REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2          ! lognormal intermediate variables
  REAL :: erf
  REAL, DIMENSION(klon) :: zfice_th
  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
  REAL, DIMENSION(klon) :: zrfl, zrfln
  REAL :: zqev, zqevt 
  REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
  REAL :: zqev0,zqevi, zqevti
  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
  REAL, DIMENSION(klon) :: zoliql, zoliqi
  REAL, DIMENSION(klon) :: zt
  REAL, DIMENSION(klon,klev) :: zrho
  REAL, DIMENSION(klon) :: zdz,iwc
  REAL :: zchau,zfroi
  REAL, DIMENSION(klon) :: zfice,zneb,znebprecip
  REAL :: zmelt,zrain,zsnow,zprecip
  REAL, DIMENSION(klon) :: dzfice
  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
  REAL :: zsolid
  REAL, DIMENSION(klon) :: qtot, qzero
  REAL, DIMENSION(klon) :: dqsl,dqsi
  REAL :: smallestreal
  !  Variables for Bergeron process
  REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
  REAL, DIMENSION(klon) :: zqpreci, zqprecl
  ! Variables precipitation energy conservation
  REAL, DIMENSION(klon) :: zmqc
  REAL :: zalpha_tr
  REAL :: zfrac_lessi
  REAL, DIMENSION(klon) :: zprec_cond
  REAL :: zmair
  REAL :: zcpair, zcpeau
  REAL, DIMENSION(klon) :: zlh_solid
  REAL, DIMENSION(klon) :: ztupnew
  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
  REAL :: zm_solid         ! for liquid -> solid conversion
  REAL, DIMENSION(klon) :: zrflclr, zrflcld
  REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
  REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr 
  REAL, DIMENSION(klon) :: ziflclr, ziflcld
  REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
  REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
  REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
  REAL, DIMENSION(klon,klev) :: velo
  REAL :: vr, ffallv
  REAL :: qlmpc, qimpc, rnebmpc
  REAL, DIMENSION(klon,klev) :: radocondi, radocondl
  REAL :: effective_zneb
  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb

  INTEGER i, k, n, kk, iter
  INTEGER, DIMENSION(klon) :: n_i
  INTEGER ncoreczq
  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
  LOGICAL iftop

  LOGICAL, DIMENSION(klon) :: lognormale
  LOGICAL, DIMENSION(klon) :: keepgoing

  CHARACTER (len = 20) :: modname = 'lscp'
  CHARACTER (len = 80) :: abort_message
  

!===============================================================================
! INITIALISATION
!===============================================================================

! Few initial checks


IF (iflag_fisrtilp_qsat < 0) THEN
    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
    CALL abort_physic(modname,abort_message,1)
ENDIF

! Few initialisations

znebprecip(:)=0.0
ctot_vol(1:klon,1:klev)=0.0
rneblsvol(1:klon,1:klev)=0.0
smallestreal=1.e-9
znebprecipclr(:)=0.0
znebprecipcld(:)=0.0
mpc_bl_points(:,:)=0

IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'

! AA for 'safety' reasons
zalpha_tr   = 0.
zfrac_lessi = 0.
beta(:,:)= 0.

! Initialisation of variables:

prfl(:,:) = 0.0
psfl(:,:) = 0.0
d_t(:,:) = 0.0
d_q(:,:) = 0.0
d_ql(:,:) = 0.0
d_qi(:,:) = 0.0
rneb(:,:) = 0.0
pfraclr(:,:)=0.0
pfracld(:,:)=0.0
cldfraliq(:,:)=0.
sigma2_icefracturb(:,:)=0.
mean_icefracturb(:,:)=0.
radocond(:,:) = 0.0
radicefrac(:,:) = 0.0
frac_nucl(:,:) = 1.0 
frac_impa(:,:) = 1.0 
rain(:) = 0.0
snow(:) = 0.0
zoliq(:)=0.0
zfice(:)=0.0
dzfice(:)=0.0
zfice_turb(:)=0.0
dzfice_turb(:)=0.0
zqprecl(:)=0.0
zqpreci(:)=0.0
zrfl(:) = 0.0
zifl(:) = 0.0
ziflprev(:)=0.0
zneb(:) = seuil_neb
zrflclr(:) = 0.0
ziflclr(:) = 0.0
zrflcld(:) = 0.0
ziflcld(:) = 0.0
tot_zneb(:) = 0.0
tot_znebn(:) = 0.0
d_tot_zneb(:) = 0.0
qzero(:) = 0.0
zdistcltop(:)=0.0
ztemp_cltop(:) = 0.0
ztupnew(:)=0.0
!--ice supersaturation
gamma_ss(:,:) = 1.0
qss(:,:) = 0.0
rnebss(:,:) = 0.0
Tcontr(:,:) = missing_val
qcontr(:,:) = missing_val
qcontr2(:,:) = missing_val
fcontrN(:,:) = 0.0
fcontrP(:,:) = 0.0
distcltop(:,:)=0.
temp_cltop(:,:)=0.

!-- poprecip
qraindiag(:,:)= 0.
qsnowdiag(:,:)= 0.
dqreva(:,:)   = 0.
dqrauto(:,:)  = 0.
dqrmelt(:,:)  = 0.
dqrfreez(:,:) = 0.
dqrcol(:,:)   = 0.
dqssub(:,:)   = 0.
dqsauto(:,:)  = 0.
dqsrim(:,:)   = 0.
dqsagg(:,:)   = 0.
dqsfreez(:,:) = 0.
dqsmelt(:,:)  = 0.
zqupnew(:)    = 0.
zqvapclr(:)   = 0.



!c_iso: variable initialisation for iso 


!===============================================================================
!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
!===============================================================================

ncoreczq=0

DO k = klev, 1, -1

    IF (k<=klev-1) THEN
       iftop=.FALSE.
    ELSE 
       iftop=.TRUE.
    ENDIF

    ! Initialisation temperature and specific humidity
    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated 
    ! d_t = temperature tendency due to lscp
    ! The temperature of the overlying layer is updated here because needed for thermalization
    DO i = 1, klon
        zt(i)=temp(i,k)
        zq(i)=qt(i,k)
        IF (.not. iftop) THEN
           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
           !--zqs(i) is the saturation specific humidity in the layer above
           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
        ENDIF
        !c_iso init of iso
    ENDDO
   
    !================================================================
    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
    IF (ok_poprecip) THEN

            CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
                              zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
                              zqvapclr, zqupnew, &
                              zrfl, zrflclr, zrflcld, &
                              zifl, ziflclr, ziflcld, &
                              dqreva(:,k),dqssub(:,k) &
                             )

    !================================================================
    ELSE

    ! --------------------------------------------------------------------
    ! P1> Thermalization of precipitation falling from the overlying layer
    ! --------------------------------------------------------------------
    ! Computes air temperature variation due to enthalpy transported by
    ! precipitation. Precipitation is then thermalized with the air in the
    ! layer. 
    ! The precipitation should remain thermalized throughout the different
    ! thermodynamical transformations. 
    ! The corresponding water mass should
    ! be added when calculating the layer's enthalpy change with 
    ! temperature
    ! See lmdzpedia page todoan
    ! todoan: check consistency with ice phase
    ! todoan: understand why several steps
    ! ---------------------------------------------------------------------
    
    IF (iftop) THEN

        DO i = 1, klon
            zmqc(i) = 0.
        ENDDO

    ELSE

        DO i = 1, klon
 
            zmair=(paprs(i,k)-paprs(i,k+1))/RG
            ! no condensed water so cp=cp(vapor+dry air)
            ! RVTMP2=rcpv/rcpd-1
            zcpair=RCPD*(1.0+RVTMP2*zq(i)) 
            zcpeau=RCPD*RVTMP2

            ! zmqc: precipitation mass that has to be thermalized with 
            ! layer's air so that precipitation at the ground has the
            ! same temperature as the lowermost layer
            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
            zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) &
             / (zcpair + zmqc(i)*zcpeau)

        ENDDO
 
    ENDIF

    ! --------------------------------------------------------------------
    ! P2> Precipitation evaporation/sublimation/melting
    ! --------------------------------------------------------------------
    ! A part of the precipitation coming from above is evaporated/sublimated/melted.
    ! --------------------------------------------------------------------
    
    IF (iflag_evap_prec>=1) THEN

        ! Calculation of saturation specific humidity
        ! depending on temperature:
        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.FALSE.,zqs(:),zdqs(:))
        ! wrt liquid water
        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.FALSE.,qsl(:),dqsl(:))
        ! wrt ice
        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.FALSE.,qsi(:),dqsi(:))

        DO i = 1, klon

            ! if precipitation
            IF (zrfl(i)+zifl(i)>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>=4) THEN
                    zrfl(i) = zrflclr(i)
                    zifl(i) = ziflclr(i)
                ENDIF

                IF (iflag_evap_prec==1) THEN
                    znebprecip(i)=zneb(i)
                ELSE
                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
                ENDIF

                IF (iflag_evap_prec>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 == 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==3) THEN
                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                ELSE IF (iflag_evap_prec>=4) THEN
                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                ELSE
                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                ENDIF

                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
                    *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 

                ! sublimation of the solid precipitation coming from above
                IF (iflag_evap_prec==3) THEN
                    zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                ELSE IF (iflag_evap_prec>=4) THEN
                     zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) &
                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                ELSE
                    zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                ENDIF

                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
                *RG*dtime/(paprs(i,k)-paprs(i,k+1))

                ! A. JAM
                ! Evaporation limit: we ensure that the layer's fraction below
                ! the cloud or the whole mesh (depending on iflag_evap_prec)
                ! does not reach saturation. In this case, we 
                ! redistribute zqev0 conserving the ratio liquid/ice 

                IF (zqevt+zqevti>zqev0) THEN
                    zqev=zqev0*zqevt/(zqevt+zqevti)
                    zqevi=zqev0*zqevti/(zqevt+zqevti)
                ELSE
                    zqev=zqevt
                    zqevi=zqevti
                ENDIF


                ! New solid and liquid precipitation fluxes after evap and sublimation
                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
                         /RG/dtime)
                zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
                         /RG/dtime)


                ! vapor, temperature, precip fluxes update
                ! vapor is updated after evaporation/sublimation (it is increased)
                zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
                ! zmqc is the total condensed water in the precip flux (it is decreased)
                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
                ! air and precip temperature (i.e., gridbox temperature)
                ! is updated due to latent heat cooling
                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
                + (zifln(i)-zifl(i))                      &
                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))

                ! New values of liquid and solid precipitation 
                zrfl(i) = zrfln(i)
                zifl(i) = zifln(i)

                ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
                !                         due to evap + sublim                                     
                                           

                IF (iflag_evap_prec>=4) THEN
                    zrflclr(i) = zrfl(i)
                    ziflclr(i) = zifl(i)
                    IF(zrflclr(i) + ziflclr(i)<=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>=4) THEN
                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
                    zrfl(i)=zrflclr(i)+zrflcld(i)
                ELSE
                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
                ENDIF


                ! c_iso: melting of isotopic precipi with zmelt (no fractionation)

                ! Latent heat of melting because of precipitation melting
                ! NB: the air + precip temperature is simultaneously updated
                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
               
                IF (iflag_evap_prec>=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)

    ENDIF ! (ok_poprecip)

    ! --------------------------------------------------------------------
    ! End precip evaporation
    ! --------------------------------------------------------------------
    
    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
    !-------------------------------------------------------

     qtot(:)=zq(:)+zmqc(:)
     CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.FALSE.,zqs(:),zdqs(:))
     DO i = 1, klon
            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
            IF (zq(i) < 1.e-15) THEN
                ncoreczq=ncoreczq+1
                zq(i)=1.e-15
            ENDIF
            ! c_iso: do something similar for isotopes

     ENDDO
    
    ! --------------------------------------------------------------------
    ! P2> Cloud formation
    !---------------------------------------------------------------------

    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach 
    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
    ! is prescribed and depends on large scale variables and boundary layer
    ! properties)
    ! The decrease in condensed part due tu latent heating is taken into
    ! account
    ! -------------------------------------------------------------------

        ! P2.1> With the PDFs (log-normal, bigaussian) 
        ! cloud properties calculation with the initial values of t and q
        ! ----------------------------------------------------------------

        ! initialise gammasat and qincloud_mpc
        gammasat(:)=1.
        qincloud_mpc(:)=0.

        IF (iflag_cld_th>=5) THEN
               ! Cloud cover and content in meshes affected by shallow convection, 
               ! are retrieved from a bi-gaussian distribution of the saturation deficit
               ! following Jam et al. 2013

               IF (iflag_cloudth_vert<=2) THEN
                  ! Old version of Arnaud Jam

                    CALL cloudth(klon,klev,k,tv,                  &
                         zq,qta,fraca,                            &
                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
                         ratqs,zqs,temp,                              &
                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)


                ELSEIF (iflag_cloudth_vert>=3 .AND. iflag_cloudth_vert<=5) THEN
                   ! Default version of Arnaud Jam
                        
                    CALL cloudth_v3(klon,klev,k,tv,                        &
                         zq,qta,fraca,                                     &
                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
                         ratqs,zqs,temp, &
                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)


                ELSEIF (iflag_cloudth_vert==6) THEN
                   ! Jean Jouhaud's version, with specific separation between surface and volume
                   ! cloud fraction Decembre 2018 

                    CALL cloudth_v6(klon,klev,k,tv,                        &
                         zq,qta,fraca,                                     &
                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
                         ratqs,zqs,temp, &
                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)

                ELSEIF (iflag_cloudth_vert == 7) THEN
                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
                   ! for boundary-layer mixed phase clouds 
                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))

                ENDIF


                DO i=1,klon
                    rneb(i,k)=ctot(i,k)
                    rneblsvol(i,k)=ctot_vol(i,k)
                    zqn(i)=qcloud(i)
                ENDDO

        ENDIF

        IF (iflag_cld_th <= 4) THEN
            
                ! lognormal 
            lognormale(:) = .TRUE.

        ELSEIF (iflag_cld_th >= 6) THEN
            
                ! lognormal distribution when no thermals 
            lognormale(:) = fraca(:,k) < min_frac_th_cld

        ELSE
                ! When iflag_cld_th=5, we always assume
                ! bi-gaussian distribution
            lognormale(:) = .FALSE.
        
        ENDIF

        DT(:) = 0.
        n_i(:)=0
        Tbef(:)=zt(:)
        qlbef(:)=0.

        ! Treatment of non-boundary layer clouds (lognormale)
        ! condensation with qsat(T) variation (adaptation)
        ! Iterative resolution to converge towards qsat
        ! with update of temperature, ice fraction and qsat at 
        ! each iteration

        ! todoan -> sensitivity to iflag_fisrtilp_qsat
        DO iter=1,iflag_fisrtilp_qsat+1

                DO i=1,klon

                    ! keepgoing = .TRUE. while convergence is not satisfied
                    keepgoing(i)=ABS(DT(i))>DDT0

                    IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN

                        ! if not convergence:

                        ! P2.2.1> cloud fraction and condensed water mass calculation
                        ! Calculated variables: 
                        ! rneb : cloud fraction
                        ! zqn : total water within the cloud
                        ! zcond: mean condensed water within the mesh
                        ! rhcl: clear-sky relative humidity
                        !---------------------------------------------------------------

                        ! new temperature that only serves in the iteration process:
                        Tbef(i)=Tbef(i)+DT(i)

                        ! Rneb, qzn and zcond for lognormal PDFs
                        qtot(i)=zq(i)+zmqc(i)

                      ENDIF

                  ENDDO

                  ! Calculation of saturation specific humidity and ice fraction
                  CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.FALSE.,zqs(:),zdqs(:))
                  CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
                  ! cloud phase determination
                  IF (iflag_t_glace>=4) THEN
                  ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
                  ENDIF

                  CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, zdistcltop(:),ztemp_cltop(:),zfice(:),dzfice(:))

                  DO i=1,klon !todoan : check if loop in i is needed

                      IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN

                        zpdf_sig(i)=ratqs(i,k)*zq(i)
                        zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
                        zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
                        zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
                        zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
                        zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
                        zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
                        zpdf_e1(i)=1.-erf(zpdf_e1(i))
                        zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
                        zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
                        zpdf_e2(i)=1.-erf(zpdf_e2(i))

                        IF ((.NOT.ok_ice_sursat).OR.(Tbef(i)>t_glace_min)) THEN

                          IF (zpdf_e1(i)<1.e-10) THEN
                              rneb(i,k)=0.
                              zqn(i)=gammasat(i)*zqs(i)
                          ELSE
                              rneb(i,k)=0.5*zpdf_e1(i)
                              zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
                          ENDIF

                          rnebss(i,k)=0.0   !--added by OB (needed because of the convergence loop on time)
                          fcontrN(i,k)=0.0  !--idem
                          fcontrP(i,k)=0.0  !--idem
                          qss(i,k)=0.0      !--idem

                       ELSE ! in case of ice supersaturation by Audran

                        !------------------------------------
                        ! ICE SUPERSATURATION
                        !------------------------------------

                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), &
                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),                 &
                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                   & 
                             Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )

                        ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))

                        ! If vertical heterogeneity, change fraction by volume as well
                        IF (iflag_cloudth_vert>=3) THEN
                            ctot_vol(i,k)=rneb(i,k)
                            rneblsvol(i,k)=ctot_vol(i,k)
                        ENDIF


                       ! P2.2.2> Approximative calculation of temperature variation DT 
                       !           due to condensation.
                       ! Calculated variables: 
                       ! dT : temperature change due to condensation
                       !---------------------------------------------------------------

                
                        IF (zfice(i)<1) THEN
                            cste=RLVTT
                        ELSE
                            cste=RLSTT
                        ENDIF
                        
                        ! LEA_R : check formule
                        qlbef(i)=max(0.,zqn(i)-zqs(i))
                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
                              *qlbef(i)*dzfice(i)
                        ! here we update a provisory temperature variable that only serves in the iteration
                        ! process
                        DT(i)=num/denom
                        n_i(i)=n_i(i)+1

                    ENDIF ! end keepgoing

                ENDDO     ! end loop on i

        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1

        ! P2.3> Final quantities calculation 
        ! Calculated variables: 
        !   rneb : cloud fraction
        !   zcond: mean condensed water in the mesh
        !   zqn  : mean water vapor in the mesh
        !   zfice: ice fraction in clouds
        !   zt   : temperature
        !   rhcl : clear-sky relative humidity
        ! ----------------------------------------------------------------


        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
        IF (iflag_t_glace>=4) THEN
           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
           distcltop(:,k)=zdistcltop(:)
           temp_cltop(:,k)=ztemp_cltop(:)
        ENDIF

        ! Partition function depending on temperature
        CALL icefrac_lscp(klon, Tbef, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)

        ! Partition function depending on tke for non shallow-convective clouds 
        IF (iflag_icefrac >= 1) THEN

           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), qice_save(:,k), ziflcld, zqn, &
           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))

        ENDIF

        ! Water vapor update, Phase determination and subsequent latent heat exchange
        DO i=1, klon
            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the 
            ! iflag_cloudth_vert=7 and specific param is activated
            IF (mpc_bl_points(i,k) > 0) THEN
                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
                ! following line is very strange and probably wrong
                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
                ! water vapor update and partition function if thermals
                zq(i) = zq(i) - zcond(i)        
                zfice(i)=zfice_th(i)
            ELSE
                ! Checks on rneb, rhcl and zqn
                IF (rneb(i,k) <= 0.0) THEN
                    zqn(i) = 0.0
                    rneb(i,k) = 0.0
                    zcond(i) = 0.0
                    rhcl(i,k)=zq(i)/zqs(i)
                ELSE IF (rneb(i,k) >= 1.0) THEN
                    zqn(i) = zq(i)
                    rneb(i,k) = 1.0
                    zcond(i) = MAX(0.0,zqn(i)-zqs(i))
                    rhcl(i,k)=1.0
                ELSE
                    zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
                    ! following line is very strange and probably wrong:
                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param) 
                    IF (iflag_icefrac >= 1) THEN
                        IF (lognormale(i)) THEN  
                           zcond(i)  = zqliq(i) + zqice(i)
                           zfice(i)=zfice_turb(i)
                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
                        ENDIF
                    ENDIF
                ENDIF

                ! water vapor update 
                zq(i) = zq(i) - zcond(i)

            ENDIF

            ! c_iso : routine that computes in-cloud supersaturation
            ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input) 

            ! temperature update due to phase change
            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & 
       * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
        ENDDO

        ! If vertical heterogeneity, change volume fraction 
        IF (iflag_cloudth_vert >= 3) THEN
          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
        ENDIF

    ! ----------------------------------------------------------------
    ! P3> Precipitation formation 
    ! ----------------------------------------------------------------

    !================================================================
    IF (ok_poprecip) THEN

      DO i = 1, klon
        zoliql(i) = zcond(i) * ( 1. - zfice(i) )
        zoliqi(i) = zcond(i) * zfice(i)
      ENDDO

      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
                            ctot_vol(:,k), ptconv(:,k), &
                            zt, zq, zoliql, zoliqi, zfice, &
                            rneb(:,k), znebprecipclr, znebprecipcld, &
                            zrfl, zrflclr, zrflcld, &
                            zifl, ziflclr, ziflcld, &
                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
                            dqsmelt(:,k), dqsfreez(:,k) &
                            )

      DO i = 1, klon
        zoliq(i) = zoliql(i) + zoliqi(i)
        IF ( zoliq(i) > 0. ) THEN
                zfice(i) = zoliqi(i)/zoliq(i) 
        ELSE  
                zfice(i) = 0.0
        ENDIF

        ! calculation of specific content of condensates seen by radiative scheme 
        IF (ok_radocond_snow) THEN
           radocond(i,k) = zoliq(i)
           radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
           radocondi(i,k)= radocond(i,k)*zfice(i)+qsnowdiag(i,k)
        ELSE
           radocond(i,k) = zoliq(i)
           radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
           radocondi(i,k)= radocond(i,k)*zfice(i)
        ENDIF
      ENDDO

    !================================================================
    ELSE

    ! LTP:
    IF (iflag_evap_prec >= 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,k),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,k),znebprecipcld(i)) 
                IF (znebprecipcld(i) > 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,k) - d_tot_zneb(i) - zneb(i)))
                IF (znebprecipclr(i) > 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 of zoliq and radocond variables

    DO i = 1, klon
            zoliq(i) = zcond(i)
            zoliqi(i)= zoliq(i)*zfice(i)
            zoliql(i)= zoliq(i)*(1.-zfice(i))
            ! c_iso : initialisation of zoliq* also for isotopes
            zrho(i,k)  = pplay(i,k) / zt(i) / RD
            zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
            iwc(i)   = 0. 
            zneb(i)  = MAX(rneb(i,k),seuil_neb)
            radocond(i,k) = zoliq(i)/REAL(niter_lscp+1)
            radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1)
            radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
    ENDDO

       
    DO n = 1, niter_lscp

        DO i=1,klon
            IF (rneb(i,k)>0.0) THEN
                iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
            ENDIF
        ENDDO

        CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))

        DO i = 1, klon

            IF (rneb(i,k)>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)==seuil_neb) THEN
                    zprecip = 0.0
                    zsnow = 0.0
                    zrain= 0.0
                ELSE

                    IF (ptconv(i,k)) 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).

                    ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
                    ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
                    !.........................................................
                    IF ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) THEN

                    ! 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).
                        effective_zneb=ctot_vol(i,k)
                    ELSE
                        effective_zneb=zneb(i)
                    ENDIF 


                    IF (iflag_autoconversion == 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 == 2) THEN
                    ! exact resolution, niter_lscp=1 is sufficient but works only
                    ! with iflag_vice=0
                       IF (zoliqi(i) > 0.) THEN
                          zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & 
                          +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
                       ELSE
                          zfroi=0.
                       ENDIF
                    ELSE
                    ! old explicit resolution with subtimesteps
                       zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k) 
                    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 >= 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

                radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1)
                radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1)
                radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1)

            ENDIF ! rneb >0

        ENDDO  ! i = 1,klon

    ENDDO      ! n = 1,niter

    ! Precipitation flux calculation

    DO i = 1, klon
            
            IF (iflag_evap_prec>=4) THEN
                ziflprev(i)=ziflcld(i)
            ELSE
                ziflprev(i)=zifl(i)*zneb(i)
            ENDIF

            IF (rneb(i,k) > 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

                    zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
                    zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
                    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
                    coef1 = rneb(i,k)*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,k)*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>=4) THEN
                    zrflcld(i) = zrflcld(i)+zqprecl(i) &
                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &
                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
                    znebprecipcld(i) = rneb(i,k)
                    zrfl(i) = zrflcld(i) + zrflclr(i)
                    zifl(i) = ziflcld(i) + ziflclr(i)
                ELSE
                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
                    zifl(i) = zifl(i)+ zqpreci(i)        &
                   *(paprs(i,k)-paprs(i,k+1))/(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>=4) THEN

        DO i=1,klon

            IF ((zrflclr(i) + ziflclr(i)) > 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)) > 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
        !IF ( ((1-zfice(i))*zoliq(i) .GT. 0.) .AND. (zt(i) .LE. 233.15) ) THEN
        !PRINT*,'WARNING LEA OLIQ A <-40°C '
        !PRINT*,'zt,Tbef,oliq,oice,cldfraliq,icefrac,rneb',zt(i),Tbef(i),(1-zfice(i))*zoliq(i),zfice(i)*zoliq(i),cldfraliq(i,k),zfice(i),rneb(i,k)
        !ENDIF 
        ENDDO

    ENDIF


    ENDIF ! ok_poprecip

    ! End of precipitation formation
    ! --------------------------------


    ! Calculation of cloud condensates amount seen by radiative scheme
    !-----------------------------------------------------------------

    ! Calculation of the concentration of condensates seen by the radiation scheme
    ! for liquid phase, we take radocondl
    ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true) 
    ! we recalculate radocondi to account for contributions from the precipitation flux
    ! TODO: for the moment, we deactivate this option when ok_poprecip=.TRUE.

    IF ((ok_radocond_snow) .AND. (k < klev) .AND. (.NOT. ok_poprecip)) THEN
        ! for the solid phase (crystals + snowflakes)
        ! we recalculate radocondi 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) /= 0.0) THEN
              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
           ELSE
              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
           ENDIF
           radocond(i,k)=radocondl(i,k)+radocondi(i,k)
        ENDDO
    ENDIF

    ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
    DO i=1,klon
        IF (radocond(i,k) > 0.) THEN
            radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
        ENDIF
    ENDDO

    ! ----------------------------------------------------------------
    ! P4> Wet scavenging
    ! ----------------------------------------------------------------

    !Scavenging through nucleation in the layer

    DO i = 1,klon
        
        IF(zcond(i)>zoliq(i)+1.e-10) THEN
            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
        ELSE
            beta(i,k) = 0.
        ENDIF

        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG

        IF (rneb(i,k)>0.0.AND.zprec_cond(i)>0.) THEN

            IF (temp(i,k) >= t_glace_min) THEN
                zalpha_tr = a_tr_sca(3)
            ELSE
                zalpha_tr = a_tr_sca(4)
            ENDIF

            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 

            ! Nucleation with a factor of -1 instead of -0.5
            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))

        ENDIF

    ENDDO

    ! Scavenging through impaction in the underlying layer

    DO kk = k-1, 1, -1

        DO i = 1, klon

            IF (rneb(i,k)>0.0.AND.zprec_cond(i)>0.) THEN

                IF (temp(i,kk) >= t_glace_min) THEN
                    zalpha_tr = a_tr_sca(1)
                ELSE
                    zalpha_tr = a_tr_sca(2)
                ENDIF

              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi

            ENDIF

        ENDDO

    ENDDO

    !--save some variables for ice supersaturation

    DO i = 1, klon
        ! for memory 
        rneb_seri(i,k) = rneb(i,k)

        ! for diagnostics
        rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)

        qvc(i,k) = zqs(i) * rneb(i,k)
        qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--added by OB because of pathologic cases with the lognormal
        qcld(i,k) = qvc(i,k) + zcond(i)
   ENDDO
   !q_sat
   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.FALSE.,qsatl(:,k),zdqs(:))
   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.FALSE.,qsats(:,k),zdqs(:))



    ! Outputs:
    !-------------------------------
    ! Precipitation fluxes at layer interfaces
    ! + precipitation fractions +
    ! temperature and water species tendencies
    DO i = 1, klon
        psfl(i,k)=zifl(i)
        prfl(i,k)=zrfl(i)
        pfraclr(i,k)=znebprecipclr(i)
        pfracld(i,k)=znebprecipcld(i)
        d_ql(i,k) = (1-zfice(i))*zoliq(i)
        d_qi(i,k) = zfice(i)*zoliq(i)
        d_q(i,k) = zq(i) - qt(i,k)
        ! c_iso: same for isotopes
        d_t(i,k) = zt(i) - temp(i,k)
    ENDDO


ENDDO


  ! Rain or snow at the surface (depending on the first layer temperature)
  DO i = 1, klon
      snow(i) = zifl(i)
      rain(i) = zrfl(i)
      ! c_iso final output
  ENDDO

  IF (ncoreczq>0) THEN
      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
  ENDIF


RETURN

END SUBROUTINE lscp
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

END MODULE lmdz_lscp
