Ignore:
Timestamp:
Feb 27, 2025, 11:35:32 AM (4 months ago)
Author:
gmilcareck
Message:

Generic PCM:
Rewriting of the rayleigh scattering module. Before the commit, the Rayleigh scattering was separated into a constant part and a wavelength-dependent part.
However, in the constant part, because everything was added up/multiplied together before being added to the module,
we lose the information about where these numbers come from and if we want to add a new molecule or update a molecule... well... good luck.
Now, each molecule depends only on 2 physical parameters: the refractive index and the King Factor.
If you want add another molecule, you just need to add these 2 parameters and the references (please).
For water, the rayleigh scattering now depends on temperature and pressure.
GM

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/blackl.F

    r959 r3654  
    2020      return
    2121      end
     22
     23      subroutine blackn(blalong,blat,blae)
     24
     25      implicit double precision (a-h,o-z)
     26
     27      ! physical constants
     28      sigma=5.67032D-8
     29      pi=datan(1.d0)*4.d0
     30      c0=2.9979d+08
     31      h=6.6262d-34
     32      cbol=1.3806d-23
     33      rind=1.d0
     34      c=c0/rind
     35      c1=h*(c**2)
     36      c2=h*c/cbol
     37
     38
     39      blae=2.d0*pi*c1*blalong**3/(dexp(c2*blalong/blat)-1.d0)
     40
     41
     42      return
     43      end

  • trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90

    r3641 r3654  
    1       subroutine calc_rayleigh
     1      subroutine calc_rayleigh(qvar,muvar,PMID,TMID,tauray)
    22
    33!==================================================================
    … …  
    77!     Average the Rayleigh scattering in each band, weighting the
    88!     average by the blackbody function at temperature tstellar.
    9 !     Works for an arbitrary mix of gases (currently CO2, N2 and
    10 !     H2 are possible).
     9!     Works for an arbitrary mix of gases.
    1110!     
    1211!     Authors
    … …  
    1413!     Robin Wordsworth (2010)
    1514!     Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995).
    16 !     Noé Clément (2022) : Additionnal comments & Methane+CO Rayleigh
    17 !
     15!     Noe Clement (2022) : Additionnal comments & Methane+CO Rayleigh
     16!     Gwenael Milcareck (2025): Rewriting the code
     17!
    1818!     Called by
    1919!     ---------
    … …  
    2626!==================================================================
    2727
    28       use radinc_h, only: L_NSPECTV
    29       use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep
    30       use gases_h, only: ngasmx, vgas, gnom, gfrac, igas_CO2, igas_H2, &
    31                            igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO, &
    32                            igas_O2
    33       use comcstfi_mod, only: g, mugaz, pi
     28      use radinc_h, only: L_NSPECTV, L_LEVELS
     29      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, scalep
     30      use gases_h, only: ngasmx, vgas, gnom, gfrac, massmol, igas_CO2, igas_H2, &
     31                           igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO, igas_Ar, igas_O2
     32      use comcstfi_mod, only: g, pi
     33      use callkeys_mod, only: strictboundrayleigh
    3434
    3535      implicit none
    3636
    37       real*8 wl
    38       integer N,Nfine,ifine,igas
     37      real, intent(in) :: qvar(L_LEVELS) ! mol/mol
     38      real, intent(in) :: muvar(L_LEVELS) ! g/mol
     39      real, intent(in) :: PMID(L_LEVELS) ! mbar
     40      real, intent(in) :: TMID(L_LEVELS) ! K
     41      real, intent(out) :: tauray(L_LEVELS,L_NSPECTV)
     42      real*8 wl,wn
     43      integer N,Nfine,ifine,igas,k
    3944      parameter(Nfine=500.0)
    40       real*8 :: P0 = 9.423D+6   ! Rayleigh scattering reference pressure in pascals. P0 = 93*101325 Pa   
    41 
    42       logical typeknown
    43       real*8 tauvar,tauvarvar,tausum,tausumvar,tauwei,tauweivar,bwidth,bstart
     45      real*8 :: Fk ! King factor for the depolarization
     46      real*8 :: ng(L_LEVELS) ! real refractive index
     47      real*8 :: P0(L_LEVELS) ! reference pressure
     48      real*8 :: T0(L_LEVELS) ! reference temperature
     49     
     50      ! Parameters for H2O
     51      real*8 :: a0, a1, a2, a3, a4, a5, a6, a7
     52      real*8 :: luv, lir
     53      real*8 :: rhor,tr,lr
     54      real*8 :: rho(L_LEVELS),rhos(L_LEVELS),ts(L_LEVELS)
     55      real*8 :: b(L_LEVELS)
     56
     57      real*8 mass_frac(ngasmx,L_LEVELS)
     58      real*8 tauvar(L_LEVELS),tausum(L_LEVELS)
     59      real*8 tauwei,bwidth,bstart
    4460      double precision df
    4561
    46       real*8 tauconsti(ngasmx)
    47       real*8 tauvari(ngasmx)
    48 
     62      real*8 tauconsti(ngasmx,L_LEVELS)
     63      real*8 tauvari(ngasmx,L_LEVELS)
     64     
     65      ! Miscellaneous :
     66      character(len=200) :: message
     67      character(len=10),parameter :: subname="rayleigh"
     68      logical, save :: firstcall=.true.
     69     
    4970      integer icantbewrong
    5071
    51       ! we multiply by scalep here (100) because plev, which is used in optcv,
    52       ! is in units of mBar, so we need to convert.
     72      ! This module calculates the Rayleigh scattering (also known as the Cabannes peak)
     73      ! Rayleigh wings, Brillouin scattering and Raman scattering are not taken into account.
     74
    5375      ! we calculate here TAURAY which is in m2/mBar
    5476
    55       ! tau0/p0=tau/p (Hansen 1974 : equation (2.31) )
    56       ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4)      (Hansen 1974 : equation (2.29) )
    57       ! where delta=n-1 and N0 is an amagat
    58 
    59       ! In exo_k : (8*pi^3)/(3*N0^2) * 4 * delta^2 ------- is written : (24*pi^3)/(N0^2) * (4/9) * delta^2
    60 
    61       ! (tauconsti * tauvari) = scalep * cross_section_molecule_exo_k /  ( gravity * molecular_mass )
    62       ! (tauconsti * tauvari) in m2/mBar because of scalep factor
    63       ! cross_section_molecule in m2/molecule
    64       ! molecular_mass is the mass of th considered molecule
    65 
    66       typeknown=.false.
     77      ! The cross section for ith particles of small size compared to the wavenumber
     78      ! and in the electric dipole approximation is:
     79      ! sigma_i = 24*pi**3*wn**4/N**2 * ((n_i(wn)**2 - 1)/(n_i(wn)**2 + 2))**2 * Fk_i(wn)
     80      ! nu is the wavenumber
     81      ! N is the number density of the gas (molecule/m3)
     82      ! n_i is the real refractive index of the ith gas
     83      ! Fk_i is the King factor of ith gas equals to (6+3*delta_i)/(6-7*delta_i)
     84      ! where delta_i is the depolarization factor of the ith gas
     85     
     86      ! The rayleigh opacity is expressed by:
     87      ! tau_r = P/(g*mu) * sum_{i=1}^Ntot [ x_i*sigma_i ]
     88      ! P is the pressure
     89      ! g is the standard gravity
     90      ! mu is the mean molecular weight
     91      ! x_i is the mass fraction of the ith gas
     92      ! The pressure P dependence is calculated in optcv.F90
     93     
     94      if(firstcall) then
     95       
     96        if ((BWNV(L_NSPECTV+1).gt.60000.).and.(strictboundrayleigh)) then
     97          message="Rayleigh scattering is unknown for wn>60000 cm-1 - some singularities are present for higher wavenumber - if you know what you are doing, use strictboundrayleigh=.false."
     98          call abort_physic(subname,message,1)
     99        elseif ((BWNV(L_NSPECTV+1).gt.60000.).and..not.(strictboundrayleigh)) then
     100          print*,'**********************************************'
     101          print*,' we allow model to continue with wn>60000 cm-1'
     102          print*,' ... we assume we know what you are doing ... '
     103          print*,' ... but do not let this happen too often ... '
     104          print*,'**********************************************'
     105        endif
     106        firstcall = .false.
     107      endif
     108     
     109
    67110      do igas=1,ngasmx
    68          if(igas.eq.vgas)then
     111         ! Convert qvar mol/mol -> kg/kg
     112         if((igas.eq.vgas).and.(maxval(QVAR(:)).ge.1.e-2))then
    69113            print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering '
     114            mass_frac(igas,:) = QVAR(:)*massmol(igas)/muvar(:)
     115         elseif((igas/=vgas).and.(gfrac(igas).ge.1.e-2))then
     116            mass_frac(igas,:) = gfrac(igas)*(1.-QVAR(:))*massmol(igas)/muvar(:)
     117         else
     118            print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
     119            'as its mixing ratio is less than 0.01.'
     120            ! ignore variable gas in Rayleigh calculation
     121            ! ignore gases of mixing ratio < 0.01 in Rayleigh calculation
     122            mass_frac(igas,:) = 0.0
    70123         endif
    71          if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
    72             print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
    73             'as its mixing ratio is less than 0.05.'
    74             ! ignore variable gas in Rayleigh calculation
    75             ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
    76             tauconsti(igas) = 0.0
    77          else
    78             if(igas.eq.igas_CO2) then
    79                ! Hansen 1974 : equation (2.32)
    80                tauconsti(igas) = (8.7/g)*1.527*scalep/P0
    81             elseif(igas.eq.igas_N2)then
    82                ! Hansen 1974 : equation (2.30) ---- (P0/93.0) is equal to 101325 Pa
    83                tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
    84             elseif(igas.eq.igas_H2O)then
    85                ! tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0
    86                tauconsti(igas) = 1.98017E-10*scalep/(g*mugaz*1E4)
    87             elseif(igas.eq.igas_H2)then
    88                tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0
    89                ! tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
    90                ! uses n=1.000132 from Optics, Fourth Edition.
    91                ! was out by a factor 2!
    92 
    93                ! below is exo_k formula : this tauconsti is consistent with commented tauvari for H2 below
    94                ! tauconsti(igas) = (1/(g*2.*1.6726*1E-27)) * 1E16 * scalep
    95             elseif(igas.eq.igas_He)then
    96                tauconsti(igas) = (10.0/g)*0.00086*scalep/101325.0
    97             elseif(igas.eq.igas_CH4)then
    98                ! For CH4 we use the formalism of exo_k (developed by Jeremy Leconte)
    99                ! the relation between exo_k formalism and LMDZ formalism is : (tauconsti*tauvari)=scalep*sigma_exok/(gravity*molecular_mass)
    100                ! Values are taken from Caldas (2019), equation 12 & appendix D
    101                ! 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 is a constant
    102                ! 4/9=(2/3)**2 is approximately ((n+1)/(n**2+2))**2
    103                ! 1E24 = (1E6)**4 because wavelength is in micron
    104                ! 16.04*1.6726*1E-27 is CH4 molecular mass
    105                tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep
    106            
    107             elseif(igas.eq.igas_CO)then
    108                ! see above for explanation
    109                ! 28.01*1.6726*1E-27 is CH4 molecular mass
    110                tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *28.01*1.6726*1e-27))  * scalep
    111             else
    112                print*,'Error in calc_rayleigh: Gas species not recognised!'
    113                print*,"igas=",igas," gnom(igas)=",trim(gnom(igas))
    114                call abort_physic("calc_rayleigh","Invalid gas",1)
    115             endif
    116 
    117             if((gfrac(igas).eq.1.0).and.(vgas.eq.0))then
    118                print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.'
    119                typeknown=.true.
    120             endif
    121 
    122          endif
     124         tauvari(igas,:) = 0.
    123125      enddo
    124 
    125       if(.not.typeknown)then
    126          print*,'Rayleigh scattering is for a mixed gas atmosphere.'
    127          typeknown=.true.
    128       endif
    129 
     126     
     127      ! ATTENTION, au delà de 60000 cm-1, pour toutes les molécules, il y a des singularités en lien avec la formule d'interpolation.
     128     
    130129   
    131130      do N=1,L_NSPECTV
     131     
     132         ! The refractive index depend on temperature and pressure
     133         ! It isn't the case here. Must be implemented in the future...
     134         ! But in the current scientific litterature (2024), it's difficult
     135         ! to find something that depends on temperature and pressure...
     136         ! except for H2O
    132137         
    133138         tausum = 0.0
    134139         tauwei = 0.0
    135          tausumvar = 0.0
    136          tauweivar = 0.0
    137140         bstart = 10000.0/BWNV(N+1) ! BWNV is in cm-1 so 10000.0/BWNV is in micron
    138141         bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1))
    139142         do ifine=1,Nfine
    140143            wl=bstart+dble(ifine)*bwidth/Nfine
    141 
    142             tauvar=0.0
    143             tauvarvar=0.0
     144            wn=BWNV(N)+dble(ifine)*(BWNV(N+1)-BWNV(N))/Nfine
     145
     146            tauvar(:)=0.0
    144147            do igas=1,ngasmx
    145                if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
    146                   ! ignore variable gas in Rayleigh calculation
    147                   ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
    148                   tauvari(igas)   = 0.0
    149                else
    150                   if(igas.eq.igas_CO2)then
    151                      tauvari(igas) = (1.0+0.013/wl**2)/wl**4
    152                   elseif(igas.eq.igas_N2)then
    153                      tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
    154                   elseif(igas.eq.igas_H2O)then
    155                      ! tauvari(igas) = 1.0/wl**4 ! to be improved...
    156                      tauvari(igas) = (5.7918E6/(2.38E2-1/wl**2)+1.679E5/(57.36E0-1/wl**2))**2/wl**4
    157                   elseif(igas.eq.igas_H2)then
    158                      tauvari(igas) = 1.0/wl**4
    159 
    160                      ! below is exo_k formula : this tauvari is consistent with commented tauconsti for H2 above
    161                      ! tauvari(igas) = (8.14E-49+1.28E-50/wl**2+1.61E-51/wl**4) / wl**4
    162                   elseif(igas.eq.igas_He)then
    163                      tauvari(igas) = 1.0/wl**4
    164                   elseif(igas.eq.igas_CH4)then
    165                      tauvari(igas) = (4.6662E-4+4.02E-6/wl**2)**2 / wl**4
    166                   elseif(igas.eq.igas_CO)then
    167                      tauvari(igas) = (2.2851E-4 + 0.456E4/(5.101816E9 - 1E8/wl**2))**2 / wl**4
    168                   else
    169                      print*,'Error in calc_rayleigh: Gas species not recognised!'
    170                      call abort_physic("calc_rayleigh","Gas species not recognised",1)
    171                   endif
    172                endif
    173 
    174                if(igas.eq.vgas) then
    175                   tauvarvar=tauvarvar+tauconsti(igas)*tauvari(igas)
    176                   tauvar=tauvar+0.0*0.0*gfrac(igas)
    177                else
    178                   tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas)
    179                endif
    180 
    181             enddo
    182 
    183             call blackl(dble(wl*1e-6),dble(tstellar),df)
    184             df=df*bwidth/Nfine
     148               if (maxval(mass_frac(igas,:)).ge.1e-2) then
     149                 
     150                 if(igas.eq.igas_CO2)then
     151                     ! Sneep et al, 2005
     152                     ! doi:10.1016/j.jqsrt.2004.07.025
     153                     T0(:) = 288.15
     154                     P0(:) = 1.01325e5
     155                     ! ng -> valid range of the measurements : 0.1807 - 1.8172 um
     156                     ng(:) = 1. + 1.1427e3*(5799.25/(128908.9**2 - wn**2) + 120.05/(89223.8**2 - wn**2) + 5.3334/(75037.5**2 - wn**2) + 4.3244/(67837.7**2 - wn**2) + 0.1218145e-4/(2418.136**2 - wn**2)) ! there is an error on the paper 1.1427e6 -> 1.1427e3
     157                     Fk = 1.1364 + 25.3e-12*wn**2
     158                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4 ! wn*100 -> cm-1 to m-1
     159                 elseif(igas.eq.igas_N2)then
     160                     ! Thalman et al, 2014
     161                     ! doi:10.1016/j.jqsrt.2014.05.030
     162                     T0(:) = 288.15
     163                     P0(:) = 1.01325e5
     164                     if(wn.gt.21360)then !between 21360 and 39370 cm-1. We extrapolate above.
     165                       ng(:) = 1. + (5677.465 + 318.81874e13/(14.4e9 - wn**2))*1e-8  !there is an error on the paper e12 -> e13
     166                     else !between 4860 and 21360 cm-1. We extrapolate below.
     167                       ng(:) = 1. + (6498.2 + 307.4335e13/(14.4e9 - wn**2))*1e-8
     168                     endif
     169                     Fk = 1.034 + 3.17e-12*wn
     170                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     171                 elseif(igas.eq.igas_H2O)then
     172                     ! Harvey et al, 1998
     173                     ! doi:10.1063/1.556029
     174                     ! doi:10.1364/AO.35.001566 for wn<4840 cm-1
     175                     a0 =  0.244257733
     176                     a1 = 9.74634476e-3
     177                     a2 = -3.73234996e-3
     178                     a3 = 2.68678472e-4
     179                     a4 = 1.58920570e-3
     180                     a5 = 2.45934259e-3
     181                     a6 = 0.900704920
     182                     a7 = -1.66626219e-2
     183                     luv = 0.2292020
     184                     lir = 5.432937
     185                     Tr = 273.15
     186                     rhor = 1000.
     187                     T0(:) = tmid(:)
     188                     P0(:) = pmid(:)*scalep
     189                     lr = 0.589
     190                     rho(:) = mass_frac(igas,:)*muvar(:)/massmol(igas)*P0(:)/(8.314462*T0(:)/muvar(:)/1000.)
     191                     rhos(:) = rho(:)/rhor
     192                     ts(:) = T0(:)/Tr
     193                     
     194                     b(:) = (a0 + a1*rhos(:) + a2*ts(:) + a3*ts(:)*(10000./wn/lr)**2 + a4/(10000./wn/lr)**2 + a5/((10000./wn/lr)**2 - luv**2) + a6/((10000./wn/lr)**2 - lir**2) + a7*rhos(:)**2)*rhos(:)
     195                     ! ng -> valid range of the measurements : 0.2 - 1.1 um
     196                     ng(:) = sqrt(2.*b(:)+1.)/sqrt(1.-b(:))
     197                     if(wn<4840.) then ! necessary to prevent a singularity at 3230 cm-1
     198                       ng(:) = 1. + 1.022e-8*(295.235 + 2.6422*(wn*1e-4)**2 - 0.032380*(wn*1e-4)**4 + 0.004028*(wn*1e-4)**6)
     199                     endif
     200                     Fk = (6.+3.*3e-4)/(6.-7.*3e-4) ! delta=3e-4 Murphy 1977 doi:10.1063/1.434794
     201                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     202                 elseif(igas.eq.igas_H2)then
     203                     ! Peck and Hung, 1977
     204                     ! doi:10.1364/JOSA.67.001550
     205                     T0(:) = 273.15
     206                     P0(:) = 1.01325e5
     207                     ! ng -> valid range of the measurements : 0.1680 - 1.6945 um
     208                     ng(:) = 1. + (14895.6/(180.7 - (wn*1e-4)**2) + 4903.7/(92.-(wn*1e-4)**2))*1e-6
     209                     Fk = (6.+3.*0.02)/(6.-7.*0.02) ! delta=0.02 Hansen 1974
     210                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     211                 elseif(igas.eq.igas_He)then
     212                     ! Thalman et al, 2014
     213                     ! doi:10.1016/j.jqsrt.2014.05.030
     214                     T0(:) = 288.15
     215                     P0(:) = 1.01325e5
     216                     ! ng -> valid range of the measurements : 0.2753 - 20.5813 um
     217                     ng(:) = 1. + (2283. + 1.8102e13/(1.5342e10 - wn**2))*1e-8
     218                     Fk = 1.
     219                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     220                 elseif(igas.eq.igas_CH4)then
     221                     ! Sneep et al, 2005
     222                     ! doi:10.1016/j.jqsrt.2004.07.025
     223                     T0(:) = 288.15
     224                     P0(:) = 1.01325e5
     225                     ! ng -> valid range of the measurements : 0.3251 - 0.6330 um
     226                     ng(:) = 1. + 46662e-8 + 4.02e-14*wn**2
     227                     Fk = 1.
     228                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     229                 elseif(igas.eq.igas_CO)then
     230                     ! Sneep et al, 2005
     231                     ! doi:10.1016/j.jqsrt.2004.07.025
     232                     T0(:) = 288.15
     233                     P0(:) = 1.01325e5
     234                     ! ng -> valid range of the measurements : 0.168 - 0.288 um
     235                     ng(:) = 1. + 22851e-8 + 0.456e4/(71427.**2 - wn**2)
     236                     Fk = 1.016
     237                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     238                 elseif(igas.eq.igas_Ar)then
     239                     ! Thalman et al, 2014
     240                     ! doi:10.1016/j.jqsrt.2014.05.030
     241                     T0(:) = 288.15
     242                     P0(:) = 1.01325e5
     243                     ! ng -> valid range of the measurements : 0.288 - 0.546 um
     244                     ng(:) = 1. + (6432.135 + 286.06021e12/(14.4e9 - wn**2))*1e-8
     245                     Fk = 1.
     246                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     247                 elseif(igas.eq.igas_O2)then
     248                     ! Thalman et al, 2014
     249                     ! doi:10.1016/j.jqsrt.2014.05.030
     250                     T0(:) = 273.15
     251                     P0(:) = 1.01325e5
     252                     ! ng -> valid range of the measurements : 0.303 - 2.0 um
     253                     ng(:) = 1. + (20564.8 + 2.480899e13/(4.09e9 - wn**2))*1e-8
     254                     Fk = 1.09 + 1.385e-11*wn**2 + 1.488e-20*wn**4
     255                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     256                 else
     257                     print*,'No rayleigh scattering for ',trim(gnom(igas)),'. No data found.'
     258                 endif
     259                 
     260                 tauconsti(igas,:) = 24.*pi**3 *6.022140857E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380648813E-23*T0(:)))**2)
     261                 ! N=P/(kB*T)
     262                 ! pmid*scalep -> mbar to Pa
     263                 ! muvar/1000 -> g/mol to kg/mol
     264               
     265                 tauvar(:)=tauvar(:)+tauconsti(igas,:)*tauvari(igas,:)
     266                 
     267               endif !greater than 0.01
     268
     269            enddo !ngasmx
     270
     271            call blackn(dble(wn*1e2),dble(tstellar),df)
     272            df=df*(BWNV(N+1)-BWNV(N))/Nfine
    185273            tauwei=tauwei+df
    186             tausum=tausum+tauvar*df
    187             tauweivar=tauweivar+df
    188             tausumvar=tausumvar+tauvarvar*df
     274            tausum(:)=tausum(:)+tauvar(:)*df
    189275         
    190          enddo
    191          TAURAY(N)=tausum/tauwei
    192          TAURAYVAR(N)=tausumvar/tauweivar
    193          ! we multiply by scalep here (100) because plev, which is used in optcv,
    194          ! is in units of mBar, so we need to convert.
    195 
    196       end do
    197 
    198       IF (L_NSPECTV > 6) THEN
    199          icantbewrong = L_NSPECTV-6
    200          print*,'At 1 atm and lambda = ',WAVEV(icantbewrong),' um'
    201          print*,'tau_R = ',TAURAY(icantbewrong)*1013.25
    202          print*,'sig_R = ',TAURAY(icantbewrong)*g*mugaz*1.67e-27*100, &
    203                'cm^2 molecule^-1'
    204       ENDIF
     276         enddo !Nfine
     277         TAURAY(:,N)=tausum(:)/tauwei
     278
     279      end do !L_NSPECTV
     280
    205281
    206282   end subroutine calc_rayleigh

  • trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

    r3259 r3654  
    2222                             glat_ig, gweight, pfgasref, pgasmax, pgasmin, &
    2323                             pgasref, tgasmax, tgasmin, tgasref, scalep, &
    24                              ubari, wnoi, stellarf, glat, dwnv, dwni, tauray
     24                             ubari, wnoi, stellarf, glat, dwnv, dwni
    2525      use watercommon_h, only: psat_water, epsi
    2626      use datafile_mod, only: datadir
    … …  
    11491149
    11501150            call optcv(dtauv,tauv,taucumv,plevrad,                 &
    1151                  qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
     1151                 qxvaer,qsvaer,gvaer,wbarv,cosbv,tauaero,          &
    11521152                 tmid,pmid,taugsurf,qvar,muvarrad,fracvari)
    11531153

  • trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90

    r3641 r3654  
    1919      logical,save :: strictboundcia                                     
    2020!$OMP THREADPRIVATE(strictboundcia)
     21      logical,save :: strictboundrayleigh                                   
     22!$OMP THREADPRIVATE(strictboundrayleigh)
    2123      logical,save :: callthermos
    2224      logical,save :: force_conduction

  • trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90

    r3641 r3654  
    438438     call getin_p("rayleigh",rayleigh)
    439439     if (is_master) write(*,*)trim(rname)//": rayleigh = ",rayleigh
     440     
     441     if (is_master) write(*,*) trim(rname)//&
     442       ": prohibit rayleigh calculations outside wavenumber grid?"
     443     strictboundrayleigh=.true. ! default value
     444     call getin_p("strictboundrayleigh",strictboundrayleigh)
     445     if (is_master) write(*,*) trim(rname)//&
     446       ": strictboundrayleigh = ",strictboundrayleigh
    440447
    441448     if (is_master) write(*,*) trim(rname)//&

  • trunk/LMDZ.GENERIC/libf/phystd/optcv.F90

    r3641 r3654  
    77SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
    88     QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
    9      TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR,FRACVAR)
     9     TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR,FRACVAR)
    1010
    1111  use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND
    … …  
    1515  use comcstfi_mod, only: g, r, mugaz
    1616  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec, &
    17                           generic_continuum_database
     17                          generic_continuum_database,rayleigh
    1818  use recombin_corrk_mod, only: corrk_recombin, gasv_recomb
    1919  use tpindex_mod, only: tpindex
    … …  
    7171  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
    7272  real*8  ANS, TAUGAS
    73   real*8,intent(in) :: TAURAY(L_NSPECTV)
     73  real*8  TAURAY(L_LEVELS,L_NSPECTV)
    7474  real*8  TRAY(L_LEVELS,L_NSPECTV)
    7575  real*8  DPR(L_LEVELS), U(L_LEVELS)
    … …  
    171171  end do
    172172 
    173   ! Rayleigh scattering
     173!=======================================================================
     174!     Set up the wavelength independent part of the Rayleigh scattering.
     175!     WAVEV is in microns.  There is no Rayleigh scattering in the IR.
     176
     177      if(rayleigh) then
     178         call calc_rayleigh(QVAR,MUVAR,PMID,TMID,TAURAY)
     179      else
     180         print*,'setspv: No Rayleigh scattering, check for NaN in output!'
     181         do NW=1,L_NSPECTV
     182            TAURAY(:,NW) = 1E-16
     183         end do
     184      endif
     185 
     186  ! Computation of pressure dependant part of Rayleigh scattering
    174187  do NW=1,L_NSPECTV
    175188     TRAY(1:4,NW)   = 1d-30
    176189     do K=5,L_LEVELS
    177         TRAY(K,NW)   = TAURAY(NW) * DPR(K)
     190        TRAY(K,NW)   = TAURAY(K,NW) * DPR(K)
    178191     end do                    ! levels
    179192  end do

  • trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90

    r2972 r3654  
    3434!                  each spectral interval.  Values are for 1 AU,
    3535!                  scaled to the planetary distance elsewhere.
    36 !     TAURAY     - Array (NSPECTV elements) of the pressure-independent
    37 !                  part of Rayleigh scattering optical depth.
    38 !     TAURAYVAR  - Array (NSPECTV elements) of the pressure-independent
    39 !                  part of Rayleigh scattering optical depth for the variable gas.
    4036!     FZEROI     - Fraction of zeros in the IR CO2 k-coefficients, for
    4137!                  each temperature, pressure, and spectral interval
    … …  
    6965      REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi
    7066      REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV) !BWNV read by master in setspv
    71       REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV), TAURAYVAR(L_NSPECTV)
     67      REAL*8 STELLARF(L_NSPECTV)
    7268!$OMP THREADPRIVATE(WNOI,DWNI,WAVEI,&
    7369        !$OMP WNOV,DWNV,WAVEV,&
    74         !$OMP STELLARF,TAURAY,TAURAYVAR)
     70        !$OMP STELLARF
    7571
    7672      REAL*8 blami(L_NSPECTI+1)

  • trunk/LMDZ.GENERIC/libf/phystd/setspv.F90

    r2831 r3654  
    55!     Purpose
    66!     -------
    7 !     Set up spectral intervals, stellar spectrum and Rayleigh
    8 !     opacity in the shortwave.
     7!     Set up spectral intervals and stellar spectrum in the shortwave.
    98!     
    109!     Authors
    … …  
    2524      use radinc_h,    only: L_NSPECTV, corrkdir, banddir
    2625      use radcommon_h, only: BWNV,BLAMV,WNOV,DWNV,WAVEV, &
    27                              STELLARF,TAURAY
     26                             STELLARF
    2827      use datafile_mod, only: datadir
    29       use callkeys_mod, only: Fat1AU,rayleigh
     28      use callkeys_mod, only: Fat1AU
    3029
    3130      implicit none
    … …  
    133132      print*,' '
    134133
    135 
    136 !=======================================================================
    137 !     Set up the wavelength independent part of the Rayleigh scattering.
    138 !     The pressure dependent part will be computed elsewhere (OPTCV).
    139 !     WAVEV is in microns.  There is no Rayleigh scattering in the IR.
    140 
    141       if(rayleigh) then
    142          call calc_rayleigh
    143       else
    144          print*,'setspv: No Rayleigh scattering, check for NaN in output!'
    145          do N=1,L_NSPECTV
    146             TAURAY(N) = 1E-16
    147          end do
    148       endif
    149 
    150134      RETURN
    151135    END subroutine setspv
Note: See TracChangeset for help on using the changeset viewer.