subroutine calc_rayleigh !================================================================== ! ! Purpose ! ------- ! Average the Rayleigh scattering in each band, weighting the ! average by the blackbody function at temperature tstellar. ! Works for an arbitrary mix of gases (currently CO2, N2 and ! H2 are possible). ! ! Authors ! ------- ! 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 ! --------- ! setspv.F ! ! Calls ! ----- ! none ! !================================================================== use radinc_h, only: L_NSPECTV 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, igas_CH4, igas_CO use comcstfi_mod, only: g, mugaz, pi implicit none real*8 wl integer N,Nfine,ifine,igas parameter(Nfine=500.0) real*8 :: P0 = 9.423D+6 ! Rayleigh scattering reference pressure in pascals. P0 = 93*101325 Pa logical typeknown real*8 tauvar,tauvarvar,tausum,tausumvar,tauwei,tauweivar,bwidth,bstart double precision df real*8 tauconsti(ngasmx) real*8 tauvari(ngasmx) integer icantbewrong ! 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. do igas=1,ngasmx 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 print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// & 'as its mixing ratio is less than 0.05.' ! ignore variable gas in Rayleigh calculation ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation tauconsti(igas) = 0.0 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*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) ! 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!' print*,"igas=",igas," gnom(igas)=",trim(gnom(igas)) call abort_physic("calc_rayleigh","Invalid gas",1) endif if((gfrac(igas).eq.1.0).and.(vgas.eq.0))then print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.' typeknown=.true. endif endif enddo if(.not.typeknown)then print*,'Rayleigh scattering is for a mixed gas atmosphere.' typeknown=.true. endif do N=1,L_NSPECTV tausum = 0.0 tauwei = 0.0 tausumvar = 0.0 tauweivar = 0.0 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 wl=bstart+dble(ifine)*bwidth/Nfine tauvar=0.0 tauvarvar=0.0 do igas=1,ngasmx if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then ! ignore variable gas in Rayleigh calculation ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation tauvari(igas) = 0.0 else if(igas.eq.igas_CO2)then tauvari(igas) = (1.0+0.013/wl**2)/wl**4 elseif(igas.eq.igas_N2)then 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) = (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!' call abort_physic("calc_rayleigh","Gas species not recognised",1) endif endif 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 tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas) endif enddo call blackl(dble(wl*1e-6),dble(tstellar),df) df=df*bwidth/Nfine tauwei=tauwei+df tausum=tausum+tauvar*df tauweivar=tauweivar+df tausumvar=tausumvar+tauvarvar*df enddo TAURAY(N)=tausum/tauwei TAURAYVAR(N)=tausumvar/tauweivar ! we multiply by scalep here (100) because plev, which is used in optcv, ! is in units of mBar, so we need to convert. end do 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, & 'cm^2 molecule^-1' ENDIF end subroutine calc_rayleigh