!$gpum horizontal klon ngrid
MODULE lmdz_lscp_main

IMPLICIT NONE

CONTAINS

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE lscp(klon, klev, dtime, missing_val,         &
     paprs, pplay, omega, temp, qt, ql_seri, qi_seri,   &
     ptconv, ratqs, sigma_qtherm,                       &
     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
     pfraclr, pfracld,                                  &
     cldfraliq, cldfraliqth,                            &
     sigma2_icefracturb,sigma2_icefracturbth,           &
     mean_icefracturb,mean_icefracturbth,               &
     radocond, radicefrac, rain, snow,                  &
     frac_impa, frac_nucl, beta,                        &
     prfl, psfl, rhcl, qta, fraca,                      &
     tv, pspsk, tla, thl, wth, iflag_cld_th,            &
     iflag_ice_thermo, distcltop, temp_cltop,           &
     tke, tke_dissip,                                   &
     entr_therm, detr_therm,                            &
     cell_area,                                         &
     cf_seri, rvc_seri, u_seri, v_seri,                 &
     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
     dcf_sub, dcf_con, dcf_mix,          &
     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
     Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi,&
     dqi_avi, dqvc_avi, flight_dist, flight_h2o,        &
     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_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
USE lmdz_lscp_condensation, ONLY : cloudth, cloudth_v3, cloudth_v6, condensation_cloudth
USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld

! Use du module lmdz_lscp_ini contenant les constantes
USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
USE lmdz_lscp_ini, ONLY : min_frac_th_cld, temp_nowater
USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
USE lmdz_lscp_ini, ONLY : ok_lscp_mergecond, gamma_mixth

IMPLICIT NONE

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

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

  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels 
  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
 

  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)   :: omega           ! vertical velocity [Pa/s]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg] 
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ql_seri         ! liquid specific content [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qi_seri         ! ice specific content [kg/kg]
  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke             ! turbulent kinetic energy [m2/s2]
  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke_dissip      ! TKE dissipation [m2/s3]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: entr_therm      ! thermal plume entrainment rate * dz [kg/s/m2] 
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: detr_therm      ! thermal plume detrainment rate * dz [kg/s/m2]


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

  !Inputs associated with thermal plumes

  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds 
  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 potential temperature within thermals [K]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: wth                 ! vertical velocity within thermals [m/s]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: sigma_qtherm        ! controls the saturation deficit distribution width in thermals [-]



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


  ! INPUT/OUTPUT condensation and ice supersaturation
  !--------------------------------------------------
  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]

  ! INPUT/OUTPUT aviation
  !--------------------------------------------------
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
  
  ! 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 fraction [-]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliqth      ! liquid fraction of cloud fraction in thermals [-]
  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)  :: sigma2_icefracturbth ! Variance of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturbth   ! Mean of the diagnostic supersaturation distribution in thermals (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+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 condensation and ice supersaturation

  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat 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)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-]  
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]      
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg]  

  ! for contrails and aviation

  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)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]


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

  ! for thermals

  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


  ! LOCAL VARIABLES:
  !----------------
  REAL, DIMENSION(klon,klev) :: ctot, rnebth, ctot_vol
  REAL, DIMENSION(klon,klev) :: wls                                 !-- large scalce vertical velocity [m/s]
  REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi
  REAL, DIMENSION(klon) :: zqsth, zqslth, zqsith, zdqsth, zdqslth, zdqsith
  REAL :: zdelta
  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,Tbefth,Tbefthm1,qlibef,DT                  ! temperature, humidity and temp. variation during condensation iteration
  REAL :: num,denom
  REAL :: cste
  REAL, DIMENSION(klon) :: qincloud
  REAL, DIMENSION(klon) :: zrfl, zifl
  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqth, zqn, zqnl
  REAL, DIMENSION(klon) :: zoliql, zoliqi
  REAL, DIMENSION(klon) :: zt, zp
  REAL, DIMENSION(klon) :: zfice, zficeth, zficeenv, zneb, zcf, zqi_ini, zsnow
  REAL, DIMENSION(klon) :: dzfice, dzficeth, dzficeenv
  REAL, DIMENSION(klon) :: qtot, zeroklon
  ! Variables precipitation energy conservation
  REAL, DIMENSION(klon) :: zmqc
  REAL :: zalpha_tr
  REAL :: zfrac_lessi
  REAL, DIMENSION(klon) :: zprec_cond
  REAL, DIMENSION(klon) :: zlh_solid
  REAL, DIMENSION(klon) :: ztupnew
  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
  REAL, DIMENSION(klon) :: zrflclr, zrflcld
  REAL, DIMENSION(klon) :: ziflclr, ziflcld
  REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
  REAL, DIMENSION(klon) :: tot_zneb
  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop, zdeltaz
  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl, zqliqth, zqiceth, zqvapclth, sursat_e, invtau_e ! for icefrac_lscp_turb

  ! for quantity of condensates seen by radiation
  REAL, DIMENSION(klon) :: zradocond, zradoice
  REAL, DIMENSION(klon) :: zrho_up, zvelo_up
  
  ! for condensation and ice supersaturation
  REAL, DIMENSION(klon) :: qvc, qvcl, shear
  REAL :: delta_z
  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
  ! Constants used for calculating ratios that are advected (using a parent-child
  ! formalism). This is not done in the dynamical core because at this moment,
  ! only isotopes can use this parent-child formalism. Note that the two constants
  ! are the same as the one use in the dynamical core, being also defined in
  ! dyn3d_common/infotrac.F90
  REAL :: min_qParent, min_ratio

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

  LOGICAL, DIMENSION(klon) :: stratiform_or_all_distrib,pticefracturb
  LOGICAL, DIMENSION(klon) :: keepgoing

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

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

! Few initial checks


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

! 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
rnebth(:,:)=0.0
pfraclr(:,:)=0.0
pfracld(:,:)=0.0
cldfraliq(:,:)=0.
sigma2_icefracturb(:,:)=0.
mean_icefracturb(:,:)=0.
cldfraliqth(:,:)=0.
sigma2_icefracturbth(:,:)=0.
mean_icefracturbth(:,:)=0.
radocond(:,:) = 0.0
radicefrac(:,:) = 0.0
frac_nucl(:,:) = 1.0 
frac_impa(:,:) = 1.0 
rain(:) = 0.0
snow(:) = 0.0
zrfl(:) = 0.0
zifl(:) = 0.0
zneb(:) = seuil_neb
zrflclr(:) = 0.0
ziflclr(:) = 0.0
zrflcld(:) = 0.0
ziflcld(:) = 0.0
tot_zneb(:) = 0.0
zeroklon(:) = 0.0
zdistcltop(:)=0.0
ztemp_cltop(:) = 0.0
ztupnew(:)=0.0
ctot_vol(:,:)=0.0
rneblsvol(:,:)=0.0
znebprecip(:)=0.0
znebprecipclr(:)=0.0
znebprecipcld(:)=0.0
distcltop(:,:)=0.
temp_cltop(:,:)=0.


!--Ice supersaturation
gamma_cond(:,:) = 1.
qissr(:,:)      = 0.
issrfra(:,:)    = 0.
dcf_sub(:,:)    = 0.
dcf_con(:,:)    = 0.
dcf_mix(:,:)    = 0.
dqi_adj(:,:)    = 0.
dqi_sub(:,:)    = 0.
dqi_con(:,:)    = 0.
dqi_mix(:,:)    = 0.
dqvc_adj(:,:)   = 0.
dqvc_sub(:,:)   = 0.
dqvc_con(:,:)   = 0.
dqvc_mix(:,:)   = 0.
fcontrN(:,:)    = 0.
fcontrP(:,:)    = 0.
Tcontr(:,:)     = missing_val
qcontr(:,:)     = missing_val
qcontr2(:,:)    = missing_val
dcf_avi(:,:)    = 0.
dqi_avi(:,:)    = 0.
dqvc_avi(:,:)   = 0.
qvc(:)          = 0.
shear(:)        = 0.
min_qParent     = 1.e-30
min_ratio       = 1.e-16

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

!-- cloud phase useful variables
wls(:,:)=-omega(:,:) / RG / (pplay(:,:)/RD/temp(:,:))
zqliq(:)=0.
zqice(:)=0.
zqvapcl(:)=0.
zqliqth(:)=0.
zqiceth(:)=0.
zqvapclth(:)=0.
sursat_e(:)=0.
invtau_e(:)=0.
pticefracturb(:)=.FALSE.


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

ncoreczq=0

DO k = klev, 1, -1

    IF (k.LE.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)
        zp(i)=pplay(i,k)
        zqi_ini(i)=qi_seri(i,k)
        zcf(i) = 0.
        zfice(i) = 1.0   ! initialized at 1 as by default we assume mpc to be at ice saturation
        dzfice(i) = 0.0
        zficeth(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
        dzficeth(i) = 0.0
        zficeenv(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
        dzficeenv(i) = 0.0
        

        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

    ! --------------------------------------------------------------------
    ! P1> Precipitation processes, before cloud formation:
    !     Thermalization of precipitation falling from the overlying layer AND
    !     Precipitation evaporation/sublimation/melting
    !---------------------------------------------------------------------
   
    !================================================================
    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
    IF ( ok_poprecip ) THEN

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

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

      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
                        zrfl, zrflclr, zrflcld, &
                        zifl, ziflclr, ziflcld, &
                        dqreva(:,k), dqssub(:,k) &
                        )

    ENDIF ! (ok_poprecip)
    
    ! Calculation of qsat,L/cp*dqsat/dT and ncoreczq counter
    !-------------------------------------------------------

     qtot(:)=zq(:)+zmqc(:)
     CALL calc_qsat_ecmwf(klon,zt,qtot,zp,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) .LT. 1.e-15) THEN
                ncoreczq=ncoreczq+1
                zq(i)=1.e-15
            ENDIF
            ! c_iso: do something similar for isotopes

     ENDDO
    
    ! -------------------------------------------------------------------------
    ! P2> Cloud formation including condensation and cloud phase determination
    !--------------------------------------------------------------------------
    !
    ! 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 stratiform_or_all_distrib
        stratiform_or_all_distrib(:)=.TRUE.
        gammasat(:)=1.

        ! part of the code that is supposed to become obsolete soon
        !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        IF (.NOT. ok_lscp_mergecond) THEN
          IF (iflag_cld_th.GE.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.LE.2) THEN
                  ! Old version of Arnaud Jam

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


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


                ELSEIF (iflag_cloudth_vert.EQ.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,                                     &
                         qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
                         ratqs,zqs,temp, &
                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)

                ENDIF


                DO i=1,klon
                    rneb(i,k)=ctot(i,k)
                    ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
                    rneblsvol(i,k)=ctot_vol(i,k)
                    zqn(i)=qincloud(i)
                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
                    qvc(i) = rneb(i,k) * zqs(i)
                ENDDO

               ! Cloud phase final determination for clouds after "old" cloudth calls
               ! for those clouds, only temperature dependent phase partitioning (eventually modulated by
               ! distance to cloud top) is available
               ! compute distance to cloud top when cloud phase is computed depending on temperature
               ! and distance to cloud top
               IF (iflag_t_glace .GE. 4) THEN
                    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)

          ENDIF

          IF (iflag_cld_th .EQ. 5) THEN
               ! When iflag_cld_th=5, we always assume
                ! bi-gaussian distribution
            stratiform_or_all_distrib(:) = .FALSE.

          ELSEIF (iflag_cld_th .GE. 6) THEN
                ! stratiform distribution only when no thermals 
            stratiform_or_all_distrib(:) = fraca(:,k) < min_frac_th_cld

          ENDIF

        ENDIF ! .not. ok_lscp_mergecond
        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        
        DT(:) = 0.
        n_i(:)=0
        Tbef(:)=zt(:)
        qlibef(:)=0.
        Tbefth(:)=tla(:,k)*pspsk(:,k)
        IF (k .GT. 1) THEN
         Tbefthm1(:)=tla(:,k-1)*pspsk(:,k-1)
        ELSE
         Tbefthm1(:)=Tbefth(:)       
        ENDIF
        zqth(:)=qta(:,k)
        zdeltaz(:)=(paprs(:,k)-paprs(:,k+1))/RG/zp(:)*RD*zt(:)

        ! Treatment of stratiform clouds (lognormale or ice-sursat) or all clouds (including cloudth
        ! in case of ok_lscp_mergecond)
        ! We iterate here to take into account the change in qsat(T) and ice fraction
        ! during the condensation process 
        ! the increment in temperature is calculated using the first principle of 
        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
        ! and a clear sky fraction)
        ! note that we assume that the vapor in the cloud is at saturation for this calculation     

        DO iter=1,iflag_fisrtilp_qsat+1

                keepgoing(:) = .FALSE.

                DO i=1,klon

                    ! keepgoing = .true. while convergence is not satisfied

                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. stratiform_or_all_distrib(i)) THEN

                        ! if not convergence:
                        ! we calculate a new iteration
                        keepgoing(i) = .TRUE.

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

                        ! total water including mass of precip
                        qtot(i)=zq(i)+zmqc(i)

                      ENDIF

                  ENDDO

                  ! Calculation of saturation specific humidity 
                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,0,.false.,zqs,zdqs)
                  ! also in thermals
                  CALL calc_qsat_ecmwf(klon,Tbefth,zqth,zp,RTT,0,.false.,zqsth,zdqsth)

                  IF (iflag_icefrac .GE. 1) THEN
                  ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction
                  ! and cloud condensed water content. Principle from Dietlitcher et al. 2018, GMD
                  ! we make this option works only for the physically-based and tke-dependent param from Raillard et al.
                  ! (iflag_icefrac>=1)
                      CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,1,.false.,zqsl,zdqsl)
                      CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,2,.false.,zqsi,zdqsi)
                      DO i=1,klon
                           zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i)
                           zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i)
                      ENDDO
                  ENDIF
                  IF (iflag_icefrac .GE. 2) THEN
                  ! same adjustment for thermals
                      CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,1,.false.,zqslth,zdqslth)
                      CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,2,.false.,zqsith,zdqsith)
                      DO i=1,klon
                           zqsth(i)=zficeth(i)*zqsith(i)+(1.-zficeth(i))*zqslth(i)
                           zdqsth(i)=zficeth(i)*zdqsith(i)+zqsith(i)*dzficeth(i)+(1.-zficeth(i))*zdqslth(i)-zqslth(i)*dzficeth(i)
                      ENDDO
                  ENDIF

                  CALL calc_gammasat(klon,Tbef,qtot,zp,gammasat,dgammasatdt)
                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)

                  ! Cloud condensation based on subgrid pdf
                  !--AB Activates a condensation scheme that allows for
                  !--ice supersaturation and contrails evolution from aviation
                  IF (ok_ice_supersat) THEN

                    !--Calculate the shear value (input for condensation and ice supersat)
                    DO i = 1, klon
                      !--Cell thickness [m]
                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
                      IF ( iftop ) THEN
                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
                      ELSEIF ( k .EQ. 1 ) THEN
                        ! surface
                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
                      ELSE
                        ! other layers
                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
                      ENDIF
                    ENDDO

                    !---------------------------------------------
                    !--   CONDENSATION AND ICE SUPERSATURATION  --
                    !---------------------------------------------

                    CALL condensation_ice_supersat( &
                        klon, dtime, missing_val, &
                        zp, paprs(:,k), paprs(:,k+1), &
                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
                        shear, tke_dissip(:,k), cell_area, &
                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
                        flight_dist(:,k), flight_h2o(:,k), &
                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))


                  ELSE
                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)

                  CALL condensation_lognormal( &
                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
                       keepgoing, rneb(:,k), zqn, qvc)


                  ENDIF ! .NOT. ok_ice_supersat

                  ! Volume cloud fraction
                  rneblsvol(:,k)=rneb(:,k)


                  IF  (ok_lscp_mergecond) THEN
                      ! in that case we overwrite rneb, rneblsvol and zqn for shallow convective clouds only 
                      CALL condensation_cloudth(klon,Tbef,zq,zqth,fraca(:,k), &
                                     pspsk(:,k),zp,tla(:,k), &
                                     ratqs(:,k),sigma_qtherm(:,k),zqsth,zqs,zqn,rneb(:,k),rnebth(:,k),rneblsvol(:,k), &
                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
                  ENDIF



                  ! Cloud phase determination


                  IF (iflag_icefrac .LE. 1) THEN
                     ! "old" phase partitioning depending on temperature and eventually distance to cloud top everywhere
                     IF (iflag_t_glace .GE. 4) THEN
                        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)
                  ENDIF

                  IF (iflag_icefrac .EQ. 1) THEN
                     ! phase partitioning depending on turbulence, vertical velocity and ice crystal microphysics 
                     ! only in stratiform clouds in the mixed phase regime (Raillard et al. 2025)
                     ! it overwrites temperature-dependent phase partitioning except for convective boundary layer clouds
                     pticefracturb(:) = (fraca(:,k) < min_frac_th_cld) .AND. (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT)
                     DO i=1,klon
                        qincloud(i)=zqn(i)
                        zcf(i)=MIN(MAX(rneb(i,k), 0.),1.)
                        sursat_e(i) = 0.
                        invtau_e(i) = 0.
                     ENDDO
                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), zqi_ini,    &
                     ziflcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, & 
                     zficeenv, dzficeenv, cldfraliq(:,k),sigma2_icefracturb(:,k),mean_icefracturb(:,k))                     
                     DO i=1,klon
                        IF (pticefracturb(i)) THEN
                          zfice(i)=zficeenv(i)
                          dzfice(i)=dzficeenv(i)
                        ENDIF
                     ENDDO
                   

                  ELSEIF (iflag_icefrac .GE. 2) THEN
                     ! compute also phase partitioning also in thermal clouds (neglecting tke in thermals as first assumption)
                     ! moreover, given the upward velocity of thermals, we assume that precipitation falls in the environment
                     ! hence we assume that no snow falls in thermals.
                     ! we activate only in the mixed phase regime not to distrub the cirrus param at very cold T
                     pticefracturb(:) = (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT)
                     !Thermals
                     DO i=1,klon
                        IF (fraca(i,k) .GT. min_frac_th_cld) THEN
                           zcf(i)=MIN(MAX(rnebth(i,k),0.), 1.)/fraca(i,k)
                           qincloud(i)=zqn(i)*fraca(i,k)
                        ELSE
                           zcf(i) = 0.
                           qincloud(i) = 0.
                        ENDIF 
                        sursat_e(i)=cloudth_senv(i,k)/zqsi(i)
                        invtau_e(i)=gamma_mixth*MAX(entr_therm(i,k)-detr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i)
                     ENDDO
                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbefth, zp, paprs(:,k), paprs(:,k+1), wth(:,k), zqi_ini,  &
                     zeroklon, qincloud, zcf, zeroklon, zeroklon, sursat_e, invtau_e, zqliqth, zqvapclth, zqiceth, &
                     zficeth, dzficeth,cldfraliqth(:,k), sigma2_icefracturbth(:,k), mean_icefracturbth(:,k))
                     !Environment
                     DO i=1,klon
                        qincloud(i)=zqn(i)*(1.-fraca(i,k))
                        zcf(i)=MIN(MAX(rneb(i,k)-rnebth(i,k), 0.),1.)/(1.-fraca(i,k))
                        IF (k .GT. 1) THEN
                           ! evaluate the mixing sursaturation using saturation deficit at level below
                           ! as air pacels detraining into clouds have not (less) seen yet entrainement from above
                           sursat_e(i)=cloudth_sth(i,k-1)/(zqsith(i)+zdqsith(i)*RCPD/RLSTT*(Tbefthm1(i)-Tbefth(i)))
                           ! mixing is assumed to scales with intensity of net detrainment/entrainment rate (D/dz-E/dz) / rho
                           invtau_e(i)=gamma_mixth*MAX(detr_therm(i,k)-entr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i)
                        ELSE
                           sursat_e(i)=0.
                           invtau_e(i)=0.
                        ENDIF
                     ENDDO
                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), zqi_ini,    &
                     ziflcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, &
                     zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k)) 
  
                    ! adjust zfice to account for condensates in thermals'fraction
                     DO i=1,klon
                        IF ( zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i) .GT. 0.) THEN
                          zfice(i)=MIN(1., (zqiceth(i)+zqice(i))/(zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i)))
                          dzfice(i)=0. ! dxice/dT=0. always when using icefrac_lscp_turb
                        ENDIF
                     ENDDO

                  ENDIF ! iflag_icefrac




                  DO i=1,klon
                      IF (keepgoing(i)) THEN

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

                
                        IF (zfice(i).LT.1) THEN
                            cste=RLVTT
                        ELSE
                            cste=RLSTT
                        ENDIF
                        
                        IF ( ok_unadjusted_clouds ) THEN
                          !--AB We relax the saturation adjustment assumption
                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
                          IF ( rneb(i,k) .GT. eps ) THEN
                            qlibef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
                          ELSE
                            qlibef(i) = 0.
                          ENDIF
                        ELSE
                          qlibef(i)=max(0.,zqn(i)-zqs(i))
                        ENDIF

                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlibef(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)      &
                              *qlibef(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.2> 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
        ! ----------------------------------------------------------------


        ! Water vapor and condensed water update, subsequent latent heat exchange for each cloud type
        !---------------------------------------------------------------------------------------------
        DO i=1, klon
                ! Checks on rneb, rhcl and zqn
                IF (rneb(i,k) .LE. 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) .GE. 1.0) THEN
                    zqn(i) = zq(i)
                    rneb(i,k) = 1.0
                    IF ( ok_unadjusted_clouds ) THEN
                      !--AB We relax the saturation adjustment assumption
                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
                      zcond(i) = MAX(0., zqn(i) - qvc(i))
                    ELSE
                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
                    ENDIF
                    rhcl(i,k)=1.0
                ELSE
                    IF ( ok_unadjusted_clouds ) THEN
                      !--AB We relax the saturation adjustment assumption
     
                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
                    ELSE
                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
                    ENDIF

                    ! following line is very strange and probably wrong
                    rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
                    ! Correct calculation of clear-sky relative humidity. To activate once
                    ! issues related to zqn>zq in bi-gaussian clouds are corrected
                    !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
                    !--This is slighly approximated, the actual formula is
                    !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
                    !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
                    !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
                    ! rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)

                ENDIF

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

                       
                ! 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

    ! ----------------------------------------------------------------
    ! P3> Precipitation processes, after cloud formation 
    !     - precipitation formation, melting/freezing
    ! ----------------------------------------------------------------

    ! Initiate the quantity of liquid and solid condensates
    ! Note that in the following, zcond is the total amount of condensates
    ! including newly formed precipitations (i.e., condensates formed by the
    ! condensation process), while zoliq is the total amount of condensates in
    ! the cloud (i.e., on which precipitation processes have applied)
    DO i = 1, klon
      zoliq(i)  = zcond(i)
      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
      zoliqi(i) = zcond(i) * zfice(i)
    ENDDO

    !================================================================
    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
    IF (ok_poprecip) THEN

      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), zp, &
                            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)
      ENDDO

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

      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
                            ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
                            rneb(:,k), znebprecipclr, znebprecipcld, &
                            zneb, tot_zneb, zrho_up, zvelo_up, &
                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
                            )

    ENDIF ! ok_poprecip

    ! End of precipitation processes after cloud formation
    ! ----------------------------------------------------

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

    DO i=1,klon

      IF (ok_poprecip) THEN
        IF (ok_radocond_snow) THEN
          radocond(i,k) = zoliq(i)
          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
        ELSE
          radocond(i,k) = zoliq(i)
          zradoice(i) = zoliqi(i)
        ENDIF
      ELSE
        radocond(i,k) = zradocond(i)
      ENDIF

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

    ! ----------------------------------------------------------------
    ! P5> Wet scavenging
    ! ----------------------------------------------------------------

    !Scavenging through nucleation in the layer

    DO i = 1,klon
        
        IF(zcond(i).GT.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).GT.0.0.AND.zprec_cond(i).GT.0.) THEN

            IF (temp(i,k) .GE. temp_nowater) THEN
                zalpha_tr = a_tr_sca(3)
            ELSE
                zalpha_tr = a_tr_sca(4)
            ENDIF

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

            ! Nucleation with a factor of -1 instead of -0.5
            zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))

        ENDIF

    ENDDO

    ! Scavenging through impaction in the underlying layer

    DO kk = k-1, 1, -1

        DO i = 1, klon

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

                IF (temp(i,kk) .GE. temp_nowater) THEN
                    zalpha_tr = a_tr_sca(1)
                ELSE
                    zalpha_tr = a_tr_sca(2)
                ENDIF

              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
              frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi

            ENDIF

        ENDDO

    ENDDO
    
    !------------------------------------------------------------
    ! P6 > write diagnostics and outputs
    !------------------------------------------------------------
   
    !--AB Write diagnostics and tracers for ice supersaturation
    IF ( ok_ice_supersat ) THEN
      CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,1,.false.,qsatl(:,k),zdqs)
      CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,2,.false.,qsati(:,k),zdqs)

      DO i = 1, klon

        IF ( zoliq(i) .LE. 0. ) THEN
          !--If everything was precipitated, the remaining empty cloud is dissipated
          !--and everything is transfered to the subsaturated clear sky region
          rneb(i,k) = 0.
          qvc(i) = 0.
        ENDIF

        cf_seri(i,k) = rneb(i,k)

        IF ( .NOT. ok_unadjusted_clouds ) THEN
          qvc(i) = zqs(i) * rneb(i,k)
        ENDIF
        IF ( zq(i) .GT. min_qParent ) THEN
          rvc_seri(i,k) = qvc(i) / zq(i)
        ELSE
          rvc_seri(i,k) = min_ratio
        ENDIF
        !--The MIN barrier is NEEDED because of:
        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
        !--The MAX barrier is a safeguard that should not be activated
        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)

        !--Diagnostics
        gamma_cond(i,k) = gammasat(i)
        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
        qcld(i,k) = qvc(i) + zoliq(i)
      ENDDO
    ENDIF

    ! 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)
        distcltop(i,k)=zdistcltop(i)
        temp_cltop(i,k)=ztemp_cltop(i)
        d_q(i,k) = zq(i) - qt(i,k)
        d_t(i,k) = zt(i) - temp(i,k)

        IF (ok_bug_phase_lscp) THEN
           d_ql(i,k) = (1-zfice(i))*zoliq(i)
           d_qi(i,k) = zfice(i)*zoliq(i)
        ELSE
           d_ql(i,k) = zoliql(i)
           d_qi(i,k) = zoliqi(i)    
        ENDIF

    ENDDO


ENDDO ! loop on k from top to bottom


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