MODULE LSCP_TOOLS_MOD

    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
    
    
    IMPLICIT NONE

    INCLUDE "nuage.h"
    INCLUDE "fisrtilp.h"

    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,cvel,dvel,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=iwc(i)*1000.    ! iwc in g/m3
        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
        dvel=0.2
        cvel=fallv_tun*65.0*(rho(i)**0.2)*(150./phpa)**0.15

    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
        dvel=0.2
        cvel=velo(i)/((iwc(i)*rho(i))**dvel)

    ELSE
        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
        velo(i) = fallv_tun*3.29/2.0 * ((iwc(i))**0.16)
        dvel=0.16
        cvel=fallv_tun*3.29/2.0*(rho(i)**0.16)
    ENDIF

    ENDDO

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

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE ICEFRAC_LSCP(klon, temp, sig, 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

    IMPLICIT none


    INCLUDE "YOMCST.h"
    INCLUDE "nuage.h"
    INCLUDE "clesphys.h"


  ! nuage.h contains:
  ! t_glace_min: if T < Tmin, the cloud is only made of water ice
  ! t_glace_max: if T > Tmax, the cloud is only made of liquid water
  ! exposant_glace: controls the sharpness of the transition

    INTEGER, INTENT(IN)                 :: klon       ! number of horizontal grid points
    REAL, INTENT(IN), DIMENSION(klon)   :: temp       ! temperature
    REAL, INTENT(IN), DIMENSION(klon)   :: sig
    REAL, INTENT(OUT), DIMENSION(klon)  :: icefrac
    REAL, INTENT(OUT), DIMENSION(klon)  :: dicefracdT


    INTEGER i
    REAL    sig0,www,tmin_tmp,liqfrac_tmp
    REAL    Dv, denomdep,beta,qsi,dqsidt
    INTEGER exposant_glace_old
    REAL t_glace_min_old
    LOGICAL ice_thermo

    sig0=0.8
    t_glace_min_old = RTT - 15.0
    ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
    IF (ice_thermo) THEN
      exposant_glace_old = 2
    ELSE
      exposant_glace_old = 6
    ENDIF


! calculation of icefrac and dicefrac/dT

    DO i=1,klon

    IF (iflag_t_glace.EQ.1) THEN
            ! Transition to ice close to surface for T<Tmax
            ! w=1 at the surface and 0 for sig < sig0
            www=(max(sig(i)-sig0,0.))/(1.-sig0)
    ELSEIF (iflag_t_glace.GE.2) THEN
            ! No convertion to ice close to surface
            www = 0.
    ENDIF
     
    tmin_tmp=www*t_glace_max+(1.-www)*t_glace_min
    liqfrac_tmp=  (temp(i)-tmin_tmp) / (t_glace_max-tmin_tmp)
    liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)
       
    IF (iflag_t_glace.GE.3) THEN
            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

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


    RETURN

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



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


    IMPLICIT none

    include "YOMCST.h"
    include "YOETHF.h"
    include "FCTTRE.h"

    REAL, INTENT(IN) :: temp     ! temperature in K
    REAL, INTENT(IN) :: qtot     ! total specific water in kg/kg
    REAL, INTENT(IN) :: 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) :: qs      ! saturation specific humidity [kg/kg]
    REAL, INTENT(OUT) :: dqs     ! derivation of saturation specific humidity wrt T


    REAL delta, cor, cvm5

   
    IF (phase .EQ. 1) THEN
        delta=0.
    ELSEIF (phase .EQ. 2) THEN
        delta=1.
    ELSE
        delta=MAX(0.,SIGN(1.,tref-temp))
    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))
    ENDIF

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



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


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE CALC_GAMMASAT(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
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    IMPLICIT none

    include "YOMCST.h"
    include "YOETHF.h"
    include "FCTTRE.h"
    include "nuage.h" 


    REAL, INTENT(IN) :: temp     ! temperature in K
    REAL, INTENT(IN) :: qtot     ! total specific water in kg/kg

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

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

    REAL qsi,qsl,fac,dqsl,dqsi,fcirrus
    REAL, PARAMETER :: acirrus=2.349
    REAL, PARAMETER :: bcirrus=259.0


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

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

        ELSEIF ((temp .LT. RTT) .AND. (temp .GT. t_glace_min)) THEN
            
            IF (iflag_gammasat .GE. 2) THEN          
               gammasat=qsl/qsi
               dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
            ELSE
               gammasat=1.
               dgammasatdt=0.
            ENDIF

        ELSE

            IF (iflag_gammasat .GE.1) THEN
            ! homogeneous freezing of aerosols, according to
            ! Koop, 2000 and Karcher 2008, QJRMS
            ! 'Cirrus regime'
               fcirrus=acirrus-temp/bcirrus
               IF (fcirrus .LT. qsl/qsi) THEN
                  gammasat=qsl/qsi
                  dgammasatdt=(dqsl*qsi-dqsi*qsl)/qsi/qsi
               ELSE
                  gammasat=fcirrus
                  dgammasatdt=-1.0/bcirrus
               ENDIF
           
            ELSE

               gammasat=1.
               dgammasatdt=0.

            ENDIF

        ENDIF
   



END SUBROUTINE CALC_GAMMASAT


!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

END MODULE LSCP_TOOLS_MOD