MODULE MP2M_CLOUDS_METHODS
    !============================================================================
    !
    !     Purpose
    !     -------
    !     Cloud model miscellaneous methods module.
    !
    !     The module contains miscellaneous methods used in the clouds of the model.
    !     The module contains four interfaces (8 methods):
    !        - mm_sigX   | Compute surface tension
    !        - mm_fshape | Compute shape factor
    !        - mm_qsatX  | Compute saturation molar mixing ratio
    !        - mm_LheatX | Compute latent heat released
    !
    !     Authors
    !     -------
    !     B. de Batz de Trenquelléon (10/2025)
    !
    !============================================================================

    USE MP2M_MPREC
    USE MP2M_GLOBALS
    USE LINT_DSET
    USE LINT_LOCATORS
    IMPLICIT NONE

    PRIVATE

    PUBLIC :: mm_sigX, mm_fshape, mm_qsatX, mm_LheatX

    !============================================================================
    !                             INTERFACES
    !============================================================================

    !! Interface to surface tension computation functions.
    !! The method computes the surface tension of a given species at given temperature(s).
    INTERFACE mm_sigX
      MODULE PROCEDURE sigX_sc,sigX_ve
    END INTERFACE mm_sigX

    !! Interface to shape factor computation functions.
    !! The method computes the shape factor for the heterogeneous nucleation.
    INTERFACE mm_fshape
      MODULE PROCEDURE fshape_sc,fshape_ve
    END INTERFACE mm_fshape

    !! Interface to saturation molar mixing ratio computation functions.
    !! The method computes the molar mixing ratio at saturation of a given species at given temperature(s)
    !! and pressure level(s).
    INTERFACE mm_qsatX
      MODULE PROCEDURE qsatX_sc,qsatX_ve
    END INTERFACE mm_qsatX
    
    !! Interface to latent heat computation functions.
    !! The method computes the latent heat released of a given species at given temperature(s).
    INTERFACE mm_LheatX
      MODULE PROCEDURE LheatX_sc,LheatX_ve
    END INTERFACE mm_LheatX

    CONTAINS

    !============================================================================
    !                   CLOUD CONDENSATION NUCLEI METHODS
    !============================================================================

    ! FUNCTION mm_sigX(temp,xESP):
    ! xESP must always be given as a scalar. If temp is given as a vector, then the method
    ! computes the result for all the temperatures and returns a vector of same size than temp.
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    FUNCTION sigX_sc(temp,xESP) RESULT(res)
        !! Get the surface tension between a given species and the air (scalar).
        !! The method computes the surface tension equation as given in Reid et al. (1986) p. 637 (eq. 12-3.6).
        !!
        REAL(kind=mm_wp), INTENT(in) :: temp ! Temperature (K).
        TYPE(mm_esp), INTENT(in)     :: xESP ! Specie properties.
        REAL(kind=mm_wp)             :: res  ! Surface tension (N.m-1).

        ! Local variables:
        REAL(kind=mm_wp) :: Tr, Tbr, sig0, sig

        Tr = MIN(temp/xESP%Tc,0.99_mm_wp)
        Tbr = xESP%Tb/xESP%Tc

        sig0 = 0.1196_mm_wp*(1._mm_wp+(Tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-Tbr))-0.279_mm_wp

        sig = xESP%pc**(2._mm_wp/3._mm_wp) * xESP%Tc**(1._mm_wp/3._mm_wp) * sig0 * (1._mm_wp-Tr)**(11._mm_wp/9._mm_wp)

        ! Convertion (dyn/cm) --> (N/m):
        res = sig*1e-3_mm_wp
        RETURN
    END FUNCTION sigX_sc

    FUNCTION sigX_ve(temp,xESP) RESULT(res)
        !! Get the surface tension between a given species and the air (vector).
        !! The method computes the surface tension equation as given in Reid et al. (1986) p. 637 (eq. 12-3.6).
        !!
        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperatures (K).
        TYPE(mm_esp), INTENT(in)                   :: xESP ! Specie properties.
        REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  ! Surface tensions (N.m-1).

        ! Local variables:
        INTEGER          :: i
        REAL(kind=mm_wp) :: Tr, Tbr, sig0, sig

        Tbr = xESP%Tb/xESP%Tc

        sig0 = 0.1196_mm_wp*(1._mm_wp+(Tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-Tbr))-0.279_mm_wp

        DO i = 1, SIZE(temp)
          Tr  = MIN(temp(i)/xESP%Tc,0.99_mm_wp)
          sig = xESP%pc**(2._mm_wp/3._mm_wp) * xESP%Tc**(1._mm_wp/3._mm_wp) * sig0 * (1._mm_wp-Tr)**(11._mm_wp/9._mm_wp)
          
          ! Convertion (dyn/cm) --> (N/m):
          res(i) = sig*1e-3_mm_wp
        ENDDO
        RETURN
    END FUNCTION sigX_ve


    ! FUNCTION mm_fshape(m,x):
    ! Where m is cosine of the contact angle and x the curvature radius. m must always be
    ! given as a scalar. If x is given as a vector, then the method compute the result for each
    ! value of x and and returns a vector of same size than x.
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    FUNCTION fshape_sc(m,x) RESULT(res)
        !! Get the shape factor of a ccn (scalar).
        !! The method computes the shape factor for the heterogeneous nucleation on a fractal particle.
        !! Details about the shape factor can be found in Fletcher et al. (1958).
        !!
        REAL(kind=mm_wp), INTENT(in) :: m   ! Cosine of the contact angle.
        REAL(kind=mm_wp), INTENT(in) :: x   ! Curvature radius (r_particle/r*).
        REAL(kind=mm_wp)             :: res ! Shape factor value.

        ! Local variables:
        REAL(kind=mm_wp) :: phi, a, b, c

        IF (x > 3000._mm_wp) THEN
          res = ((2._mm_wp+m) * (1._mm_wp-m)**2) / 4._mm_wp
        ELSE
          phi = dsqrt(1._mm_wp-2._mm_wp*m*x+x**2)
          a = 1._mm_wp + ((1._mm_wp-m*x)/phi)**3
          b = (x**3) * (2._mm_wp - 3._mm_wp*(x-m)/phi + ((x-m)/phi)**3)
          c = 3._mm_wp*m*(x**2) * ((x-m)/phi-1._mm_wp)
          res = 0.5_mm_wp*(a + b + c)
        ENDIF
        RETURN
    END FUNCTION fshape_sc

    FUNCTION fshape_ve(m,x) RESULT(res)
        !! Get the shape factor of a ccn (vector).
        !! The method computes the shape factor for the heterogeneous nucleation on a fractal particle.
        !! Details about the shape factor can be found in Fletcher et al. (1958).
        !!
        REAL(kind=mm_wp), INTENT(in)               :: m   ! Cosine of the contact angle.
        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: x   ! Curvature radii (r_particle/r*).
        REAL(kind=mm_wp), DIMENSION(SIZE(x))       :: res ! Shape factor value.

        ! Local variables:
        REAL(kind=mm_wp), DIMENSION(SIZE(x)) :: phi, a, b, c

        WHERE(x > 3000._mm_wp)
          res = ((2._mm_wp+m) * (1._mm_wp-m)**2) / 4._mm_wp
        ELSEWHERE
          phi = dsqrt(1._mm_wp-2._mm_wp*m*x+x**2)
          a = 1._mm_wp + ((1._mm_wp-m*x)/phi)**3
          b = (x**3) * (2._mm_wp - 3._mm_wp*(x-m)/phi + ((x-m)/phi)**3)
          c = 3._mm_wp*m*(x**2) * ((x-m)/phi-1._mm_wp)
          res = 0.5_mm_wp*(a + b + c)
        ENDWHERE
        RETURN
    END FUNCTION fshape_ve

    !============================================================================
    !                       CONDENSATION METHODS
    !============================================================================

    ! FUNCTION mm_qsatX(temp,pres,xESP)
    ! xESP must always be given as a scalar. If temp and pres  are given as a vector,
    ! then the method computes the result for each couple of (temperature, pressure) and returns
    ! a vector of same size than temp.
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    FUNCTION qsatX_sc(temp,pres,xESP) RESULT(res)
        !! Get the molar mixing ratio of a given species at saturation (scalar).
        !! Compute saturation molar mixing ratio for condensable tracers.
        !!
        !! @warning:
        !! The formula depends on the species (and the reference)!
        !!
        REAL(kind=mm_wp), INTENT(in) :: temp ! Temperature (K).
        REAL(kind=mm_wp), INTENT(in) :: pres ! Pressure level (Pa).
        TYPE(mm_esp), INTENT(in)     :: xESP ! Specie properties.
        REAL(kind=mm_wp)             :: res  ! Saturation molar mixing ratio (mol/mol).

        ! Local variables:
        REAL(kind=mm_wp) :: fp, fsat

        fp = (1.0e5 / pres)

        ! C2H2, C6H6, HCN: Fray & Schmitt (2009)
        !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if ((xESP%name == 'C2H2').OR.(xESP%name == 'C6H6').OR.(xESP%name == 'HCN')) then
            fsat = xESP%a0_sat + xESP%a1_sat/temp + xESP%a2_sat/temp**2 + xESP%a3_sat/temp**3 + &
                   xESP%a4_sat/temp**4 + xESP%a5_sat/temp**5 + xESP%a6_sat/temp**6

            fsat = exp(fsat)

        ! C2H6: Dykyj et al. (1999)
        ! C4H2: Orton et al. (2014)
        !~~~~~~~~~~~~~~~~~~~~~~~~~~
        else if ((xESP%name == 'C2H6').OR.(xESP%name == 'C4H2')) then
            fsat = xESP%a0_sat + xESP%a1_sat/(temp+xESP%a2_sat)

            fsat = 10.**(fsat)
        
        ! Otherwise: error
        !~~~~~~~~~~~~~~~~~
        else
            write(*,'(a)') "[qsatX_sc] This is a fatal error..."
            write(*,'(a)') "Species no found for condensation!"
            call exit(111)
        endif

        res = fp * fsat
        RETURN
    END FUNCTION qsatX_sc

    FUNCTION qsatX_ve(temp,pres,xESP) RESULT(res)
        !! Get the molar mixing ratio of a given species at saturation (vector).
        !! Compute saturation molar mixing ratio for condensable tracers.
        !!
        !! @warning:
        !! The formula depends on the species (and the reference)!
        !!
        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperatures (K).
        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pres ! Pressure levels (Pa).
        TYPE(mm_esp), INTENT(in)                   :: xESP ! Specie properties.
        REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  ! Saturation molar mixing ratio (mol/mol).

        REAL(kind=mm_wp), DIMENSION(SIZE(temp)) :: fp, fsat

        fp = (1.0e5 / pres)

        ! C2H2, C6H6, HCN: Fray & Schmitt (2009)
        !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        if ((xESP%name == 'C2H2').OR.(xESP%name == 'C6H6').OR.(xESP%name == 'HCN')) then
            fsat = xESP%a0_sat + xESP%a1_sat/temp + xESP%a2_sat/temp**2 + xESP%a3_sat/temp**3 + &
                   xESP%a4_sat/temp**4 + xESP%a5_sat/temp**5 + xESP%a6_sat/temp**6

            fsat = exp(fsat)
        
        ! C2H6: Dykyj et al. (1999)
        ! C4H2: Orton et al. (2014)
        !~~~~~~~~~~~~~~~~~~~~~~~~~~
        else if ((xESP%name == 'C2H6').OR.(xESP%name == 'C4H2')) then
            fsat = xESP%a0_sat + xESP%a1_sat/(temp+xESP%a2_sat)

            fsat = 10.**(fsat)
        
        ! Otherwise: error
        !~~~~~~~~~~~~~~~~~
        else
            write(*,'(a)') "[qsatX_ve] This is a fatal error..."
            write(*,'(a)') "Species no found for condensation!"
            call exit(111)
        endif

        res = fp * fsat
        RETURN
    END FUNCTION qsatX_ve


    ! FUNCTION mm_LheatX(temp,xESP):
    ! xESP must always be given as a scalar. If temp is given as a vector, then the method
    ! computes the result for all the temperatures and returns a vector of same size than temp.
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    FUNCTION LheatX_sc(temp,xESP) RESULT(res)
        !! Compute latent heat of a given species at given temperature (scalar).
        !! The method computes the latent heat equation as given in Reid et al. (1986) p. 220 (eq. 7-9.4).
        !!
        REAL(kind=mm_wp), INTENT(in) :: temp ! Temperature (K).
        TYPE(mm_esp), INTENT(in)     :: xESP ! Specie properties.
        REAL(kind=mm_wp)             :: res  ! Latent heat of given species at given temperature (J.kg-1).
    
        ! Local variables:
        REAL(kind=mm_wp) :: Tr, ftm

        Tr = temp / xESP%Tc
        ftm = MAX(1._mm_wp - Tr,1.e-3_mm_wp)

        res = (7.08_mm_wp*ftm**0.354_mm_wp + 10.95_mm_wp*xESP%w*ftm**0.456_mm_wp) * mm_rgas * xESP%Tc / xESP%masmol
        RETURN
    END FUNCTION LheatX_sc

    FUNCTION LheatX_ve(temp,xESP) RESULT(res)
        !! Compute latent heat of a given species at given temperature (vector).
        !! The method computes the latent heat equation as given in Reid et al. (1986) p. 220 (eq. 7-9.4).
        !!
        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp ! Temperatures (K).
        TYPE(mm_esp), INTENT(in)                   :: xESP ! Specie properties.
        REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  ! Latent heat of given species at given temperatures (J.kg-1).

        ! Local variables:
        INTEGER :: i
        REAL(kind=mm_wp) :: Tr, ftm

        DO i=1,SIZE(temp)
          Tr = temp(i) / xESP%Tc
          ftm = MAX(1._mm_wp - Tr,1.e-3_mm_wp)
          
          res(i) = (7.08_mm_wp*ftm**0.354_mm_wp + 10.95_mm_wp*xESP%w*ftm**0.456_mm_wp) * &
                    mm_rgas * xESP%Tc / xESP%masmol
        ENDDO
        RETURN
    END FUNCTION LheatX_ve

END MODULE MP2M_CLOUDS_METHODS