module generic_cloud_common_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 = 140.6931
            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 = 78.04454 !Na2S
            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.2E3
            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 = 87.003049
            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 = 97.445
            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 specie_parameters_table(specname)
        
        implicit none
    !============================================================================
    !   Load the adequate set of parameters for specname
    !   From a table of traceurs
    !
    !   Also set up few others useful parameters, such as epsi=m/mugaz, RLVTT and
    !   the metallicity_coeff.
    !   Authors 
    !   --------
    !   Noé Clément (2022) 
    !============================================================================

        character(*), intent(in) :: specname
        integer :: ios
        character(len=500):: table_traceurs_line ! table_traceurs_line lines with parameters

        open(117,file='table_tracers_condensable',form='formatted',status='old')

        read(117,'(A)') table_traceurs_line

        do
            read(117,'(A)', iostat=ios) table_traceurs_line
            
            if (index(table_traceurs_line,specname) /= 0) then
                
                if (index(table_traceurs_line,'delta_vapH='   ) /= 0) then
                    read(table_traceurs_line(index(table_traceurs_line,'delta_vapH=')+len('delta_vapH='):),*) delta_vapH
                    print*, 'delta_vapH ', delta_vapH
                end if
                if (index(table_traceurs_line,'Tref='   ) /= 0) then
                    read(table_traceurs_line(index(table_traceurs_line,'Tref=')+len('Tref='):),*) Tref
                    print*, 'Tref', Tref
                end if 
                if (index(table_traceurs_line,'Pref='   ) /= 0) then
                    read(table_traceurs_line(index(table_traceurs_line,'Pref=')+len('Pref='):),*) Pref
                    print*, 'Pref', Pref
                end if
                if (index(table_traceurs_line,'mass='   ) /= 0) then
                    read(table_traceurs_line(index(table_traceurs_line,'mass=')+len('mass='):),*) m
                    print*, 'mass', m
                end if
                if (index(table_traceurs_line,'metallicity_coeff='   ) /= 0) then
                   read(table_traceurs_line(index(table_traceurs_line,'metallicity_coeff=')+len('metallicity_coeff='):),*) metallicity_coeff
                   print*, 'metallicity_coeff', metallicity_coeff
                end if

                ios=1
            end if

            if (ios /= 0) exit
            
        end do
        ! RLVTT

        !if (is_master)
        close(117)

        RLVTT=delta_vapH/(m/1000.)

        write(*,*) 'RLVTT', RLVTT

        epsi = m/mugaz

    end subroutine specie_parameters_table

    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 generic_cloud_common_h