subroutine calc_rayleigh

!==================================================================
!     
!     Purpose
!     -------
!     Average the Rayleigh scattering in each band, weighting the 
!     average by the blackbody function at temperature tstellar.
!     
!     Authors
!     ------- 
!     Robin Wordsworth (2010)
!     Benjamin Charnay (2010)
!
!     Called by
!     ---------
!     setspv.F
!     
!     Calls
!     -----
!     none
!     
!==================================================================

      use radinc_h, only: L_NSPECTV
      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, gastype, tauray, scalep

      implicit none

#include "comcstfi.h"

      real*8 wl
      integer N,Nfine,ifine
      parameter(Nfine=500.0)
      real*8 :: P0 = 9.423D+6   ! Rayleigh scattering reference pressure in pascals.    

      logical waverage
      real*8 tauconst,tauvar,tausum,tauwei,bwidth,bstart
      double precision df

      waverage=.true. ! do we perform a blackbody weighted average by band?

      if(waverage)then
         if(gastype(1).eq.'CO2')then
            print*,'Rayleigh scattering is for a CO2 atmosphere.'
            tauconst = (8.7/g)*1.527*scalep/P0
         elseif(gastype(1).eq.'AIR')then ! AIR = Earth air
            print*,'Rayleigh scattering is for an Earth-like atmosphere.'
            tauconst = (9.81/g)*8.569E-3*scalep/(P0/93.0)
         else
            print*,'Rayleigh scattering not defined for gas type ',gastype(n)
            print*,'exiting.'
            call abort
         endif
      endif
 
      do N=1,L_NSPECTV
         
         if(waverage)then
            tausum = 0.0
            tauwei = 0.0
            bstart = 10000.0/BWNV(N+1)
            bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1))
            do ifine=1,Nfine
                wl=bstart+dble(ifine)*bwidth/Nfine

                if(gastype(1).eq.'CO2')then
                   tauvar = (1.0+0.013/wl**2)/wl**4
                elseif(gastype(1).eq.'AIR')then
                   tauvar = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
                endif

                call blackl(dble(wl*1e-6),dble(tstellar),df)
                df=df*bwidth/Nfine
                tauwei=tauwei+df
                tausum=tausum+tauvar*df
            enddo
            TAURAY(N)=tauconst*tausum/tauwei
         else
            wl = WAVEV(N)
            if(gastype(1).eq.'CO2')then
               if(N.eq.1) print*,'Rayleigh scattering is for a CO2 atmosphere.'
               TAURAY(N) = (8.7/g)*(1.527*(1.0+0.013/wl**2)/wl**4)*scalep/P0
            elseif(gastype(1).eq.'AIR')then ! AIR = Earth air
               if(N.eq.1) print*,'Rayleigh scattering is for an Earth-like atmosphere.'
               TAURAY(N) = (9.81/g)*(8.569E-3*(1.0+0.0113/wl**2+0.00013/wl**4)/wl**4)*scalep/(P0/93.0)
            else
               print*,'Rayleigh scattering not defined for gas type ',gastype(n)
               print*,'exiting.'
               call abort
            endif
         endif       

         ! 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

    end subroutine calc_rayleigh