!==============================================================================

      subroutine photolysis_online(nlayer, alt, press, temp,
     $           zmmean, rm,
     $           tau, sza, dist_sol, v_phot, e_phot, ig, ngrid,
     $           nreact)

!***********************************************************************
!
!   subject:
!   --------
!
!   photolysis online
!
!   VERSION: Extracted from LMDZ.MARS work of Franck Lefevre (Yassin Jaziri)
!   April 2019 - Yassin Jaziri add updates generic input (Yassin Jaziri)
!
!***********************************************************************

      use photolysis_mod
      use tracer_h
      use chimiedata_h, only: indexchim
      use types_asis,   only: nb_phot_hv_max, nb_phot_max, jlabel,
     $                        reactions

      implicit none

!     input 

      integer, intent(in) :: nlayer ! number of atmospheric layers
      integer, intent(in) :: ngrid  ! number of atmospheric columns

      real,    intent(in), dimension(nlayer) :: press, temp, zmmean ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
      real,    intent(in), dimension(nlayer) :: alt                 ! altitude (km)
      real,    intent(in), dimension(nlayer,nesp) :: rm             ! mixing ratios 
      real,    intent(in) :: tau                                    ! integrated aerosol optical depth at the surface
      real,    intent(in) :: sza                                    ! solar zenith angle (degrees)
      real,    intent(in) :: dist_sol                               ! solar distance (au)
      integer, intent(in) :: ig                                     ! grid point index
      integer, intent(in) :: nreact                                 ! number of reactions in reactions files

!     output

      real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot      ! photolysis rates (s-1)
      real (kind = 8), dimension(nlayer,nb_phot_max) :: e_phot      ! photolysis rates by energie (J.mol-1.s-1)
 
!     solar flux at mars

      real, dimension(nw) :: fmars                                  ! solar flux (w.m-2.nm-1)
      real :: factor

!     atmosphere

      real, dimension(nw)             :: albedo_chim                ! Surface albedo calculated on chemistry wavelenght grid
      real, dimension(nlayer)         :: colinc                     ! air column increment (molecule.cm-2) 
      real, dimension(nlayer)         :: airlev                     ! air density at each specified altitude (molec/cc)
      real, dimension(nlayer)         :: edir, edn, eup             ! normalised irradiances
      real, dimension(nlayer)         :: fdir, fdn, fup             ! normalised actinic fluxes
      real, dimension(nlayer)         :: saflux                     ! total actinic flux
      real, dimension(nlayer,nw)      :: dtrl                       ! rayleigh optical depth 
      real, dimension(nlayer,nw)      :: dtaer                      ! aerosol optical depth 
      real, dimension(nlayer,nw)      :: omaer                      ! aerosol single scattering albedo 
      real, dimension(nlayer,nw)      :: gaer                       ! aerosol asymmetry parameter
      real, dimension(nlayer,nw)      :: dtcld                      ! cloud optical depth 
      real, dimension(nlayer,nw)      :: omcld                      ! cloud single scattering albedo 
      real, dimension(nlayer,nw)      :: gcld                       ! cloud asymmetry parameter
      real, dimension(nlayer,nw)      :: dagas                      ! total gas optical depth
      real, dimension(nlayer,nw,nabs) :: dtgas                      ! optical depth for each gas 

      integer, dimension(0:nlayer)        :: nid
      real,    dimension(0:nlayer,nlayer) :: dsdh

      integer :: i, ilay, iw, ialt, iphot, ispe, ij, igas, ireact
      real    :: deltaj

!==== air column increments and rayleigh optical depth

      call setair(nlayer, nw, wl, wc, press, temp, zmmean, colinc, dtrl,
     $            airlev)

!====  set surface albedo

      call setalb(nw,wl,ig,ngrid,albedo_chim)

!==== set temperature-dependent cross-sections and optical depths

      iphot  = 0
      ij     = 0
      igas   = 0
      ireact = 0

      do while(iphot<nb_phot_hv_max)
        ij     = ij     + 1
        iphot  = iphot  + 1
        ireact = ireact + 1

        if (tdim(ij).eq.1) then
          ! Avoid to calculate several times dtgas for a same specie
          if (jlabelbis(iphot)) then
            ispe = indexchim(trim(jlabel(iphot,2)))
            do iw = 1,nw-1
              do ilay = 1,nlayer
                dtgas(ilay,iw,ij-igas) = colinc(ilay)*rm(ilay,ispe)
     $                                *xs(1,iw,ij)
              end do
            end do
          else
            igas = igas + 1
          end if
          do while(reactions(ireact)%rtype/=0)
            ireact = ireact + 1
          end do
          if (reactions(ireact)%three_prod) iphot = iphot + 1
        else
          call setsj(nb_phot_hv_max,nlayer,nw,temp,tdim(ij),
     $               xs(:tdim(ij),:,ij),
     $               xs_temp(:,ij),sj(:,:,iphot))
  
          ! Avoid to calculate several times dtgas for a same specie
          if (jlabelbis(iphot)) then
            ispe = indexchim(trim(jlabel(iphot,2)))
            do iw = 1,nw-1
              do ilay = 1,nlayer
                dtgas(ilay,iw,ij-igas) = colinc(ilay)*rm(ilay,ispe)
     $                                *sj(ilay,iw,iphot)
              end do
            end do
          else
            igas = igas + 1
          end if
  
          do iw = 1,nw-1
            do ilay = 1,nlayer
              sj(ilay,iw,iphot) = sj(ilay,iw,iphot)*qy(iw,ij)
            end do
          end do

          do while(reactions(ireact)%rtype/=0)
            ireact = ireact + 1
          end do
          if (reactions(ireact)%three_prod) then
            iphot = iphot + 1
            do iw = 1,nw-1
              do ilay = 1,nlayer
                sj(ilay,iw,iphot) = sj(ilay,iw,iphot-1)
              end do
            end do
          end if

        endif ! end if tdim .eq. 1

      end do ! end while

!     total gas optical depth

      dagas(:,:) = 0.

      do i = 1,nabs
         do iw = 1,nw-1
            do ilay = 1,nlayer
               dagas(ilay,iw) = dagas(ilay,iw) + dtgas(ilay,iw,i)
            end do
         end do
      end do

!==== set aerosol properties and optical depth

      call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)

!==== set cloud properties and optical depth

      call setcld(nlayer,nw,dtcld,omcld,gcld)

!==== slant path lengths in spherical geometry

      call sphers(nlayer,alt,sza,dsdh,nid)

!==== solar flux at mars

      factor = (1./dist_sol)**2.
      do iw = 1,nw-1
         fmars(iw) = f(iw)*factor
      end do

!==== initialise photolysis rates

      v_phot(:,1:nb_phot_hv_max) = 0.
      e_phot(:,1:nb_phot_hv_max) = 0.

!==== start of wavelength lopp

      do iw = 1,nw-1

!     monochromatic radiative transfer. outputs are:
!     normalized irradiances     edir(nlayer), edn(nlayer), eup(nlayer) 
!     normalized actinic fluxes  fdir(nlayer), fdn(nlayer), fup(nlayer)
!     where 
!     dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse

         call rtlink(nlayer, nw, iw, albedo_chim(iw),
     $               sza, dsdh, nid, dtrl,
     $               dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
     $               edir, edn, eup, fdir, fdn, fup)


!     spherical actinic flux

         do ilay = 1,nlayer
           saflux(ilay) = fmars(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay))
         end do

!     photolysis rate integration
!     (0.12/(wc(iw)*1e-9)) E(wc) en J.mol-1

         do i = 1,nb_phot_hv_max
            do ilay = 1,nlayer
               deltaj = saflux(ilay)*sj(ilay,iw,i)
               v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw))
               if (wc(iw).le.photoheat_lmax) then
                  e_phot(ilay,i) = e_phot(ilay,i)
     $                    + deltaj*(wu(iw)-wl(iw))*(0.12/(wc(iw)*1e-9))
               end if
            end do
         end do

!     eliminate small values

         where (v_phot(:,1:nb_phot_hv_max) < 1.e-30)
            v_phot(:,1:nb_phot_hv_max) = 0.
         end where
         where (e_phot(:,1:nb_phot_hv_max) < 1.e-30)
            e_phot(:,1:nb_phot_hv_max) = 0.
         end where

      end do ! iw

      contains

!==============================================================================

      subroutine setair(nlev, nw, wl, wc, press, temp, zmmean, 
     $                  colinc, dtrl, airlev)

*-----------------------------------------------------------------------------*
*=  PURPOSE:                                                                 =*
*=  computes air column increments and rayleigh optical depth                =*
*-----------------------------------------------------------------------------*

      use comcstfi_mod, only:  g, avocado

      implicit none

!     input:

      integer :: nlev, nw

      real, dimension(nw)   :: wl, wc              ! lower and central wavelength grid (nm)
      real, dimension(nlev) :: press, temp, zmmean ! pressure (hpa), temperature (k), molecular mass (g.mol-1)

!     output:

      real, dimension(nlev)    :: colinc ! air column increments (molecule.cm-2)
      real, dimension(nlev)    :: airlev ! air density at each specified altitude (molec/cc)
      real, dimension(nlev,nw) :: dtrl   ! rayleigh optical depth

!     local:

      real :: dp, nu
      real, dimension(nw) :: srayl
      integer :: ilev, iw

!     compute column increments 

      do ilev = 1, nlev-1
         dp = (press(ilev) - press(ilev+1))*100.
         colinc(ilev) = avocado*0.1*dp/(zmmean(ilev)*g)
      end do
      colinc(nlev) = avocado*0.1*press(nlev)*100./(zmmean(nlev)*g)

      do iw = 1, nw - 1

!        co2 rayleigh cross-section
!        ityaksov et al., chem. phys. lett., 462, 31-34, 2008

!         nu = 1./(wc(iw)*1.e-7)
!         srayl(iw) = 1.78e-26*nu**(4. + 0.625)
!         srayl(iw) = srayl(iw)*1.e-20 ! cm2

!        calcul Ityaksov et al., 2008 pour N2

         nu = 1./(wc(iw)*1.e-7)                     ! cm-1
         srayl(iw) = 1.8e-26*nu**(4. + 0.534)
         srayl(iw) = srayl(iw)*1.e-20

         do ilev = 1, nlev
            dtrl(ilev,iw) = colinc(ilev)*srayl(iw)   ! cm2
         end do
      end do

!     compute density

      do ilev = 1, nlev
         airlev(ilev) = press(ilev)/(1.38e-19*temp(ilev))
      end do

      end subroutine setair

!==============================================================================

      subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)

!-----------------------------------------------------------------------------*
!=  PURPOSE:                                                                 =*
!=  Set aerosol properties for each specified altitude layer.  Properties    =*
!=  may be wavelength dependent.                                             =*
!-----------------------------------------------------------------------------*

      implicit none

!     input

      integer :: nlayer                              ! number of vertical layers
      integer :: nw                                  ! number of wavelength grid points
      real, dimension(nlayer) :: alt                 ! altitude (km)
      real :: tau                                    ! integrated aerosol optical depth at the surface

!     output 

      real, dimension(nlayer,nw) :: dtaer            ! aerosol optical depth 
      real, dimension(nlayer,nw) :: omaer            ! aerosol single scattering albedo 
      real, dimension(nlayer,nw) :: gaer             ! aerosol asymmetry parameter

!     local

      integer :: ilay, iw
      real, dimension(nlayer) :: aer                                 ! dust extinction
      real :: omega, g, scaleh, gamma
      real :: dz, tautot, q0

      omega  = 0.622 ! single scattering albedo : wolff et al.(2010) at 258 nm
      g      = 0.88  ! asymmetry factor : mateshvili et al. (2007) at 210 nm
      scaleh = 10.   ! scale height (km)
      gamma  = 0.03  ! conrath parameter

      dtaer(:,:) = 0.
      omaer(:,:) = 0.
      gaer(:,:)  = 0.

!     optical depth profile:

      tautot = 0.
      do ilay = 1, nlayer-1
         dz = alt(ilay+1) - alt(ilay)
         tautot = tautot + exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
      end do
      
      q0 = tau/tautot
      do ilay = 1, nlayer-1
         dz = alt(ilay+1) - alt(ilay)
         dtaer(ilay,:) = q0*exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
         omaer(ilay,:) = omega
         gaer(ilay,:)  = g
      end do

      end subroutine setaer

!==============================================================================

      subroutine setcld(nlayer,nw,dtcld,omcld,gcld)

!-----------------------------------------------------------------------------*
!=  PURPOSE:                                                                 =*
!=  Set cloud properties for each specified altitude layer.  Properties      =*
!=  may be wavelength dependent.                                             =*
!-----------------------------------------------------------------------------*

      implicit none

!     input

      integer :: nlayer                              ! number of vertical layers
      integer :: nw                                  ! number of wavelength grid points

!     output 

      real, dimension(nlayer,nw) :: dtcld            ! cloud optical depth 
      real, dimension(nlayer,nw) :: omcld            ! cloud single scattering albedo 
      real, dimension(nlayer,nw) :: gcld             ! cloud asymmetry parameter

!     local

      integer :: ilay, iw

!     dtcld : optical depth
!     omcld : single scattering albedo
!     gcld  : asymmetry factor

      do ilay = 1, nlayer - 1
         do iw = 1, nw - 1
            dtcld(ilay,iw) = 0.       ! no clouds for the moment
            omcld(ilay,iw) = 0.99
            gcld(ilay,iw)  = 0.85
         end do
      end do

      end subroutine setcld

!==============================================================================

      subroutine sphers(nlev, z, zen, dsdh, nid)

!-----------------------------------------------------------------------------*
!=  PURPOSE:                                                                 =*
!=  Calculate slant path over vertical depth ds/dh in spherical geometry.    =*
!=  Calculation is based on:  A.Dahlback, and K.Stamnes, A new spheric model =*
!=  for computing the radiation field available for photolysis and heating   =*
!=  at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B)  =*
!-----------------------------------------------------------------------------*
!=  PARAMETERS:                                                              =*
!=  NZ      - INTEGER, number of specified altitude levels in the working (I)=*
!=            grid                                                           =*
!=  Z       - REAL, specified altitude working grid (km)                  (I)=*
!=  ZEN     - REAL, solar zenith angle (degrees)                          (I)=*
!=  DSDH    - REAL, slant path of direct beam through each layer crossed  (O)=*
!=            when travelling from the top of the atmosphere to layer i;     =*
!=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                            =*
!=  NID     - INTEGER, number of layers crossed by the direct beam when   (O)=*
!=            travelling from the top of the atmosphere to layer i;          =*
!=            NID(i), i = 0..NZ-1                                            =*
!-----------------------------------------------------------------------------*

      use comcstfi_mod, only: rad, pi

      implicit none

! input

      integer, intent(in) :: nlev
      real, dimension(nlev), intent(in) :: z
      real, intent(in) :: zen

! output

      INTEGER nid(0:nlev)
      REAL dsdh(0:nlev,nlev)

! more program constants

      REAL re, ze(nlev)
      REAL  dr
      real radius ! km

! local 

      real :: zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
      integer :: i, j, k, id, nlay

      REAL zd(0:nlev-1)

!-----------------------------------------------------------------------------

      radius = rad*1.e-3 ! rad [m] -> radius [km]
      dr = pi/180.
      zenrad = zen*dr

! number of layers:

      nlay = nlev - 1

! include the elevation above sea level to the radius of Mars:

      re = radius + z(1)

! correspondingly z changed to the elevation above Mars surface:

      DO k = 1, nlev
         ze(k) = z(k) - z(1)
      END DO

! inverse coordinate of z

      zd(0) = ze(nlev)
      DO k = 1, nlay
        zd(k) = ze(nlev - k)
      END DO

! initialise dsdh(i,j), nid(i)

      nid(:) = 0.
      dsdh(:,:) = 0.

! calculate ds/dh of every layer

      do i = 0,nlay
         rpsinz = (re + zd(i))*sin(zenrad)
 
        IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN
           nid(i) = -1
        ELSE

! Find index of layer in which the screening height lies

           id = i 
           if (zen > 90.) then
              do j = 1,nlay
                 IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND.
     $               (rpsinz .GE. ( zd(j) + re )) ) id = j
              end do
           end if
 
           do j = 1,id
              sm = 1.0
              IF (j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0)
     $            sm = -1.0
 
              rj = re + zd(j-1)
              rjp1 = re + zd(j)
 
              dhj = zd(j-1) - zd(j)
 
              ga = rj*rj - rpsinz*rpsinz
              gb = rjp1*rjp1 - rpsinz*rpsinz

              ga = max(ga, 0.)
              gb = max(gb, 0.)
 
              IF (id.GT.i .AND. j.EQ.id) THEN
                 dsj = sqrt(ga)
              ELSE
                 dsj = sqrt(ga) - sm*sqrt(gb)
              END IF
              dsdh(i,j) = dsj/dhj
            end do
            nid(i) = id
         end if
      end do ! i = 0,nlay

      end subroutine sphers

!==============================================================================

      SUBROUTINE rtlink(nlev, nw, iw, ag, zen, dsdh, nid, dtrl, 
     $                  dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
     $                  edir, edn, eup, fdir, fdn, fup)

      implicit none

! input

      integer, intent(in) :: nlev, nw, iw                       ! number of wavelength grid points
      REAL ag
      REAL zen
      REAL dsdh(0:nlev,nlev)
      INTEGER nid(0:nlev)

      REAL dtrl(nlev,nw)
      REAL dagas(nlev,nw)
      REAL dtcld(nlev,nw), omcld(nlev,nw), gcld(nlev,nw)
      REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw)

! output

      REAL edir(nlev), edn(nlev), eup(nlev)
      REAL fdir(nlev), fdn(nlev), fup(nlev)

! local:

      REAL dt(nlev), om(nlev), g(nlev)
      REAL dtabs,dtsct,dscld,dsaer,dacld,daaer
      INTEGER i, ii
      real, parameter :: largest = 1.e+36

! specific two ps2str

      REAL ediri(nlev), edni(nlev), eupi(nlev)
      REAL fdiri(nlev), fdni(nlev), fupi(nlev)

      logical, save :: delta = .true.

!_______________________________________________________________________

! initialize:

      do i = 1, nlev
         fdir(i) = 0.
         fup(i)  = 0.
         fdn(i)  = 0.
         edir(i) = 0.
         eup(i)  = 0.
         edn(i)  = 0.
      end do

      do i = 1, nlev - 1
         dscld = dtcld(i,iw)*omcld(i,iw)
         dacld = dtcld(i,iw)*(1.-omcld(i,iw))

         dsaer = dtaer(i,iw)*omaer(i,iw)
         daaer = dtaer(i,iw)*(1.-omaer(i,iw))

         dtsct = dtrl(i,iw) + dscld + dsaer
         dtabs = dagas(i,iw) + dacld + daaer

         dtabs = amax1(dtabs,1./largest)
         dtsct = amax1(dtsct,1./largest)

!     invert z-coordinate:

         ii = nlev - i
         dt(ii) = dtsct + dtabs
         om(ii) = dtsct/(dtsct + dtabs)
         IF(dtsct .EQ. 1./largest) om(ii) = 1./largest
         g(ii) = (gcld(i,iw)*dscld + 
     $            gaer(i,iw)*dsaer)/dtsct
      end do

!     call rt routine:

      call ps2str(nlev, zen, ag, dt, om, g,
     $            dsdh, nid, delta,
     $            fdiri, fupi, fdni, ediri, eupi, edni)

!     output (invert z-coordinate)

      do i = 1, nlev
         ii = nlev - i + 1
         fdir(i) = fdiri(ii)
         fup(i) = fupi(ii)
         fdn(i) = fdni(ii)
         edir(i) = ediri(ii)
         eup(i) = eupi(ii)
         edn(i) = edni(ii)
      end do

      end subroutine rtlink

*=============================================================================*

      subroutine ps2str(nlev,zen,rsfc,tauu,omu,gu,
     $                  dsdh, nid, delta,
     $                  fdr, fup, fdn, edr, eup, edn)

!-----------------------------------------------------------------------------*
!=  PURPOSE:                                                                 =*
!=  Solve two-stream equations for multiple layers.  The subroutine is based =*
!=  on equations from:  Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=*
!=  It contains 9 two-stream methods to choose from.  A pseudo-spherical     =*
!=  correction has also been added.                                          =*
!-----------------------------------------------------------------------------*
!=  PARAMETERS:                                                              =*
!=  NLEVEL  - INTEGER, number of specified altitude levels in the working (I)=*
!=            grid                                                           =*
!=  ZEN     - REAL, solar zenith angle (degrees)                          (I)=*
!=  RSFC    - REAL, surface albedo at current wavelength                  (I)=*
!=  TAUU    - REAL, unscaled optical depth of each layer                  (I)=*
!=  OMU     - REAL, unscaled single scattering albedo of each layer       (I)=*
!=  GU      - REAL, unscaled asymmetry parameter of each layer            (I)=*
!=  DSDH    - REAL, slant path of direct beam through each layer crossed  (I)=*
!=            when travelling from the top of the atmosphere to layer i;     =*
!=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                            =*
!=  NID     - INTEGER, number of layers crossed by the direct beam when   (I)=*
!=            travelling from the top of the atmosphere to layer i;          =*
!=            NID(i), i = 0..NZ-1                                            =*
!=  DELTA   - LOGICAL, switch to use delta-scaling                        (I)=*
!=            .TRUE. -> apply delta-scaling                                  =*
!=            .FALSE.-> do not apply delta-scaling                           =*
!=  FDR     - REAL, contribution of the direct component to the total     (O)=*
!=            actinic flux at each altitude level                            =*
!=  FUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
!=            the total actinic flux at each altitude level                  =*
!=  FDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
!=            the total actinic flux at each altitude level                  =*
!=  EDR     - REAL, contribution of the direct component to the total     (O)=*
!=            spectral irradiance at each altitude level                     =*
!=  EUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
!=            the total spectral irradiance at each altitude level           =*
!=  EDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
*=            the total spectral irradiance at each altitude level           =*
!-----------------------------------------------------------------------------*

      implicit none

! input:

      INTEGER nlev
      REAL zen, rsfc
      REAL tauu(nlev), omu(nlev), gu(nlev)
      REAL dsdh(0:nlev,nlev)
      INTEGER nid(0:nlev)
      LOGICAL delta

! output:

      REAL fup(nlev),fdn(nlev),fdr(nlev)
      REAL eup(nlev),edn(nlev),edr(nlev)

! local:

      REAL tausla(0:nlev), tauc(0:nlev)
      REAL mu2(0:nlev), mu, sum

! internal coefficients and matrix

      REAL lam(nlev),taun(nlev),bgam(nlev)
      REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev)
      REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev)
      REAL mu1(nlev)
      INTEGER row
      REAL a(2*nlev),b(2*nlev),d(2*nlev),e(2*nlev),y(2*nlev)

! other:

      REAL pifs, fdn0
      REAL gi(nlev), omi(nlev), tempg
      REAL f, g, om
      REAL gam1, gam2, gam3, gam4
      real, parameter :: largest = 1.e+36
      real, parameter :: precis = 1.e-7

! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
! in delta-function, modified quadrature, hemispheric constant,
! Hybrid modified Eddington-delta function metods, p633,Table1.
! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440, 
! uncomment the following two lines and the appropriate statements further
! down.
!     REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
!    >     BETA1, BETAn, amu1, subd

      REAL expon, expon0, expon1, divisr, temp, up, dn
      REAL ssfc
      INTEGER nlayer, mrows, lev

      INTEGER i, j

! Some additional program constants:

      real pi, dr
      REAL eps
      PARAMETER (eps = 1.E-3)
!_______________________________________________________________________

! MU = cosine of solar zenith angle
! RSFC = surface albedo
! TAUU =  unscaled optical depth of each layer
! OMU  =  unscaled single scattering albedo
! GU   =  unscaled asymmetry factor
! KLEV = max dimension of number of layers in atmosphere
! NLAYER = number of layers in the atmosphere
! NLEVEL = nlayer + 1 = number of levels

! initial conditions:  pi*solar flux = 1;  diffuse incidence = 0

      pifs = 1.      
      fdn0 = 0.

      nlayer = nlev - 1

      pi = acos(-1.)
      dr = pi/180.
      mu = COS(zen*dr)

!************* compute coefficients for each layer:
! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
! EXPON0 = calculation of e when TAU is zero
! EXPON1 = calculation of e when TAU is TAUN
! CUP and CDN = calculation when TAU is zero
! CUPTN and CDNTN = calc. when TAU is TAUN
! DIVISR = prevents division by zero

      do j = 0, nlev
         tauc(j) = 0.
         tausla(j) = 0.
         mu2(j) = 1./SQRT(largest)
      end do

      IF (.NOT. delta) THEN
         DO i = 1, nlayer
            gi(i) = gu(i)
            omi(i) = omu(i)
            taun(i) = tauu(i)
         END DO
      ELSE 

! delta-scaling. Have to be done for delta-Eddington approximation, 
! delta discrete ordinate, Practical Improved Flux Method, delta function,
! and Hybrid modified Eddington-delta function methods approximations

         DO i = 1, nlayer
            f = gu(i)*gu(i)
            gi(i) = (gu(i) - f)/(1 - f)
            omi(i) = (1 - f)*omu(i)/(1 - omu(i)*f)       
            taun(i) = (1 - omu(i)*f)*tauu(i)
         END DO
      END IF

! calculate slant optical depth at the top of the atmosphere when zen>90.
! in this case, higher altitude of the top layer is recommended.

      IF (zen .GT. 90.0) THEN
         IF (nid(0) .LT. 0) THEN
            tausla(0) = largest
         ELSE
            sum = 0.0
            DO j = 1, nid(0)
               sum = sum + 2.*taun(j)*dsdh(0,j)
            END DO
            tausla(0) = sum 
         END IF
      END IF

      DO 11, i = 1, nlayer
          g = gi(i)
          om = omi(i)
          tauc(i) = tauc(i-1) + taun(i)

! stay away from 1 by precision.  For g, also stay away from -1

          tempg = AMIN1(abs(g),1. - precis)
          g = SIGN(tempg,g)
          om = AMIN1(om,1.-precis)

! calculate slant optical depth
              
          IF (nid(i) .LT. 0) THEN
             tausla(i) = largest
          ELSE
             sum = 0.0
             DO j = 1, MIN(nid(i),i)
                sum = sum + taun(j)*dsdh(i,j)
             END DO
             DO j = MIN(nid(i),i)+1,nid(i)
                sum = sum + 2.*taun(j)*dsdh(i,j)
             END DO
             tausla(i) = sum 
             IF (tausla(i) .EQ. tausla(i-1)) THEN
                mu2(i) = SQRT(largest)
             ELSE
                mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
                mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)),
     $                     mu2(i) )
             END IF
          END IF

!** the following gamma equations are from pg 16,289, Table 1
!** save mu1 for each approx. for use in converting irradiance to actinic flux

! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):

c       gam1 =  (7. - om*(4. + 3.*g))/4.
c       gam2 = -(1. - om*(4. - 3.*g))/4.
c       gam3 = (2. - 3.*g*mu)/4.
c       gam4 = 1. - gam3
c       mu1(i) = 0.5

* quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):

c          gam1 = 1.7320508*(2. - om*(1. + g))/2.
c          gam2 = 1.7320508*om*(1. - g)/2.
c          gam3 = (1. - 1.7320508*g*mu)/2.
c          gam4 = 1. - gam3
c          mu1(i) = 1./sqrt(3.)
         
* hemispheric mean (Toon et al., 1089, JGR, 94, 16287):

           gam1 = 2. - om*(1. + g)
           gam2 = om*(1. - g)
           gam3 = (2. - g*mu)/4.
           gam4 = 1. - gam3
           mu1(i) = 0.5

* PIFM  (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
c         GAM1 = 0.25*(8. - OM*(5. + 3.*G))
c         GAM2 = 0.75*OM*(1.-G)
c         GAM3 = 0.25*(2.-3.*G*MU)
c         GAM4 = 1. - GAM3
c         mu1(i) = 0.5

* delta discrete ordinates  (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
c         GAM1 = 0.5*1.7320508*(2. - OM*(1. + G))
c         GAM2 = 0.5*1.7320508*OM*(1.-G)
c         GAM3 = 0.5*(1.-1.7320508*G*MU)
c         GAM4 = 1. - GAM3
c         mu1(i) = 1./sqrt(3.)

* Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
* in delta-function, modified quadrature, hemispheric constant,
* Hybrid modified Eddington-delta function metods, p633,Table1.
* W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
* W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
c      YLM0 = 2.
c      YLM2 = -3.*G*MU
c      YLM4 = 0.875*G**3*MU*(5.*MU**2-3.)
c      YLM6=-0.171875*G**5*MU*(15.-70.*MU**2+63.*MU**4)
c     YLM8=+0.073242*G**7*MU*(-35.+315.*MU**2-693.*MU**4
c    *+429.*MU**6)
c     YLM10=-0.008118*G**9*MU*(315.-4620.*MU**2+18018.*MU**4
c    *-25740.*MU**6+12155.*MU**8)
c     YLM12=0.003685*G**11*MU*(-693.+15015.*MU**2-90090.*MU**4
c    *+218790.*MU**6-230945.*MU**8+88179.*MU**10)
c      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
c      YLMS=0.25*YLMS
c      BETA0 = YLMS
c
c         amu1=1./1.7320508
c      YLM0 = 2.
c      YLM2 = -3.*G*amu1
c      YLM4 = 0.875*G**3*amu1*(5.*amu1**2-3.)
c      YLM6=-0.171875*G**5*amu1*(15.-70.*amu1**2+63.*amu1**4)
c     YLM8=+0.073242*G**7*amu1*(-35.+315.*amu1**2-693.*amu1**4
c    *+429.*amu1**6)
c     YLM10=-0.008118*G**9*amu1*(315.-4620.*amu1**2+18018.*amu1**4
c    *-25740.*amu1**6+12155.*amu1**8)
c     YLM12=0.003685*G**11*amu1*(-693.+15015.*amu1**2-90090.*amu1**4
c    *+218790.*amu1**6-230945.*amu1**8+88179.*amu1**10)
c      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
c      YLMS=0.25*YLMS
c      BETA1 = YLMS
c
c         BETAn = 0.25*(2. - 1.5*G-0.21875*G**3-0.085938*G**5
c    *-0.045776*G**7)


* Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
c         subd=4.*(1.-G*G*(1.-MU))
c         GAM1 = (7.-3.*G*G-OM*(4.+3.*G)+OM*G*G*(4.*BETA0+3.*G))/subd
c         GAM2 =-(1.-G*G-OM*(4.-3.*G)-OM*G*G*(4.*BETA0+3.*G-4.))/subd
c         GAM3 = BETA0
c         GAM4 = 1. - GAM3
c         mu1(i) = (1. - g*g*(1.- mu) )/(2. - g*g)

*****
* delta function  (Meador, and Weaver, 1980, JAS, 37, 630):
c         GAM1 = (1. - OM*(1. - beta0))/MU
c         GAM2 = OM*BETA0/MU
c         GAM3 = BETA0
c         GAM4 = 1. - GAM3
c         mu1(i) = mu
*****
* modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
c         GAM1 = 1.7320508*(1. - OM*(1. - beta1))
c         GAM2 = 1.7320508*OM*beta1
c         GAM3 = BETA0
c         GAM4 = 1. - GAM3
c         mu1(i) = 1./sqrt(3.)

* hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
c         GAM1 = 2.*(1. - OM*(1. - betan))
c         GAM2 = 2.*OM*BETAn
c         GAM3 = BETA0
c         GAM4 = 1. - GAM3
c         mu1(i) = 0.5

*****

* lambda = pg 16,290 equation 21
* big gamma = pg 16,290 equation 22
* if gam2 = 0., then bgam = 0. 

         lam(i) = sqrt(gam1*gam1 - gam2*gam2)

         IF (gam2 .NE. 0.) THEN
            bgam(i) = (gam1 - lam(i))/gam2
         ELSE
            bgam(i) = 0.
         END IF

         expon = EXP(-lam(i)*taun(i))

* e1 - e4 = pg 16,292 equation 44
         
         e1(i) = 1. + bgam(i)*expon
         e2(i) = 1. - bgam(i)*expon
         e3(i) = bgam(i) + expon
         e4(i) = bgam(i) - expon

* the following sets up for the C equations 23, and 24
* found on page 16,290
* prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
* which is approx equiv to shifting MU by 0.5*EPS* (MU)**3

         expon0 = EXP(-tausla(i-1))
         expon1 = EXP(-tausla(i))
          
         divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i))
         temp = AMAX1(eps,abs(divisr))
         divisr = SIGN(temp,divisr)

         up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr
         dn = om*pifs*((gam1 + 1./mu2(i))*gam4 + gam2*gam3)/divisr
         
* cup and cdn are when tau is equal to zero
* cuptn and cdntn are when tau is equal to taun

         cup(i) = up*expon0
         cdn(i) = dn*expon0
         cuptn(i) = up*expon1
         cdntn(i) = dn*expon1
 
   11 CONTINUE

***************** set up matrix ******
* ssfc = pg 16,292 equation 37  where pi Fs is one (unity).

      ssfc = rsfc*mu*EXP(-tausla(nlayer))*pifs

* MROWS = the number of rows in the matrix

      mrows = 2*nlayer     
      
* the following are from pg 16,292  equations 39 - 43.
* set up first row of matrix:

      i = 1
      a(1) = 0.
      b(1) = e1(i)
      d(1) = -e2(i)
      e(1) = fdn0 - cdn(i)

      row=1

* set up odd rows 3 thru (MROWS - 1):

      i = 0
      DO 20, row = 3, mrows - 1, 2
         i = i + 1
         a(row) = e2(i)*e3(i) - e4(i)*e1(i)
         b(row) = e1(i)*e1(i + 1) - e3(i)*e3(i + 1)
         d(row) = e3(i)*e4(i + 1) - e1(i)*e2(i + 1)
         e(row) = e3(i)*(cup(i + 1) - cuptn(i)) + 
     $        e1(i)*(cdntn(i) - cdn(i + 1))
   20 CONTINUE

* set up even rows 2 thru (MROWS - 2): 

      i = 0
      DO 30, row = 2, mrows - 2, 2
         i = i + 1
         a(row) = e2(i + 1)*e1(i) - e3(i)*e4(i + 1)
         b(row) = e2(i)*e2(i + 1) - e4(i)*e4(i + 1)
         d(row) = e1(i + 1)*e4(i + 1) - e2(i + 1)*e3(i + 1)
         e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) - 
     $        (cdn(i + 1) - cdntn(i))*e4(i + 1)
   30 CONTINUE

* set up last row of matrix at MROWS:

      row = mrows
      i = nlayer
      
      a(row) = e1(i) - rsfc*e3(i)
      b(row) = e2(i) - rsfc*e4(i)
      d(row) = 0.
      e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)

* solve tri-diagonal matrix:

      CALL tridiag(a, b, d, e, y, mrows)

**** unfold solution of matrix, compute output fluxes:

      row = 1 
      lev = 1
      j = 1
      
* the following equations are from pg 16,291  equations 31 & 32

      fdr(lev) = EXP( -tausla(0) )
      edr(lev) = mu * fdr(lev)
      edn(lev) = fdn0
      eup(lev) =  y(row)*e3(j) - y(row + 1)*e4(j) + cup(j)
      fdn(lev) = edn(lev)/mu1(lev)
      fup(lev) = eup(lev)/mu1(lev)

      DO 60, lev = 2, nlayer + 1
         fdr(lev) = EXP(-tausla(lev-1))
         edr(lev) =  mu *fdr(lev)
         edn(lev) =  y(row)*e3(j) + y(row + 1)*e4(j) + cdntn(j)
         eup(lev) =  y(row)*e1(j) + y(row + 1)*e2(j) + cuptn(j)
         fdn(lev) = edn(lev)/mu1(j)
         fup(lev) = eup(lev)/mu1(j)

         row = row + 2
         j = j + 1
   60 CONTINUE

      end subroutine ps2str

*=============================================================================*

      subroutine tridiag(a,b,c,r,u,n)

!_______________________________________________________________________
! solves tridiagonal system.  From Numerical Recipies, p. 40
!_______________________________________________________________________

      IMPLICIT NONE

! input:

      INTEGER n
      REAL a, b, c, r
      DIMENSION a(n),b(n),c(n),r(n)

! output:

      REAL u
      DIMENSION u(n)

! local:

      INTEGER j

      REAL bet, gam
      DIMENSION gam(n)
!_______________________________________________________________________

      IF (b(1) .EQ. 0.) STOP 1001
      bet   = b(1)
      u(1) = r(1)/bet
      DO 11, j = 2, n   
         gam(j) = c(j - 1)/bet
         bet = b(j) - a(j)*gam(j)
         IF (bet .EQ. 0.) STOP 2002 
         u(j) = (r(j) - a(j)*u(j - 1))/bet
   11 CONTINUE
      DO 12, j = n - 1, 1, -1  
         u(j) = u(j) - gam(j + 1)*u(j + 1)
   12 CONTINUE
!_______________________________________________________________________

      end subroutine tridiag

!==============================================================================
 
      subroutine setalb(nw,wl,ig,ngrid,albedo_chim)
 
!-----------------------------------------------------------------------------*
!=  PURPOSE:                                                                 =*
!=  Set the albedo of the surface.  The albedo is assumed to be Lambertian,  =*
!=  i.e., the reflected light is isotropic, and idependt of the direction    =*
!=  of incidence of light.  Albedo can be chosen to be wavelength dependent. =*
!-----------------------------------------------------------------------------*
!=  PARAMETERS:                                                              =*
!=  NW      - INTEGER, number of specified intervals + 1 in working       (I)=*
!=            wavelength grid                                                =*
!=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
!=           working wavelength grid                                         =*
!=  ALBEDO  - REAL, surface albedo at each specified wavelength           (O)=*
!-----------------------------------------------------------------------------*

      use chimiedata_h,  only: albedo_snow_chim, albedo_co2_ice_chim
      use slab_ice_h,    only: h_alb_ice, alb_ice_min, alb_ice_max
      use tracer_h,      only: igcm_h2o_ice, igcm_co2_ice
      use callkeys_mod,  only: ok_slab_ocean, co2cond, alb_ocean
      use phys_state_var_mod, only: albedo_bareground, 
     &                              rnat, qsurf, sea_ice,
     &                              pctsrf_sic, tsurf, capcal
      use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater

      implicit none

! input: (wavelength working grid data)

      integer, intent(in) :: ngrid  ! number of atmospheric columns
      INTEGER nw
      INTEGER ig             ! grid point index
      REAL wl(nw)

! output:

      REAL albedo_chim(nw)

! local:

      INTEGER iw
!      REAL alb
      real zfra, alb_ice, twater, hice
      real snowlayer
      parameter (snowlayer=33.0)        ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm

!     0.015: mean value from clancy et al., icarus, 49-63, 1999.

!      alb = 0.015
!      alb = albedo_phys

!      do iw = 1, nw - 1
!         albedo_chim(iw) = alb
!      end do

!     See hydrol.F90 for taking into account new tendencies

      if(nint(rnat(ig)).eq.0)then

        if(ok_slab_ocean) then
        
          zfra = MAX(0.0,MIN(1.0,qsurf(ig,igcm_h2o_ice)/45.0))     ! Snow Fraction (Critical height 45kg/m2~15cm)
          alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min)            ! Ice Albedo
     &                        *exp(-sea_ice(ig)/h_alb_ice)
          ! Albedo final calculation :
          do iw=1,nw - 1
             albedo_chim(iw) = pctsrf_sic(ig)*
     &                        (albedo_snow_chim(iw)*zfra
     &                        + alb_ice*(1.0-zfra))
     &                        + (1.-pctsrf_sic(ig))*alb_ocean
          enddo  

        else !ok_slab_ocean
            
            
!     calculate oceanic ice height including the latent heat of ice formation
!     hice is the height of oceanic ice with a maximum of maxicethick.
          hice    = qsurf(ig,igcm_h2o_ice)/rhowater ! update hice to include recent snowfall
          twater  = tsurf(ig) - hice*RLFTT*rhowater/capcal(ig)
          ! this is temperature water would have if we melted entire ocean ice layer

          if(twater .lt. T_h2O_ice_liq)then

            do iw=1,nw - 1
              albedo_chim(iw) = albedo_snow_chim(iw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
            enddo

          else

            DO iw=1,nw - 1
              albedo_chim(iw) = alb_ocean
              if(ngrid.eq.1) then
                albedo_chim(iw) = albedo_bareground(ig)
              endif
            ENDDO

          endif

        endif!(ok_slab_ocean)


!     Continent
!     ---------
      elseif (nint(rnat(ig)).eq.1) then

!     re-calculate continental albedo
        if (qsurf(ig,igcm_h2o_ice).ge.snowlayer) then
          DO iw=1,nw - 1
            albedo_chim(iw) = albedo_snow_chim(iw)
          ENDDO
        else
          DO iw=1,nw - 1
            albedo_chim(iw) = albedo_bareground(ig)
     &                      + (albedo_snow_chim(iw)
     &                         - albedo_bareground(ig))
     &                      *qsurf(ig,igcm_h2o_ice)/snowlayer
          ENDDO
        endif

      else

        print*,'Surface type not recognised in photolysis_online.F!'
        print*,'Exiting...'
        call abort

      endif

!     Re-add the albedo effects of CO2 ice if necessary
!     -------------------------------------------------
      if(co2cond)then

          if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015
            DO iw=1,nw - 1
              albedo_chim(iw) = albedo_co2_ice_chim(iw)
            ENDDO
          endif
            
      endif ! co2cond

      end subroutine setalb

!==============================================================================
 
      subroutine setsj(nd,nlayer,nw,tlay,tdiml,xsl,xs_templ,sjl)


      implicit none

! input: (wavelength working grid data)

      integer :: nd                                      ! number of photolysis rates
      integer :: nlayer                                  ! number of vertical layers
      integer :: nw                                      ! number of wavelength grid points
      integer :: tdiml                                   ! number of different temperature cross section file
      real    :: xsl(tdiml,nw)                           ! cross section (cm2) from reading files
      real    :: xs_templ(tdiml)                         ! temperature of the cross section (cm2)
      real    :: tlay(nlayer)                            ! temperature (K)

! output:

      real    :: sjl(nlayer,nw)                          ! output cross section (cm2)


! local:

      INTEGER ilay,iw,tpos

        do iw = 1,nw-1
          do ilay = 1,nlayer
            tpos = 1
            do while(tlay(ilay)>xs_templ(tpos) .and. tpos<tdiml)
              tpos = tpos + 1
            end do
            if (tpos.eq.1 .or. tpos.eq.tdiml) then
              sjl(ilay,iw) = xsl(tpos,iw)
            else
              sjl(ilay,iw) = ( (xsl(tpos,iw)-xsl(tpos-1,iw))*
     $                       (tlay(ilay)-xs_templ(tpos-1))
     $                       /(xs_templ(tpos)-xs_templ(tpos-1))
     $                       + xsl(tpos-1,iw) )
            end if
          end do
        end do

      end subroutine setsj

      end subroutine photolysis_online