module genericcommon_h
    use comcstfi_mod, only: r, cpp, mugaz
    implicit none 
    real, save :: m  ! molecular mass of the specie (g/mol)
    real, save :: delta_vapH ! Enthalpy of vaporization (J/mol)
    real,save :: Tref ! Ref temperature for Clausis-Clapeyron (K)
    real,save :: Pref ! Reference pressure for Clausius Clapeyron (Pa)
    real, save :: epsi ! fractionnal molecular mass (m/mugaz)
    real,save :: RLVTT !Latent heat of vaporization (J/kg)
    real,save :: metallicity_coeff ! Coefficient to take into account the metallicity
    contains 

    subroutine specie_parameters(specname)
        
        implicit none

    !============================================================================
    !   Load the adequate set of parameters for specname
    !   We use Clausius Clapeyron.
    !   Values are based on Vissher (2010) and Morley (2012).
    !
    !   Also set up few others useful parameters, such as epsi=m/mugaz, RLVTT and
    !   the metallicity_coeff.
    !   Authors 
    !   --------
    !   Lucas Teinturier (2022) 
    !============================================================================
        ! Inputs
        character(*), intent(in) :: specname
        print*,"entree dans specie_parameters avec specname = ",specname
        if (trim(specname) .eq. "Mg") then
            print*,"Loading data for Mg"
            delta_vapH = 521.7E3 !J/mol
            Tref = 2303.0 !K
            Pref =   1.0E5 !in Pa
            m = 24.3055
            RLVTT = delta_vapH/(m/1000.)
            metallicity_coeff=1.0*log(10.0)
            
        else if (trim(specname) .eq. "Na") then
            print*,"Loading data for Na"
            delta_vapH = 265.9E3
            Tref = 1624.0
            Pref =  1.0E5
            m = 22.9898
            RLVTT = delta_vapH/(m/1000.)
            metallicity_coeff=0.5*log(10.0)
            
        else if (trim(specname) .eq. "Fe") then 
            print*,"Loading data for Fe"
            delta_vapH = 401.9E3 
            Tref = 2903.0
            Pref  = 1.0E5
            m = 55.8450
            RLVTT = delta_vapH/(m/1000.)
            metallicity_coeff=0.0*log(10.0)
            
        else if (trim(specname) .eq. "Cr") then
            print*,"Loading data for Cr"
            delta_vapH = 394.2 
            Tref = 2749.0
            Pref  = 1.0E5
            m = 51.9961
            RLVTT = delta_vapH/(m/1000.)
            metallicity_coeff=0.0*log(10.0)
            
        else if (trim(specname) .eq. "KCl") then 
            print*,"Loading data for KCl"
            delta_vapH = 217.9E3
            Tref = 1495.0
            Pref = 1.0E5
            metallicity_coeff=0.0*log(10.0)
            m = 74.5498
            RLVTT = delta_vapH/(m/1000.)
        else if (trim(specname) .eq. "Mn") then 
            print*,"Loading data for Mn"
            delta_vapH = 455.8E3
            Tref = 2064.0
            Pref = 1.0E5
            metallicity_coeff=1.0*log(10.0)
            m = 54.9380
            RLVTT = delta_vapH/(m/1000.)
        else if (trim(specname) .eq. "Zn") then 
            print*,"Loading data for Zn"
            delta_vapH = 303.9E3
            Tref = 1238.0
            Pref = 1.0E5
            metallicity_coeff=1.0*log(10.0)
            m = 65.38
            RLVTT = delta_vapH/(m/1000.)
        else
            print*,"Unknow species (not in Mg, Fe, Na, KCl, Cr, Mn or Zn)"
        endif
        epsi = m/mugaz
    end subroutine specie_parameters

    subroutine Psat_generic(T,p,metallicity,psat,qsat)
        IMPLICIT NONE 
    !============================================================================
    !   Clausius-Clapeyron relation : 
    !   d(Ln(psat))/dT = delta_vapH/(RT^2)
    !   ->Psat = Pref * exp(-delta_vapH/R*(1/T-1/Tref))
    !
    !   Authors 
    !   --------
    !   Lucas Teinturier (2022) adapted from Psat_water of Jeremy Leconte (2012)
    !============================================================================
        !Inputs
        real, intent(in) :: T, p ! Temperature and pressure of the layer (in K and Pas)
        real, intent(in) :: metallicity ! metallycity (log10)
        !Outputs 
        real, intent(out) :: psat,qsat !saturation vapor pressure (Pa) and mass mixing ratio at saturation (kg/kg) of the layer 

        psat = pref*exp(-delta_vapH/(r*mugaz/1000.)*(1/T-1/Tref)-metallicity_coeff*metallicity) ! in Pa (because pref in Pa)

        if (psat .gt. p) then 
            qsat = 1.
        else 
            qsat = epsi *psat/(p-(1-epsi)*psat)
        endif
    end subroutine Psat_generic 

    subroutine Lcpdqsat_generic(T,p,psat,qsat,dqsat,dlnpsat)
        implicit none 

    !===============================================================================
    !   Compute dqsat = L/cp* d(qsat)/dT and d(Ln(psat)) = L/cp*d(Ln(psat))/dT
    !   we have d(ln(psat))/dT =  delta_vapH/R*(1/T^2) for clausius-clapeyron
    !   r*mugaz/1000. is the perfect gaz constant in the computation of "dummy"
    !   Authors 
    !   --------
    !   Lucas Teinturier (2022) adapted from Lcpdqsat_water of Jeremy Leconte (2012)
    !===============================================================================
        ! Inputs 
        real, intent(in) :: T ! Temperature (K)
        real, intent(in) :: p ! Pressure (Pa)
        real, intent(in) :: psat ! Saturation vapor pressure (Pa)
        real, intent(in) :: qsat ! Mass mixing ratio at saturation (kg/kg)

        ! Outputs 
        real, intent(out) :: dqsat,dlnpsat 

        ! Variables 
        real :: dummy ! used to store d(ln(psat))/dT

        dummy = delta_vapH/((r*mugaz/1000.)*(T**2)) 

        if (psat .gt. p) then 
            dqsat =0.
        else 
            dqsat = (RLVTT/cpp) *qsat*(p/(p-(1-epsi)*psat))*dummy
            dlnpsat = (RLVTT/cpp) * dummy 
        endif

    end subroutine Lcpdqsat_generic

end module genericcommon_h