source: trunk/LMDZ.GENERIC/libf/phystd/generic_cloud_common_h.F90 @ 3523

module generic_cloud_common_h

    !============================================================================
    ! Purpose
    ! -------
    ! Set up relevant constants and parameters for generic condensable tracers cycles, 
    ! and generic condensable tracers cloud properties
    !
    ! Authors
    ! -------
    ! Adapted from watercommon_h.F90
    ! by Lucas Teinturier & Noé Clément (2022)
    !============================================================================

    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_generic ! fractionnal molecular mass (m/mugaz)
    real,save :: RLVTT_generic !Latent heat of vaporization (J/kg)
    real,save :: metallicity_coeff ! Coefficient to take into account the metallicity
    real,save :: RCPV_generic ! specific heat capacity of the tracer vapor at Tref
    
    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_generic=m/mugaz, RLVTT_generic 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_generic = 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_generic = 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_generic = 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_generic = 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_generic = 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_generic = 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_generic = delta_vapH/(m/1000.)
        else
            print*,"Unknow species (not in Mg, Fe, Na, KCl, Cr, Mn or Zn)"
        endif
        epsi_generic = m/mugaz
    end subroutine specie_parameters


    subroutine specie_parameters_table(specname)
        use datafile_mod, only: datadir
        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_generic=m/mugaz, RLVTT_generic 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=trim(datadir)//'/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
                if (index(table_traceurs_line,'RCPV_generic='   ) /= 0) then
                    read(table_traceurs_line(index(table_traceurs_line,'RCPV_generic=')+len('RCPV_generic='):),*) RCPV_generic
                    print*, 'RCPV_generic', RCPV_generic
                end if
                         
                ios=1
            end if

            if (ios /= 0) exit
            
        end do

        if (ios .eq. -1) then ! If the file has been read but there is no data on the speci, then we have ios equal to -1

            write(*,*), "ERROR : the following traceur is not referenced in the table_tracers_condensable : ", specname

        end if
        
        close(117)

        RLVTT_generic=delta_vapH/(m/1000.)

        write(*,*) 'RLVTT_generic', RLVTT_generic

        epsi_generic = 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. ! is the maximum amount of vapor that a parcel can hold without condensation, in specific concentration.
        else 
            qsat = epsi_generic *psat/(p-(1-epsi_generic)*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_generic/cpp) *qsat*(p/(p-(1-epsi_generic)*psat))*dummy
            dlnpsat = (RLVTT_generic/cpp) * dummy 
        endif

    end subroutine Lcpdqsat_generic

    !==================================================================
    subroutine Tsat_generic(p,Tsat)

        implicit none

    !==================================================================
    !     Purpose
    !     -------
    !     Compute the saturation temperature
    !     for a given pressure (Pa)
    !
    !     Authors
    !     -------
    !     Noé Clément (2022)
    !
    !==================================================================

        ! input
        real p
    
        ! output
        real Tsat


        Tsat = 1/(1/Tref - (r*mugaz/1000.)/delta_vapH*Log(p/Pref))

    end subroutine Tsat_generic
    !==================================================================

end module generic_cloud_common_h