source: trunk/LMDZ.GENERIC/libf/phystd/calc_rayleigh.F90 @ 735

Last change on this file since 735 was 716, checked in by rwordsworth, 13 years ago

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

Note to Martians: I've changed the default values of 'sedimentation' and
'co2cond' in inifis to false. Both these are defined in the standard deftank
callphys.def file, so there shouldn't be any compatibility problems.

File size: 5.0 KB
RevLine 
[135]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.
[253]9!     Works for an arbitrary mix of gases (currently CO2, N2 and
10!     H2 are possible).
[135]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
[253]27      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, scalep
[471]28      use gases_h
[135]29
30      implicit none
31
32#include "comcstfi.h"
[253]33#include "callkeys.h"
[135]34
35      real*8 wl
[253]36      integer N,Nfine,ifine,igas
[135]37      parameter(Nfine=500.0)
38      real*8 :: P0 = 9.423D+6   ! Rayleigh scattering reference pressure in pascals.   
39
[253]40      logical typeknown
41      real*8 tauvar,tausum,tauwei,bwidth,bstart
[135]42      double precision df
43
[253]44      real*8 tauconsti(ngasmx)
45      real*8 tauvari(ngasmx)
[135]46
[253]47      ! tau0/p0=tau/p (Hansen 1974)
[305]48      ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4)
49      ! where delta=n-1 and N0 is an amagat
[253]50
51      typeknown=.false.
52      do igas=1,ngasmx
53         if(igas.eq.vgas)then
54            print*,'Ignoring ',gnom(igas),' in Rayleigh scattering '// &
55            'as it is variable.'
56         elseif(gfrac(igas).lt.5.e-2)then
57            print*,'Ignoring ',gnom(igas),' in Rayleigh scattering '// &
58            'as its mixing ratio is less than 0.05.'
59            ! ignore variable gas in Rayleigh calculation
[305]60            ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
[253]61            tauconsti(igas) = 0.0
[135]62         else
[253]63            if(gnom(igas).eq.'CO2')then
64               tauconsti(igas) = (8.7/g)*1.527*scalep/P0
65            elseif(gnom(igas).eq.'N2_')then
66               tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
[305]67            elseif(gnom(igas).eq.'H2O')then
68               tauconsti(igas) = (10.0/g)*9.22E-3*scalep/101325.0
[253]69            elseif(gnom(igas).eq.'H2_')then
[305]70               tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0
71               !tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
[253]72               ! uses n=1.000132 from Optics, Fourth Edition.
[305]73               ! was out by a factor 2!
[253]74            elseif(gnom(igas).eq.'He_')then
75               print*,'Helium not ready yet!'
[716]76               tauconsti(igas) = 0.0
77               call abort
[253]78            else
79               print*,'Error in calc_rayleigh: Gas species not recognised!'
80               call abort
81            endif
82
83            if(gfrac(igas).eq.1.0)then
84               print*,'Rayleigh scattering is for a pure ',gnom(igas),' atmosphere.'
85               typeknown=.true.
86            endif
87
[135]88         endif
[253]89      enddo
90
91      if(.not.typeknown)then
92         print*,'Rayleigh scattering is for a mixed gas atmosphere.'
93         typeknown=.true.
[135]94      endif
[253]95
[135]96 
97      do N=1,L_NSPECTV
98         
[253]99         tausum = 0.0
100         tauwei = 0.0
101         bstart = 10000.0/BWNV(N+1)
102         bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1))
103         do ifine=1,Nfine
104            wl=bstart+dble(ifine)*bwidth/Nfine
[135]105
[253]106            tauvar=0.0
107            do igas=1,ngasmx
108               if((igas.eq.vgas).or.(gfrac(igas).lt.5.e-2))then
109                  ! ignore variable gas in Rayleigh calculation
[305]110                  ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
[253]111                  tauvari(igas)   = 0.0
112               else
113                  if(gnom(igas).eq.'CO2')then
114                     tauvari(igas) = (1.0+0.013/wl**2)/wl**4
115                  elseif(gnom(igas).eq.'N2_')then
116                     tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
[305]117                  elseif(gnom(igas).eq.'H2O')then
118                     tauvari(igas) = 1.0/wl**4 ! to be improved...
[253]119                  elseif(gnom(igas).eq.'H2_')then
[305]120                     tauvari(igas) = 1.0/wl**4
[716]121                  elseif(gnom(igas).eq.'He_')then
122                     print*,'Helium not ready yet!'
123                     tauvari(igas) = 0.0
[253]124                  else
125                     print*,'Error in calc_rayleigh: Gas species not recognised!'
126                     call abort
127                  endif
128               endif
[135]129
[253]130               tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas)
131
[135]132            enddo
133
[253]134            call blackl(dble(wl*1e-6),dble(tstellar),df)
135            df=df*bwidth/Nfine
136            tauwei=tauwei+df
137            tausum=tausum+tauvar*df
138         
139         enddo
140         TAURAY(N)=tausum/tauwei
[135]141         ! we multiply by scalep here (100) because plev, which is used in optcv,
142         ! is in units of mBar, so we need to convert.
[253]143
[135]144      end do
145
[305]146      print*,'At 1 atm and lambda = ',WAVEV(L_NSPECTV-6),' um'
147      print*,'tau_R = ',TAURAY(L_NSPECTV-6)*1013.25
148      print*,'sig_R = ',TAURAY(L_NSPECTV-6)*g*mugaz*1.67e-27*100, &
149             'cm^2 molecule^-1'
[253]150
[135]151    end subroutine calc_rayleigh
Note: See TracBrowser for help on using the repository browser.