module calc_rayleigh_mod

implicit none

contains

      subroutine calc_rayleigh(qvar,muvar,PMID,TMID,tauray)

!==================================================================
!     
!     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.
!     
!     Authors
!     ------- 
!     Robin Wordsworth (2010)
!     Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995).
!     Noe Clement (2022) : Additionnal comments & Methane+CO Rayleigh
!     Gwenael Milcareck (2025): Rewriting the code
! 
!     Called by
!     ---------
!     setspv.F
!     
!     Calls
!     -----
!     none
!     
!==================================================================

      use radinc_h, only: L_NSPECTV, L_LEVELS
      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, scalep
      use gases_h, only: ngasmx, vgas, gnom, gfrac, massmol, igas_CO2, igas_H2, &
                           igas_H2O, igas_He, igas_N2, igas_CH4, igas_CO, igas_Ar, igas_O2
      use comcstfi_mod, only: g, pi
      use callkeys_mod, only: strictboundrayleigh

      implicit none

      real, intent(in) :: qvar(L_LEVELS) ! mol/mol
      real, intent(in) :: muvar(L_LEVELS) ! g/mol
      real, intent(in) :: PMID(L_LEVELS) ! mbar
      real, intent(in) :: TMID(L_LEVELS) ! K
      real, intent(out) :: tauray(L_LEVELS,L_NSPECTV)
      real*8 wl,wn
      integer N,Nfine,ifine,igas,k
      parameter(Nfine=500.0)
      real*8 :: Fk ! King factor for the depolarization
      real*8 :: ng(L_LEVELS) ! real refractive index 
      real*8 :: P0(L_LEVELS) ! reference pressure
      real*8 :: T0(L_LEVELS) ! reference temperature
      
      ! Parameters for H2O
      real*8 :: a0, a1, a2, a3, a4, a5, a6, a7
      real*8 :: luv, lir
      real*8 :: rhor,tr,lr
      real*8 :: rho(L_LEVELS),rhos(L_LEVELS),ts(L_LEVELS)
      real*8 :: b(L_LEVELS)

      real*8 mass_frac(ngasmx,L_LEVELS)
      real*8 tauvar(L_LEVELS),tausum(L_LEVELS)
      real*8 tauwei,bwidth,bstart
      double precision df

      real*8 tauconsti(ngasmx,L_LEVELS)
      real*8 tauvari(ngasmx,L_LEVELS)
      
      ! Miscellaneous :
      character(len=200) :: message
      character(len=10),parameter :: subname="rayleigh"
      logical, save :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)      

      integer icantbewrong

      ! This module calculates the Rayleigh scattering (also known as the Cabannes peak)
      ! Rayleigh wings, Brillouin scattering and Raman scattering are not taken into account.

      ! we calculate here TAURAY which is in m2/mBar

      ! The cross section for ith particles of small size compared to the wavenumber
      ! and in the electric dipole approximation is:
      ! sigma_i = 24*pi**3*wn**4/N**2 * ((n_i(wn)**2 - 1)/(n_i(wn)**2 + 2))**2 * Fk_i(wn)
      ! nu is the wavenumber
      ! N is the number density of the gas (molecule/m3)
      ! n_i is the real refractive index of the ith gas
      ! Fk_i is the King factor of ith gas equals to (6+3*delta_i)/(6-7*delta_i)
      ! where delta_i is the depolarization factor of the ith gas
      
      ! The rayleigh opacity is expressed by:
      ! tau_r = P/(g*mu) * sum_{i=1}^Ntot [ x_i*sigma_i ]
      ! P is the pressure
      ! g is the standard gravity
      ! mu is the mean molecular weight
      ! x_i is the mass fraction of the ith gas
      ! The pressure P dependence is calculated in optcv.F90
      
      if(firstcall) then
        
        if ((BWNV(L_NSPECTV+1).gt.60000.).and.(strictboundrayleigh)) then
          message="Rayleigh scattering is unknown for wn>60000 cm-1 - some singularities are present for higher wavenumber - if you know what you are doing, use strictboundrayleigh=.false."
          call abort_physic(subname,message,1)
        elseif ((BWNV(L_NSPECTV+1).gt.60000.).and..not.(strictboundrayleigh)) then
          print*,'**********************************************'
          print*,' we allow model to continue with wn>60000 cm-1' 
          print*,' ... we assume we know what you are doing ... '
          print*,' ... but do not let this happen too often ... '
          print*,'**********************************************'
        endif
        firstcall = .false.
      endif
      

      do igas=1,ngasmx
         ! Convert qvar mol/mol -> kg/kg
         if((igas.eq.vgas).and.(maxval(QVAR(:)).ge.1.e-2))then
           ! print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering '
            mass_frac(igas,:) = QVAR(:)*massmol(igas)/muvar(:)
         elseif((igas/=vgas).and.(gfrac(igas).ge.1.e-2))then
            mass_frac(igas,:) = gfrac(igas)*(1.-QVAR(:))*massmol(igas)/muvar(:)
         else 
           ! 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
            mass_frac(igas,:) = 0.0
         endif
         tauvari(igas,:) = 0.
      enddo
      
      ! WARNING, beyond 60000 cm-1, for all molecules, there are singularities due to the interpolation formula.
      
   
      do N=1,L_NSPECTV
      
         ! The refractive index depend on temperature and pressure
         ! It isn't the case here. Must be implemented in the future...
         ! But in the current scientific litterature (2024), it's difficult 
         ! to find something that depends on temperature and pressure...
         ! except for H2O
         
         tausum = 0.0
         tauwei = 0.0
         bstart = 10000.0/BWNV(N+1) ! BWNV is in cm-1 so 10000.0/BWNV is in micron
         bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1))
         do ifine=1,Nfine
            wl=bstart+dble(ifine)*bwidth/Nfine
            wn=BWNV(N)+dble(ifine)*(BWNV(N+1)-BWNV(N))/Nfine

            tauvar(:)=0.0
            do igas=1,ngasmx
               if (maxval(mass_frac(igas,:)).ge.1e-2) then
                 
                 if(igas.eq.igas_CO2)then
                     ! Sneep et al, 2005
                     ! doi:10.1016/j.jqsrt.2004.07.025
                     T0(:) = 288.15
                     P0(:) = 1.01325e5
                     ! ng -> valid range of the measurements : 0.1807 - 1.8172 um
                     ng(:) = 1. + 1.1427e3*(5799.25/(128908.9**2 - wn**2) + 120.05/(89223.8**2 - wn**2) + 5.3334/(75037.5**2 - wn**2) + 4.3244/(67837.7**2 - wn**2) + 0.1218145e-4/(2418.136**2 - wn**2)) ! there is an error on the paper 1.1427e6 -> 1.1427e3
                     Fk = 1.1364 + 25.3e-12*wn**2
                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4 ! wn*100 -> cm-1 to m-1
                 elseif(igas.eq.igas_N2)then
                     ! Thalman et al, 2014
                     ! doi:10.1016/j.jqsrt.2014.05.030
                     T0(:) = 288.15
                     P0(:) = 1.01325e5
                     if(wn.gt.21360)then !between 21360 and 39370 cm-1. We extrapolate above.
                       ng(:) = 1. + (5677.465 + 318.81874e13/(14.4e9 - wn**2))*1e-8  !there is an error on the paper e12 -> e13
                     else !between 4860 and 21360 cm-1. We extrapolate below.
                       ng(:) = 1. + (6498.2 + 307.4335e13/(14.4e9 - wn**2))*1e-8
                     endif
                     Fk = 1.034 + 3.17e-12*wn
                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
                 elseif(igas.eq.igas_H2O)then
                     ! Harvey et al, 1998
                     ! doi:10.1063/1.556029
                     ! doi:10.1364/AO.35.001566 for wn<4840 cm-1
                     a0 =  0.244257733
                     a1 = 9.74634476e-3
                     a2 = -3.73234996e-3
                     a3 = 2.68678472e-4
                     a4 = 1.58920570e-3
                     a5 = 2.45934259e-3
                     a6 = 0.900704920
                     a7 = -1.66626219e-2
                     luv = 0.2292020
                     lir = 5.432937
                     Tr = 273.15
                     rhor = 1000.
                     T0(:) = tmid(:)
                     P0(:) = pmid(:)*scalep
                     lr = 0.589
                     rho(:) = mass_frac(igas,:)*muvar(:)/massmol(igas)*P0(:)/(8.314463*T0(:)/(muvar(:)/1000.))
                     rhos(:) = rho(:)/rhor
                     ts(:) = T0(:)/Tr
                     
                     b(:) = (a0 + a1*rhos(:) + a2*ts(:) + a3*ts(:)*(10000./wn/lr)**2 + a4/(10000./wn/lr)**2 + a5/((10000./wn/lr)**2 - luv**2) + a6/((10000./wn/lr)**2 - lir**2) + a7*rhos(:)**2)*rhos(:)
                     ! ng -> valid range of the measurements : 0.2 - 1.1 um
                     ng(:) = sqrt(2.*b(:)+1.)/sqrt(1.-b(:))
                     if(wn<4840.) then ! necessary to prevent a singularity at 3230 cm-1
                       ng(:) = 1. + 1.022e-8*(295.235 + 2.6422*(wn*1e-4)**2 - 0.032380*(wn*1e-4)**4 + 0.004028*(wn*1e-4)**6)
                     endif
                     Fk = (6.+3.*3e-4)/(6.-7.*3e-4) ! delta=3e-4 Murphy 1977 doi:10.1063/1.434794
                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
                 elseif(igas.eq.igas_H2)then
                     ! Peck and Hung, 1977
                     ! doi:10.1364/JOSA.67.001550
                     T0(:) = 273.15
                     P0(:) = 1.01325e5
                     ! ng -> valid range of the measurements : 0.1680 - 1.6945 um
                     ng(:) = 1. + (14895.6/(180.7 - (wn*1e-4)**2) + 4903.7/(92.-(wn*1e-4)**2))*1e-6
                     Fk = (6.+3.*0.02)/(6.-7.*0.02) ! delta=0.02 Hansen 1974
                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
                 elseif(igas.eq.igas_He)then
                     ! Thalman et al, 2014
                     ! doi:10.1016/j.jqsrt.2014.05.030
                     T0(:) = 288.15
                     P0(:) = 1.01325e5
                     ! ng -> valid range of the measurements : 0.2753 - 20.5813 um
                     ng(:) = 1. + (2283. + 1.8102e13/(1.5342e10 - wn**2))*1e-8
                     Fk = 1.
                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
                 elseif(igas.eq.igas_CH4)then
                     ! Sneep et al, 2005
                     ! doi:10.1016/j.jqsrt.2004.07.025
                     T0(:) = 288.15
                     P0(:) = 1.01325e5
                     ! ng -> valid range of the measurements : 0.3251 - 0.6330 um
                     ng(:) = 1. + 46662e-8 + 4.02e-14*wn**2
                     Fk = 1. 
                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
                 elseif(igas.eq.igas_CO)then
                     ! Sneep et al, 2005
                     ! doi:10.1016/j.jqsrt.2004.07.025
                     T0(:) = 288.15
                     P0(:) = 1.01325e5
                     ! ng -> valid range of the measurements : 0.168 - 0.288 um
                     ng(:) = 1. + 22851e-8 + 0.456e4/(71427.**2 - wn**2)
                     Fk = 1.016
                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
                 elseif(igas.eq.igas_Ar)then
                     ! Thalman et al, 2014
                     ! doi:10.1016/j.jqsrt.2014.05.030
                     T0(:) = 288.15
                     P0(:) = 1.01325e5
                     ! ng -> valid range of the measurements : 0.288 - 0.546 um
                     ng(:) = 1. + (6432.135 + 286.06021e12/(14.4e9 - wn**2))*1e-8
                     Fk = 1.
                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
                 elseif(igas.eq.igas_O2)then
                     ! Thalman et al, 2014
                     ! doi:10.1016/j.jqsrt.2014.05.030
                     T0(:) = 273.15
                     P0(:) = 1.01325e5
                     ! ng -> valid range of the measurements : 0.303 - 2.0 um
                     ng(:) = 1. + (20564.8 + 2.480899e13/(4.09e9 - wn**2))*1e-8
                     Fk = 1.09 + 1.385e-11*wn**2 + 1.488e-20*wn**4
                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
                 else
                     print*,'No rayleigh scattering for ',trim(gnom(igas)),'. No data found.'
                 endif
                 
                 tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                 ! N=P/(kB*T)
                 ! pmid*scalep -> mbar to Pa
                 ! muvar/1000 -> g/mol to kg/mol
               
                 tauvar(:)=tauvar(:)+tauconsti(igas,:)*tauvari(igas,:)
                 
               endif !greater than 0.01

            enddo !ngasmx

            call blackn(dble(wn*1e2),dble(tstellar),df)
            df=df*(BWNV(N+1)-BWNV(N))/Nfine
            tauwei=tauwei+df
            tausum(:)=tausum(:)+tauvar(:)*df
         
         enddo !Nfine
         TAURAY(:,N)=tausum(:)/tauwei

      end do !L_NSPECTV


   end subroutine calc_rayleigh

end module calc_rayleigh_mod