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

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Contacts: JB Madeleine jbm@lmd.ipsl.fr ;  E Vignon etienne.vignon@lmd.ipsl.fr
!--------------------------------------------------------------------------------
! First version: 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)
!
!--------------------------------------------------------------------------------
! References:
!
! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
! - Borella et al. 2025, JAMES, doi:10.1029/2024MS004918
! - Raillard et al. 2026, JAMES, doi:10.22541/essoar.175096287.71557703/v1
! - Vignon et al. 2026, ACP, doi:0.5194/egusphere-2025-4641
! -------------------------------------------------------------------------------
! Code structure:
!
! the formation of cloud and precipitation is treated as a top-down process,
! starting from the top of the atmosphere down to the ground
!
! P1>     Precipitation processes before cloud formation
! P2>     Cloud formation and cloud phase determination
! P3>     Precipitation formation and interaction with clouds
! P4>     Computation of cloud water seen by radiation
! P5>     Wet scavenging
! P6>     Diagnostics and outputs
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!

! USE de modules contenant des fonctions.
      USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat, distance_to_cloud_top
      USE lmdz_lscp_phase, ONLY: icefrac_lscp, icefrac_lscp_turb
      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: iflag_cld_th, iflag_ice_thermo
      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

      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

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

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

         ! 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))
            END IF
            !c_iso init of iso
         END DO

         ! --------------------------------------------------------------------
         ! 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 (poprecip) 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) &
                                   )

         END IF ! (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
            END IF
            ! c_iso: do something similar for isotopes

         END DO

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

         ! initialise gammasat and stratiform_or_all_distrib
         stratiform_or_all_distrib(:) = .TRUE.
         gammasat(:) = 1.

        !!!! part of the code that is obsolete and should be deleted 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)

               END IF

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

               ! 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)
               END IF
               CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)

            END IF

            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

            END IF

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

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

               END IF

            END DO

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

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

               ! P2.1> 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)

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

            ! P2.2> 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)
               END IF
               CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
            END IF

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

            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(:) = .TRUE.
               !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.
                  END IF
                  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)
               END DO
               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.
                  END IF
               END DO
               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
                  END IF
               END DO

            END IF ! iflag_icefrac

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

                  ! Approximative calculation of temperature variation DT due to condensation.
                  !----------------------------------------------------------------------------

                  IF (zfice(i) .LT. 1) THEN
                     cste = RLVTT
                  ELSE
                     cste = RLSTT
                  END IF

                  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.
                     END IF
                  ELSE
                     qlibef(i) = max(0., zqn(i) - zqs(i))
                  END IF

                  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

               END IF ! end keepgoing

            END DO     ! end loop on i

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

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

            END IF

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

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

         !================================================================
         ! Flag for the new and more microphysical treatment of precipitation (poprecip) 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)
            END DO

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

         END IF ! 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)
               END IF
            ELSE
               ! During the autoconversion to precipitation with the historical precipitation scheme (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
               radocond(i, k) = zradocond(i)
            END IF

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

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

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

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

            END IF

         END DO

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

                  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

               END IF

            END DO

         END DO

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

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

               IF (.NOT. ok_unadjusted_clouds) THEN
                  qvc(i) = zqs(i)*rneb(i, k)
               END IF
               IF (zq(i) .GT. min_qParent) THEN
                  rvc_seri(i, k) = qvc(i)/zq(i)
               ELSE
                  rvc_seri(i, k) = min_ratio
               END IF
               !--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)
            END DO
         END IF

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

         END DO

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

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

      RETURN

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

END MODULE lmdz_lscp_main
