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 N2, H2
!     and CH4 are possible).
!     
!     Authors
!     ------- 
!     Robin Wordsworth (2010)
!     Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995).
!     B. de Batz de Trenquelleon (2023) : Added CH4 rayleigh.
!
!     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
      use comcstfi_mod, only: g, mugaz, pi

      implicit none

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

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

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

      integer icantbewrong

       ! we multiply by scalep here (100) because plev, which is used in optcv,
       ! is in units of mBar, so we need to convert.
       ! we calculate here TAURAY which is in m2/mBar

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

       ! In exo_k : (8*pi^3)/(3*N0^2) * 4 * delta^2 ------- is written : (24*pi^3)/(N0^2) * (4/9) * delta^2

       ! (tauconsti * tauvari) = scalep * cross_section_molecule_exo_k /  ( gravity * molecular_mass )
       ! (tauconsti * tauvari) in m2/mBar because of scalep factor
       ! cross_section_molecule in m2/molecule
       ! molecular_mass is the mass of th considered molecule

      typeknown=.false.
      do igas=1,ngasmx
         if(gfrac(igas,nivref).lt.1.e-2)then
            print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
               'as its mixing ratio is less than 0.01.' 
            ! ignore variable gas in Rayleigh calculation
            ! ignore gases of mixing ratio < 0.01 in Rayleigh calculation
            tauconsti(igas) = 0.0
         else
            if(igas.eq.igas_N2)then
               tauconsti(igas) = (9.81/g)*8.569E-3*scalep/(P0/93.0)
            elseif(igas.eq.igas_H2)then
               tauconsti(igas) = (10.0/g)*0.0241*scalep/101325.0
               !tauconsti(igas) = (10.0/g)*0.0487*scalep/(101325.0)
               ! uses n=1.000132 from Optics, Fourth Edition.
               ! was out by a factor 2!
            elseif(igas.eq.igas_CH4)then
               ! For CH4 we use the formalism of exo_k (developed by Jeremy Leconte)
               ! the relation between exo_k formalism and LMDZ formalism is : (tauconsti*tauvari)=scalep*sigma_exok/(gravity*molecular_mass)
               ! Values are taken from Caldas (2019), equation 12 & appendix D
               ! 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 is a constant
               ! 4/9=(2/3)**2 is approximately ((n+1)/(n**2+2))**2
               ! 1E24 = (1E6)**4 because wavelength is in micron
               ! 16.04*1.6726*1E-27 is CH4 molecular mass
               tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep             
            else
               print*,'Error in calc_rayleigh: Gas species not recognised!'
               call abort
            endif

            if(gfrac(igas,nivref).eq.1.0)then
               print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.'
               typeknown=.true.
            endif

         endif
      enddo

      if(.not.typeknown)then
         print*,'Rayleigh scattering is for a mixed gas atmosphere.'
         typeknown=.true.
      endif

 
      do N=1,L_NSPECTV
         
         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

            tauvar=0.0
            do igas=1,ngasmx
               if(gfrac(igas,nivref).lt.1.e-2)then
                  ! ignore variable gas in Rayleigh calculation
                  ! ignore gases of mixing ratio < 0.01 in Rayleigh calculation
                  tauvari(igas)   = 0.0
               else
                  if(igas.eq.igas_N2)then
                     tauvari(igas) = (1.0+0.0113/wl**2+0.00013/wl**4)/wl**4
                  elseif(igas.eq.igas_H2)then
                     tauvari(igas) = 1.0/wl**4 
		  elseif(igas.eq.igas_CH4)then
                      tauvari(igas) = (4.6662E-4+4.02E-6/wl**2)**2 / wl**4
                  else
                     print*,'Error in calc_rayleigh: Gas species not recognised!'
                     call abort
                  endif
               endif

               tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas,nivref)

            enddo

            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)*g*mugaz*1.67e-27*100, &
               'cm^2 molecule^-1'
      ENDIF

    end subroutine calc_rayleigh