! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later

  !*****************************************************************************
  !                                                                            *
  !    This module contains variables and subroutines for  calculating         *
  !    chemical loss of species                                                *
  !                                                                            *
  !*****************************************************************************

module chemistry_mod

  use netcdf
  use par_mod
  use com_mod
  use date_mod
  use windfields_mod,    only: rho,nxmax,nymax,nzmax
  use totals_mod,        only: chem_loss
  use netcdf_output_mod, only: nf90_err

  ! reagent field variables 

  implicit none

  integer, allocatable, dimension(:)      :: nxCL, nyCL, nzCL
  integer                                 :: nxjr, nyjr
  real, allocatable, dimension(:,:)       :: lonCL, latCL, altCL
  real, allocatable, dimension(:)         :: dxCL, dyCL
  real                                    :: dxjr, dyjr
  real, allocatable, dimension(:,:,:,:,:) :: CL_field           ! chemical loss fields at input resolution
  real, allocatable, dimension(:,:,:,:)   :: clfield_cur        ! chemical loss fields current hour
  integer, dimension(2)                   :: memCLtime          ! time of fields in memory (sec)  
  integer                                 :: curCLhour          ! current hour since start of simulation
  real(kind=dp), dimension(2)             :: CL_time            ! julian date of fields in memory
  real, allocatable, dimension(:,:,:)     :: jrate_average      ! monthly average actinic flux 
  real, allocatable, dimension(:)         :: lonjr, latjr

  private :: zenithangle, photo_O1D

  contains

  subroutine readreagents()

  !*****************************************************************************
  !                                                                            *
  !   This routine reads names and input paths of chemical reagents            *
  !                                                                            *
  !   Author: R. Thompson, Sep 2023                                            *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! preagent              Reagent name                                         *
  ! preag_path            Path to reagent fields                               *
  ! phourly               Logical for hourly interpolate (0 = no, 1 = yes)     *
  !                                                                            *
  !*****************************************************************************

    implicit none

    character(len=256) :: preag_path
    character(len=16)  :: preagent, punit
    integer :: phourly
    integer,parameter :: unitreagents=2, unitreagentsout=3
    integer :: readerror
    integer :: j

    ! declare namelist
    namelist /reagent_params/ &
         preagent, preag_path, phourly, punit

    preagent="" ! read failure indicator value
    preag_path=""
    punit=""
    phourly=0
    reag_hourly(:)=0

    open(unitreagents,file=path(1)(1:length(1))//'REAGENTS',status='old',form='formatted',iostat=readerror)
    if (readerror.ne.0) then
      ! no REAGENTS file
      nreagent=0
      go to 999
    endif

    ! prepare namelist output if requested
    if (nmlout.and.lroot) then
      open(unitreagentsout,file=path(2)(1:length(2))//'REAGENTS.namelist',access='append',status='replace',iostat=readerror)
      if (readerror.ne.0) then
        write(*,*) '#### FLEXPART MODEL ERROR CANNOT CREATE FILE:   ####'
        write(*,*) '#### ',path(2)(1:length(2))//'REAGENTS.namelist ####'
      endif
    endif

    ! try namelist input
    read(unitreagents,reagent_params,iostat=readerror)
    close(unitreagents)
    if (readerror.ne.0) then
      ! not namelist format
      nreagent=0
      go to 999
    endif

    ! namelist input
    open(unitreagents,file=path(1)(1:length(1))//'REAGENTS',status='old',form='formatted',iostat=readerror)
    j=0
    do while (readerror.eq.0)
      j=j+1
      if (j.gt.maxreagent) then
        write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY REAGENTS #### '
        write(*,*) ' #### MAXIMUM NUMBER IS ',maxreagent,'        #### '
        write(*,*) ' #### PLEASE MAKE CHANGES IN FILE REAGENTS    #### '
        stop
      endif
      read(unitreagents,reagent_params,iostat=readerror)
      reagents(j)=preagent
      reag_path(j)=preag_path
      reag_hourly(j)=phourly
      reag_unit(j)=punit
      ! namelist output
      if (nmlout.and.lroot) then
        write(unitreagentsout,nml=reagent_params)
      endif
    end do
    nreagent=j-1
    if (lroot) then
      write(*,*) 'Number of reagents: ',nreagent
      write(*,*) 'Reagent names: ',reagents(1:nreagent)
      write(*,*) 'Reagent units: ',reag_unit(1:nreagent)
    endif
    close(unitreagents)
    close(unitreagentsout)

999 continue ! no reagents file

  end subroutine readreagents

  subroutine getchemfield(itime)

  !*****************************************************************************
  !                                                                            *
  !    This routine reads the chemical reagent fields into memory              *
  !                                                                            *
  !    Author: Rona Thompson, Sep 2023                                         *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! itime            time since start of simulation in sec                     *
  !                                                                            *
  !*****************************************************************************

    implicit none

    integer          :: itime
    integer          :: jjjjmmdd, hhmmss, mm, eomday
    integer          :: nr, memid
    real(kind=dp)    :: jd, jdmid
    character(len=2) :: amonth
    character(len=256):: CL_name, jr_name
    logical          :: lexist

    print*, 'getchemfield: ldirect*memCLtime(1) = ',ldirect*memCLtime(1)
    print*, 'getchemfield: ldirect*memCLtime(2) = ',ldirect*memCLtime(2)

    ! Check fields are available for the current time step
    !*****************************************************

    if ((ldirect*memCLtime(1).le.ldirect*itime).and. &
         (ldirect*memCLtime(2).gt.ldirect*itime)) then

      ! The rightfields are already in memory -> don't do anything
      return 

    else if ((ldirect*memCLtime(2).le.ldirect*itime).and. &
         (memCLtime(2).ne.0.)) then

      ! Current time is after 2nd chem field
      !*************************************

      write(*,*) 'Updating CL fields... '

      memCLtime(1)=memCLtime(2)  ! time in sec
      CL_time(1)=CL_time(2)      ! julian date
      memid=2

      ! Get date/time of next chem field 
      !*********************************
      ! assumes fields are monthly
     
      call caldate(CL_time(1), jjjjmmdd, hhmmss)
      eomday=calceomday(jjjjmmdd/100)
      memCLtime(2)=memCLtime(1)+ldirect*eomday*24*3600   ! time in sec
      CL_time(2)=CL_time(1)+real(ldirect*eomday,kind=dp) ! julian date 
      write(*,*) 'getchemfield: memCLtime = ',memCLtime(1), memCLtime(2)
      write(*,*) 'getchemfield: CL_time = ',CL_time(1), CL_time(2)
      call caldate(CL_time(2), jjjjmmdd,hhmmss)
      mm=int((jjjjmmdd-(jjjjmmdd/10000)*10000)/100)
      write(amonth,'(I2.2)') mm

      do nr=1,nreagent

        write(*,*) 'Updating CL field for: ',trim(reagents(nr))

        CL_field(:,:,:,nr,1)=CL_field(:,:,:,nr,2)

        ! Read new chem field and store in 2nd position
        !**********************************************

        write(CL_name,'(A)') trim(reag_path(nr))//trim(reagents(nr))//'_'//amonth//'.nc'
        inquire(file=CL_name,exist=lexist)
        if (lexist) then
          write(*,*) 'Reading CL field: ',CL_name
          call readchemfield(CL_name, memid, nr)
        else
          write(*,*) '#### FLEXPART ERROR                ####'
          write(*,*) '#### CHEMISTRY FIELD NOT FOUND     ####'
          write(*,*) '#### '//trim(CL_name)//' ####'
          error stop
        endif

        if (reag_hourly(nr).gt.0) then

          ! Read average jrates and store in 2nd position
          !**********************************************

          write(jr_name,'(A)') trim(reag_path(nr))//'jrate_average.nc'
          inquire(file=jr_name,exist=lexist)
          if (lexist) then
            write(*,*) 'Reading jrate field: ',jr_name
            call readjrate(jr_name, memid, mm)
          else
            write(*,*) '#### FLEXPART ERROR                ####'
            write(*,*) '#### JRATE_AVERAGE NOT FOUND       ####'
            write(*,*) '#### '//trim(jr_name)//' ####'
            error stop
          endif

        endif

      end do ! nreagent

    else

      ! No chem field in memory that can be used
      !******************************************  

      write(*,*) 'Reading two new CL fields...'

      ! Get date/time of both chem fields 
      !**********************************
      ! assumes fields are monthly

      do memid=1,2
        if (memid.eq.1) then
          jd=bdate+real(ldirect*itime,kind=dp)/86400._dp
          call caldate(jd, jjjjmmdd, hhmmss)
          ! middle of month day
          jdmid=juldate(int(jjjjmmdd/100)*100+15,0) 
          !! testing
          print*, 'getchemfield: jjjjmmdd, jjjjmmdd_mid = ',jjjjmmdd, int(jjjjmmdd/100)*100+15
          ! julian date of fields
          if (ldirect.gt.0) then
            if (jd.ge.jdmid) then
              ! use current month
              CL_time(memid)=jdmid 
            else
              ! use last month
              eomday=calceomday(jjjjmmdd/100)
              CL_time(memid)=jdmid-real(eomday,kind=dp)
            endif
          else
            if (jd.ge.jdmid) then
              ! use next month
              CL_time(memid)=jdmid+real(eomday,kind=dp)
            else
              ! use current month
              CL_time(memid)=jdmid
            endif
          endif
          !! testing
          print*, 'getchemfield: memid, jd, jdmid, CL_time = ',memid,jd,jdmid,CL_time(1)
        else
          call caldate(jd, jjjjmmdd, hhmmss)
          eomday=calceomday(jjjjmmdd/100)
          CL_time(memid)=CL_time(memid-1)+real(ldirect*eomday,kind=dp)
        endif
        ! time of field in seconds from start
        memCLtime(memid)=int((CL_time(memid)-bdate)*86400._dp)

        write(*,*) 'getchemfield: memid, memCLtime = ',memCLtime(memid)
        write(*,*) 'getchemfield: memid, CL_time = ',CL_time(memid)

        call caldate(CL_time(memid), jjjjmmdd, hhmmss)
        mm=int((jjjjmmdd-(jjjjmmdd/10000)*10000)/100)
        write(amonth,'(I2.2)') mm

        do nr=1,nreagent

          write(*,*) 'Reading two new CL fields for: ',trim(reagents(nr))

          ! Read new chem field
          !********************

          write(CL_name,'(A)') trim(reag_path(nr))//trim(reagents(nr))//'_'//amonth//'.nc'
          inquire(file=CL_name,exist=lexist)
          if (lexist) then
            write(*,*) 'Reading CL field: ',CL_name
            call readchemfield(CL_name, memid, nr)
          else
            write(*,*) '#### FLEXPART ERROR                ####'
            write(*,*) '#### CHEMISTRY FIELD NOT FOUND     ####'
            write(*,*) '#### '//trim(CL_name)//' ####'
            error stop
          endif
 
          ! Read average jrates
          !********************

          if (reag_hourly(nr).gt.0) then
            ! Read average jrates into memory 
            write(jr_name,'(A)') trim(reag_path(nr))//'jrate_average.nc'
            inquire(file=jr_name,exist=lexist)
            if (lexist) then
              write(*,*) 'Reading jrate field: ',jr_name
              call readjrate(jr_name, memid, mm)
            else
              write(*,*) '#### FLEXPART ERROR                ####'
              write(*,*) '#### JRATE_AVERAGE NOT FOUND       ####'
              write(*,*) '#### '//trim(jr_name)//' ####'
              error stop
            endif
          endif ! reag_hourly

        end do ! nreagent

      end do ! memid
      
    endif ! update hourly fields


  end subroutine getchemfield


  subroutine readchemfield(CL_name,memid,nr)

  !*****************************************************************************
  !                                                                            *
  !    Reads the chemical reagent fields into memory                           *
  !                                                                            *
  !    Author: Rona Thompson, Sep 2023                                         *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! CL_name           name of chemical reagent file                            *
  ! memid             time index to chemical field variable                    *
  ! nr                reagent index to chemical field variable                 *
  !                                                                            * 
  !*****************************************************************************

    implicit none

    character(len=256) :: CL_name
    integer            :: memid,nr,len
    integer            :: ncid,dimid,varid

    ! Read netcdf file
    !******************

    ! open file
    call nf90_err( nf90_open(trim(CL_name),nf90_NOWRITE,ncid) )

    ! longitude
    call nf90_err( nf90_inq_dimid(ncid,'lon',dimid) )
    call nf90_err( nf90_inquire_dimension(ncid,dimid,len=len) )
    if (.not.allocated(lonCL)) then
      allocate(nxCL(nreagent))
      allocate(dxCL(nreagent))
      allocate(lonCL(nxmax,nreagent))
    endif
    nxCL(nr)=len
    call nf90_err( nf90_inq_varid(ncid,'lon',varid) )
    call nf90_err( nf90_get_var(ncid,varid,lonCL(1:nxCL(nr),nr)) )
    dxCL(nr)=abs(lonCL(2,nr)-lonCL(1,nr))

    ! latitude
    call nf90_err( nf90_inq_dimid(ncid,'lat',dimid) )
    call nf90_err( nf90_inquire_dimension(ncid,dimid,len=len) )
    if (.not.allocated(latCL)) then
      allocate(nyCL(nreagent))
      allocate(dyCL(nreagent))
      allocate(latCL(nymax,nreagent))
    endif
    nyCL(nr)=len
    call nf90_err( nf90_inq_varid(ncid,'lat',varid) )
    call nf90_err( nf90_get_var(ncid,varid,latCL(1:nyCL(nr),nr)) )
    dyCL(nr)=abs(latCL(2,nr)-latCL(1,nr))

    ! altitude
    call nf90_err( nf90_inq_dimid(ncid,'lev',dimid) )
    call nf90_err( nf90_inquire_dimension(ncid,dimid,len=len) )
    if (.not.allocated(altCL)) then
      allocate(nzCL(nreagent))
      allocate(altCL(nzmax,nreagent))
    endif
    nzCL(nr)=len
    call nf90_err( nf90_inq_varid(ncid,'lev',varid) )
    call nf90_err( nf90_get_var(ncid,varid,altCL(1:nzCL(nr),nr)) )

    ! chemical field
    call nf90_err( nf90_inq_varid(ncid,trim(reagents(nr)),varid) )
    if (.not.allocated(CL_field)) then
      allocate(CL_field(nxmax,nymax,nzmax,nreagent,2))
      allocate(clfield_cur(nxmax,nymax,nzmax,nreagent))
    endif
    call nf90_err( nf90_get_var(ncid,varid,CL_field(1:nxCL(nr),1:nyCL(nr),1:nzCL(nr),nr,memid)) )

    ! close file
    call nf90_err( nf90_close(ncid) )

    !! testing
!    print*, 'readchemfield: nxCL, nyCL, nzCL = ',nxCL(nr), nyCL(nr), nzCL(nr)
!    print*, 'readchemfield: lonCL = ',lonCL(1:nxCL(nr),nr)
!    print*, 'readchemfield: latCL = ',latCL(1:nyCL(nr),nr)

    return

  end subroutine readchemfield

  subroutine getchemhourly(itime)

  !*****************************************************************************
  !                                                                            *
  !    This routine interpolates the chemistry fields to current hour and      *
  !    if required using information on solar zenith angle                     *
  !                                                                            *
  !    Author: Rona Thompson, Mar 2024                                         *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! itime [s]            actual simulation time [s]                            *
  !                                                                            *
  !*****************************************************************************

    use point_mod,      only: xlon0, ylat0, dx, dy
    use windfields_mod, only : height, nz

    implicit none

    integer            :: itime, curhour, interp_time
    integer            :: nr, kz, ix, jy, ixm, jym, ixp, jyp, indz, indzm, ii
    real               :: dt1, dt2, dtt, sza, jrate, jrate_cur, r
    real, dimension(2) :: r1
    real               :: dz1, dz2, dz, ddx, ddy, rddx, rddy, p1, p2, p3, p4
    integer            :: jjjjmmdd, hhmmss
    integer            :: jrx, jry
    real(kind=dp)      :: jul
    real, parameter    :: avog = 6.02214e23 ! Avogadro constant (1/mol)
    !! testing
!    character(len=4) :: atime
!    character(len=20) :: frmt
!    real, dimension(nxjr,nyjr) :: jscalar, sza_grid
    
    !! testing
!    jscalar(:,:)=1.
!    sza_grid(:,:)=0.

    ! current hour of simulation
    curhour=itime/3600
    write(*,*) 'getchemhourly: curhour, curCLhour = ',curhour, curCLhour

    jul=bdate+real(itime,kind=dp)/86400._dp
    call caldate(jul,jjjjmmdd,hhmmss)

    if ((ldirect*curCLhour.eq.ldirect*curhour).and.(ldirect*itime.gt.0)) then
      ! chemical field is for current hour -> don't do anything
      return
    else
      ! interpolate to middle of hour
      curCLhour=curhour
      interp_time=curhour*3600+1800
      dt1=float(interp_time-memCLtime(1))
      dt2=float(memCLtime(2)-interp_time)
      dtt=1./(dt1+dt2)
      !! testing
      print*, 'getchemhourly: dt1, dt2, dtt = ',dt1,dt2,dtt 
      do nr=1,nreagent
        write(*,*) 'Interpolating fields for reagent: ',reagents(nr)
        clfield_cur(:,:,:,nr)=(dt2*CL_field(:,:,:,nr,1) + dt1*CL_field(:,:,:,nr,2))*dtt
        if (reag_unit(nr).eq.'mol/mol') then
          ! convert to molecule/cm3
!$OMP PARALLEL &
!$OMP PRIVATE(kz,ix,jy,indz,indzm,ixm,jym,ixp,jyp,ddx,ddy,rddx,rddy,&
!$OMP p1,p2,p3,p4,r1,ii,dz1,dz2,dz,r)
!$OMP DO
          do kz=1,nzCL(nr)
            ! assume chem fields vertical coordinate is in meters
            indzm=nz-1
            do indz=1,nz
              if (height(indz).gt.altCL(kz,nr)) then
                indzm=indz-1
                exit
              endif
            end do
            dz1=altCL(kz,nr)-height(indzm)
            dz2=height(indz)-altCL(kz,nr)
            dz=1./(dz1+dz2)
            do jy=1,nyCL(nr)
              jym=int((latCL(jy,nr)-ylat0)/dy)
              jyp=jym+1
              ddy=(latCL(jy,nr)-ylat0)/dy-real(jym)
              rddy=1.-ddy
              do ix=1,nxCL(nr)
                ixm=int((lonCL(ix,nr)-xlon0)/dx)
                ixp=ixm+1
                ddx=(lonCL(ix,nr)-xlon0)/dx-real(ixm)
                rddx=1.-ddx
                p1=rddx*rddy
                p2=ddx*rddy
                p3=rddx*ddy
                p4=ddx*ddy
                !! testing
!                if ((ix.lt.5).and.(jy.gt.10.and.jy.lt.20).and.kz.eq.1) then
!                  print*, 'getchemhourly: lonCL, xlon, latCL, ylat = ',lonCL(ix,nr), xlon0+ixm*dx, latCL(jy,nr), ylat0+jyp*dy
!                  print*, 'getchemhourly: altCL, height = ',altCL(kz,nr), (dz2*height(indzm)+dz1*height(indzm+1))*dz
!                endif
                ! take density from first field (accurate enough)
                do ii=1,2
                  indz=indzm+ii-1
                  r1(ii)=p1*rho(ixm,jym,indz,1)+&
                         p2*rho(ixp,jym,indz,1)+&
                         p3*rho(ixm,jyp,indz,1)+&
                         p4*rho(ixp,jyp,indz,1)
                end do
                r=(dz2*r1(1)+dz1*r1(2))*dz
                ! vmr*Avog*P/(RT)/1e6
                ! using P/T = rho*rair
                clfield_cur(ix,jy,kz,nr)=clfield_cur(ix,jy,kz,nr)*avog*r*r_air/rgas/1.e6 
              end do
            end do
          end do
!$OMP END DO
!$OMP END PARALLEL
        endif
        if (reag_hourly(nr).gt.0) then
          ! use actinic flux (jrate) for interpolation
!$OMP PARALLEL &
!$OMP PRIVATE(kz,jy,jry,ix,jrx,sza,jrate,jrate_cur) 
!$OMP DO
          do kz=1,nzCL(nr)
            do jy=1,nyCL(nr)
              ! jrate_average dimensions given as grid midpoints
              jry=int((latCL(jy,nr)-(latjr(1)-0.5*dyjr))/dyjr)+1
              jry=min(jry,nyCL(nr))
              !! testing
!              if (kz.eq.1.and.jy.lt.10) print*, 'latCL, latjr = ',latCL(jy), latjr(jry)
              do ix=1,nxCL(nr)
                ! jrate_average dimensions given as grid midpoints
                jrx=int((lonCL(ix,nr)-(lonjr(1)-0.5*dxjr))/dxjr)+1
                jrx=min(jrx,nxCL(nr))
                !! testing
!                if (kz.eq.1.and.jy.lt.10.and.ix.lt.10) print*, 'lonCL, lonjr = ',lonCL(ix), lonjr(jrx)
                ! solar zenith angle
                sza=zenithangle(latjr(jry),lonjr(jrx),jjjjmmdd,hhmmss)
                ! calculate J(O1D) (jrate)
                jrate=photo_O1D(sza)
                jrate_cur=(dt2*jrate_average(jrx,jry,1) + &
                           dt1*jrate_average(jrx,jry,2))*dtt
                ! apply hourly correction to chem field
                if(jrate_cur.gt.0.) then
                  clfield_cur(ix,jy,kz,nr)=clfield_cur(ix,jy,kz,nr)*jrate/jrate_cur
                  !! testing
!                  if (kz.eq.1) then
!                    jscalar(ix,jy)=jrate/jrate_cur
!                    sza_grid(ix,jy)=sza
!                  endif
                endif
              end do
            end do
          end do
!$OMP END DO
!$OMP END PARALLEL
        endif ! reag_hourly
      end do ! nreagent
    endif ! curhour

    !! testing
!    write(frmt,fmt='(A,I4,A)') '(',nxCL,'(E14.6))'
!    write(atime,fmt='(I4.4)') curhour
!    open(600,file=path(2)(1:length(2))//'clfield_'//atime//'.txt',action='write',status='replace')
!    do kz=1,nzCL
!      do jy=1,nyCL
!        write(600,fmt=frmt) clfield_cur(:,jy,kz,1)
!      end do
!    end do
!    close(600)
!    write(frmt,fmt='(A,I4,A)') '(',nxjr,'(E14.6))'
!    open(600,file=path(2)(1:length(2))//'jscalar_'//atime//'.txt',action='write',status='replace')
!    do jy=1,nyjr
!      write(600,fmt=frmt) jscalar(:,jy)
!    end do
!    close(600)
!    write(frmt,fmt='(A,I4,A)') '(',nxjr,'(E14.6))'
!    open(600,file=path(2)(1:length(2))//'sza_'//atime//'.txt',action='write',status='replace')
!    do jy=1,nyjr
!      write(600,fmt=frmt) sza_grid(:,jy)
!    end do
!    close(600)

  end subroutine getchemhourly

  subroutine chemreaction(itime)

  !*****************************************************************************
  !                                                                            *
  !    This routine computes the chemical loss of each species and updates     *
  !    the particle mass                                                       *
  !                                                                            *
  !    Author: Rona Thompson, Sep 2023                                         *
  !    updated for v11, Jan 2024                                               *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! itime [s]            actual simulation time [s]                            *
  !                                                                            *
  !*****************************************************************************

    use particle_mod,   only: count, part, mass
    use point_mod,      only: xlon0, ylat0, dx, dy
    use windfields_mod, only: xresoln, yresoln, xln, yln, xrn, yrn, height, tt, nz
    use omp_lib

    implicit none

    integer               :: jpart,itime
    integer               :: ii,i,j,ks,ix,jy,kz,jrx,jry,nr
    integer               :: ngrid,interp_time,n,indz
    integer               :: jjjjmmdd,hhmmss,mm,hh,m1,m2
    integer               :: clx,cly,clz,clzm,ithread
    real, dimension(nzmax) :: altCLtop
    real, dimension(2)    :: cl_tmp
    real                  :: xlon,ylat
    real                  :: xtn,ytn
    real                  :: dt1,dt2,dtt,dtt1,dtt2,dttt,dz1,dz2,dzz
    real                  :: restmass,clreacted,cl_cur
    real                  :: jrate,jrate_cur,sza
    real                  :: clrate,temp
    real, parameter       :: smallnum = tiny(0.0) ! smallest number that can be handled
    real(kind=dp)         :: jul
    real                  :: lonjrx,latjry
#ifdef _OPENMP
    real(kind=dp), allocatable, dimension(:,:,:) :: chem_loss_tmp
#endif

    ! use middle of synchronisation time step
    interp_time=nint(itime+0.5*lsynctime)
    dtt1=float(interp_time-memtime(1))
    dtt2=float(memtime(2)-interp_time)
    dttt=1./(dtt1+dtt2)

    ! initialization
    chem_loss(:,:)=0d0

#ifdef _OPENMP
    allocate( chem_loss_tmp(nreagent,nspec,numthreads) )
    chem_loss_tmp(:,:,:) = 0d0
#endif

    ! Loop over particles
    !*****************************************

!$OMP PARALLEL &
!$OMP PRIVATE(ii,jpart,ngrid,j,xtn,ytn,ix,jy, &
!$OMP   xlon,ylat,clx,cly,clz,clzm,kz,altCLtop,dz1,dz2,dzz,nr,i, &
!$OMP   cl_cur,indz,temp,ks,clrate,restmass,clreacted,ithread) 

#ifdef _OPENMP
    ithread = OMP_GET_THREAD_NUM()+1 ! Starts with 1
#else
    ithread = 1
#endif

!$OMP DO
    do ii=1,count%alive

      jpart=count%ialive(ii)
      ! Determine which nesting level to be used
      ngrid=0
      do j=numbnests,1,-1 ! check if need +/- dxn below
        if ((part(jpart)%xlon.gt.xln(j)).and.(part(jpart)%xlon.lt.xrn(j)).and. &
             (part(jpart)%ylat.gt.yln(j)).and.(part(jpart)%ylat.lt.yrn(j))) then
          ngrid=j
          exit
        endif
      end do

      ! Determine nested grid coordinates
      if (ngrid.gt.0) then
        xtn=(real(part(jpart)%xlon)-xln(ngrid))*xresoln(ngrid)
        ytn=(real(part(jpart)%ylat)-yln(ngrid))*yresoln(ngrid)
        ix=int(xtn)
        jy=int(ytn)
      else
        ix=int(part(jpart)%xlon)
        jy=int(part(jpart)%ylat)
      endif

      ! Get CL from nearest grid-cell for current hour
      !***********************************************

      ! world coordinates
      xlon=real(part(jpart)%xlon)*dx+xlon0
      if (xlon.gt.180.) then
        xlon=xlon-360.
      endif
      ylat=real(part(jpart)%ylat)*dy+ylat0
      if (ylat.gt.90.) then
        ylat=90.
      endif

      do nr=1,nreagent

        ! get position in the chem field
        ! assumes chem field dimensions given as grid midpoints
        clx=min(nxCL(nr),int((xlon-(lonCL(1,nr)-0.5*dxCL(nr)))/dxCL(nr))+1)
        cly=min(nyCL(nr),int((ylat-(latCL(1,nr)-0.5*dyCL(nr)))/dyCL(nr))+1)   

        ! get the level of the chem field for the particle
        ! z is the z-coord of the trajectory above model orography in metres
        ! altCL is the height of the centre of the level in the chem field above orography
        do kz=2,nzCL(nr)
          altCLtop(kz-1)=altCL(kz-1,nr)+0.5*(altCL(kz,nr)-altCL(kz-1,nr))
        end do
        altCLtop(nzCL(nr))=altCL(nzCL(nr),nr)+0.5*(altCL(nzCL(nr),nr)-altCL(nzCL(nr)-1,nr))
        clzm=nzCL(nr)-1
        do clz=1,nzCL(nr)
          if (real(part(jpart)%z).lt.altCLtop(clz)) then
            clzm=clz-1
            exit
          endif
        end do
        clz=min(nzCL(nr),clz)
        if (clzm.eq.0 ) then
          dz1=1.
          dz2=1.
          clzm=clz
        else
          dz1=real(part(jpart)%z)-altCL(clzm,nr)
          dz2=altCL(clz,nr)-real(part(jpart)%z)
        endif
        if (dz1.eq.(-1.*dz2)) then
          dz1=1.
          dz2=1.
        endif
        dzz=1./(dz1+dz2)

        !! testing
!        if(ii.lt.100) print*, 'chemreaction: nr, clx, cly, clz, clzm = ',nr, clx, cly, clz, clzm

        ! chem reagent for particle time and location 
        cl_cur=(dz2*clfield_cur(clx,cly,clzm,nr) + &
                dz1*clfield_cur(clx,cly,clz,nr))*dzz

        ! Compute chemical loss
        !**********************

        if (cl_cur.gt.smallnum) then
          indz=nz
          do kz=2,nz
            if (height(kz).gt.part(jpart)%z) then
              indz=kz-1
              exit
            endif
          end do
          temp=(tt(ix,jy,indz,1)*dtt2 &
                 + tt(ix,jy,indz,2)*dtt1)*dttt
          do ks=1,nspec
            if (reaccconst(nr,ks).gt.0.) then
              ! k = CT^Nexp(-D/T)[reagent]
              clrate=reaccconst(nr,ks)*(temp**reacnconst(nr,ks))*exp(-1.*reacdconst(nr,ks)/temp)*cl_cur
              ! new particle mass
              restmass=mass(jpart,ks)*exp(-1.*clrate*abs(lsynctime))
              if (restmass.gt.smallnum) then
                clreacted=mass(jpart,ks)-restmass
                mass(jpart,ks)=restmass
              else
                clreacted=mass(jpart,ks)
                mass(jpart,ks)=0.
              endif
#ifdef _OPENMP
              chem_loss_tmp(nr,ks,ithread)=chem_loss_tmp(nr,ks,ithread)+real(clreacted,kind=dp)
#else
              chem_loss(nr,ks)=chem_loss(nr,ks)+real(clreacted,kind=dp)
#endif
            endif
          end do  ! nspec
        endif   

      end do  ! nreagent

    end do  ! loop over all particles
!$OMP END DO

!$OMP END PARALLEL

#ifdef _OPENMP
  do ithread=1,numthreads
    chem_loss(:,:) = chem_loss(:,:)+chem_loss_tmp(:,:,ithread)
  end do
  deallocate( chem_loss_tmp )
#endif

  end subroutine chemreaction


  subroutine readjrate(jr_name,memid,mm)

  !*****************************************************************************
  !                                                                            *
  !    Reads the average actinic flux rates into memory                        *
  !                                                                            *
  !    Author: Rona Thompson, Sep 2023                                         *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! jr_name           name of the file                                         *
  ! memid             time index to chemical field variable                    *
  ! mm                month to read                                            *
  !                                                                            * 
  !*****************************************************************************

    implicit none

    character(len=256) :: jr_name
    integer            :: memid,mm
    integer            :: ncid,dimid,varid

    ! Read netcdf file
    !******************

    ! open file
    call nf90_err( nf90_open(trim(jr_name),nf90_NOWRITE,ncid) )

    ! longitude
    call nf90_err( nf90_inq_dimid(ncid,'longitude',dimid) )
    call nf90_err( nf90_inquire_dimension(ncid,dimid,len=nxjr) )
    if (.not.allocated(lonjr)) allocate(lonjr(nxjr))
    call nf90_err( nf90_inq_varid(ncid,'longitude',varid) )
    call nf90_err( nf90_get_var(ncid,varid,lonjr) )
    dxjr=abs(lonjr(2)-lonjr(1))

    ! latitude
    call nf90_err( nf90_inq_dimid(ncid,'latitude',dimid) )
    call nf90_err( nf90_inquire_dimension(ncid,dimid,len=nyjr) )
    if (.not.allocated(latjr)) allocate(latjr(nyjr))
    call nf90_err( nf90_inq_varid(ncid,'latitude',varid) )
    call nf90_err( nf90_get_var(ncid,varid,latjr) )
    dyjr=abs(latjr(2)-latjr(1))

    ! jrate field
    call nf90_err( nf90_inq_varid(ncid,'jrate',varid) )
    if (.not.allocated(jrate_average)) then
      allocate(jrate_average(nxjr,nyjr,2))
    endif
    call nf90_err( nf90_get_var(ncid,varid,jrate_average(:,:,memid),start=(/1,1,mm/),count=(/nxjr,nyjr,1/)) )

    ! close file
    call nf90_err( nf90_close(ncid) )

    return

  end subroutine readjrate

  function photo_O1D(sza) result(jrate)

  !*****************************************************************************
  !                                                                            *
  !    Calculates J(O1D) photolysis rate based on solar zenith angle           *
  !                                                                            *
  !    Author: A. Stohl, Nov 2014                                              *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  !    INPUT:                                                                  *
  !    sza        solar zenith angle (degrees)                                 *
  !                                                                            *
  !    OUTPUT:                                                                 *
  !    photo_O1D  J(O1D) photoylsis rate                                       *
  !                                                                            *
  !*****************************************************************************

    implicit none

    real, intent(in)  :: sza
    real :: jrate
    integer :: iz,ik
    real :: z1,z2,zg,f1,f2,dummy
    real :: photo_NO2
    integer, parameter :: nzenith=11
    real, parameter :: pi=3.1415927
    real, dimension(nzenith) :: zangle,fact_photo

    ! zangle: zenith angles for which fact_photo is tabulated
    ! fact_photo: conversion of photolysis rate of NO2 to photolysis 
    !     rate of O3 into O1D as a function of solar zenith angle

    zangle=(/0.,10.,20.,30.,40.,50.,60.,70.,78.,86.,90.0001/)
    fact_photo=(/0.4616E-02,0.4478E-02,0.4131E-02,0.3583E-02,0.2867E-02,&
      &0.2081E-02,0.1235E-02,0.5392E-03,0.2200E-03,0.1302E-03,0.0902E-03/)

    if (sza.lt.90.) then
      do iz=1,nzenith-1
        if(sza.ge.zangle(iz)) ik=iz
      end do
      z1=1./cos(zangle(ik)*pi/180.)
      z2=1./cos(zangle(ik+1)*pi/180.)
      zg=1./cos(sza*pi/180.)
      dummy=(zg-z1)/(z2-z1)
      f1=alog(fact_photo(ik))
      f2=alog(fact_photo(ik+1))
      photo_NO2=1.45e-2*exp(-0.4/cos(sza*pi/180.))
      jrate=photo_NO2*exp(f1+(f2-f1)*dummy)
    else
      jrate=0.
    endif

  end function photo_O1D


  function zenithangle(ylat,xlon,jjjjmmdd,hhmmss) result(sza)

  !*****************************************************************************
  !                                                                            *
  !    This function returns the sinus of solar elevation as a function        *
  !    of geographic longitde, latitude and GMT-Time.                          *
  !                                                                            *
  !    Author: G. Wotawa, Nov 1993                                             *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  !     INPUT:                                                                 *
  !                                                                            *
  !            ylat          geographical latitude  [DEG]                      *
  !            xlon          geographical longitude [DEG]                      *
  !            jjjj          Year                                              *
  !            mm            Month                                             *
  !            dd            Day                                               *
  !            hh            Hour                                              *
  !            minute        Minute                                            *
  !                                                                            *
  !*****************************************************************************

    implicit none

    integer,       intent(in)  :: jjjjmmdd,hhmmss
    real,          intent(in)  :: ylat,xlon
    real    :: sza
    integer :: jjjj,mm,id,iu,minute 
    integer :: ndaynum
    real :: sinsol,solelev
    real :: rnum,rylat,ttime,dekl,rdekl,eq
    real,parameter :: pi=3.1415927

    jjjj=jjjjmmdd/10000
    mm=jjjjmmdd/100-jjjj*100
    id=jjjjmmdd-jjjj*10000-mm*100
    iu=hhmmss/10000
    minute=hhmmss/100-100*iu

    ndaynum=31*(mm-1)+id
    if(mm.gt.2) ndaynum=ndaynum-int(0.4*mm+2.3)
    if((mm.gt.2).and.(jjjj/4*4.eq.jjjj)) ndaynum=ndaynum+1

    rnum=2.*pi*ndaynum/365.
    rylat=pi*ylat/180.
    ttime=real(iu)+real(minute)/60.

    dekl=0.396+3.631*sin(rnum)+0.038*sin(2.*rnum)+0.077*sin(3.*rnum)- &
         22.97*cos(rnum)-0.389*cos(2.*rnum)-0.158*cos(3.*rnum)
    rdekl=pi*dekl/180.

    eq=(0.003-7.343*sin(rnum)-9.47*sin(2.*rnum)- &
         0.329*sin(3.*rnum)-0.196*sin(4.*rnum)+ &
         0.552*cos(rnum)-3.020*cos(2.*rnum)- &
         0.076*cos(3.*rnum)-0.125*cos(4.*rnum))/60.
    sinsol=sin(rylat)*sin(rdekl)+cos(rylat)*cos(rdekl)* &
         cos((ttime-12.+xlon/15.+eq)*pi/12.)
    ! Calculate the maximum solar elevation on that day
    !sinsol=sin(rylat)*sin(rdekl)+cos(rylat)*cos(rdekl)*
    !    &       cos((eq)*pi/12.)
    solelev=asin(sinsol)*180./pi
    sza=90.-solelev

  end function zenithangle


end module chemistry_mod

