Ignore:
Timestamp:
Oct 12, 2023, 10:30:22 AM (15 months ago)
Author:
slebonnois
Message:

BBT : Update for the titan microphysics and cloud model

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/muphytitan/mm_methods.f90

    r1897 r3083  
    1 ! Copyright 2013-2015,2017 Université de Reims Champagne-Ardenne
    2 ! Contributor: J. Burgalat (GSMA, URCA)
     1! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne
     2! Contributors: J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
    33! email of the author : jeremie.burgalat@univ-reims.fr
    4 ! 
     4!
    55! This software is a computer program whose purpose is to compute
    66! microphysics processes using a two-moments scheme.
    7 ! 
     7!
    88! This library is governed by the CeCILL-B license under French law and
    9 ! abiding by the rules of distribution of free software.  You can  use, 
     9! abiding by the rules of distribution of free software.  You can  use,
    1010! modify and/ or redistribute the software under the terms of the CeCILL-B
    1111! license as circulated by CEA, CNRS and INRIA at the following URL
    12 ! "http://www.cecill.info". 
    13 ! 
     12! "http://www.cecill.info".
     13!
    1414! As a counterpart to the access to the source code and  rights to copy,
    1515! modify and redistribute granted by the license, users are provided only
    1616! with a limited warranty  and the software's author,  the holder of the
    1717! economic rights,  and the successive licensors  have only  limited
    18 ! liability. 
    19 ! 
     18! liability.
     19!
    2020! In this respect, the user's attention is drawn to the risks associated
    2121! with loading,  using,  modifying and/or developing or reproducing the
    … …  
    2525! professionals having in-depth computer knowledge. Users are therefore
    2626! encouraged to load and test the software's suitability as regards their
    27 ! requirements in conditions enabling the security of their systems and/or 
    28 ! data to be ensured and,  more generally, to use and operate it in the 
    29 ! same conditions as regards security. 
    30 ! 
     27! requirements in conditions enabling the security of their systems and/or
     28! data to be ensured and,  more generally, to use and operate it in the
     29! same conditions as regards security.
     30!
    3131! The fact that you are presently reading this means that you have had
    3232! knowledge of the CeCILL-B license and that you accept its terms.
    … …  
    3535!! summary: Model miscellaneous methods module.
    3636!! author: J. Burgalat
    37 !! date: 2013-2015,2017
     37!! date: 2013-2015,2017,2022
     38!! corrections: B. de Batz de Trenquelléon (2023)
    3839
    3940MODULE MM_METHODS
    … …  
    4344  !!
    4445  !! All thermodynamic functions related to cloud microphysics (i.e. [[mm_methods(module):mm_lHeatX(interface)]],
    45   !! [[mm_methods(module):mm_sigX(interface)]] and [[mm_methods(module):mm_psatX(interface)]]) compute related equations 
     46  !! [[mm_methods(module):mm_sigX(interface)]] and [[mm_methods(module):mm_psatX(interface)]]) compute related equations
    4647  !! from \cite{reid1986}. A version of the book is freely available [here](http://f3.tiera.ru/3/Chemistry/References/Poling%20B.E.,%20Prausnitz%20J.M.,%20O'Connell%20J.P.%20The%20Properties%20of%20Gases%20and%20Liquids%20(5ed.,%20MGH,%202000)(ISBN%200070116822)(803s).pdf).
    4748  !!
    4849  !! The module defines the following functions/subroutines/interfaces:
    4950  !!
    50   !! | name        | description 
     51  !! | name        | description
    5152  !! | :---------: | :-------------------------------------------------------------------------------------
    5253  !! | mm_lheatx   | Compute latent heat released
    … …  
    6465  IMPLICIT NONE
    6566
    66   PRIVATE 
    67 
    68   PUBLIC  :: mm_sigX, mm_LheatX, mm_psatX, mm_qsatx, mm_fshape, &
    69              mm_get_kco, mm_get_kfm, mm_eta_g, mm_lambda_g
     67  PRIVATE
     68
     69  PUBLIC  :: mm_sigX, mm_LheatX, mm_psatX, mm_qsatx, mm_ysatX, mm_fshape, &
     70    mm_get_kco, mm_get_kfm, mm_eta_g, mm_lambda_g
    7071
    7172  ! ---- INTERFACES
    … …  
    7879  !! FUNCTION mm_sigX(temp,xESP)
    7980  !! ```
    80   !! 
    81   !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method 
     81  !!
     82  !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method
    8283  !! computes the result for all the temperatures and returns a vector of same size than __temp__.
    8384  INTERFACE mm_sigX
    8485    MODULE PROCEDURE sigx_sc,sigx_ve
    85   END INTERFACE
     86  END INTERFACE mm_sigX
    8687
    8788  !> Interface to Latent heat computation functions.
    88   !! 
     89  !!
    8990  !! The method computes the latent heat released of a given specie at given temperature(s).
    9091  !!
    … …  
    9394  !! ```
    9495  !!
    95   !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method 
     96  !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method
    9697  !! computes the result for all the temperatures and returns a vector of same size than __temp__.
    9798  INTERFACE mm_LheatX
    9899    MODULE PROCEDURE lheatx_sc,lheatx_ve
    99   END INTERFACE
     100  END INTERFACE mm_LheatX
    100101
    101102  !> Interface to saturation vapor pressure computation functions.
    … …  
    104105  !! FUNCTION mm_psatX(temp,xESP)
    105106  !! ```
    106   !! 
     107  !!
    107108  !! The method computes the saturation vapor pressure of a given specie at given temperature(s).
    108109  !!
    109   !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method 
     110  !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method
    110111  !! computes the result for all the temperatures and returns a vector of same size than __temp__.
    111112  INTERFACE mm_psatX
    112113    MODULE PROCEDURE psatx_sc,psatx_ve
    113   END INTERFACE
    114 
    115   !! Interface to saturation mass mixing ratio computaiton functions.
    116   !!
    117   !! The method computes the mass mixing ratio at saturation of a given specie at given temperature(s) 
     114  END INTERFACE mm_psatX
     115
     116  !> Interface to saturation mass mixing ratio computation functions.
     117  !!
     118  !! The method computes the mass mixing ratio at saturation of a given specie at given temperature(s)
    118119  !! and pressure level(s).
    119120  !!
    120121  !! ```fortran
    121   !! FUNCTION mm_qsatX(temp,pres,xESP) 
    122   !! ```
    123   !!
    124   !! __xESP__ must always be given as a scalar. If __temp__ and __pres__  are given as a vector (of same 
     122  !! FUNCTION mm_qsatX(temp,pres,xESP)
     123  !! ```
     124  !!
     125  !! __xESP__ must always be given as a scalar. If __temp__ and __pres__  are given as a vector (of same
    125126  !! size !), then the method computes the result for each couple of (temperature, pressure) and returns
    126127  !! a vector of same size than __temp__.
    127128  INTERFACE mm_qsatx
    128129    MODULE PROCEDURE qsatx_sc,qsatx_ve
    129   END INTERFACE
     130  END INTERFACE mm_qsatx
     131
     132  !> Interface to saturation molar mixing ratio computation functions.
     133  !!
     134  !! The method computes the molar mixing ratio at saturation of a given specie at given temperature(s)
     135  !! and pressure level(s) [Fray and Schmidt (2009)].
     136  !!
     137  !! ```fortran
     138  !! FUNCTION mm_ysatX(temp,pres,xESP)
     139  !! ```
     140  !!
     141  !! __xESP__ must always be given as a scalar. If __temp__ and __pres__  are given as a vector (of same
     142  !! size !), then the method computes the result for each couple of (temperature, pressure) and returns
     143  !! a vector of same size than __temp__.
     144  INTERFACE mm_ysatX
     145    MODULE PROCEDURE ysatX_sc,ysatX_ve
     146  END INTERFACE mm_ysatX
    130147
    131148  !> Interface to shape factor computation functions.
    … …  
    137154  !! ```
    138155  !!
    139   !! Where __m__ is cosine of the contact angle and __x__ the curvature radius. __m__ must always be 
    140   !! given as a scalar. If __x__ is given as a vector, then the method compute the result for each 
     156  !! Where __m__ is cosine of the contact angle and __x__ the curvature radius. __m__ must always be
     157  !! given as a scalar. If __x__ is given as a vector, then the method compute the result for each
    141158  !! value of __x__ and and returns a vector of same size than __x__.
    142159  INTERFACE mm_fshape
    143160    MODULE PROCEDURE fshape_sc,fshape_ve
    144   END INTERFACE
    145 
    146   CONTAINS
     161  END INTERFACE mm_fshape
     162
     163CONTAINS
    147164
    148165  FUNCTION fshape_sc(cost,rap) RESULT(res)
    149166    !! Get the shape factor of a ccn (scalar).
    150167    !!
    151     !! The method computes the shape factor for the heterogeneous nucleation on a fractal particle. 
     168    !! The method computes the shape factor for the heterogeneous nucleation on a fractal particle.
    152169    !! Details about the shape factor can be found in \cite{prup1978}.
    153170    REAL(kind=mm_wp), INTENT(in) :: cost !! Cosine of the contact angle.
    … …  
    160177      phi = dsqrt(1._mm_wp-2._mm_wp*cost*rap+rap**2)
    161178      a = 1._mm_wp + ( (1._mm_wp-cost*rap)/phi )**3
    162       b = (rap**3) * (2._mm_wp-3._mm_wp*(rap-cost)/phi+((rap-cost)/phi)**3)
     179      b = (rap**3) * (2._mm_wp - 3._mm_wp*(rap-cost)/phi + ((rap-cost)/phi)**3)
    163180      c = 3._mm_wp * cost * (rap**2) * ((rap-cost)/phi-1._mm_wp)
    164181      res = 0.5_mm_wp*(a+b+c)
    … …  
    177194    WHERE(rap > 3000._mm_wp)
    178195      res = ((2._mm_wp+cost)*(1._mm_wp-cost)**2)/4._mm_wp
    179     ELSEWHERE 
     196    ELSEWHERE
    180197      phi = dsqrt(1._mm_wp-2._mm_wp*cost*rap+rap**2)
    181198      a = 1._mm_wp + ((1._mm_wp-cost*rap)/phi )**3
    … …  
    192209    !! The method computes the latent heat equation as given in \cite{reid1986} p. 220 (eq. 7-9.4).
    193210    IMPLICIT NONE
    194     ! - DUMMY 
     211    ! - DUMMY
    195212    REAL(kind=mm_wp), INTENT(in) :: temp !! temperature (K).
    196213    TYPE(mm_esp), INTENT(in)     :: xESP !! Specie properties.
    197214    REAL(kind=mm_wp) :: res              !! Latent heat of given specie at given temperature (\(J.kg^{-1}\)).
    198215    REAL(kind=mm_wp) :: ftm
    199     ftm=MIN(1._mm_wp-temp/xESP%tc,1.e-3_mm_wp)
     216    ftm=MAX(1._mm_wp-temp/xESP%tc,1.e-3_mm_wp)
    200217    res = mm_rgas*xESP%tc*(7.08_mm_wp*ftm**0.354_mm_wp+10.95_mm_wp*xESP%w*ftm**0.456_mm_wp)/xESP%masmol
    201218  END FUNCTION LHeatX_sc
    202    
     219
    203220  FUNCTION LHeatX_ve(temp,xESP) RESULT(res)
    204221    !! Compute latent heat of a given specie at given temperature (vector).
    … …  
    208225    TYPE(mm_esp), INTENT(in)                   :: xESP !! Specie properties.
    209226    REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  !! Latent heat of given specie at given temperatures (\(J.kg^{-1}\)).
    210     REAL(kind=mm_wp) :: ftm 
     227    REAL(kind=mm_wp) :: ftm
    211228    INTEGER      :: i
    212229    DO i=1,SIZE(temp)
    213       ftm=MIN(1._mm_wp-temp(i)/xESP%tc,1.e-3_mm_wp)
     230      ftm=MAX(1._mm_wp-temp(i)/xESP%tc,1.e-3_mm_wp)
    214231      res(i) = mm_rgas*xESP%tc*(7.08_mm_wp*ftm**0.354_mm_wp+10.95_mm_wp*xESP%w*ftm**0.456_mm_wp) / &
    215                xESP%masmol
     232        xESP%masmol
    216233    ENDDO
    217234  END FUNCTION LHeatX_ve
    … …  
    219236  FUNCTION sigX_sc(temp,xESP) RESULT(res)
    220237    !! Get the surface tension between a given specie and the air (scalar).
    221     !! 
     238    !!
    222239    !! The method computes the surface tension equation as given in \cite{reid1986} p. 637 (eq. 12-3.6).
    223240    REAL(kind=mm_wp), INTENT(in) :: temp !! temperature (K).
    … …  
    229246    sig = 0.1196_mm_wp*(1._mm_wp+(tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-tbr))-0.279_mm_wp
    230247    sig = xESP%pc**(2._mm_wp/3._mm_wp)*xESP%tc**(1._mm_wp/3._mm_wp)*sig*(1._mm_wp-tr)**(11._mm_wp/9._mm_wp)
    231     res = sig*1e3_mm_wp ! dyn/cm2 -> N/m
     248    res = sig*1e-3_mm_wp ! dyn/cm -> N/m
    232249  END FUNCTION sigX_sc
    233  
     250
    234251  FUNCTION sigX_ve(temp,xESP) RESULT(res)
    235252    !! Get the surface tension between a given specie and the air (vector).
    … …  
    240257    REAL(kind=mm_wp), DIMENSION(SIZE(temp)) :: res     !! Surface tensions (\(N.m^{-1}\)).
    241258    INTEGER      :: i
    242     REAL(kind=mm_wp) :: tr,tbr,sig
     259    REAL(kind=mm_wp) :: tr,tbr,sig0,sig
    243260    tbr = xESP%tb/xESP%tc
    244     sig = 0.1196_mm_wp*(1._mm_wp+(tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-tbr))-0.279_mm_wp
     261    sig0 = 0.1196_mm_wp*(1._mm_wp+(tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-tbr))-0.279_mm_wp
    245262    DO i=1,SIZE(temp)
    246263      tr     = MIN(temp(i)/xESP%tc,0.99_mm_wp)
    247       sig    = xESP%pc**(2._mm_wp/3._mm_wp)*xESP%tc**(1._mm_wp/3._mm_wp)*sig*(1._mm_wp-tr)**(11._mm_wp/9._mm_wp)
    248       res(i) = sig*1e3_mm_wp ! dyn/cm2 -> N/m
     264      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)
     265      res(i) = sig*1e-3_mm_wp ! dyn/cm -> N/m
    249266    ENDDO
    250267  END FUNCTION sigX_ve
    … …  
    252269  FUNCTION psatX_sc(temp,xESP) RESULT(res)
    253270    !! Get saturation vapor pressure for a given specie at given temperature (scalar).
    254     !! 
     271    !!
    255272    !! The method computes the saturation vapor pressure equation given in \cite{reid1986} p. 657 (eq. 1).
    256     !!
    257     !! @warning
    258     !! This subroutine accounts for a specific Titan feature:
    259     !! If __xESP__ corresponds to \(CH_{4}\), the saturation vapor presure is multiplied by 0.85
    260     !! to take into account its dissolution in \(N_{2}\).
    261273    REAL(kind=mm_wp), INTENT(in) :: temp !! Temperature (K).
    262274    TYPE(mm_esp), INTENT(in)     :: xESP !! Specie properties.
    … …  
    266278    IF (x > 0._mm_wp) THEN
    267279      qsat = (1._mm_wp-x)**(-1) * &
    268       (xESP%a_sat*x + xESP%b_sat*x**1.5_mm_wp + xESP%c_sat*x**3 + xESP%d_sat*x**6)
     280        (xESP%a_sat*x + xESP%b_sat*x**1.5_mm_wp + xESP%c_sat*x**2.5_mm_wp + xESP%d_sat*x**5_mm_wp)
    269281    ELSE
    270       qsat = XESP%a_sat*x/abs(1._mm_wp-x)     ! approx for  t > tc
    271     ENDIF
    272     !  Special case : ch4 : x0.85 (dissolution in N2)
    273     IF (xESP%name == "ch4") THEN
    274       res = xESP%pc*dexp(qsat)*0.85_mm_wp
    275     ELSE
    276       res = xESP%pc*dexp(qsat)
    277     ENDIF
     282      qsat = XESP%a_sat*x/abs(1._mm_wp-x) ! approx for t > tc
     283    ENDIF
     284    res = xESP%pc*exp(qsat)
    278285    ! now convert bar to Pa
    279286    res = res * 1e5_mm_wp
    … …  
    282289  FUNCTION psatX_ve(temp,xESP) RESULT(res)
    283290    !! Get saturation vapor pressure for a given specie at given temperature (vector).
    284     !! 
     291    !!
    285292    !! See [[mm_methods(module):psatX_sc(function)]].
    286293    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp !! Temperatures (K).
    … …  
    292299      x = 1._mm_wp-temp(i)/xESP%tc
    293300      IF (x > 0._mm_wp) THEN
    294         qsat = (1._mm_wp-x)**(-1) * & 
    295         (xESP%a_sat*x + xESP%b_sat*x**1.5_mm_wp + xESP%c_sat*x**3 + xESP%d_sat*x**6)
     301        qsat = (1._mm_wp-x)**(-1) * &
     302          (xESP%a_sat*x + xESP%b_sat*x**1.5_mm_wp + xESP%c_sat*x**2.5_mm_wp + xESP%d_sat*x**5_mm_wp)
    296303      ELSE
    297         qsat = XESP%a_sat*x/abs(1._mm_wp-x)     ! approx for t > tc
     304        qsat = XESP%a_sat*x/abs(1._mm_wp-x) ! approx for t > tc
    298305      ENDIF
    299       res(i) = xESP%pc*dexp(qsat)
    300       !  Peculiar case : ch4 : x0.85 (dissolution in N2)
    301       IF (xESP%name == "ch4") res(i) = res(i)* 0.85_mm_wp
     306      res(i) = xESP%pc*exp(qsat)
    302307    ENDDO
    303     res = res * 1e5_mm_wp  ! bar -> Pa
     308    ! now convert bar to Pa
     309    res = res * 1e5_mm_wp
    304310  END FUNCTION psatX_ve
    305311
    306312  FUNCTION qsatX_sc(temp,pres,xESP) RESULT(res)
    307313    !! Get the mass mixing ratio of a given specie at saturation (scalar).
     314    !!
     315    !! @warning
     316    !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
     317    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
    308318    REAL(kind=mm_wp), INTENT(in) :: temp !! Temperature (K).
    309319    REAL(kind=mm_wp), INTENT(in) :: pres !! Pressure level (Pa).
    310320    TYPE(mm_esp), INTENT(in)    :: xESP  !! Specie properties.
    311321    REAL(kind=mm_wp) :: res              !! Mass mixing ratio of the specie.
    312     REAL(kind=mm_wp) :: x,psat
     322    REAL(kind=mm_wp) :: psat
    313323    psat = mm_psatX(temp,xESP)
    314324    res = (psat / pres) * xESP%fmol2fmas
     325    ! Peculiar case : CH4 : x0.85 (dissolution in N2)
     326    IF (xESP%name == "CH4") THEN
     327      res = res * 0.85_mm_wp
     328      IF (mm_debug) WRITE(*,'(a)') "[DEBUG] mm_qsat: applying .85 factor to qsat for CH4 specie (N2 dissolution)"
     329    ENDIF
    315330  END FUNCTION qsatX_sc
    316331
    317332  FUNCTION qsatX_ve(temp,pres,xESP) RESULT(res)
    318333    !! Get the mass mixing ratio of a given specie at saturation (vector).
     334    !!
     335    !! @warning
     336    !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
     337    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
    319338    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp !! Temperatures (K).
    320339    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pres !! Pressure levels (Pa).
    … …  
    324343    psat = mm_psatX(temp,xESP)
    325344    res = (psat / pres) * xESP%fmol2fmas
     345    ! Peculiar case : CH4 : x0.85 (dissolution in N2)
     346    IF (xESP%name == "CH4") THEN
     347      res = res * 0.85_mm_wp
     348      IF (mm_debug) WRITE(*,'(a)') "[DEBUG] mm_qsat: applying .85 factor to qsat for CH4 specie (N2 dissolution)"
     349    ENDIF
    326350  END FUNCTION qsatX_ve
     351
     352  FUNCTION ysatX_sc(temp,pres,xESP) RESULT(res)
     353    !! Get the molar mixing ratio of a given specie at saturation (scalar).
     354    !!
     355    !! @warning
     356    !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
     357    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
     358    REAL(kind=mm_wp), INTENT(in) :: temp !! Temperature (K).
     359    REAL(kind=mm_wp), INTENT(in) :: pres !! Pressure level (Pa).
     360    TYPE(mm_esp), INTENT(in)     :: xESP !! Specie properties.
     361    REAL(kind=mm_wp)             :: res  !! Molar mixing ratio of the specie.
     362
     363    ! Fray and Schmidt (2009)
     364    IF(xESP%name == "C2H2") THEN
     365      res = (1.0e5 / pres) * exp(1.340e1 - 2.536e3/temp)
     366   
     367    ELSE IF(xESP%name == "C2H6") THEN
     368      res = (1.0e5 / pres) * exp(1.511e1 - 2.207e3/temp - 2.411e4/temp**2 + 7.744e5/temp**3 - 1.161e7/temp**4 + 6.763e7/temp**5)
     369   
     370    ELSE IF(xESP%name == "HCN") THEN
     371      res = (1.0e5 / pres) * exp(1.393e1 - 3.624e3/temp - 1.325e5/temp**2 + 6.314e6/temp**3 - 1.128e8/temp**4)
     372   
     373    ELSE IF (xESP%name == "CH4") THEN
     374      res = (1.0e5 / pres) * exp(1.051e1 - 1.110e3/temp - 4.341e3/temp**2 + 1.035e5/temp**3 - 7.910e5/temp**4)
     375      res = res * 0.85_mm_wp
     376      !IF (res < 0.014) THEN
     377      !  res = 0.014
     378      !ENDIF
     379    ENDIF
     380  END FUNCTION ysatX_sc
     381
     382  FUNCTION ysatX_ve(temp,pres,xESP) RESULT(res)
     383    !! Get the molar mixing ratio of a given specie at saturation (vector).
     384    !!
     385    !! @warning
     386    !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
     387    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
     388    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp !! Temperatures (K).
     389    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pres !! Pressure levels (Pa).
     390    TYPE(mm_esp), INTENT(in)                   :: xESP !! Specie properties.
     391    REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  !! Molar mixing ratios of the specie.
     392
     393    ! Fray and Schmidt (2009)
     394    IF(xESP%name == "C2H2") THEN
     395      res = (1.0e5 / pres) * exp(1.340e1 - 2.536e3/temp)
     396   
     397    ELSE IF(xESP%name == "C2H6") THEN
     398      res = (1.0e5 / pres) * exp(1.511e1 - 2.207e3/temp - 2.411e4/temp**2 + 7.744e5/temp**3 - 1.161e7/temp**4 + 6.763e7/temp**5)
     399   
     400    ELSE IF(xESP%name == "HCN") THEN
     401      res = (1.0e5 / pres) * exp(1.393e1 - 3.624e3/temp - 1.325e5/temp**2 + 6.314e6/temp**3 - 1.128e8/temp**4)
     402   
     403    ! Peculiar case : CH4 : x0.85 (dissolution in N2)
     404    ELSE IF (xESP%name == "CH4") THEN
     405      res = (1.0e5 / pres) * exp(1.051e1 - 1.110e3/temp - 4.341e3/temp**2 + 1.035e5/temp**3 - 7.910e5/temp**4)
     406      res = res * 0.85_mm_wp
     407      !WHERE (res(:) < 0.014) res(:) = 0.014
     408    ENDIF
     409  END FUNCTION ysatX_ve
    327410
    328411  ELEMENTAL FUNCTION mm_get_kco(t) RESULT(res)
Note: See TracChangeset for help on using the changeset viewer.