!$gpum horizontal klon
MODULE lmdz_lscp_phase

   IMPLICIT NONE
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! This module contains routines that determine
! the phase of the clouds computed in lscp
! routines
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

CONTAINS

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      ! Compute the ice fraction (see e.g. Doutriaux-Boucher & Quaas 2004, section 2.2.)
      ! as a function of temperature
      ! see also Fig 3 of Madeleine et al. 2020, JAMES, 10.1029/2020MS002046
      ! An option is also available to make the cloud phase depend on distance from cloud
      ! top. See Lea Raillard's PhD manuscript
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      USE print_control_mod, ONLY: lunout, prt_level
      USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
      USE lmdz_lscp_ini, ONLY: RTT, dist_liq, temp_nowater

      IMPLICIT NONE

      INTEGER, INTENT(IN)                 :: klon              !-- number of horizontal grid points
      INTEGER, INTENT(IN)                 :: iflag_ice_thermo  !-- option for cloud phase determination
      REAL, INTENT(IN), DIMENSION(klon)   :: temp              !-- temperature [K]
      REAL, INTENT(IN), DIMENSION(klon)   :: distcltop         !-- distance to cloud top [m]
      REAL, INTENT(IN), DIMENSION(klon)   :: temp_cltop        !-- temperature of cloud top [K]
      REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac           !-- ice fraction in clouds [0-1]
      REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT        !-- dicefraction/dT [K-1]

      INTEGER i
      REAL liqfrac_tmp, dicefrac_tmp
      REAL Dv, denomdep, beta, qsi, dqsidt
      LOGICAL ice_thermo

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

      IF ((iflag_t_glace .LT. 2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN
         abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'
         CALL abort_physic(modname, abort_message, 1)
      END IF

      IF (.NOT. ((iflag_ice_thermo .EQ. 1) .OR. (iflag_ice_thermo .GE. 3))) THEN
         abort_message = 'lscp cannot be used without ice thermodynamics'
         CALL abort_physic(modname, abort_message, 1)
      END IF

      DO i = 1, klon

         ! old function with sole dependence upon temperature
         IF (iflag_t_glace .EQ. 2) THEN
            liqfrac_tmp = (temp(i) - t_glace_min)/(t_glace_max - t_glace_min)
            liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0)
            icefrac(i) = (1.0 - liqfrac_tmp)**exposant_glace
            IF (icefrac(i) .GT. 0.) THEN
               dicefracdT(i) = exposant_glace*(icefrac(i)**(exposant_glace - 1.)) &
                               /(t_glace_min - t_glace_max)
            END IF

            IF ((icefrac(i) .EQ. 0) .OR. (icefrac(i) .EQ. 1)) THEN
               dicefracdT(i) = 0.
            END IF

         END IF

         ! function of temperature used in CMIP6 physics
         IF (iflag_t_glace .EQ. 3) THEN
            liqfrac_tmp = (temp(i) - t_glace_min)/(t_glace_max - t_glace_min)
            liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0)
            icefrac(i) = 1.0 - liqfrac_tmp**exposant_glace
            IF ((icefrac(i) .GT. 0.) .AND. (liqfrac_tmp .GT. 0.)) THEN
               dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace - 1.)) &
                               /(t_glace_min - t_glace_max)
            ELSE
               dicefracdT(i) = 0.
            END IF
         END IF

         ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top
         ! and then decreases with decreasing height

         !with linear function of temperature at cloud top
         IF (iflag_t_glace .EQ. 4) THEN
            liqfrac_tmp = (temp(i) - t_glace_min)/(t_glace_max - t_glace_min)
            liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0)
            icefrac(i) = MAX(MIN(1., 1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)), 0.)
            dicefrac_tmp = -temp(i)/(t_glace_max - t_glace_min)
            dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)
            IF ((liqfrac_tmp .LE. 0) .OR. (liqfrac_tmp .GE. 1)) THEN
               dicefracdT(i) = 0.
            END IF
         END IF

         ! with CMIP6 function of temperature at cloud top
         IF ((iflag_t_glace .EQ. 5) .OR. (iflag_t_glace .EQ. 7)) THEN
            liqfrac_tmp = (temp(i) - t_glace_min)/(t_glace_max - t_glace_min)
            liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0)
            liqfrac_tmp = liqfrac_tmp**exposant_glace
            icefrac(i) = MAX(MIN(1., 1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)), 0.)
            IF ((liqfrac_tmp .LE. 0) .OR. (liqfrac_tmp .GE. 1)) THEN
               dicefracdT(i) = 0.
            ELSE
               dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace - 1.))/(t_glace_min - t_glace_max) &
                               *exp(-distcltop(i)/dist_liq)
            END IF
         END IF

         ! with modified function of temperature at cloud top
         ! to get largere values around 260 K, works well with t_glace_min = 241K
         IF (iflag_t_glace .EQ. 6) THEN
            IF (temp(i) .GT. t_glace_max) THEN
               liqfrac_tmp = 1.
            ELSE
               liqfrac_tmp = -((temp(i) - t_glace_max)/(t_glace_max - t_glace_min))**2 + 1.
            END IF
            liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0)
            icefrac(i) = MAX(MIN(1., 1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)), 0.)
            IF ((liqfrac_tmp .LE. 0) .OR. (liqfrac_tmp .GE. 1)) THEN
               dicefracdT(i) = 0.
            ELSE
               dicefracdT(i) = 2*((temp(i) - t_glace_max)/(t_glace_max - t_glace_min))/(t_glace_max - t_glace_min) &
                               *exp(-distcltop(i)/dist_liq)
            END IF
         END IF

         ! if temperature or temperature of cloud top <-40°C,
         IF (iflag_t_glace .GE. 4) THEN
            IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN
               icefrac(i) = 1.
               dicefracdT(i) = 0.
            END IF
         END IF

      END DO ! klon

      RETURN

   END SUBROUTINE ICEFRAC_LSCP
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, pticefracturb, temp, pplay, paprsdn, paprsup, wvel, qice_ini, snowcld, qtot_incl, cldfra, tke,   &
         tke_dissip, sursat_e, invtau_e, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      ! Compute the liquid, ice and vapour content (+ice fraction) based
      ! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)
      ! L.Raillard (23/09/24)
      ! E.Vignon (03/2025) : additional elements for treatment of convective
      !                      boundary layer clouds
      ! References:
      ! - Raillard et al. 2026, JAMES, doi:10.22541/essoar.175096287.71557703/v1
      ! - Vignon et al. 2026, ACP, doi:0.5194/egusphere-2025-4641
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      USE lmdz_lscp_ini, ONLY: prt_level, lunout
      USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI
      USE lmdz_lscp_ini, ONLY: seuil_neb, temp_nowater
      USE lmdz_lscp_ini, ONLY: naero5, gamma_snwretro, gamma_taud, capa_crystal, rho_ice
      USE lmdz_lscp_ini, ONLY: eps, snow_fallspeed, fix_icnc, icnc0
      USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf

      IMPLICIT NONE

      INTEGER, INTENT(IN)                           :: klon                !--number of horizontal grid points
      REAL, INTENT(IN)                           :: dtime               !--time step [s]
      LOGICAL, INTENT(IN), DIMENSION(klon)    :: pticefracturb       !--grid points concerned by this routine
      REAL, INTENT(IN), DIMENSION(klon)    :: temp                !--temperature
      REAL, INTENT(IN), DIMENSION(klon)    :: pplay               !--pressure in the middle of the layer           [Pa]
      REAL, INTENT(IN), DIMENSION(klon)    :: paprsdn             !--pressure at the bottom interface of the layer [Pa]
      REAL, INTENT(IN), DIMENSION(klon)    :: paprsup             !--pressure at the top interface of the layer    [Pa]
      REAL, INTENT(IN), DIMENSION(klon)    :: wvel                !--vertical velocity                             [m/s]
      REAL, INTENT(IN), DIMENSION(klon)    :: qtot_incl           !--specific total cloud water in-cloud content   [kg/kg]
      REAL, INTENT(IN), DIMENSION(klon)    :: cldfra              !--cloud fraction in gridbox                     [-]
      REAL, INTENT(IN), DIMENSION(klon)    :: tke                 !--turbulent kinetic energy                      [m2/s2]
      REAL, INTENT(IN), DIMENSION(klon)    :: tke_dissip          !--TKE dissipation                               [m2/s3]

      REAL, INTENT(IN), DIMENSION(klon)    :: qice_ini            !--initial specific ice content gridbox-mean     [kg/kg]
      REAL, INTENT(IN), DIMENSION(klon)    :: snowcld             !--in-cloud snowfall flux                        [kg/m2/s]
      REAL, INTENT(IN), DIMENSION(klon)    :: sursat_e            !--environment supersaturation                   [-]
      REAL, INTENT(IN), DIMENSION(klon)    :: invtau_e            !--inverse time-scale of mixing with environment [s-1]
      REAL, INTENT(OUT), DIMENSION(klon)    :: qliq                !--specific liquid content gridbox-mean          [kg/kg]
      REAL, INTENT(OUT), DIMENSION(klon)    :: qvap_cld            !--specific cloud vapor content, gridbox-mean    [kg/kg]
      REAL, INTENT(OUT), DIMENSION(klon)    :: qice                !--specific ice content gridbox-mean             [kg/kg]

      REAL, INTENT(INOUT), DIMENSION(klon)    :: icefrac             !--fraction of ice in condensed water            [-]
      REAL, INTENT(INOUT), DIMENSION(klon)    :: dicefracdT

      REAL, INTENT(OUT), DIMENSION(klon)    :: cldfraliq           !--fraction of cldfra where liquid               [-]
      REAL, INTENT(OUT), DIMENSION(klon)    :: sigma2_icefracturb  !--Sigma2 of the ice supersaturation PDF         [-]
      REAL, INTENT(OUT), DIMENSION(klon)    :: mean_icefracturb    !--Mean of the ice supersaturation PDF           [-]

      REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati           !--specific humidity saturation values
      INTEGER :: i

      REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl                  !--In-cloud specific quantities                  [kg/kg]
      REAL :: water_vapor_diff                                               !--Water-vapour diffusion coeff in air (f(T,P))  [m2/s]
      REAL :: air_thermal_conduct                                            !--Thermal conductivity of air (f(T))            [J/m/K/s]
      REAL :: C0                                                             !--Lagrangian structure function                 [-]
      REAL :: tau_dissipturb
      REAL :: invtau_phaserelax
      REAL :: sigma2_pdf
      REAL :: ai, bi, B0
      REAL :: sursat_iceliq
      REAL :: sursat_equ
      REAL :: liqfra_max
      REAL :: sursat_iceext
      REAL :: nb_crystals                                                    !--number concentration of ice crystals          [#/m3]
      REAL :: moment1_PSD                                                    !--1st moment of ice PSD
      REAL :: N0_PSD, lambda_PSD                                             !--parameters of the exponential PSD

      REAL :: cldfra1D
      REAL :: rho_air
      REAL :: psati                                                          !--saturation vapor pressure wrt ice             [Pa]
      REAL :: sigmaw2                                                        !--variance of vertical turbulent velocity       [m2/s2]

      REAL :: tempvig1, tempvig2

      tempvig1 = -21.06 + RTT
      tempvig2 = -30.35 + RTT
      C0 = 10.                                                    !--value assumed in Field2014
      sursat_iceext = -0.1
      qzero(:) = 0.
      cldfraliq(:) = 0.
      qliq(:) = 0.
      qice(:) = 0.
      qvap_cld(:) = 0.
      sigma2_icefracturb(:) = 0.
      mean_icefracturb(:) = 0.

      !--wrt liquid
      CALL calc_qsat_ecmwf(klon, temp(:), qzero(:), pplay(:), RTT, 1, .false., qsatl(:), dqsatl(:))
      !--wrt ice
      CALL calc_qsat_ecmwf(klon, temp(:), qzero(:), pplay(:), RTT, 2, .false., qsati(:), dqsati(:))

      DO i = 1, klon
         rho_air = pplay(i)/temp(i)/RD
         ! assuming turbulence isotropy, tke=3/2*sigmaw2
         sigmaw2 = 2./3*tke(i)
         ! because cldfra is intent in, but can be locally modified due to test
         cldfra1D = MIN(cldfra(i),1.)
         ! T>0°C, no ice allowed
         IF (temp(i) .GE. RTT) THEN
            qvap_cld(i) = qsatl(i)*cldfra1D
            qliq(i) = MAX(0.0, qtot_incl(i) - qsatl(i))*cldfra1D
            qice(i) = 0.
            cldfraliq(i) = 1.
            icefrac(i) = 0.
            dicefracdT(i) = 0.

            ! T<-38°C, no liquid allowed
         ELSE IF (temp(i) .LE. temp_nowater) THEN
            qvap_cld(i) = qsati(i)*cldfra1D
            qliq(i) = 0.
            qice(i) = MAX(0.0, qtot_incl(i) - qsati(i))*cldfra1D
            cldfraliq(i) = 0.
            icefrac(i) = 1.
            dicefracdT(i) = 0.
         ELSE 
            !---------------------------------------------------------
            !--             MIXED PHASE TEMPERATURE REGIME
            !---------------------------------------------------------
            !--In the mixed phase regime (-38°C< T <0°C) we distinguish
            !--3 possible subcases.
            !--1.  No pre-existing ice
            !--2A. Pre-existing ice and no turbulence
            !--2B. Pre-existing ice and turbulence

            ! gamma_snwretro controls the contribution of snowflakes to the negative feedback
            ! note that for reasons related to inetarctions with the condensation iteration in lscp_main
            ! we consider here the mean snowflake concentration in the mesh (not the in-cloud concentration)
            ! when poprecip is active, it will be worth testing considering the incloud fraction, dividing
            ! by znebprecipcld
            ! qiceini_incl  = qice_ini(i) / cldfra1D + &
            !              gamma_snwretro * snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) )
            ! assuming constant snowfall velocity

            ! activate param for concerned grid points and for cloudy conditions
            IF ((pticefracturb(i)) .AND. (cldfra(i) .GT. 0.)) THEN
               
               qiceini_incl = qice_ini(i)/cldfra1D + gamma_snwretro*snowcld(i)/pplay(i)*RD*temp(i)/snow_fallspeed

               !--1. No preexisting ice and no mixing with environment: if vertical motion, fully liquid
               !--cloud else fully iced cloud
               IF ((qiceini_incl .LT. eps) .AND. (invtau_e(i) .LT. eps)) THEN
                  IF ((wvel(i) + sqrt(sigmaw2) .GT. eps) .OR. (tke(i) .GT. eps)) THEN
                     qvap_cld(i) = qsatl(i)*cldfra1D
                     qliq(i) = MAX(0., qtot_incl(i) - qsatl(i))*cldfra1D
                     qice(i) = 0.
                     cldfraliq(i) = 1.
                     icefrac(i) = 0.
                     dicefracdT(i) = 0.
                  ELSE
                     qvap_cld(i) = qsati(i)*cldfra1D
                     qliq(i) = 0.
                     qice(i) = MAX(0., qtot_incl(i) - qsati(i))*cldfra1D
                     cldfraliq(i) = 0.
                     icefrac(i) = 1.
                     dicefracdT(i) = 0.
                  END IF

                  !--2. Pre-existing ice and/or mixing with environment:computation of ice properties for
                  !--feedback
               ELSE

                  sursat_iceliq = qsatl(i)/qsati(i) - 1.
                  psati = qsati(i)*pplay(i)/(RD/RV)

                  !--We assume an exponential ice PSD whose parameters
                  !--are computed following Morrison&Gettelman 2008
                  !--Ice number density is assumed equals to a constant value
                  !--if fix_icnc = True
                  !--or to the INP density
                  !--which is for naero5>0 a function of temperature (DeMott 2010)
                  !--bi and B0 are microphysical function characterizing
                  !--vapor/ice interactions
                  !--tau_phase_relax is the typical time of vapor deposition
                  !--onto ice crystals

                  !--For naero5<=0 INP density is derived from the empirical fit
                  !--from MARCUS campaign from Vignon 2021
                  !--/!\ Note that option is very specific and should be use for
                  !--the Southern Ocean and the Antarctic

                  IF (fix_icnc) THEN 
                  
                     nb_crystals = icnc0 * 1000.! from L-1 to m-3
                     
                  ELSE

                     IF (naero5 .LE. 0) THEN
                        IF (temp(i) .GT. tempvig1) THEN
                           nb_crystals = 1.e3*10**(-0.14*(temp(i) - tempvig1) - 2.88)
                        ELSE IF (temp(i) .GT. tempvig2) THEN
                           nb_crystals = 1.e3*10**(-0.31*(temp(i) - tempvig1) - 2.88)
                        ELSE
                           nb_crystals = 1.e3*10**(0.)
                        END IF
                     ELSE
                        nb_crystals = 1.e3*5.94e-5*(RTT - temp(i))**3.33*naero5**(0.0264*(RTT - temp(i)) + 0.0033)
                     END IF
                  
                  ENDIF

                  lambda_PSD = ((RPI*rho_ice*nb_crystals)/(rho_air*MAX(qiceini_incl, eps)))**(1./3.)
                  N0_PSD = nb_crystals*lambda_PSD
                  moment1_PSD = N0_PSD/lambda_PSD**2

                  !--Formulae for air thermal conductivity and water vapor diffusivity
                  !--comes respectively from Beard and Pruppacher (1971)
                  !--and  Hall and Pruppacher (1976)

                  air_thermal_conduct = (5.69 + 0.017*(temp(i) - RTT))*1.e-3*4.184
                  water_vapor_diff = 2.11*1e-5*(temp(i)/RTT)**1.94*(101325/pplay(i))

                  bi = 1./((qsati(i) + qsatl(i))/2.) + RLSTT**2/RCPD/RV/temp(i)**2
                  B0 = 4.*RPI*capa_crystal*1./(RLSTT**2/air_thermal_conduct/RV/temp(i)**2 &
                                               + RV*temp(i)/psati/water_vapor_diff)
                  invtau_phaserelax = bi*B0*moment1_PSD

                  ai = RG/RD/temp(i)*(RD*RLSTT/RCPD/RV/temp(i) - 1.)
                  sursat_equ = (ai*wvel(i) + sursat_e(i)*invtau_e(i))/(invtau_phaserelax + invtau_e(i))
                  ! as sursaturation is by definition lower than -1 and
                  ! because local supersaturation > 1 are never found in the atmosphere

                  !--2A. No TKE : stationnary binary solution depending on vertical velocity and mixing with env.
                  ! If Sequ > Siw liquid cloud, else ice cloud
                  IF (tke_dissip(i) .LE. eps) THEN
                     sigma2_icefracturb(i) = 0.
                     mean_icefracturb(i) = sursat_equ
                     IF (sursat_equ .GT. sursat_iceliq) THEN
                        qvap_cld(i) = qsatl(i)*cldfra1D
                        qliq(i) = MAX(0., qtot_incl(i) - qsatl(i))*cldfra1D
                        qice(i) = 0.
                        cldfraliq(i) = 1.
                        icefrac(i) = 0.
                        dicefracdT(i) = 0.
                     ELSE
                        qvap_cld(i) = qsati(i)*cldfra1D
                        qliq(i) = 0.
                        qice(i) = MAX(0., qtot_incl(i) - qsati(i))*cldfra1D
                        cldfraliq(i) = 0.
                        icefrac(i) = 1.
                        dicefracdT(i) = 0.
                     END IF

                     !--2B. TKE and ice : ice supersaturation PDF
                     !--we compute the cloud liquid properties with a Gaussian PDF
                     !--of ice supersaturation F(Si) (Field2014, Furtado2016).
                     !--Parameters of the PDF are function of turbulence and
                     !--microphysics/existing ice.
                  ELSE

                     !--Tau_dissipturb is the time needed for turbulence to decay
                     !--due to viscosity
                     tau_dissipturb = gamma_taud*2.*sigmaw2/tke_dissip(i)/C0

                     !--------------------- PDF COMPUTATIONS ---------------------
                     !--Formulae for sigma2_pdf (variance), mean of PDF in Raillard2025
                     !--cloud liquid fraction and in-cloud liquid content are given
                     !--by integrating resp. F(Si) and Si*F(Si)
                     !--Liquid is limited by the available water vapor trough a
                     !--maximal liquid fraction
                     !--qice_ini(i) / cldfra1D = qiceincld without precip

       liqfra_max = MAX(0., (MIN(1., (qtot_incl(i) - (qice_ini(i)/cldfra1D) - qsati(i)*(1 + sursat_iceext))/(qsatl(i) - qsati(i)))))
                     sigma2_pdf = 1./2.*(ai**2)*sigmaw2*tau_dissipturb/(invtau_phaserelax + invtau_e(i))
                     ! sursat ranges between -1 and 1, so we prevent sigma2 so exceed 1
                     cldfraliq(i) = 0.5*(1.-erf((sursat_iceliq - sursat_equ)/(SQRT(2.*sigma2_pdf))))
                     IF (cldfraliq(i) .GT. liqfra_max) THEN
                        cldfraliq(i) = liqfra_max
                     END IF

                     qliq_incl = qsati(i)*SQRT(sigma2_pdf)/SQRT(2.*RPI)*EXP(-1.*(sursat_iceliq - sursat_equ)**2./(2.*sigma2_pdf)) &
                                 - qsati(i)*cldfraliq(i)*(sursat_iceliq - sursat_equ)

                     sigma2_icefracturb(i) = sigma2_pdf
                     mean_icefracturb(i) = sursat_equ

                     !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION  ------------

                     IF ((qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps)) THEN
                        qliq_incl = 0.
                        cldfraliq(i) = 0.
                     END IF

                     !--Specific humidity is the max between qsati and the weighted mean between
                     !--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are
                     !--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1)
                     !--The whole cloud can therefore be supersaturated but never subsaturated.

                     qvap_incl = MAX(qsati(i), (1.-cldfraliq(i))*(sursat_iceext + 1.)*qsati(i) + cldfraliq(i)*qsatl(i))

                     IF (qvap_incl .GE. qtot_incl(i)) THEN
                        qvap_incl = qsati(i)
                        qliq_incl = MAX(0.0, qtot_incl(i) - qvap_incl)
                        qice_incl = 0.

                     ELSEIF ((qvap_incl + qliq_incl) .GE. qtot_incl(i)) THEN
                        qliq_incl = MAX(0.0, qtot_incl(i) - qvap_incl)
                        qice_incl = 0.
                     ELSE
                        qice_incl = qtot_incl(i) - qvap_incl - qliq_incl
                     END IF

                     qvap_cld(i) = qvap_incl*cldfra1D
                     qliq(i) = qliq_incl*cldfra1D
                     qice(i) = qice_incl*cldfra1D
                     IF ((qice(i) + qliq(i)) .GT. 0.) THEN
                        icefrac(i) = qice(i)/(qice(i) + qliq(i))
                     ELSE
                        icefrac(i) = 1. ! to keep computation of qsat wrt ice in condensation loop in lmdz_lscp_main
                     END IF
                     dicefracdT(i) = 0.

                  END IF ! Enough TKE

               END IF ! End qini

            END IF ! pticefracturb and cldfra

         END IF ! MPC temp

      END DO ! klon

   END SUBROUTINE ICEFRAC_LSCP_TURB
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

END MODULE lmdz_lscp_phase
