Changeset 2889 for trunk/LMDZ.GENERIC


Ignore:
Timestamp:
Feb 7, 2023, 4:13:20 PM (23 months ago)
Author:
jleconte
Message:

calc_rayleigh with more molecules & more comments

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90

    r2529 r2889  
    1414!     Robin Wordsworth (2010)
    1515!     Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995).
     16!     Noé Clément (2022) : Additionnal comments & Methane+CO Rayleigh
    1617!
    1718!     Called by
     
    2829      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep
    2930      use gases_h, only: ngasmx, vgas, gnom, gfrac, igas_CO2, igas_H2, &
    30                          igas_H2O, igas_He, igas_N2
    31       use comcstfi_mod, only: g, mugaz
     31                           igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO
     32      use comcstfi_mod, only: g, mugaz, pi
    3233
    3334      implicit none
     
    3637      integer N,Nfine,ifine,igas
    3738      parameter(Nfine=500.0)
    38       real*8 :: P0 = 9.423D+6   ! Rayleigh scattering reference pressure in pascals.    
     39      real*8 :: P0 = 9.423D+6   ! Rayleigh scattering reference pressure in pascals. P0 = 93*101325 Pa   
    3940
    4041      logical typeknown
     
    4748      integer icantbewrong
    4849
    49       ! tau0/p0=tau/p (Hansen 1974)
    50       ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4)
     50      ! we multiply by scalep here (100) because plev, which is used in optcv,
     51      ! is in units of mBar, so we need to convert.
     52      ! we calculate here TAURAY which is in m2/mBar
     53
     54      ! tau0/p0=tau/p (Hansen 1974 : equation (2.31) )
     55      ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4)      (Hansen 1974 : equation (2.29) )
    5156      ! where delta=n-1 and N0 is an amagat
     57
     58      ! In exo_k : (8*pi^3)/(3*N0^2) * 4 * delta^2 ------- is written : (24*pi^3)/(N0^2) * (4/9) * delta^2
     59
     60      ! (tauconsti * tauvari) = scalep * cross_section_molecule_exo_k /  ( gravity * molecular_mass )
     61      ! (tauconsti * tauvari) in m2/mBar because of scalep factor
     62      ! cross_section_molecule in m2/molecule
     63      ! molecular_mass is the mass of th considered molecule
    5264
    5365      typeknown=.false.
     
    5567         if(igas.eq.vgas)then
    5668            print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering '
    57         endif
    58         if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
     69        endif
     70        if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
    5971            print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
    6072            'as its mixing ratio is less than 0.05.'
     
    6476         else
    6577            if(igas.eq.igas_CO2) then
     78               ! Hansen 1974 : equation (2.32)
    6679               tauconsti(igas) = (8.7/g)*1.527*scalep/P0
    6780            elseif(igas.eq.igas_N2)then
     81               ! Hansen 1974 : equation (2.30) ---- (P0/93.0) is equal to 101325 Pa
    6882               tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
    6983            elseif(igas.eq.igas_H2O)then
    70 !              tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0
    71                tauconsti(igas) = 1.98017E-10/(g*mugaz*100.)
     84               ! tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0
     85               tauconsti(igas) = 1.98017E-10*scalep/(g*mugaz*1E4)
    7286            elseif(igas.eq.igas_H2)then
    7387               tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0
    74                !tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
     88               ! tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
    7589               ! uses n=1.000132 from Optics, Fourth Edition.
    7690               ! was out by a factor 2!
     91
     92               ! below is exo_k formula : this tauconsti is consistent with commented tauvari for H2 below
     93               ! tauconsti(igas) = (1/(g*2.*1.6726*1E-27)) * 1E16 * scalep
    7794            elseif(igas.eq.igas_He)then
    7895               tauconsti(igas) = (10.0/g)*0.00086*scalep/101325.0
     96            elseif(igas.eq.igas_CH4)then
     97               ! For CH4 we use the formalism of exo_k (developed by Jeremy Leconte)
     98               ! the relation between exo_k formalism and LMDZ formalism is : (tauconsti*tauvari)=scalep*sigma_exok/(gravity*molecular_mass)
     99               ! Values are taken from Caldas (2019), equation 12 & appendix D
     100               ! 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 is a constant
     101               ! 4/9=(2/3)**2 is approximately ((n+1)/(n**2+2))**2
     102               ! 1E24 = (1E6)**4 because wavelength is in micron
     103               ! 16.04*1.6726*1E-27 is CH4 molecular mass
     104               tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep
     105           
     106            elseif(igas.eq.igas_CO)then
     107               ! see above for explanation
     108               ! 28.01*1.6726*1E-27 is CH4 molecular mass
     109               tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *28.01*1.6726*1e-27))  * scalep
    79110            else
    80111               print*,'Error in calc_rayleigh: Gas species not recognised!'
     
    91122      enddo
    92123
     124      !!!!!!! methane rayleigh verification
     125      tauconsti(igas_CH4) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep
     126                     !* (4.6662E-4+4.02E-6/0.55**2)**2 / 0.55**4
     127
     128      write(*,*) 'methane rayleigh verification', 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  &
     129               *  (4.6662E-4+4.02E-6/0.55**2)**2 / 0.55**4
     130     
     131      !!!!!!! methane rayleigh verification
     132
    93133      if(.not.typeknown)then
    94134         print*,'Rayleigh scattering is for a mixed gas atmosphere.'
     
    96136      endif
    97137
    98  
     138  
    99139      do N=1,L_NSPECTV
    100140         
     
    103143         tausumvar = 0.0
    104144         tauweivar = 0.0
    105          bstart = 10000.0/BWNV(N+1)
     145         bstart = 10000.0/BWNV(N+1) ! BWNV is in cm-1 so 10000.0/BWNV is in micron
    106146         bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1))
    107147         do ifine=1,Nfine
     
    121161                     tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
    122162                  elseif(igas.eq.igas_H2O)then
    123 !                    tauvari(igas) = 1.0/wl**4 ! to be improved...
     163                     ! tauvari(igas) = 1.0/wl**4 ! to be improved...
    124164                     tauvari(igas) = (5.7918E6/(2.38E2-1/wl**2)+1.679E5/(57.36E0-1/wl**2))**2/wl**4
    125165                  elseif(igas.eq.igas_H2)then
    126166                     tauvari(igas) = 1.0/wl**4
     167
     168                     ! below is exo_k formula : this tauvari is consistent with commented tauconsti for H2 above
     169                     ! tauvari(igas) = (8.14E-49+1.28E-50/wl**2+1.61E-51/wl**4) / wl**4
    127170                  elseif(igas.eq.igas_He)then
    128171                     tauvari(igas) = 1.0/wl**4
     172                  elseif(igas.eq.igas_CH4)then
     173                     tauvari(igas) = (4.6662E-4+4.02E-6/wl**2)**2 / wl**4
     174                  elseif(igas.eq.igas_CO)then
     175                     tauvari(igas) = (2.2851E-4 + 0.456E4/(5.101816E9 - 1E8/wl**2))**2 / wl**4
    129176                  else
    130177                     print*,'Error in calc_rayleigh: Gas species not recognised!'
     
    135182               if(igas.eq.vgas) then
    136183                  tauvarvar=tauvarvar+tauconsti(igas)*tauvari(igas)
     184                  tauvar=tauvar+0.0*0.0*gfrac(igas)
     185               elseif(igas.eq.igas_CH4) then !!!!!!! methane rayleigh verification
     186                  tauvari(igas_CH4) = (4.6662E-4+4.02E-6/wl**2)**2 / wl**4
     187                  tauvarvar=tauvarvar+tauconsti(igas_CH4)*tauvari(igas_CH4)*gfrac(igas_CH4)
    137188                  tauvar=tauvar+0.0*0.0*gfrac(igas)
    138189               else
     
    158209
    159210      IF (L_NSPECTV > 6) THEN
    160         icantbewrong = L_NSPECTV-6
    161         print*,'At 1 atm and lambda = ',WAVEV(icantbewrong),' um'
    162         print*,'tau_R = ',TAURAY(icantbewrong)*1013.25
    163         print*,'sig_R = ',TAURAY(icantbewrong)*g*mugaz*1.67e-27*100, &
     211         icantbewrong = L_NSPECTV-6
     212         print*,'At 1 atm and lambda = ',WAVEV(icantbewrong),' um'
     213         print*,'tau_R = ',TAURAY(icantbewrong)*1013.25
     214         print*,'sig_R = ',TAURAY(icantbewrong)*g*mugaz*1.67e-27*100, &
    164215               'cm^2 molecule^-1'
     216         !!!!!!! methane rayleigh verification
     217         print*,'methane rayleigh verification : tau_R_CH4 = ',TAURAYVAR(icantbewrong)*1013.25
     218         !!!!!!! methane rayleigh verification
    165219      ENDIF
    166220
    167     end subroutine calc_rayleigh
     221   end subroutine calc_rayleigh
Note: See TracChangeset for help on using the changeset viewer.