MODULE lmdz_lscp_tools

   IMPLICIT NONE

CONTAINS

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE FALLICE_VELOCITY(klon, iwc, temp, rho, pres, ptconv, velo)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! routine that calculates the ice crystals fall velocity
! References:
! - Heymsfield 1977, doi: 10.1175/1520-0469(1977)034<0367:PDISIC>2.0.CO;2
! - Heymsfield and Donner 1990, doi: 10.1175/1520-0469(1990)047<1865:ASFPIC>2.0.CO;2
! - Stubenrauch et al. 2019, doi: 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
         END IF

         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
            end if

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

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

   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 [K]
      REAL, INTENT(IN), DIMENSION(klon) :: qtot     !-- total specific water [kg/kg]
      REAL, INTENT(IN), DIMENSION(klon) :: pressure !-- pressure [Pa]
      REAL, INTENT(IN)                  :: tref     ! reference temperature [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     ! dqs/dT [kg/kg/K]

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

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

         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)

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      ! routine that calculates the gammasat parameter that determines the
      ! homogeneous condensation thresholds for cold (<0oC) clouds
      ! condensation at q>gammasat*qsat
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      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 [K]
      REAL, INTENT(IN), DIMENSION(klon) :: qtot         !-- total specific water [kg/kg]

      REAL, INTENT(IN), DIMENSION(klon) :: pressure     !-- pressure [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 [K-1]

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

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

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

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

            END IF

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

         END IF

      END DO

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

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon, klev, k, temp, pplay, paprs, rneb, distcltop1D, temp_cltop)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      ! routine that calculates the distance to cloud top for cloud phase
      ! determination (when cloud phase is assumed to depend on both
      ! temperature and distance to cloud top, see Lea Raillard's PhD
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      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 [K]
      REAL, INTENT(IN), DIMENSION(klon, klev) :: pplay !-- pressure middle layer [Pa]
      REAL, INTENT(IN), DIMENSION(klon, klev + 1) :: paprs !-- pressure interfaces [Pa]
      REAL, INTENT(IN), DIMENSION(klon, klev) :: rneb    !-- cloud fraction [-]

      REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D !-- distance from cloud top [m]
      REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop  !-- temperature of cloud top [K]

      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))
            END IF
            dzinterf = rd*temp(i, kvert)/rg*log(pplay(i, kvert - 1)/pplay(i, kvert))
            zlay(i, kvert) = zlay(i, kvert - 1) + dzinterf
         END DO
      END DO

      DO i = 1, klon
         k_top = k
         IF (rneb(i, k) .LE. tresh_cl) THEN
            bool_cl = .FALSE.
         ELSE
            bool_cl = .TRUE.
         END IF

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