| 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 | ! |
|---|
| 16 | ! Called by |
|---|
| 17 | ! --------- |
|---|
| 18 | ! setspv.F |
|---|
| 19 | ! |
|---|
| 20 | ! Calls |
|---|
| 21 | ! ----- |
|---|
| 22 | ! none |
|---|
| 23 | ! |
|---|
| 24 | !================================================================== |
|---|
| 25 | |
|---|
| 26 | use radinc_h, only: L_NSPECTV |
|---|
| 27 | use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, scalep |
|---|
| 28 | |
|---|
| 29 | implicit none |
|---|
| 30 | |
|---|
| 31 | #include "comcstfi.h" |
|---|
| 32 | #include "callkeys.h" |
|---|
| 33 | #include "gases.h" |
|---|
| 34 | |
|---|
| 35 | real*8 wl |
|---|
| 36 | integer N,Nfine,ifine,igas |
|---|
| 37 | parameter(Nfine=500.0) |
|---|
| 38 | real*8 :: P0 = 9.423D+6 ! Rayleigh scattering reference pressure in pascals. |
|---|
| 39 | |
|---|
| 40 | logical typeknown |
|---|
| 41 | real*8 tauvar,tausum,tauwei,bwidth,bstart |
|---|
| 42 | double precision df |
|---|
| 43 | |
|---|
| 44 | real*8 tauconsti(ngasmx) |
|---|
| 45 | real*8 tauvari(ngasmx) |
|---|
| 46 | |
|---|
| 47 | ! tau0/p0=tau/p (Hansen 1974) |
|---|
| 48 | ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0) * 4*delta^2/(g*mugaz*lambda^4) |
|---|
| 49 | ! Then calculate tau0 = 1.16e-20 * 4*delta^2/(g*mugaz*lambda^4 [in um]) |
|---|
| 50 | ! where delta=n-1 |
|---|
| 51 | |
|---|
| 52 | typeknown=.false. |
|---|
| 53 | do igas=1,ngasmx |
|---|
| 54 | !if((igas.eq.vgas).or.(gfrac(igas).lt.1.e-4))then |
|---|
| 55 | if(igas.eq.vgas)then |
|---|
| 56 | print*,'Ignoring ',gnom(igas),' in Rayleigh scattering '// & |
|---|
| 57 | 'as it is variable.' |
|---|
| 58 | elseif(gfrac(igas).lt.5.e-2)then |
|---|
| 59 | print*,'Ignoring ',gnom(igas),' in Rayleigh scattering '// & |
|---|
| 60 | 'as its mixing ratio is less than 0.05.' |
|---|
| 61 | ! ignore variable gas in Rayleigh calculation |
|---|
| 62 | ! ignore gases of mixing ratio < 1e-4 in Rayleigh calculation |
|---|
| 63 | tauconsti(igas) = 0.0 |
|---|
| 64 | else |
|---|
| 65 | if(gnom(igas).eq.'CO2')then |
|---|
| 66 | tauconsti(igas) = (8.7/g)*1.527*scalep/P0 |
|---|
| 67 | elseif(gnom(igas).eq.'N2_')then |
|---|
| 68 | tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0) |
|---|
| 69 | elseif(gnom(igas).eq.'H2_')then |
|---|
| 70 | tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0) |
|---|
| 71 | ! uses n=1.000132 from Optics, Fourth Edition. |
|---|
| 72 | ! around four times more opaque than CO2 |
|---|
| 73 | ! around 5.5 times more opaque than air |
|---|
| 74 | elseif(gnom(igas).eq.'He_')then |
|---|
| 75 | print*,'Helium not ready yet!' |
|---|
| 76 | else |
|---|
| 77 | print*,'Error in calc_rayleigh: Gas species not recognised!' |
|---|
| 78 | call abort |
|---|
| 79 | endif |
|---|
| 80 | |
|---|
| 81 | if(gfrac(igas).eq.1.0)then |
|---|
| 82 | print*,'Rayleigh scattering is for a pure ',gnom(igas),' atmosphere.' |
|---|
| 83 | typeknown=.true. |
|---|
| 84 | endif |
|---|
| 85 | |
|---|
| 86 | endif |
|---|
| 87 | enddo |
|---|
| 88 | |
|---|
| 89 | if(.not.typeknown)then |
|---|
| 90 | print*,'Rayleigh scattering is for a mixed gas atmosphere.' |
|---|
| 91 | typeknown=.true. |
|---|
| 92 | endif |
|---|
| 93 | |
|---|
| 94 | |
|---|
| 95 | do N=1,L_NSPECTV |
|---|
| 96 | |
|---|
| 97 | tausum = 0.0 |
|---|
| 98 | tauwei = 0.0 |
|---|
| 99 | bstart = 10000.0/BWNV(N+1) |
|---|
| 100 | bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1)) |
|---|
| 101 | do ifine=1,Nfine |
|---|
| 102 | wl=bstart+dble(ifine)*bwidth/Nfine |
|---|
| 103 | |
|---|
| 104 | tauvar=0.0 |
|---|
| 105 | do igas=1,ngasmx |
|---|
| 106 | !if((igas.eq.vgas).or.(gfrac(igas).lt.1.e-4))then |
|---|
| 107 | if((igas.eq.vgas).or.(gfrac(igas).lt.5.e-2))then |
|---|
| 108 | ! ignore variable gas in Rayleigh calculation |
|---|
| 109 | ! ignore gases of mixing ratio < 1e-4 in Rayleigh calculation |
|---|
| 110 | tauvari(igas) = 0.0 |
|---|
| 111 | else |
|---|
| 112 | if(gnom(igas).eq.'CO2')then |
|---|
| 113 | tauvari(igas) = (1.0+0.013/wl**2)/wl**4 |
|---|
| 114 | elseif(gnom(igas).eq.'N2_')then |
|---|
| 115 | tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4 |
|---|
| 116 | elseif(gnom(igas).eq.'H2_')then |
|---|
| 117 | tauvari(igas) = 1.0/wl**4 ! no wvl dependence of ref. index here - improve? |
|---|
| 118 | else |
|---|
| 119 | print*,'Error in calc_rayleigh: Gas species not recognised!' |
|---|
| 120 | call abort |
|---|
| 121 | endif |
|---|
| 122 | endif |
|---|
| 123 | |
|---|
| 124 | tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas) |
|---|
| 125 | |
|---|
| 126 | enddo |
|---|
| 127 | |
|---|
| 128 | call blackl(dble(wl*1e-6),dble(tstellar),df) |
|---|
| 129 | df=df*bwidth/Nfine |
|---|
| 130 | tauwei=tauwei+df |
|---|
| 131 | tausum=tausum+tauvar*df |
|---|
| 132 | |
|---|
| 133 | enddo |
|---|
| 134 | TAURAY(N)=tausum/tauwei |
|---|
| 135 | ! we multiply by scalep here (100) because plev, which is used in optcv, |
|---|
| 136 | ! is in units of mBar, so we need to convert. |
|---|
| 137 | |
|---|
| 138 | end do |
|---|
| 139 | |
|---|
| 140 | print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-5), & |
|---|
| 141 | ' um, tau_R = ',TAURAY(L_NSPECTV-5)*1013.25 |
|---|
| 142 | print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-6), & |
|---|
| 143 | ' um, tau_R = ',TAURAY(L_NSPECTV-6)*1013.25 |
|---|
| 144 | |
|---|
| 145 | end subroutine calc_rayleigh |
|---|
