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).
!
!     Called by
!     ---------
!     setspv.F
!     
!     Calls
!     -----
!     none
!     
!==================================================================

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

      implicit none
#include "YOMCST.h"

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

      real*8 tauvar,tausum,tausumvar,tauwei,tauweivar,bwidth,bstart
      double precision df

      real*8 tauconsti(ngasmx)
      real*8 tauvari(ngasmx)

      integer icantbewrong

      ! tau0/p0=tau/p (Hansen 1974)
      ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(RG*RMD*lambda^4)
      ! where delta=n-1 and N0 is an amagat

!! ADAPTED FOR FIXED VENUS ATMOSPHERE: 96.5% CO2, 3.5% N2
!! using CO2 = 1 ; N2 = 2
      tauconsti(1) = (8.7/RG)*1.527*scalep/P0
      tauconsti(2) = (9.81/RG)*8.569E-3*scalep/(P0/93.0)
 
      do N=1,L_NSPECTV
         
         tausum = 0.0
         tauwei = 0.0
         tausumvar = 0.0
         tauweivar = 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

            tauvar=0.0
            tauvari(1) = (1.0+0.013/wl**2)/wl**4
            tauvari(2) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
            tauvar=tauconsti(1)*tauvari(1)*0.965 + tauconsti(2)*tauvari(2)*0.035

            call blackl(dble(wl*1e-6),dble(tstellar),df)
            df=df*bwidth/Nfine
            tauwei=tauwei+df
            tausum=tausum+tauvar*df
         
         enddo
         TAURAY(N)=tausum/tauwei
         ! 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)*RG*RMD*1.67e-27*100, &
               'cm^2 molecule^-1'
      ENDIF

    end subroutine calc_rayleigh