Changeset 2889 for trunk/LMDZ.GENERIC/libf

Timestamp:

Feb 7, 2023, 4:13:20 PM (21 months ago)

Author:

jleconte

Message:

calc_rayleigh with more molecules & more comments

File:

• trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90 (modified) (12 diffs)

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

@@ -14 +14 @@
 !     Robin Wordsworth (2010)
 !     Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995).
+!     Noé Clément (2022) : Additionnal comments & Methane+CO Rayleigh
 !
 !     Called by
@@ -28 +29 @@
       use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep
       use gases_h, only: ngasmx, vgas, gnom, gfrac, igas_CO2, igas_H2, &
-                         igas_H2O, igas_He, igas_N2
-      use comcstfi_mod, only: g, mugaz
+                           igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO
+      use comcstfi_mod, only: g, mugaz, pi
 
       implicit none
@@ -36 +37 @@
       integer N,Nfine,ifine,igas
       parameter(Nfine=500.0)
-      real*8 :: P0 = 9.423D+6   ! Rayleigh scattering reference pressure in pascals.    
+      real*8 :: P0 = 9.423D+6   ! Rayleigh scattering reference pressure in pascals. P0 = 93*101325 Pa   
 
       logical typeknown
@@ -47 +48 @@
       integer icantbewrong
 
-      ! tau0/p0=tau/p (Hansen 1974)
-      ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4)
+      ! we multiply by scalep here (100) because plev, which is used in optcv,
+      ! is in units of mBar, so we need to convert.
+      ! we calculate here TAURAY which is in m2/mBar
+
+      ! tau0/p0=tau/p (Hansen 1974 : equation (2.31) )
+      ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4)      (Hansen 1974 : equation (2.29) )
       ! where delta=n-1 and N0 is an amagat
+
+      ! In exo_k : (8*pi^3)/(3*N0^2) * 4 * delta^2 ------- is written : (24*pi^3)/(N0^2) * (4/9) * delta^2
+
+      ! (tauconsti * tauvari) = scalep * cross_section_molecule_exo_k /  ( gravity * molecular_mass )
+      ! (tauconsti * tauvari) in m2/mBar because of scalep factor
+      ! cross_section_molecule in m2/molecule
+      ! molecular_mass is the mass of th considered molecule
 
       typeknown=.false.
@@ -55 +67 @@
          if(igas.eq.vgas)then
             print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering '
-         endif
-         if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
+         endif
+         if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
             print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
             'as its mixing ratio is less than 0.05.' 
@@ -64 +76 @@
          else
             if(igas.eq.igas_CO2) then
+               ! Hansen 1974 : equation (2.32)
                tauconsti(igas) = (8.7/g)*1.527*scalep/P0
             elseif(igas.eq.igas_N2)then
+               ! Hansen 1974 : equation (2.30) ---- (P0/93.0) is equal to 101325 Pa
                tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
             elseif(igas.eq.igas_H2O)then
-!               tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 
-               tauconsti(igas) = 1.98017E-10/(g*mugaz*100.) 
+               ! tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 
+               tauconsti(igas) = 1.98017E-10*scalep/(g*mugaz*1E4)
             elseif(igas.eq.igas_H2)then
                tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0
-               !tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
+               ! tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
                ! uses n=1.000132 from Optics, Fourth Edition.
                ! was out by a factor 2!
+
+               ! below is exo_k formula : this tauconsti is consistent with commented tauvari for H2 below
+               ! tauconsti(igas) = (1/(g*2.*1.6726*1E-27)) * 1E16 * scalep
             elseif(igas.eq.igas_He)then
                tauconsti(igas) = (10.0/g)*0.00086*scalep/101325.0
+            elseif(igas.eq.igas_CH4)then
+               ! For CH4 we use the formalism of exo_k (developed by Jeremy Leconte)
+               ! the relation between exo_k formalism and LMDZ formalism is : (tauconsti*tauvari)=scalep*sigma_exok/(gravity*molecular_mass)
+               ! Values are taken from Caldas (2019), equation 12 & appendix D
+               ! 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 is a constant
+               ! 4/9=(2/3)**2 is approximately ((n+1)/(n**2+2))**2
+               ! 1E24 = (1E6)**4 because wavelength is in micron
+               ! 16.04*1.6726*1E-27 is CH4 molecular mass
+               tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep
+            
+            elseif(igas.eq.igas_CO)then
+               ! see above for explanation
+               ! 28.01*1.6726*1E-27 is CH4 molecular mass
+               tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *28.01*1.6726*1e-27))  * scalep
             else
                print*,'Error in calc_rayleigh: Gas species not recognised!'
@@ -91 +122 @@
       enddo
 
+      !!!!!!! methane rayleigh verification 
+      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
+                     !* (4.6662E-4+4.02E-6/0.55**2)**2 / 0.55**4
+
+      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))  &
+               *  (4.6662E-4+4.02E-6/0.55**2)**2 / 0.55**4
+      
+      !!!!!!! methane rayleigh verification 
+
       if(.not.typeknown)then
          print*,'Rayleigh scattering is for a mixed gas atmosphere.'
@@ -96 +136 @@
       endif
 
- 
+   
       do N=1,L_NSPECTV
          
@@ -103 +143 @@
          tausumvar = 0.0
          tauweivar = 0.0
-         bstart = 10000.0/BWNV(N+1)
+         bstart = 10000.0/BWNV(N+1) ! BWNV is in cm-1 so 10000.0/BWNV is in micron
          bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1))
          do ifine=1,Nfine
@@ -121 +161 @@
                      tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
                   elseif(igas.eq.igas_H2O)then
-!                     tauvari(igas) = 1.0/wl**4 ! to be improved...
+                     ! tauvari(igas) = 1.0/wl**4 ! to be improved...
                      tauvari(igas) = (5.7918E6/(2.38E2-1/wl**2)+1.679E5/(57.36E0-1/wl**2))**2/wl**4 
                   elseif(igas.eq.igas_H2)then
                      tauvari(igas) = 1.0/wl**4 
+
+                     ! below is exo_k formula : this tauvari is consistent with commented tauconsti for H2 above
+                     ! tauvari(igas) = (8.14E-49+1.28E-50/wl**2+1.61E-51/wl**4) / wl**4
                   elseif(igas.eq.igas_He)then
                      tauvari(igas) = 1.0/wl**4 
+                  elseif(igas.eq.igas_CH4)then
+                     tauvari(igas) = (4.6662E-4+4.02E-6/wl**2)**2 / wl**4
+                  elseif(igas.eq.igas_CO)then
+                     tauvari(igas) = (2.2851E-4 + 0.456E4/(5.101816E9 - 1E8/wl**2))**2 / wl**4
                   else
                      print*,'Error in calc_rayleigh: Gas species not recognised!'
@@ -135 +182 @@
                if(igas.eq.vgas) then
                   tauvarvar=tauvarvar+tauconsti(igas)*tauvari(igas)
+                  tauvar=tauvar+0.0*0.0*gfrac(igas)
+               elseif(igas.eq.igas_CH4) then !!!!!!! methane rayleigh verification 
+                  tauvari(igas_CH4) = (4.6662E-4+4.02E-6/wl**2)**2 / wl**4
+                  tauvarvar=tauvarvar+tauconsti(igas_CH4)*tauvari(igas_CH4)*gfrac(igas_CH4)
                   tauvar=tauvar+0.0*0.0*gfrac(igas)
                else
@@ -158 +209 @@
 
       IF (L_NSPECTV > 6) THEN
-        icantbewrong = L_NSPECTV-6
-        print*,'At 1 atm and lambda = ',WAVEV(icantbewrong),' um'
-        print*,'tau_R = ',TAURAY(icantbewrong)*1013.25
-        print*,'sig_R = ',TAURAY(icantbewrong)*g*mugaz*1.67e-27*100, &
+         icantbewrong = L_NSPECTV-6
+         print*,'At 1 atm and lambda = ',WAVEV(icantbewrong),' um'
+         print*,'tau_R = ',TAURAY(icantbewrong)*1013.25
+         print*,'sig_R = ',TAURAY(icantbewrong)*g*mugaz*1.67e-27*100, &
                'cm^2 molecule^-1'
+         !!!!!!! methane rayleigh verification 
+         print*,'methane rayleigh verification : tau_R_CH4 = ',TAURAYVAR(icantbewrong)*1013.25
+         !!!!!!! methane rayleigh verification 
       ENDIF
 
-    end subroutine calc_rayleigh
+   end subroutine calc_rayleigh