| 1 | subroutine calc_rayleigh |
|---|
| 2 | |
|---|
| 3 | !================================================================== |
|---|
| 4 | ! |
|---|
| 5 | ! Purpose |
|---|
| 6 | ! ------- |
|---|
| 7 | ! Average the Rayleigh scattering in each band, weighting the |
|---|
| 8 | ! average by the blackbody function at temperature tstellar. |
|---|
| 9 | ! Works for an arbitrary mix of gases (currently CO2, N2 and |
|---|
| 10 | ! H2 are possible). |
|---|
| 11 | ! |
|---|
| 12 | ! Authors |
|---|
| 13 | ! ------- |
|---|
| 14 | ! Robin Wordsworth (2010) |
|---|
| 15 | ! Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995). |
|---|
| 16 | ! Noé Clément (2022) : Additionnal comments & Methane+CO Rayleigh |
|---|
| 17 | ! |
|---|
| 18 | ! Called by |
|---|
| 19 | ! --------- |
|---|
| 20 | ! setspv.F |
|---|
| 21 | ! |
|---|
| 22 | ! Calls |
|---|
| 23 | ! ----- |
|---|
| 24 | ! none |
|---|
| 25 | ! |
|---|
| 26 | !================================================================== |
|---|
| 27 | |
|---|
| 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 | use comcstfi_mod, only: g, mugaz, pi |
|---|
| 33 | |
|---|
| 34 | implicit none |
|---|
| 35 | |
|---|
| 36 | real*8 wl |
|---|
| 37 | integer N,Nfine,ifine,igas |
|---|
| 38 | parameter(Nfine=500.0) |
|---|
| 39 | real*8 :: P0 = 9.423D+6 ! Rayleigh scattering reference pressure in pascals. P0 = 93*101325 Pa |
|---|
| 40 | |
|---|
| 41 | logical typeknown |
|---|
| 42 | real*8 tauvar,tauvarvar,tausum,tausumvar,tauwei,tauweivar,bwidth,bstart |
|---|
| 43 | double precision df |
|---|
| 44 | |
|---|
| 45 | real*8 tauconsti(ngasmx) |
|---|
| 46 | real*8 tauvari(ngasmx) |
|---|
| 47 | |
|---|
| 48 | integer icantbewrong |
|---|
| 49 | |
|---|
| 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) ) |
|---|
| 56 | ! 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 |
|---|
| 64 | |
|---|
| 65 | typeknown=.false. |
|---|
| 66 | do igas=1,ngasmx |
|---|
| 67 | if(igas.eq.vgas)then |
|---|
| 68 | print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering ' |
|---|
| 69 | endif |
|---|
| 70 | if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then |
|---|
| 71 | print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// & |
|---|
| 72 | 'as its mixing ratio is less than 0.05.' |
|---|
| 73 | ! ignore variable gas in Rayleigh calculation |
|---|
| 74 | ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation |
|---|
| 75 | tauconsti(igas) = 0.0 |
|---|
| 76 | else |
|---|
| 77 | if(igas.eq.igas_CO2) then |
|---|
| 78 | ! Hansen 1974 : equation (2.32) |
|---|
| 79 | tauconsti(igas) = (8.7/g)*1.527*scalep/P0 |
|---|
| 80 | elseif(igas.eq.igas_N2)then |
|---|
| 81 | ! Hansen 1974 : equation (2.30) ---- (P0/93.0) is equal to 101325 Pa |
|---|
| 82 | tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0) |
|---|
| 83 | elseif(igas.eq.igas_H2O)then |
|---|
| 84 | ! tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0 |
|---|
| 85 | tauconsti(igas) = 1.98017E-10*scalep/(g*mugaz*1E4) |
|---|
| 86 | elseif(igas.eq.igas_H2)then |
|---|
| 87 | tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0 |
|---|
| 88 | ! tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0) |
|---|
| 89 | ! uses n=1.000132 from Optics, Fourth Edition. |
|---|
| 90 | ! 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 |
|---|
| 94 | elseif(igas.eq.igas_He)then |
|---|
| 95 | 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 |
|---|
| 110 | else |
|---|
| 111 | print*,'Error in calc_rayleigh: Gas species not recognised!' |
|---|
| 112 | print*,"igas=",igas," gnom(igas)=",trim(gnom(igas)) |
|---|
| 113 | call abort_physic("calc_rayleigh","Invalid gas",1) |
|---|
| 114 | endif |
|---|
| 115 | |
|---|
| 116 | if((gfrac(igas).eq.1.0).and.(vgas.eq.0))then |
|---|
| 117 | print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.' |
|---|
| 118 | typeknown=.true. |
|---|
| 119 | endif |
|---|
| 120 | |
|---|
| 121 | endif |
|---|
| 122 | enddo |
|---|
| 123 | |
|---|
| 124 | if(.not.typeknown)then |
|---|
| 125 | print*,'Rayleigh scattering is for a mixed gas atmosphere.' |
|---|
| 126 | typeknown=.true. |
|---|
| 127 | endif |
|---|
| 128 | |
|---|
| 129 | |
|---|
| 130 | do N=1,L_NSPECTV |
|---|
| 131 | |
|---|
| 132 | tausum = 0.0 |
|---|
| 133 | tauwei = 0.0 |
|---|
| 134 | tausumvar = 0.0 |
|---|
| 135 | tauweivar = 0.0 |
|---|
| 136 | bstart = 10000.0/BWNV(N+1) ! BWNV is in cm-1 so 10000.0/BWNV is in micron |
|---|
| 137 | bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1)) |
|---|
| 138 | do ifine=1,Nfine |
|---|
| 139 | wl=bstart+dble(ifine)*bwidth/Nfine |
|---|
| 140 | |
|---|
| 141 | tauvar=0.0 |
|---|
| 142 | tauvarvar=0.0 |
|---|
| 143 | do igas=1,ngasmx |
|---|
| 144 | if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then |
|---|
| 145 | ! ignore variable gas in Rayleigh calculation |
|---|
| 146 | ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation |
|---|
| 147 | tauvari(igas) = 0.0 |
|---|
| 148 | else |
|---|
| 149 | if(igas.eq.igas_CO2)then |
|---|
| 150 | tauvari(igas) = (1.0+0.013/wl**2)/wl**4 |
|---|
| 151 | elseif(igas.eq.igas_N2)then |
|---|
| 152 | tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4 |
|---|
| 153 | elseif(igas.eq.igas_H2O)then |
|---|
| 154 | ! tauvari(igas) = 1.0/wl**4 ! to be improved... |
|---|
| 155 | tauvari(igas) = (5.7918E6/(2.38E2-1/wl**2)+1.679E5/(57.36E0-1/wl**2))**2/wl**4 |
|---|
| 156 | elseif(igas.eq.igas_H2)then |
|---|
| 157 | tauvari(igas) = 1.0/wl**4 |
|---|
| 158 | |
|---|
| 159 | ! below is exo_k formula : this tauvari is consistent with commented tauconsti for H2 above |
|---|
| 160 | ! tauvari(igas) = (8.14E-49+1.28E-50/wl**2+1.61E-51/wl**4) / wl**4 |
|---|
| 161 | elseif(igas.eq.igas_He)then |
|---|
| 162 | tauvari(igas) = 1.0/wl**4 |
|---|
| 163 | elseif(igas.eq.igas_CH4)then |
|---|
| 164 | tauvari(igas) = (4.6662E-4+4.02E-6/wl**2)**2 / wl**4 |
|---|
| 165 | elseif(igas.eq.igas_CO)then |
|---|
| 166 | tauvari(igas) = (2.2851E-4 + 0.456E4/(5.101816E9 - 1E8/wl**2))**2 / wl**4 |
|---|
| 167 | else |
|---|
| 168 | print*,'Error in calc_rayleigh: Gas species not recognised!' |
|---|
| 169 | call abort_physic("calc_rayleigh","Gas species not recognised",1) |
|---|
| 170 | endif |
|---|
| 171 | endif |
|---|
| 172 | |
|---|
| 173 | if(igas.eq.vgas) then |
|---|
| 174 | tauvarvar=tauvarvar+tauconsti(igas)*tauvari(igas) |
|---|
| 175 | tauvar=tauvar+0.0*0.0*gfrac(igas) |
|---|
| 176 | else |
|---|
| 177 | tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas) |
|---|
| 178 | endif |
|---|
| 179 | |
|---|
| 180 | enddo |
|---|
| 181 | |
|---|
| 182 | call blackl(dble(wl*1e-6),dble(tstellar),df) |
|---|
| 183 | df=df*bwidth/Nfine |
|---|
| 184 | tauwei=tauwei+df |
|---|
| 185 | tausum=tausum+tauvar*df |
|---|
| 186 | tauweivar=tauweivar+df |
|---|
| 187 | tausumvar=tausumvar+tauvarvar*df |
|---|
| 188 | |
|---|
| 189 | enddo |
|---|
| 190 | TAURAY(N)=tausum/tauwei |
|---|
| 191 | TAURAYVAR(N)=tausumvar/tauweivar |
|---|
| 192 | ! we multiply by scalep here (100) because plev, which is used in optcv, |
|---|
| 193 | ! is in units of mBar, so we need to convert. |
|---|
| 194 | |
|---|
| 195 | end do |
|---|
| 196 | |
|---|
| 197 | IF (L_NSPECTV > 6) THEN |
|---|
| 198 | icantbewrong = L_NSPECTV-6 |
|---|
| 199 | print*,'At 1 atm and lambda = ',WAVEV(icantbewrong),' um' |
|---|
| 200 | print*,'tau_R = ',TAURAY(icantbewrong)*1013.25 |
|---|
| 201 | print*,'sig_R = ',TAURAY(icantbewrong)*g*mugaz*1.67e-27*100, & |
|---|
| 202 | 'cm^2 molecule^-1' |
|---|
| 203 | ENDIF |
|---|
| 204 | |
|---|
| 205 | end subroutine calc_rayleigh |
|---|
