MODULE lmdz_lscp_tools

    IMPLICIT NONE

CONTAINS

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)

    ! Ref:
    ! Stubenrauch, C. J., Bonazzola, M.,
    ! Protopapadaki, S. E., & Musat, I. (2019).
    ! New cloud system metrics to assess bulk
    ! ice cloud schemes in a GCM. Journal of
    ! Advances in Modeling Earth Systems, 11,
    ! 3212–3234. https://doi.org/10.1029/2019MS001642
    
    use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc
    use lmdz_lscp_ini, only: cice_velo, dice_velo 

    IMPLICIT NONE

    INTEGER, INTENT(IN) :: klon
    REAL, INTENT(IN), DIMENSION(klon) :: iwc       ! specific ice water content [kg/m3]
    REAL, INTENT(IN), DIMENSION(klon) :: temp      ! temperature [K]
    REAL, INTENT(IN), DIMENSION(klon) :: rho       ! dry air density [kg/m3]
    REAL, INTENT(IN), DIMENSION(klon) :: pres      ! air pressure [Pa]
    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv    ! convective point  [-]

    REAL, INTENT(OUT), DIMENSION(klon) :: velo    ! fallspeed velocity of crystals [m/s]


    INTEGER i
    REAL logvm,iwcg,tempc,phpa,fallv_tun
    REAL m2ice, m2snow, vmice, vmsnow
    REAL aice, bice, asnow, bsnow
    

    DO i=1,klon

        IF (ptconv(i)) THEN
            fallv_tun=ffallv_con
        ELSE
            fallv_tun=ffallv_lsc
        ENDIF

        tempc=temp(i)-273.15 ! celcius temp
        iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0
        phpa=pres(i)/100.    ! pressure in hPa

    IF (iflag_vice .EQ. 1) THEN
        ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
        if (tempc .GE. -60.0) then

            logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) &
                    -0.0516344*log(iwcg)+0.00216078*tempc + 1.9714    
            velo(i)=exp(logvm)
        else
            velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.15
        endif
        
        velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s

    ELSE IF (iflag_vice .EQ. 2) THEN
        ! so called  PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019
        aice=0.587
        bice=2.45
        asnow=0.0444
        bsnow=2.1
        
        m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* & 
                exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc)))   &
                **(1./(0.807+bice*0.00581+0.0457*bice**2))

        vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ &
                 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) &
                 *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice)

       
        vmice=vmice*((1000./phpa)**0.2)
      
        m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* &
               exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc)))         &
               **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2))


        vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)&
                  *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)&
                  *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow)

        vmsnow=vmsnow*((1000./phpa)**0.35)
        velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s

    ELSE
        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
        velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo)
    ENDIF
    ENDDO

END SUBROUTINE FALLICE_VELOCITY
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  
  ! Compute the ice fraction 1-xliq (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
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    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
    REAL, INTENT(IN), DIMENSION(klon)   :: temp              ! temperature
    REAL, INTENT(IN), DIMENSION(klon)   :: distcltop         ! distance to cloud top
    REAL, INTENT(IN), DIMENSION(klon)   :: temp_cltop        ! temperature of cloud top
    INTEGER, INTENT(IN)                 :: iflag_ice_thermo
    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT


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

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


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

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

        ENDIF

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

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

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

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

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

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


   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

   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 :: 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
     ! because cldfra is intent in, but can be locally modified due to test 
     cldfra1D = cldfra(i)
     ! activate param for concerned grid points and for cloudy conditions
     IF ((pticefracturb(i)) .AND. (cldfra(i) .GT. 0.)) THEN
        IF (cldfra(i) .GE. 1.0) THEN
           cldfra1D = 1.0
        END IF
        
        ! 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.


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

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

           !--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 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 (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.)
                 ENDIF
              ELSE
                 nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033)
              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.
                 ENDIF
                 
              !--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. * 2./3. * tke(i) / 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 ) *  2./3. * tke(i) * 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
                 ENDIF 
                 
                 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
                 ENDIF
                 dicefracdT(i) = 0.

              END IF ! Enough TKE
           
           END IF ! End qini

        END IF ! ! MPC temperature

     END IF ! pticefracturb and cldfra
   
   ENDDO ! klon
END SUBROUTINE ICEFRAC_LSCP_TURB
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ! Calculate qsat following ECMWF method
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


USE yoethf_mod_h
        USE yomcst_mod_h
IMPLICIT NONE


    include "FCTTRE.h"

    INTEGER, INTENT(IN) :: klon  ! number of horizontal grid points
    REAL, INTENT(IN), DIMENSION(klon) :: temp     ! temperature in K
    REAL, INTENT(IN), DIMENSION(klon) :: qtot     ! total specific water in kg/kg
    REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa
    REAL, INTENT(IN)                  :: tref     ! reference temperature in K
    LOGICAL, INTENT(IN) :: flagth     ! flag for qsat calculation for thermals
    INTEGER, INTENT(IN) :: phase 
    ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid)
    !        1=liquid
    !        2=solid

    REAL, INTENT(OUT), DIMENSION(klon) :: qs      ! saturation specific humidity [kg/kg]
    REAL, INTENT(OUT), DIMENSION(klon) :: dqs     ! derivation of saturation specific humidity wrt T

    REAL delta, cor, cvm5
    INTEGER i

    DO i=1,klon

    IF (phase .EQ. 1) THEN
        delta=0.
    ELSEIF (phase .EQ. 2) THEN
        delta=1.
    ELSE
        delta=MAX(0.,SIGN(1.,tref-temp(i)))
    ENDIF

    IF (flagth) THEN
    cvm5=R5LES*(1.-delta) + R5IES*delta
    ELSE
    cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta
    cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))
    ENDIF

    qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)
    qs(i)=MIN(0.5,qs(i))
    cor=1./(1.-RETV*qs(i))
    qs(i)=qs(i)*cor
    dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)

    END DO

END SUBROUTINE CALC_QSAT_ECMWF
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ! programme that calculates the gammasat parameter that determines the
    ! homogeneous condensation thresholds for cold (<0oC) clouds
    ! condensation at q>gammasat*qsat
    ! Etienne Vignon, March 2021
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT
    use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez

    IMPLICIT NONE


    INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
    REAL, INTENT(IN), DIMENSION(klon) :: temp         ! temperature in K
    REAL, INTENT(IN), DIMENSION(klon) :: qtot         ! total specific water in kg/kg

    REAL, INTENT(IN), DIMENSION(klon) :: pressure     ! pressure in Pa

    REAL, INTENT(OUT), DIMENSION(klon) :: gammasat    ! coefficient to multiply qsat with to calculate saturation
    REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature

    REAL, DIMENSION(klon) ::  qsi,qsl,dqsl,dqsi
    REAL  f_homofreez, fac

    INTEGER i
    
        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl)
        CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi)

    DO i = 1, klon

        IF ( temp(i) .GE. RTT ) THEN
            ! warm clouds: condensation at saturation wrt liquid
            gammasat(i) = 1.
            dgammasatdt(i) = 0.

        ELSE
            ! cold clouds: qsi > qsl
            
            ! homogeneous freezing of aerosols, according to
            ! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS)
            ! 'Cirrus regime'
            ! if f_homofreez > qsl / qsi, liquid nucleation
            ! if f_homofreez < qsl / qsi, homogeneous freezing of aerosols
            ! Note: f_homofreez = qsl / qsi for temp ~= -38degC
            f_homofreez = a_homofreez - temp(i) / b_homofreez
            
            IF ( iflag_gammasat .GE. 3 ) THEN
              ! condensation at homogeneous freezing threshold for temp < -38 degC
              ! condensation at liquid saturation for temp > -38 degC
              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
                gammasat(i) = f_homofreez
                dgammasatdt(i) = - 1. / b_homofreez
              ELSE
                gammasat(i) = qsl(i) / qsi(i)
                dgammasatdt(i) = ( dqsl(i) * qsi(i) - dqsi(i) * qsl(i) ) / qsi(i) / qsi(i)
              ENDIF

            ELSEIF ( iflag_gammasat .EQ. 2 ) THEN
              ! condensation at homogeneous freezing threshold for temp < -38 degC
              ! condensation at a threshold linearly decreasing between homogeneous
              ! freezing and ice saturation for -38 degC < temp < temp_nowater
              ! condensation at ice saturation for temp > temp_nowater
              ! If temp_nowater = 235.15 K, this is equivalent to iflag_gammasat = 1
              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
                gammasat(i) = f_homofreez
                dgammasatdt(i) = - 1. / b_homofreez
              ELSEIF ( temp(i) .LE. temp_nowater ) THEN
                ! Here, we assume that f_homofreez = qsl / qsi for temp = -38 degC = 235.15 K
                dgammasatdt(i) = ( a_homofreez - 235.15 / b_homofreez - 1. ) &
                               / ( 235.15 - temp_nowater )
                gammasat(i) = dgammasatdt(i) * ( temp(i) - temp_nowater ) + 1.
              ELSE
                gammasat(i) = 1.
                dgammasatdt(i) = 0.
              ENDIF

            ELSEIF ( iflag_gammasat .EQ. 1 ) THEN
              ! condensation at homogeneous freezing threshold for temp < -38 degC
              ! condensation at ice saturation for temp > -38 degC
              IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN
                gammasat(i) = f_homofreez
                dgammasatdt(i) = - 1. / b_homofreez
              ELSE
                gammasat(i) = 1.
                dgammasatdt(i) = 0.
              ENDIF

            ELSE
              ! condensation at ice saturation for temp < -38 degC
              ! condensation at ice saturation for temp > -38 degC
              gammasat(i) = 1.
              dgammasatdt(i) = 0.

            ENDIF

            ! Note that the delta_hetfreez parameter allows to linearly decrease the
            ! condensation threshold between the calculated threshold and the ice saturation
            ! for delta_hetfreez = 1, the threshold is the calculated condensation threshold
            ! for delta_hetfreez = 0, the threshold is the ice saturation
            gammasat(i) = ( 1. - delta_hetfreez ) + delta_hetfreez * gammasat(i)
            dgammasatdt(i) = delta_hetfreez * dgammasatdt(i)

        ENDIF
   
    END DO


END SUBROUTINE CALC_GAMMASAT
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   
   USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl

   IMPLICIT NONE
   
   INTEGER, INTENT(IN) :: klon,klev                !number of horizontal and vertical grid points
   INTEGER, INTENT(IN) :: k                        ! vertical index
   REAL, INTENT(IN), DIMENSION(klon,klev) :: temp  ! temperature in K
   REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa
   REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa
   REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb  ! cloud fraction

   REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D  ! distance from cloud top
   REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop     ! temperature of cloud top
   
   REAL dzlay(klon,klev)
   REAL zlay(klon,klev)
   REAL dzinterf
   INTEGER i,k_top, kvert
   LOGICAL bool_cl


   DO i=1,klon
         ! Initialization height middle of first layer
          dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2))
          zlay(i,1) = dzlay(i,1)/2

          DO kvert=2,klev
                 IF (kvert.EQ.klev) THEN
                       dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert)))
                 ELSE
                       dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1))
                 ENDIF
                       dzinterf       = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert))
                       zlay(i,kvert)  = zlay(i,kvert-1) + dzinterf
           ENDDO
   ENDDO
   
   DO i=1,klon
          k_top = k
          IF (rneb(i,k) .LE. tresh_cl) THEN
                 bool_cl = .FALSE.
          ELSE
                 bool_cl = .TRUE.
          ENDIF

          DO WHILE ((bool_cl) .AND. (k_top .LE. klev))
          ! find cloud top
                IF (rneb(i,k_top) .GT. tresh_cl) THEN
                      k_top = k_top + 1
                ELSE
                      bool_cl = .FALSE.
                      k_top   = k_top - 1
                ENDIF
          ENDDO
          k_top=min(k_top,klev)

          !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to
          !interf for layer of cloud top
          distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2
          temp_cltop(i)  = temp(i,k_top)
   ENDDO ! klon

END SUBROUTINE DISTANCE_TO_CLOUD_TOP
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION GAMMAINC ( p, x )

!*****************************************************************************80
!
!! GAMMAINC computes the regularized lower incomplete Gamma Integral
!
!  Modified:
!
!    20 January 2008
!
!  Author:
!
!    Original FORTRAN77 version by B Shea.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    B Shea,
!    Algorithm AS 239:
!    Chi-squared and Incomplete Gamma Integral,
!    Applied Statistics,
!    Volume 37, Number 3, 1988, pages 466-473.
!
!  Parameters:
!
!    Input, real X, P, the parameters of the incomplete
!    gamma ratio.  0 <= X, and 0 < P.
!
!    Output, real GAMMAINC, the value of the incomplete
!    Gamma integral.
!
  IMPLICIT NONE

  REAL A
  REAL AN
  REAL ARG
  REAL B
  REAL C
  REAL, PARAMETER :: ELIMIT = - 88.0E+00
  REAL GAMMAINC
  REAL, PARAMETER :: OFLO = 1.0E+37
  REAL P
  REAL, PARAMETER :: PLIMIT = 1000.0E+00
  REAL PN1
  REAL PN2
  REAL PN3
  REAL PN4
  REAL PN5
  REAL PN6
  REAL RN
  REAL, PARAMETER :: TOL = 1.0E-14
  REAL X
  REAL, PARAMETER :: XBIG = 1.0E+08

  GAMMAINC = 0.0E+00

  IF ( X == 0.0E+00 ) THEN
    GAMMAINC = 0.0E+00
    RETURN
  END IF
!
!  IF P IS LARGE, USE A NORMAL APPROXIMATION.
!
  IF ( PLIMIT < P ) THEN

    PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) &
    + 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 )

    GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) )
    RETURN

  END IF
!
!  IF X IS LARGE SET GAMMAD = 1.
!
  IF ( XBIG < X ) THEN
    GAMMAINC = 1.0E+00
    RETURN
  END IF
!
!  USE PEARSON'S SERIES EXPANSION.
!  (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM).
!
  IF ( X <= 1.0E+00 .OR. X < P ) THEN

    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 )
    C = 1.0E+00
    GAMMAINC = 1.0E+00
    A = P

    DO

      A = A + 1.0E+00
      C = C * X / A
      GAMMAINC = GAMMAINC + C

      IF ( C <= TOL ) THEN
        EXIT
      END IF

    END DO

    ARG = ARG + LOG ( GAMMAINC )

    IF ( ELIMIT <= ARG ) THEN
      GAMMAINC = EXP ( ARG )
    ELSE
      GAMMAINC = 0.0E+00
    END IF
!
!  USE A CONTINUED FRACTION EXPANSION.
!
  ELSE

    ARG = P * LOG ( X ) - X - LOG_GAMMA ( P )
    A = 1.0E+00 - P
    B = A + X + 1.0E+00
    C = 0.0E+00
    PN1 = 1.0E+00
    PN2 = X
    PN3 = X + 1.0E+00
    PN4 = X * B
    GAMMAINC = PN3 / PN4

    DO

      A = A + 1.0E+00
      B = B + 2.0E+00
      C = C + 1.0E+00
      AN = A * C
      PN5 = B * PN3 - AN * PN1
      PN6 = B * PN4 - AN * PN2

      IF ( PN6 /= 0.0E+00 ) THEN

        RN = PN5 / PN6

        IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN
          EXIT
        END IF

        GAMMAINC = RN

      END IF

      PN1 = PN3
      PN2 = PN4
      PN3 = PN5
      PN4 = PN6
!
!  RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE.
!
      IF ( OFLO <= ABS ( PN5 ) ) THEN
        PN1 = PN1 / OFLO
        PN2 = PN2 / OFLO
        PN3 = PN3 / OFLO
        PN4 = PN4 / OFLO
      END IF

    END DO

    ARG = ARG + LOG ( GAMMAINC )

    IF ( ELIMIT <= ARG ) THEN
      GAMMAINC = 1.0E+00 - EXP ( ARG )
    ELSE
      GAMMAINC = 1.0E+00
    END IF

  END IF

  RETURN
END FUNCTION GAMMAINC
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

END MODULE lmdz_lscp_tools


