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
    
    use lscp_ini_mod, only: iflag_vice, ffallv_con, ffallv_lsc

    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,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=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
        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)/((iwcg/1000.*rho(i))**dvel)

    ELSE
        ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
        velo(i) = fallv_tun*3.29/2.0*((iwcg/1000.)**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, iflag_ice_thermo, 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
    USE lscp_ini_mod, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace
    USE lscp_ini_mod, ONLY : RTT

    IMPLICIT NONE


    INTEGER, INTENT(IN)                 :: klon       ! number of horizontal grid points
    REAL, INTENT(IN), DIMENSION(klon)   :: temp       ! temperature
    REAL, INTENT(IN), DIMENSION(klon)   :: sig
    INTEGER, INTENT(IN)                 :: iflag_ice_thermo
    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(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    ! Calculate qsat following ECMWF method
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    IMPLICIT NONE

    include "YOMCST.h"
    include "YOETHF.h"
    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 lscp_ini_mod, only: iflag_gammasat, t_glace_min, RTT

    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  fcirrus, fac
    REAL, PARAMETER :: acirrus=2.349
    REAL, PARAMETER :: bcirrus=259.0

    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.

        ELSEIF ((temp(i) .LT. RTT) .AND. (temp(i) .GT. t_glace_min)) THEN
            
            IF (iflag_gammasat .GE. 2) THEN          
               gammasat(i)=qsl(i)/qsi(i)
               dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
            ELSE
               gammasat(i)=1.
               dgammasatdt(i)=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(i)/bcirrus
               IF (fcirrus .LT. qsl(i)/qsi(i)) THEN
                  gammasat(i)=qsl(i)/qsi(i)
                  dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)
               ELSE
                  gammasat(i)=fcirrus
                  dgammasatdt(i)=-1.0/bcirrus
               ENDIF
           
            ELSE

               gammasat(i)=1.
               dgammasatdt(i)=0.

            ENDIF

        ENDIF
   
    END DO


END SUBROUTINE CALC_GAMMASAT


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

END MODULE LSCP_TOOLS_MOD