module photolysis_mod

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

    use types_asis
    implicit none
  
    ! photolysis
  
    integer, save :: nabs                      ! number of absorbing gases
  
    ! spectral grid
  
    integer, save :: nw                        ! number of spectral intervals (low-res)
    integer, save :: mopt                      ! high-res/low-res switch
  
    real, allocatable, save :: wl(:)           ! lower  wavelength for each interval (nm)
    real, allocatable, save :: wc(:)           ! center wavelength for each interval (nm)
    real, allocatable, save :: wu(:)           ! upper  wavelength for each interval (nm)
    real, save :: photoheat_lmax               ! maximum wavelength until photochemical heating is calculated
  
    ! solar flux
  
    real, allocatable, save :: fstar1AU(:)     ! stellar flux (photon.s-1.nm-1.cm-2) at 1 au
  
    ! cross-sections and quantum yields
  
    real,    allocatable, save :: qy(:,:)      ! photodissociation yield
    real,    allocatable, save :: xs(:,:,:)    ! absorption cross-section (cm2)
    real,    allocatable, save :: sj(:,:,:)    ! general cross-section array by photodissociation yield (cm2)
    real,    allocatable, save :: xs_temp(:,:) ! absorption cross-section temperature (K)
    integer, allocatable, save :: tdim(:)      ! absorption cross-section temperature dimension
    logical, allocatable, save :: jlabelbis(:) ! check in jlabel if the species is included in more than 1 photolysis
  
    contains
  
    subroutine init_photolysis(nlayer,nreact)
  
    use types_asis, only: reactions
    use tracer_h
    use chimiedata_h, only: indexchim
    use ioipsl_getin_p_mod, only: getin_p
  
    integer, intent(in) :: nlayer              ! number of atmospheric layers
    integer, intent(in) :: nreact              ! number of reactions in reactions files
    integer :: iphot,tdim_max,i3prod,ij,iw,ilay,ireact
    integer :: specheck(nesp)
  
    ! initialise on-line photolysis
  
    ! mopt = 1 high-resolution
    ! mopt = 2 low-resolution (recommended for gcm use)
    ! mopt = 3 input file reading grid
  
    mopt = 3
  
    ! set wavelength grid
  
    call gridw(nw,wl,wc,wu,mopt)
  
    allocate(fstar1AU(nw))
  
    ! read and grid solar flux data
   
    call rdsolarflux(nw,wl,wc,fstar1AU)
  
    ! read maximum wavelength until photochemical heating is calculated
    write(*,*) "Maximum wavelength until photochemical heating is calculated"
    photoheat_lmax=wc(nw-1)         ! default value last point of wavelength grid
    call getin_p("photoheat_lmax",photoheat_lmax)
    write(*,*) "photoheat_lmax = ",photoheat_lmax
  
    ! calculate nabs number of absorbing gases
    allocate(jlabelbis(nb_phot_hv_max))
    nabs = 0
    specheck(:) = 0
    jlabelbis(:) = .true.
    do iphot=1,nb_phot_hv_max
      if (specheck(indexchim(trim(jlabel(iphot,2)))) .eq. 0) then
        specheck(indexchim(trim(jlabel(iphot,2)))) = 1
        nabs = nabs + 1
      else
        jlabelbis(iphot) = .false.
      endif
    end do
  
    ! Get temperature dimension and allocate
    allocate(tdim(nb_hv_max))
    tdim(:)  = 1
    tdim_max = 1
    do iphot=1,nb_hv_max
      read(jparamline(iphot),*) tdim(iphot)
      if (tdim(iphot) > tdim_max) tdim_max = tdim(iphot)
    end do
  
    ! allocation
    allocate(qy(nw,nb_hv_max))
    allocate(xs(tdim_max,nw,nb_hv_max))
    allocate(sj(nlayer,nw,nb_phot_hv_max))
    allocate(xs_temp(tdim_max,nb_hv_max))
    xs(:,:,:)    = 0.
    xs_temp(:,:) = 0.
  
    print*, 'WARNING: for all branching photolysis of one specie'
    print*, '         the cross section files should be the same'
    print*, '         they are used for optical depths calculation'
    print*, '         only the quantum yield should be different'
  
    i3prod = 0
    ireact = 0
  
    ! read and grid cross-sections
    do iphot=1,nb_hv_max
      ireact = ireact + 1
      call rdxsi(nw,wl,xs(:tdim(iphot),:,iphot),jparamline(iphot),xs_temp(:tdim(iphot),iphot),qy(:,iphot),tdim(iphot),jlabel(iphot+i3prod,2))
      do while(reactions(ireact)%rtype/=0)
        ireact = ireact + 1
      end do
      if (reactions(ireact)%three_prod) i3prod = i3prod + 1
    end do
  
    ! init sj for no temperature dependent cross sections (tdim.eq.1)
    iphot  = 0
    ij     = 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
         do iw = 1,nw-1
           do ilay = 1,nlayer
             sj(ilay,iw,iphot) = xs(1,iw,ij)*qy(iw,ij)
           end do
         end do
       endif
       do while(reactions(ireact)%rtype/=0)
         ireact = ireact + 1
       end do
       if (reactions(ireact)%three_prod) then
         iphot = iphot + 1
         if (tdim(ij).eq.1) then
           do iw = 1,nw-1
             do ilay = 1,nlayer
               sj(ilay,iw,iphot) = sj(ilay,iw,iphot-1)
             end do
           end do
         endif
       end if
    enddo

    end subroutine init_photolysis
  
    !==============================================================================
  
    subroutine gridw(nw,wl,wc,wu,mopt)
  
        use datafile_mod, only: datadir
        use ioipsl_getin_p_mod, only: getin_p
  
        !==========================================================================
        ! Create the wavelength grid for all interpolations and radiative transfer 
        ! calculations.  Grid may be irregularly spaced.  Wavelengths are in nm.   
        ! No gaps are allowed within the wavelength grid.                          
        !==========================================================================
  
        implicit none
  
        !     input
   
        integer :: mopt    ! high-res/low-res switch
  
        !     output
  
        integer :: nw      ! number of wavelength grid points
        real, allocatable :: wl(:), wc(:), wu(:)   ! lower, center, upper wavelength for each interval
  
        !     local
  
        real :: wincr    ! wavelength increment
        integer :: iw, kw, ierr
        character(len=200) :: fil
        character(len=150) :: gridwfile
  
        !     mopt = 1    high-resolution mode (3789 intervals)
        !
        !                   0-108 nm :  1.0  nm
        !                 108-124 nm :  0.1  nm
        !                 124-175 nm :  0.5  nm
        !                 175-205 nm :  0.01 nm
        !                 205-365 nm :  0.5  nm
        !                 365-850 nm :  5.0  nm  
        !
        !     mopt = 2    low-resolution mode
        !
        !                    0-60 nm :  6.0 nm
        !                   60-80 nm :  2.0 nm
        !                   80-85 nm :  5.0 nm
        !                  85-117 nm :  2.0 nm
        !                 117-120 nm :  5.0 nm
        !                 120-123 nm :  0.2 nm
        !                 123-163 nm :  5.0 nm
        !                 163-175 nm :  2.0 nm
        !                 175-205 nm :  0.5 nm
        !                 205-245 nm :  5.0 nm
        !                 245-415 nm : 10.0 nm
        !                 415-815 nm : 50.0 nm
        !
        !     mopt = 3    input file reading grid
        !
        !                 read "gridw.dat" in datadir/cross_sections/
  
        if (mopt == 1) then   ! high-res
  
        nw = 3789
        allocate(wl(nw))
        allocate(wu(nw))
        allocate(wc(nw))
        ! define wavelength intervals of width 1.0 nm from 0 to 108 nm:
  
        kw = 0
        wincr = 1.0
        do iw = 0, 107
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        end do
  
        ! define wavelength intervals of width 0.1 nm from 108 to 124 nm:
  
        wincr = 0.1
        do iw = 1080, 1239, 1
          kw = kw + 1
          wl(kw) = real(iw)/10.
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        end do
  
        ! define wavelength intervals of width 0.5 nm from 124 to 175 nm:
  
        wincr = 0.5
        do iw = 1240, 1745, 5
          kw = kw + 1
          wl(kw) = real(iw)/10.
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        end do
  
        ! define wavelength intervals of width 0.01 nm from 175 to 205 nm:
  
        wincr = 0.01
        do iw = 17500, 20499, 1
          kw = kw + 1
          wl(kw) = real(iw)/100.
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        end do
  
        ! define wavelength intervals of width 0.5 nm from 205 to 365 nm:
  
        wincr = 0.5
        do iw = 2050, 3645, 5
          kw = kw + 1
          wl(kw) = real(iw)/10.
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        end do
  
        ! define wavelength intervals of width 5.0 nm from 365 to 855 nm:
  
        wincr = 5.0
        do iw = 365, 850, 5
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        end do
        wl(kw+1) = wu(kw)
  
        !============================================================
  
        else if (mopt == 2) then   ! low-res
  
        nw = 162
        allocate(wl(nw))
        allocate(wu(nw))
        allocate(wc(nw))
  
        ! define wavelength intervals of width 6.0 nm from 0 to 60 nm:
  
        kw = 0
        wincr = 6.0
        DO iw = 0, 54, 6
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        END DO
  
        ! define wavelength intervals of width 2.0 nm from 60 to 80 nm:
  
        wincr = 2.0
        DO iw = 60, 78, 2
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        END DO
  
        ! define wavelength intervals of width 5.0 nm from 80 to 85 nm:
  
        wincr = 5.0
        DO iw = 80, 80
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        END DO
  
        ! define wavelength intervals of width 2.0 nm from 85 to 117 nm:
  
        wincr = 2.0
        DO iw = 85, 115, 2
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        END DO
  
        ! define wavelength intervals of width 3.0 nm from 117 to 120 nm:
  
        wincr = 3.0
        DO iw = 117, 117
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        END DO
  
        ! define wavelength intervals of width 0.2 nm from 120 to 123 nm:
  
        wincr = 0.2
        DO iw = 1200, 1228, 2
          kw = kw + 1
          wl(kw) = real(iw)/10.
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        ENDDO
  
        ! define wavelength intervals of width 5.0 nm from 123 to 163 nm:
  
        wincr = 5.0
        DO iw = 123, 158, 5
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        ENDDO
  
        ! define wavelength intervals of width 2.0 nm from 163 to 175 nm:
  
        wincr = 2.0
        DO iw = 163, 173, 2
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        ENDDO
  
        ! define wavelength intervals of width 0.5 nm from 175 to 205 nm:
  
        wincr = 0.5
        DO iw = 1750, 2045, 5
          kw = kw + 1
          wl(kw) = real(iw)/10.
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        ENDDO
  
        ! define wavelength intervals of width 5.0 nm from 205 to 245 nm:
  
        wincr = 5.
        DO iw = 205, 240, 5
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        ENDDO
  
        ! define wavelength intervals of width 10.0 nm from 245 to 415 nm:
  
        wincr = 10.0
        DO iw = 245, 405, 10
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        ENDDO
  
        ! define wavelength intervals of width 50.0 nm from 415 to 865 nm:
  
        wincr = 50.0
        DO iw = 415, 815, 50
          kw = kw + 1
          wl(kw) = real(iw)
          wu(kw) = wl(kw) + wincr
          wc(kw) = (wl(kw) + wu(kw))/2.
        ENDDO
  
        wl(kw+1) = wu(kw)
  
        !============================================================

        else if (mopt == 3) then   ! input file
  
           ! look for a " gridwfile= ..." option in def files
           write(*,*) "Input wavelenght grid files for photolysis online is:"
           gridwfile = "gridw.dat" ! default
           call getin_p("gridwfile",gridwfile) ! default path
           write(*,*) " gridwfile = ",trim(gridwfile)
           write(*,*) 'Please use 1 and only 1 skipping line in ',trim(gridwfile)
  
           ! Opening file
           fil = trim(datadir)//'/cross_sections/'//gridwfile
           print*, 'wavelenght grid : ', fil
           OPEN(UNIT=10,FILE=fil,STATUS='old',iostat=ierr)
  
           if (ierr /= 0) THEN
             write(*,*)'Error : cannot open wavelenght grid file ', trim(gridwfile)
             write(*,*)'It should be in :',trim(datadir),'/cross_sections/'
             write(*,*)'1) You can change the directory in callphys.def'
             write(*,*)'   with:'
             write(*,*)'   datadir=/path/to/the/directory'
             write(*,*)'2) You can change the input wavelenght grid file name in'
             write(*,*)'   callphys.def with:'
             write(*,*)'   gridwfile=filename'
             stop
           end if
  
           READ(10,*)
  
           READ(10,*) nw
           allocate(wl(nw))
           allocate(wu(nw))
           allocate(wc(nw))
  
           kw = 0
           READ(10,*) wl(1)
           DO iw = 2, nw
              READ(10,*) wl(iw)
              wu(iw-1) = wl(iw)
              wc(iw-1) = (wl(iw-1) + wu(iw-1))/2.
              kw = kw + 1
           ENDDO
  
           close(10)
        end if  ! mopt
  
        print*, 'number of spectral intervals : ', kw+1
        
    end subroutine gridw
  
    !==============================================================================
   
    subroutine rdsolarflux(nw,wl,wc,fstar1AU)
  
        !     Read and re-grid solar flux data.            
  
        use datafile_mod, only: datadir
        use ioipsl_getin_p_mod, only: getin_p
  
        implicit none
  
        !     input
   
        integer :: nw      ! number of wavelength grid points
        real, dimension(nw) :: wl, wc   ! lower and central wavelength for each interval
  
        !     output
  
        real, dimension(nw) :: fstar1AU  ! stellar flux (photon.s-1.nm-1.cm-2)
  
        !     local 
  
        integer :: iw, nhead, n, i, ierr, kin, naddpnt
  
        real, parameter   :: deltax = 1.e-4
        real, allocatable :: x1(:), y1(:)      ! input solar flux
        real, dimension(nw)    :: yg1          ! gridded solar flux
  
        character(len=200) :: fil
        character(len=150) :: stellarflux
  
        kin     = 10  ! input logical unit
        naddpnt = 4   ! number of adding point
        nhead   = 1   ! number of skipping line at the beggining of the file
  
        ! look for a " stellarflux= ..." option in def files
        write(*,*) "Input stellar spectra files for photolysis online is:"
        stellarflux = "Sun_Claire2012_0.0Ga.txt" ! default
        call getin_p("stellarflux",stellarflux) ! default path
        write(*,*) " stellarflux = ",trim(stellarflux)
        write(*,*) 'Please use ',nhead,' and only ',nhead,' skipping line in ',trim(stellarflux)
  
        ! Opening file
        fil = trim(datadir)//'/stellar_spectra/photochem_stellar_spectra/'//stellarflux
        print*, 'solar flux : ', fil
        OPEN(UNIT=kin,FILE=fil,STATUS='old',iostat=ierr)
  
        if (ierr /= 0) THEN
          write(*,*)'Error : cannot open stellarflux file ', trim(stellarflux)
          write(*,*)'It should be in :',trim(datadir),'/stellar_spectra/photochem_stellar_spectra/'
          write(*,*)'1) You can change the data directory in callphys.def'
          write(*,*)'   with:'
          write(*,*)'   datadir=/path/to/the/directory'
          write(*,*)'2) You can change the input stellarflux file name in'
          write(*,*)'   callphys.def with:'
          write(*,*)'   stellarflux=filename'
          stop
        end if
  
        ! Get number of line in the file
        DO i = 1, nhead
           READ(kin,*)
        ENDDO
        n = 0
        do
          read(kin,*,iostat=ierr)
          if (ierr<0) exit
          n = n + 1
        end do
        rewind(kin)
        DO i = 1, nhead
           READ(kin,*)
        ENDDO
        allocate(x1(n+naddpnt))
        allocate(y1(n+naddpnt))
  
        ! Read stellar flux
        DO i = 1, n
           READ(kin,*) x1(i), y1(i)  !->  x1 nm and y1 in photon.s-1.nm-1.cm-2
        ENDDO
        CLOSE (kin)
  
        ! Interpolation on the grid
        CALL addpnt(x1,y1,n+naddpnt,n,x1(1)*(1.-deltax),0.)
        CALL addpnt(x1,y1,n+naddpnt-1,n,          0.,0.)
        CALL addpnt(x1,y1,n+naddpnt-2,n,x1(n)*(1.+deltax),0.)
        CALL addpnt(x1,y1,n+naddpnt-3,n,      1.e+38,0.)
        CALL inter2(nw,wl,yg1,n,x1,y1,ierr)
        IF (ierr .NE. 0) THEN
           WRITE(*,*) ierr, fil
           STOP
        ENDIF
  
        ! factor to convert to photon.s-1.nm-1.cm-2 :
        ! 5.039e11 = 1.e-4*1e-9/(hc = 6.62e-34*2.998e8)
  
        ! fstar1AU need to be in photon.s-1.nm-1.cm-2
        DO iw = 1, nw-1
           !f(iw) = yg1(iw)*wc(iw)*5.039e11  ! If yg1 in w.m-2.nm-1
           fstar1AU(iw) = yg1(iw)
        ENDDO
  
    end subroutine rdsolarflux
  
    !==============================================================================
  
    subroutine addpnt ( x, y, ld, n, xnew, ynew )
  
    !-----------------------------------------------------------------------------*
    !=  PURPOSE:                                                                 =*
    !=  Add a point <xnew,ynew> to a set of data pairs <x,y>.  x must be in      =*
    !=  ascending order                                                          =*
    !-----------------------------------------------------------------------------*
    !=  PARAMETERS:                                                              =*
    !=  X    - REAL vector of length LD, x-coordinates                       (IO)=*
    !=  Y    - REAL vector of length LD, y-values                            (IO)=*
    !=  LD   - INTEGER, dimension of X, Y exactly as declared in the calling  (I)=*
    !=         program                                                           =*
    !=  N    - INTEGER, number of elements in X, Y.  On entry, it must be:   (IO)=*
    !=         N < LD.  On exit, N is incremented by 1.                          =*
    !=  XNEW - REAL, x-coordinate at which point is to be added               (I)=*
    !=  YNEW - REAL, y-value of point to be added                             (I)=*
    !-----------------------------------------------------------------------------*
  
    IMPLICIT NONE
  
    ! calling parameters
  
    INTEGER ld, n
    REAL x(ld), y(ld)
    REAL xnew, ynew
    INTEGER ierr
  
    ! local variables
  
    INTEGER insert
    INTEGER i
  
    !-----------------------------------------------------------------------
  
    ! initialize error flag
  
    ierr = 0
  
    ! check n<ld to make sure x will hold another point
  
    IF (n .GE. ld) THEN
       WRITE(0,*) '>>> ERROR (ADDPNT) <<<  Cannot expand array '
       WRITE(0,*) '                        All elements used.'
       STOP
    ENDIF
  
    insert = 1
    i = 2
  
    ! check, whether x is already sorted.
    ! also, use this loop to find the point at which xnew needs to be inserted
    ! into vector x, if x is sorted.
  
 10   CONTINUE
      IF (i .LT. n) THEN
        IF (x(i) .LT. x(i-1)) THEN
           print*, x(i-1), x(i)
           WRITE(0,*) '>>> ERROR (ADDPNT) <<<  x-data must be in ascending order!'
           STOP
        ELSE
           IF (xnew .GT. x(i)) insert = i + 1
        ENDIF
        i = i+1
        GOTO 10
      ENDIF
  
    ! if <xnew,ynew> needs to be appended at the end, just do so,
    ! otherwise, insert <xnew,ynew> at position INSERT
  
    IF ( xnew .GT. x(n) ) THEN
   
        x(n+1) = xnew
        y(n+1) = ynew
    
    ELSE
  
    ! shift all existing points one index up
  
        DO i = n, insert, -1
          x(i+1) = x(i)
          y(i+1) = y(i)
        ENDDO
  
    ! insert new point
  
        x(insert) = xnew
        y(insert) = ynew
    
    ENDIF
  
    ! increase total number of elements in x, y
  
    n = n+1
  
    end subroutine addpnt
  
    !==============================================================================
  
    subroutine inter2(ng,xg,yg,n,x,y,ierr)
  
    !-----------------------------------------------------------------------------*
    !=  PURPOSE:                                                                 =*
    !=  Map input data given on single, discrete points onto a set of target     =*
    !=  bins.                                                                    =*
    !=  The original input data are given on single, discrete points of an       =*
    !=  arbitrary grid and are being linearly interpolated onto a specified set  =*
    !=  of target bins.  In general, this is the case for most of the weighting  =*
    !=  functions (action spectra, molecular cross section, and quantum yield    =*
    !=  data), which have to be matched onto the specified wavelength intervals. =*
    !=  The average value in each target bin is found by averaging the trapezoi- =*
    !=  dal area underneath the input data curve (constructed by linearly connec-=*
    !=  ting the discrete input values).                                         =*
    !=  Some caution should be used near the endpoints of the grids.  If the     =*
    !=  input data set does not span the range of the target grid, an error      =*
    !=  message is printed and the execution is stopped, as extrapolation of the =*
    !=  data is not permitted.                                                   =*
    !=  If the input data does not encompass the target grid, use ADDPNT to      =*
    !=  expand the input array.                                                  =*
    !-----------------------------------------------------------------------------*
    !=  PARAMETERS:                                                              =*
    !=  NG  - INTEGER, number of bins + 1 in the target grid                  (I)=*
    !=  XG  - REAL, target grid (e.g., wavelength grid);  bin i is defined    (I)=*
    !=        as [XG(i),XG(i+1)] (i = 1..NG-1)                                   =*
    !=  YG  - REAL, y-data re-gridded onto XG, YG(i) specifies the value for  (O)=*
    !=        bin i (i = 1..NG-1)                                                =*
    !=  N   - INTEGER, number of points in input grid                         (I)=*
    !=  X   - REAL, grid on which input data are defined                      (I)=*
    !=  Y   - REAL, input y-data                                              (I)=*
    !-----------------------------------------------------------------------------*
  
    IMPLICIT NONE
  
    ! input:
    INTEGER ng, n
    REAL x(n), y(n), xg(ng)
  
    ! output:
    REAL yg(ng)
  
    ! local:
    REAL area, xgl, xgu
    REAL darea, slope
    REAL a1, a2, b1, b2
    INTEGER ngintv
    INTEGER i, k, jstart
    INTEGER ierr
    !_______________________________________________________________________
  
    ierr = 0
  
    !  test for correct ordering of data, by increasing value of x
  
    DO 10, i = 2, n
       IF (x(i) .LE. x(i-1)) THEN
          ierr = 1
          WRITE(*,*)'data not sorted'
          WRITE(*,*) x(i), x(i-1)
          RETURN
       ENDIF
 10 CONTINUE     

    DO i = 2, ng
      IF (xg(i) .LE. xg(i-1)) THEN
         ierr = 2
        WRITE(0,*) '>>> ERROR (inter2) <<<  xg-grid not sorted!'
        RETURN
      ENDIF
    ENDDO
  
    ! check for xg-values outside the x-range
  
    IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN
        WRITE(0,*) '>>> ERROR (inter2) <<<  Data do not span grid.  '
        WRITE(0,*) '                        Use ADDPNT to expand data and re-run.'
        STOP
    ENDIF
  
    !  find the integral of each grid interval and use this to 
    !  calculate the average y value for the interval      
    !  xgl and xgu are the lower and upper limits of the grid interval
  
    jstart = 1
    ngintv = ng - 1
    DO 50, i = 1,ngintv
  
    ! initialize:
  
        area = 0.0
        xgl = xg(i)
        xgu = xg(i+1)
  
    !  discard data before the first grid interval and after the 
    !  last grid interval
    !  for internal grid intervals, start calculating area by interpolating
    !  between the last point which lies in the previous interval and the
    !  first point inside the current interval
  
        k = jstart
        IF (k .LE. n-1) THEN
  
    !  if both points are before the first grid, go to the next point
 30         CONTINUE
              IF (x(k+1) .LE. xgl) THEN
                 jstart = k - 1
                 k = k+1
                 IF (k .LE. n-1) GO TO 30
              ENDIF
  
  
    !  if the last point is beyond the end of the grid, complete and go to the next
    !  grid
 40         CONTINUE
               IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN          

                  jstart = k-1
  
    ! compute x-coordinates of increment
  
                  a1 = MAX(x(k),xgl)
                  a2 = MIN(x(k+1),xgu)
  
    ! if points coincide, contribution is zero
  
                  IF (x(k+1).EQ.x(k)) THEN
                     darea = 0.e0
                  ELSE
                     slope = (y(k+1) - y(k))/(x(k+1) - x(k))
                     b1 = y(k) + slope*(a1 - x(k))
                     b2 = y(k) + slope*(a2 - x(k))
                     darea = (a2 - a1)*(b2 + b1)/2.
                  ENDIF
  
    !  find the area under the trapezoid from a1 to a2
  
                  area = area + darea
  
    ! go to next point
                
                  k = k+1
                  GO TO 40

              ENDIF
          ENDIF
  
    !  calculate the average y after summing the areas in the interval
  
            yg(i) = area/(xgu - xgl)
  
 50 CONTINUE
  
    end subroutine inter2
  
    !==============================================================================

    subroutine rdxsi(nw,wl,xsi,jparamlinei,xs_tempi,qyi,tdimi,species)
  
    !-----------------------------------------------------------------------------*
    !=  PURPOSE:                                                                 =*
    !=  Read molecular absorption cross section.  Re-grid data to match          =*
    !=  specified wavelength working grid.                                       =*
    !-----------------------------------------------------------------------------*
    !=  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                                         =*
    !=  XS     - REAL, molecular absoprtion cross section (cm^2) at each      (O)=*
    !=           specified wavelength                                            =*
    !-----------------------------------------------------------------------------*
  
    use datafile_mod, only: datadir
  
    IMPLICIT NONE
  
    !     input
  
    integer, intent(in) :: nw       ! number of wavelength grid points
    integer, intent(in) :: tdimi    ! temperature dimension
    real,    intent(in) :: wl(nw)   ! lower and central wavelength for each interval (nm)
    character(len = 20) :: species  ! species which is photolysed
    character(len = 300), intent(in) :: jparamlinei ! line of jonline parameters
  
    !     output
  
    real, intent(inout) :: qyi(nw)         ! photodissociation yield
    real, intent(inout) :: xsi(tdimi,nw)   ! absorption cross-section (cm2)
    real, intent(inout) :: xs_tempi(tdimi) !  absorption cross-section temperature (K)
  
    !     local
  
    character(len = 50)  :: sigfile(tdimi) ! cross section file name
    character(len = 50)  :: qyfile         ! quantum yield file name
    CHARACTER(len = 120) :: fil            ! path files
    integer :: tdimdummy, ifile, itemp, ierr, nheadxs, nheadqy, i, n, naddpnt
    integer :: kin                         ! input logical unit
    real, allocatable :: wlf(:), xsf(:)    ! input cross section
    real, allocatable :: wlqy(:), qyf(:)   ! input photodissociation yield
    real, parameter   :: deltax = 1.e-4
  
    kin     = 10     ! input logical unit
    naddpnt = 4      ! number of adding point (careful: same for xs and qy)
    nheadxs = 1      ! number of skipping line at the beggining of the cross section files
    nheadqy = 1      ! number of skipping line at the beggining of the quantum yield files
    sigfile(:) = ''
    qyfile = ''
  
    read(jparamlinei,*) tdimdummy, (xs_tempi(itemp), itemp=1,tdimi), (sigfile(ifile), ifile=1,tdimi), qyfile
  
    ! Check if xs_tempi is sorted from low to high values
    do i = 1,tdimi-1
      if (xs_tempi(i)>xs_tempi(i+1)) then
        print*, 'ERROR: temperature cross section file'
        print*, '       has to be sorted from low to high values'
        print*, '       Check reactfile file'
        stop
      end if
    end do
  
    ! Cross section
    do itemp=1,tdimi

      fil = trim(datadir)//'/cross_sections/'//trim(sigfile(itemp))
      print*, 'section efficace '//trim(species)//': ', fil

      OPEN(UNIT=kin,FILE=fil,STATUS='old',iostat=ierr)

      if (ierr /= 0) THEN
         write(*,*)'Error : cannot open cross_sections file ', trim(sigfile(itemp))
         write(*,*)'It should be in :',trim(datadir),'/cross_sections/'
         write(*,*)'1) You can change the datadir directory in callphys.def'
         write(*,*)'   with:'
         write(*,*)'   datadir=/path/to/the/directory'
         write(*,*)'2) You can check if the file is in datadir/cross_sections/'
         write(*,*)'3) You can change the input cross_sections file name in'
         write(*,*)'   chemestry/reactfile file for the specie:',trim(species)
         stop
      end if

      ! Get number of line in the file
      DO i = 1, nheadxs
         READ(kin,*)
      ENDDO
      n = 0
      do
        read(kin,*,iostat=ierr)
        if (ierr<0) exit
        n = n + 1
      end do
      rewind(kin)
      DO i = 1, nheadxs
         READ(kin,*)
      ENDDO

      allocate(wlf(n+naddpnt))
      allocate(xsf(n+naddpnt))

      ! Read cross section
      DO i = 1, n
         READ(kin,*) wlf(i), xsf(i)  !-> xsf in cm2 and wlf in nm
      ENDDO
      CLOSE (kin)

      CALL addpnt(wlf,xsf,n+naddpnt,n,wlf(1)*(1.-deltax),0.)
      CALL addpnt(wlf,xsf,n+naddpnt-1,n,               0.,0.)
      CALL addpnt(wlf,xsf,n+naddpnt-2,n,wlf(n)*(1.+deltax),0.)
      CALL addpnt(wlf,xsf,n+naddpnt-3,n,           1.e+38,0.)
      CALL inter2(nw,wl,xsi(itemp,:),n,wlf,xsf,ierr)
      IF (ierr .NE. 0) THEN
         WRITE(*,*) ierr, fil
         STOP
      ENDIF

      deallocate(wlf)
      deallocate(xsf)

    end do

    ! Photodissociation yield
    fil = trim(datadir)//'/cross_sections/'//trim(qyfile)
    print*, 'photodissociation yield '//trim(species)//': ', fil

    OPEN(UNIT=kin,FILE=fil,STATUS='old',iostat=ierr)

    if (ierr /= 0) THEN
       write(*,*)'Error : cannot open photodissociation yield file ', trim(qyfile)
       write(*,*)'It should be in :',trim(datadir),'/cross_sections/'
       write(*,*)'1) You can change the datadir directory in callphys.def'
       write(*,*)'   with:'
       write(*,*)'   datadir=/path/to/the/directory'
       write(*,*)'2) You can check if the file is in datadir/cross_sections/'
       write(*,*)'3) You can change the input photodissociation yield file name in'
       write(*,*)'   chemestry/reactfile file for the specie:',trim(species)
       stop
    end if

    ! Get number of line in the file
    DO i = 1, nheadqy
       READ(kin,*)
    ENDDO
    n = 0
    do
      read(kin,*,iostat=ierr)
      if (ierr<0) exit
      n = n + 1
    end do
    rewind(kin)
    DO i = 1, nheadqy
       READ(kin,*)
    ENDDO
    allocate(wlqy(n+naddpnt))
    allocate(qyf(n+naddpnt))

    ! Read photodissociation yield
    DO i = 1, n
       READ(kin,*) wlqy(i), qyf(i)  !-> wlqy in nm
    ENDDO
    CLOSE (kin)

    CALL addpnt(wlqy,qyf,n+naddpnt,n,wlqy(1)*(1.-deltax),0.)
    CALL addpnt(wlqy,qyf,n+naddpnt-1,n,                 0.,0.)
    CALL addpnt(wlqy,qyf,n+naddpnt-2,n,wlqy(n)*(1.+deltax),0.)
    CALL addpnt(wlqy,qyf,n+naddpnt-3,n,             1.e+38,0.)
    CALL inter2(nw,wl,qyi,n,wlqy,qyf,ierr)
    IF (ierr .NE. 0) THEN
       WRITE(*,*) ierr, fil
       STOP
    ENDIF

    end subroutine rdxsi
   
end module photolysis_mod