!*****************************************************************
!
!     Photochemical routine 
!
!     Author: Franck Lefevre
!             Benjamin Charnay
!             Yassin Jaziri
!     ------
!
!     Version: 27/05/2016
!     update 06/03/2021 generic tracer/network + photolysis online (Yassin Jaziri)
!
!*****************************************************************

subroutine photochemistry_asis(nlayer, ngrid,                                   &
                          ig, lswitch, zycol, sza, fractcol, ptimestep, press,  &
                          alt, temp, dens, zmmean, dist_sol, surfdust1d,        &
                          surfice1d, tau, iter,zdtchim)

      use callkeys_mod
      use comcstfi_mod, only:  r,cpp,avocado,pi
      use tracer_h
      use types_asis
      use chimiedata_h
      use photolysis_mod

implicit none

!===================================================================
!     inputs:
!===================================================================

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

integer :: ig                 ! grid point index
      
real :: sza                   ! solar zenith angle (deg)
real :: fractcol              ! day fraction
real :: ptimestep             ! physics timestep (s)
real :: press(nlayer)         ! pressure (hpa)
real :: alt(nlayer)           ! altitude (km)
real :: temp(nlayer)          ! temperature (k)
real :: dens(nlayer)          ! density (cm-3)
real :: zmmean(nlayer)        ! mean molar mass (g/mole)
real :: dist_sol              ! sun distance (au) 
real :: surfdust1d(nlayer)    ! dust surface area (cm2/cm3)
real :: surfice1d(nlayer)     ! ice surface area (cm2/cm3)
real :: tau                   ! optical depth at 7 hpa

!===================================================================
!     input/output:
!===================================================================
      
real :: zycol(nlayer,nesp)    ! chemical species volume mixing ratio

!===================================================================
!     output:
!===================================================================
      
integer :: iter(nlayer)       ! iteration counter
real    :: zdtchim(nlayer)    ! temperature modification by ozone

!===================================================================
!     local:
!===================================================================

integer :: phychemrat         ! (physical timestep)/(nominal chemical timestep)
integer :: ilev, iesp, iphot, iw
integer :: lswitch

logical, save :: firstcall = .true.

real :: stimestep           ! standard timestep for the chemistry (s) 
real :: ctimestep           ! real timestep for the chemistry (s) 
real :: dt_guess            ! first-guess timestep (s) 
real :: dt_corrected        ! corrected timestep (s) 
real :: dt_min = 1.         ! minimum allowed timestep (s) 
real :: dtg                 ! correction factor for the timestep (s) 
real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
real :: time                ! internal time (between 0 and ptimestep, in s)
real :: rho(nlayer)         ! mass density (kg/m-3)

real, dimension(nlayer,nesp)            :: rm   ! mixing ratios 
real (kind = 8), dimension(nesp)        :: cold ! number densities at previous timestep (molecule.cm-3) 
real (kind = 8), dimension(nlayer,nesp) :: c    ! number densities at current timestep (molecule.cm-3) 
real (kind = 8), dimension(nesp)        :: cnew ! number densities at next timestep (molecule.cm-3) 
 
! reaction rates
 
real (kind = 8), allocatable, save     :: v_phot(:,:)
real (kind = 8), allocatable, save     :: v_3(:,:)
real (kind = 8), allocatable, save     :: v_4(:,:)
real (kind = 8), allocatable, save     :: e_phot(:,:)    ! photolysis rates by energie (J.mol-1.s-1)

integer, save                          :: nreact         ! number of reactions in reactions files

! matrix

real (kind = 8), dimension(nesp,nesp)  :: mat, mat1
integer, dimension(nesp)               :: indx
integer                                :: code

! production and loss terms (for first-guess solution only)

real (kind = 8), dimension(nesp)       :: prod, loss

!===================================================================
!     initialisation of the reaction indexes
!===================================================================

if (firstcall) then
   print*,'photochemistry: initialize indexes'
   call indice(nreact)
   allocate(v_phot(nlayer,nb_phot_max))
   allocate(v_3(nlayer,nb_reaction_3_max))
   allocate(v_4(nlayer,nb_reaction_4_max))
   allocate(v_phot_3d(ngrid,nlayer,nb_phot_hv_max))
   allocate(e_phot(nlayer,nb_phot_max))
   v_phot(:,:)      = 0.
   v_3(:,:)         = 0.
   v_4(:,:)         = 0.
   v_phot_3d(:,:,:) = 0.
   e_phot(:,:)      = 0.

   ! Need to be done after subroutine indice
   if (jonline) then
     print*,'calchim: Read UV absorption cross-sections'
     ! Jonline cannot run without photolysis
     if (nb_phot_hv_max == 0) then
       print*,'Jonline cannot run without photolysis'
       print*,'set jonline=.false. in callphys.def'
       print*,'or set photolysis reactions in monoreact file'
       call abort
     end if
     call init_photolysis(nlayer,nreact)     ! on-line photolysis
     allocate(albedo_snow_chim(nw))
     allocate(albedo_co2_ice_chim(nw))
    ! Step 1 : Initialisation of the Spectral Albedos.
     DO iw=1,nw
      albedo_snow_chim(iw)=albedosnow
      albedo_co2_ice_chim(iw)=albedoco2ice

      ! Spectral Snow Albedo Calculation.
      call albedo_snow_calc(wc(iw)*1.0e-3,        &
                     albedo_snow_chim(iw),  & 
                     albedosnow)
      
     ENDDO
   end if

   firstcall = .false.
end if

!===================================================================
!     initialisation of mixing ratios and densities       
!===================================================================

call gcmtochim(nlayer, zycol, lswitch, nesp, dens, rm, c) 

!===================================================================
!     interpolation of photolysis rates in the lookup table      
!===================================================================

if (jonline) then
   if (sza <= 95.) then
      call photolysis_online(nlayer, alt, press, temp, zmmean, rm,   &
                             tau, sza, dist_sol, v_phot, e_phot, ig, ngrid, nreact)
      if (ngrid.eq.1) then
        do iphot = 1,nb_phot_hv_max
          v_phot(:,iphot) = v_phot(:,iphot)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
          e_phot(:,iphot) = e_phot(:,iphot)* 0.25 / cos(sza*pi/180.) ! globally averaged = divide by 4
        end do
      elseif(diurnal .eqv. .false.) then
        do iphot = 1,nb_phot_hv_max
          v_phot(:,iphot) = v_phot(:,iphot) * fractcol
          e_phot(:,iphot) = e_phot(:,iphot) * fractcol
        end do
      endif
   else ! night
      v_phot(:,:) = 0.
      e_phot(:,:) = 0.
   end if
   ! save photolysis for output
   v_phot_3d(ig,:,:) = v_phot(:,1:nb_phot_hv_max)
else if (nb_phot_hv_max /= 0) then
   call photolysis_asis(nlayer, ngrid, lswitch, press, temp, sza, fractcol,tau, zmmean, dist_sol, & 
                     rm(:,indexchim('co2')), rm(:,indexchim('o3')), rm(:,indexchim('ch4')), v_phot)
   ! save photolysis for output
   v_phot_3d(ig,:,:) = v_phot(:,1:nb_phot_hv_max)
end if

!===================================================================
!     reaction rates                                     
!===================================================================

call reactionrates(nlayer, nreact, c, lswitch, dens, &
                   press, temp, zmmean, sza, surfdust1d, surfice1d, v_phot, v_3, v_4)

zdtchim(:) = 0.
rho(:)     = (press(:)*1e2)/(r*temp(:))

!===================================================================
!     temperature modification by photolysis reaction                                    
!===================================================================

if (photoheat .and. jonline) then ! for now just with jonline

  do iphot = 1,nb_phot_hv_max
    zdtchim(:nlayer-1) = zdtchim(:nlayer-1) + e_phot(:nlayer-1,iphot)*c(:nlayer-1,indexchim(trim(jlabel(iphot,2))))/(cpp*(rho(:nlayer-1)*1e-6)*avocado)
  end do

else

!! o + o2 -> o3
!!dE = 24 ! kcal.mol-1
!!zdtchim(:) = zdtchim(:) + 4.184*24e3*v_4(:,1)*c(:,indexchim('o'))*c(:,indexchim('o2'))*1e6/(cpp*rho*avocado)
!
!! o3 -> o2 + o1d
!! E(250nm) = 480 kJ.mol-1
!j_o3_o1d   = 3
!zdtchim(:) = zdtchim(:) + 480e3*v_phot(:,j_o3_o1d)*c(:,indexchim('o3'))/(cpp*(rho*1e-6)*avocado)
!
!! o3 -> o2 + o
!! E(350nm) = 343 kJ.mol-1
!j_o3_o     = 4
!zdtchim(:) = zdtchim(:) + 343e3*v_phot(:,j_o3_o)*c(:,indexchim('o3'))/(cpp*(rho*1e-6)*avocado)

zdtchim(:) = 0.

end if

!===================================================================
!     stimestep : standard chemical timestep (s)                           
!     ctimestep : real chemical timestep (s),
!                 taking into account the physical timestep                           
!===================================================================

stimestep = 600. ! standard value : 10 mn

phychemrat = nint(ptimestep/stimestep)
phychemrat = 1

ctimestep = ptimestep/real(phychemrat)

!===================================================================
!     loop over levels         
!===================================================================

do ilev = 1,lswitch - 1

!  initializations

   time = 0.
   iter(ilev) = 0
   dt_guess = ctimestep
   cold(:) = c(ilev,:)

!  internal loop for the chemistry

   do while (time < ptimestep)

   iter(ilev) = iter(ilev) + 1    ! iteration counter
  
!  first-guess: fill matrix

   call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)

!  adaptative evaluation of the sub time step

   call define_dt(nesp, dt_corrected, dt_guess, ctimestep, cold(:), c(ilev,:),  &
                  mat1, prod, loss, dens(ilev))

   if (time + dt_corrected > ptimestep) then
      dt_corrected = ptimestep - time
   end if

!  if (dt_corrected /= dt_guess) then  ! the timestep has been modified

!  form the matrix identity + mat*dt_corrected

   mat(:,:) = mat1(:,:)*dt_corrected
   do iesp = 1,nesp
      mat(iesp,iesp) = 1. + mat(iesp,iesp)
   end do

!  solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)

   cnew(:) = c(ilev,:)

#ifdef LAPACK
   call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
#else
   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
   call abort
#endif

!  end if

!  eliminate small values

   where (cnew(:)/dens(ilev) < 1.e-30)
      cnew(:) = 0.
   end where

!  update concentrations

   cold(:)   = c(ilev,:)
   c(ilev,:) = cnew(:)
   cnew(:)   = 0.

!  increment internal time

   time = time + dt_corrected
   dt_guess = dt_corrected     ! first-guess timestep for next iteration

   end do ! while (time < ptimestep)

end do ! ilev

!===================================================================
!     save chemical species for the gcm       
!===================================================================

call chimtogcm(nlayer, zycol, lswitch, nesp, dens, c) 

contains

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

 subroutine define_dt(nesp, dtnew, dtold, ctimestep, cold, ccur, mat1, &
                      prod, loss, dens)

!================================================================
! iterative evaluation of the appropriate time step dtnew
! according to curvature criterion based on
! e = 2 Rtol [r Cn+1 -(1-r) Cn + Cn-1 ]/[(1+r) Cn]
! with r = (tn - tn-1)/(tn+1 - tn)
!================================================================

implicit none

! input

integer :: nesp  ! number of species in the chemistry

real :: dtold, ctimestep
real (kind = 8), dimension(nesp)      :: cold, ccur
real (kind = 8), dimension(nesp,nesp) :: mat1
real (kind = 8), dimension(nesp)      :: prod, loss
real                                  :: dens

! output

real :: dtnew

! local

real (kind = 8), dimension(nesp)      :: cnew
real (kind = 8), dimension(nesp,nesp) :: mat
real (kind = 8) :: atol, ratio, e, es, coef

integer                  :: code, iesp, iter
integer, dimension(nesp) :: indx

real :: dttest

! parameters

real (kind = 8), parameter :: dtmin   = 10.      ! minimum time step (s)
real (kind = 8), parameter :: vmrtol  = 1.e-11   ! absolute tolerance on vmr
real (kind = 8), parameter :: rtol    = 1./0.05  ! 1/rtol recommended value : 0.1-0.02
integer,         parameter :: niter   = 3        ! number of iterations
real (kind = 8), parameter :: coefmax = 2.
real (kind = 8), parameter :: coefmin = 0.1
logical                    :: fast_guess = .true.

dttest = dtold   ! dttest = dtold = dt_guess

atol = vmrtol*dens ! absolute tolerance in molecule.cm-3

do iter = 1,niter

if (fast_guess) then

! first guess : fast semi-implicit method

   do iesp = 1, nesp
      cnew(iesp) = (ccur(iesp) + prod(iesp)*dttest)/(1. + loss(iesp)*dttest)
   end do

else

! first guess : form the matrix identity + mat*dt_guess

   mat(:,:) = mat1(:,:)*dttest
   do iesp = 1,nesp
      mat(iesp,iesp) = 1. + mat(iesp,iesp)
   end do

! form right-hand side (RHS) of the system

   cnew(:) = ccur(:)

! solve the linear system  M*Cn+1 = Cn (RHS in cnew, then replaced by solution)

#ifdef LAPACK
      call dgesv(nesp,1,mat,nesp,indx,cnew,nesp,code)
#else
   write(*,*) "photochemistry_asis error, missing LAPACK routine dgesv"
   call abort
#endif

end if

! ratio old/new subtimestep

ratio = dtold/dttest

! e : local error indicatocitr

e = 0.

do iesp = 1,nesp
   es = 2.*abs((ratio*cnew(iesp) - (1. + ratio)*ccur(iesp) + cold(iesp))   &
         /(1. + ratio)/max(ccur(iesp),atol))

   if (es > e) then
      e = es
   end if
end do
e = rtol*e

! timestep correction

if (e <= 0.) then
  coef = coefmax
else
  coef = max(coefmin, min(coefmax,0.8/sqrt(e)))
end if

dttest = max(dtmin,dttest*coef)
dttest = min(ctimestep,dttest)

end do ! iter

! new timestep

dtnew = dttest

end subroutine define_dt


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

 subroutine reactionrates(nlayer, nreact, c,                    &
                          lswitch, dens, press, t, zmmean, sza, &
                          surfdust1d, surfice1d,                &
                          v_phot, v_3, v_4)
 
!================================================================
! compute reaction rates                                        !
!----------------------------------------------------------------
! reaction               type                array              !
!----------------------------------------------------------------
! A + B    --> C + D     bimolecular         v_4                !
! A + A    --> B + C     quadratic           v_3                !
! A + C    --> B + C     quenching           v_phot             !
! A + ice  --> B + C     heterogeneous       v_phot             !
!================================================================

use comcstfi_mod
use types_asis
use pfunc
use tracer_h
use chimiedata_h

implicit none

!----------------------------------------------------------------------
!     input
!----------------------------------------------------------------------

integer, intent(in)     :: nlayer             ! number of atmospheric layers
integer, intent(in)     :: nreact             ! number of reactions in reactions files
integer                 :: lswitch            ! interface level between lower
                                              ! atmosphere and thermosphere chemistries
real, dimension(nlayer) :: dens               ! total number density (molecule.cm-3)
real, dimension(nlayer) :: press              ! pressure (hPa)
real, dimension(nlayer) :: t                  ! temperature (K)
real, dimension(nlayer) :: zmmean             ! mean molar mass (g/mole)
real                    :: sza                ! solar zenith angle (deg)
real, dimension(nlayer) :: surfdust1d         ! dust surface area (cm2.cm-3)
real, dimension(nlayer) :: surfice1d          ! ice surface area (cm2.cm-3)
real (kind = 8), dimension(nlayer,nesp) :: c  ! species number density (molecule.cm-3)

!----------------------------------------------------------------------
!     output
!----------------------------------------------------------------------

real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4

!----------------------------------------------------------------------
!     local
!----------------------------------------------------------------------

logical,save     :: firstcall = .true.
integer          :: ilev, ireact
integer          :: nb_phot, nb_reaction_3, nb_reaction_4
integer          :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
real, dimension(nlayer) :: ratec ! rate coefficient

!----------------------------------------------------------------------
!     initialisation
!----------------------------------------------------------------------

      nb_phot       = nb_phot_hv_max
      nb_reaction_3 = 0
      nb_reaction_4 = 0

      nb_hv         = 0
      nb_pfunc1     = 0
      nb_pfunc2     = 0
      nb_pfunc3     = 0

!----------------------------------------------------------------------
!     reactions 
!----------------------------------------------------------------------

do ireact = 1,nreact
  if (reactions(ireact)%rtype==0) then
    nb_hv = nb_hv + 1
    cycle
  end if

  ! get right coefficient rate function
  if (reactions(ireact)%rtype==1) then
    nb_pfunc1 = nb_pfunc1 + 1
    if (reactions(ireact)%third_body==0) then            !! third_body check
      ratec = pfunc1(nlayer,t,dens,pfunc1_param(nb_pfunc1))
    else
      ratec = pfunc1(nlayer,t,c(:,reactions(ireact)%third_body),pfunc1_param(nb_pfunc1))
    end if
  else if (reactions(ireact)%rtype==2) then
    nb_pfunc2 = nb_pfunc2 + 1
    if (reactions(ireact)%third_body==0) then            !! third_body check
      ratec = pfunc2(nlayer,t,dens,pfunc2_param(nb_pfunc2))
    else
      ratec = pfunc2(nlayer,t,c(:,reactions(ireact)%third_body),pfunc2_param(nb_pfunc2))
    end if
  else if (reactions(ireact)%rtype==3) then
    nb_pfunc3 = nb_pfunc3 + 1
    if (reactions(ireact)%third_body==0) then            !! third_body check
      ratec = pfunc3(nlayer,t,dens,pfunc3_param(nb_pfunc3))
    else
      ratec = pfunc3(nlayer,t,c(:,reactions(ireact)%third_body),pfunc3_param(nb_pfunc3))
    end if
  else
    print*, 'Error in reactionrates: wrong index coefficient rate vphot'
    print*, 'rtype(ireact) = ',reactions(ireact)%rtype
    print*, 'It should be 1 or 2 '
    print*, 'Please check the reaction ',ireact
    if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
    call abort
  end if

  ! fill the right type index
  if (reactions(ireact)%vtype=="v4") then
    nb_reaction_4 = nb_reaction_4 + 1
    v_4(:,nb_reaction_4) = ratec(:)
    if (reactions(ireact)%three_prod) then             ! three_prod check
      nb_reaction_4 = nb_reaction_4 + 1
      v_4(:,nb_reaction_4) = ratec(:)
    end if
  else if (reactions(ireact)%vtype=="v3") then
    nb_reaction_3 = nb_reaction_3 + 1
    v_3(:,nb_reaction_3) = ratec(:)
    if (reactions(ireact)%three_prod) then             ! three_prod check
      nb_reaction_3 = nb_reaction_3 + 1
      v_3(:,nb_reaction_3) = ratec(:)
    end if
  else if (reactions(ireact)%vtype=="vphot") then
    nb_phot = nb_phot + 1
    v_phot(:,nb_phot) = ratec(:)
    if (reactions(ireact)%three_prod) then             ! three_prod check
      nb_phot = nb_phot + 1
      v_phot(:,nb_phot) = ratec(:)
    end if
  else
    print*, 'Error in reactionrates: wrong type coefficient rate'
    print*, 'vtype(ireact) = ',reactions(ireact)%vtype
    print*, 'It should be v4, v3 or vphot '
    print*, 'Please check the reaction ',ireact
    if (ireact>nreact) print*, 'Please check the reaction count into the code beacause ireact > nreact is not possible'
    call abort
  end if

end do

call reactionrates_HC(nlayer,nesp,dens,t,press,zmmean,sza,c,v_phot,v_4,v_3,nb_phot,nb_reaction_4,nb_reaction_3)

!===========================================================
!  check dimensions 
!===========================================================

if (firstcall) then
 if ((nb_phot /= nb_phot_max)             .or.  &
     (nb_reaction_3 /= nb_reaction_3_max) .or.  &
     (nb_reaction_4 /= nb_reaction_4_max)) then
    print*, 'nb_phot       = ', nb_phot
    print*, 'nb_reaction_4 = ', nb_reaction_4
    print*, 'nb_reaction_3 = ', nb_reaction_3
    print*, 'wrong dimensions in reactionrates' 
    print*, 'nb_phot_max       = ', nb_phot_max
    print*, 'nb_reaction_4_max = ', nb_reaction_4_max
    print*, 'nb_reaction_3_max = ', nb_reaction_3_max
    print*, '------ hard coding reaction ------'
    print*, 'nphot_hard_coding = ', nphot_hard_coding
    print*, 'n4_hard_coding    = ', n4_hard_coding
    print*, 'n3_hard_coding    = ', n3_hard_coding
    call abort
 end if
 if ((nb_hv /= nb_hv_max)         .or.  &
     (nb_pfunc1 /= nb_pfunc1_max) .or.  &
     (nb_pfunc2 /= nb_pfunc2_max) .or.  &
     (nb_pfunc3 /= nb_pfunc3_max)) then
    print*, 'nb_hv         = ', nb_hv
    print*, 'nb_pfunc1     = ', nb_pfunc1
    print*, 'nb_pfunc2     = ', nb_pfunc2
    print*, 'nb_pfunc3     = ', nb_pfunc3
    print*, 'wrong dimensions in reactionrates'
    print*, 'nb_hv_max     = ', nb_hv_max
    print*, 'nb_pfunc1_max = ', nb_pfunc1_max
    print*, 'nb_pfunc2_max = ', nb_pfunc2_max
    print*, 'nb_pfunc3_max = ', nb_pfunc3_max
    call abort
 end if
firstcall = .false.
end if 

end subroutine reactionrates

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

 subroutine fill_matrix(ilev, mat, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4)

!======================================================================
! filling of the jacobian matrix
!======================================================================

use types_asis
use chimiedata_h

implicit none

! input

integer             :: ilev    ! level index
integer             :: nesp    ! number of species in the chemistry
integer, intent(in) :: nlayer  ! number of atmospheric layers

real (kind = 8), dimension(nlayer,nesp)              :: c    ! number densities
real (kind = 8), dimension(nlayer,      nb_phot_max) :: v_phot
real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3
real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4

! output

real (kind = 8), dimension(nesp,nesp), intent(out) :: mat  ! matrix
real (kind = 8), dimension(nesp), intent(out)      :: prod, loss

! local

integer :: iesp
integer :: ind_phot_2,ind_phot_4,ind_phot_6
integer :: ind_3_2,ind_3_4,ind_3_6
integer :: ind_4_2,ind_4_4,ind_4_6,ind_4_8
integer :: iphot,i3,i4

real(kind = jprb) :: eps, eps_4  ! implicit/explicit coefficient

! initialisations 

mat(:,:) = 0.
prod(:)  = 0.
loss(:)  = 0.

! photodissociations
! or reactions a + c -> b + c
! or reactions a + ice -> b + c

do iphot = 1,nb_phot_max

  ind_phot_2 = indice_phot(iphot)%z2
  ind_phot_4 = indice_phot(iphot)%z4
  ind_phot_6 = indice_phot(iphot)%z6

  mat(ind_phot_2,ind_phot_2) = mat(ind_phot_2,ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
  mat(ind_phot_4,ind_phot_2) = mat(ind_phot_4,ind_phot_2) - indice_phot(iphot)%z3*v_phot(ilev,iphot)
  mat(ind_phot_6,ind_phot_2) = mat(ind_phot_6,ind_phot_2) - indice_phot(iphot)%z5*v_phot(ilev,iphot)

  loss(ind_phot_2) = loss(ind_phot_2) + indice_phot(iphot)%z1*v_phot(ilev,iphot)
  prod(ind_phot_4) = prod(ind_phot_4) + indice_phot(iphot)%z3*v_phot(ilev,iphot)*c(ilev,ind_phot_2)
  prod(ind_phot_6) = prod(ind_phot_6) + indice_phot(iphot)%z5*v_phot(ilev,iphot)*c(ilev,ind_phot_2)

end do

! reactions a + a -> b + c 

do i3 = 1,nb_reaction_3_max

  ind_3_2 = indice_3(i3)%z2
  ind_3_4 = indice_3(i3)%z4
  ind_3_6 = indice_3(i3)%z6

  mat(ind_3_2,ind_3_2) = mat(ind_3_2,ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
  mat(ind_3_4,ind_3_2) = mat(ind_3_4,ind_3_2) - indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)
  mat(ind_3_6,ind_3_2) = mat(ind_3_6,ind_3_2) - indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)

  loss(ind_3_2) = loss(ind_3_2) + indice_3(i3)%z1*v_3(ilev,i3)*c(ilev,ind_3_2)
  prod(ind_3_4) = prod(ind_3_4) + indice_3(i3)%z3*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)
  prod(ind_3_6) = prod(ind_3_6) + indice_3(i3)%z5*v_3(ilev,i3)*c(ilev,ind_3_2)*c(ilev,ind_3_2)

end do

! reactions a + b -> c + d 

eps = 1.d-10

do i4 = 1,nb_reaction_4_max

  ind_4_2 = indice_4(i4)%z2
  ind_4_4 = indice_4(i4)%z4
  ind_4_6 = indice_4(i4)%z6
  ind_4_8 = indice_4(i4)%z8

  eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps)
  eps_4 = min(eps_4,1.0_jprb)

  mat(ind_4_2,ind_4_2) = mat(ind_4_2,ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4) 
  mat(ind_4_2,ind_4_4) = mat(ind_4_2,ind_4_4) + indice_4(i4)%z1*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
  mat(ind_4_4,ind_4_2) = mat(ind_4_4,ind_4_2) + indice_4(i4)%z3*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
  mat(ind_4_4,ind_4_4) = mat(ind_4_4,ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)   
  mat(ind_4_6,ind_4_2) = mat(ind_4_6,ind_4_2) - indice_4(i4)%z5*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
  mat(ind_4_6,ind_4_4) = mat(ind_4_6,ind_4_4) - indice_4(i4)%z5*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)
  mat(ind_4_8,ind_4_2) = mat(ind_4_8,ind_4_2) - indice_4(i4)%z7*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
  mat(ind_4_8,ind_4_4) = mat(ind_4_8,ind_4_4) - indice_4(i4)%z7*v_4(ilev,i4)*eps_4*c(ilev,ind_4_2)

  loss(ind_4_2) = loss(ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*c(ilev,ind_4_4)
  loss(ind_4_4) = loss(ind_4_4) + indice_4(i4)%z3*v_4(ilev,i4)*c(ilev,ind_4_2)
  prod(ind_4_6) = prod(ind_4_6) + indice_4(i4)%z5*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)
  prod(ind_4_8) = prod(ind_4_8) + indice_4(i4)%z7*v_4(ilev,i4)*c(ilev,ind_4_2)*c(ilev,ind_4_4)

end do

end subroutine fill_matrix

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

 subroutine indice(nreact)

!================================================================
! set the "indice" arrays used to fill the jacobian matrix      !
!----------------------------------------------------------------
! reaction                                   type               !
!----------------------------------------------------------------
! A + hv   --> B + C     photolysis          indice_phot        ! 
! A + B    --> C + D     bimolecular         indice_4           !
! A + A    --> B + C     quadratic           indice_3           !
! A + C    --> B + C     quenching           indice_phot        !
! A + ice  --> B + C     heterogeneous       indice_phot        !
!================================================================

use types_asis
use datafile_mod
use ioipsl_getin_p_mod, only: getin_p
use tracer_h, only: nesp
use chimiedata_h
use callkeys_mod, only: jonline

implicit none

! output

integer, intent(out)              :: nreact          ! number of reactions in reactions files

! local

integer :: nb_phot, nb_reaction_3, nb_reaction_4
integer :: nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3
integer :: iq, ireact

character(len = 128)              :: reactfile       ! reactions table file name
integer                           :: nbq             ! number of species in reactions file
integer                           :: nlines          ! number of lines in reactions file
integer                           :: pnreact         ! number of reaction in photochemical reactions file
integer                           :: bnreact         ! number of reaction in bimolecular reactions file
integer                           :: qnreact         ! number of reaction in quadratic reactions file
integer                           :: specheck(nesp)  ! to count the species of traceur.def in reactions file
integer                           :: specheckr(nesp) ! to count the reactant species of traceur.def in reactions file
integer                           :: specheckp(nesp) ! to count the product species of traceur.def in reactions file

nb_phot        = 0
nb_reaction_3  = 0
nb_reaction_4  = 0
nb_phot_hv_max = 0

nb_hv          = 0
nb_pfunc1      = 0
nb_pfunc2      = 0
nb_pfunc3      = 0

!===========================================================
!  Read chemical reactions
!===========================================================

! Get number of reactions
nlines         = 0
nreact         = 0
pnreact        = 0
bnreact        = 0
qnreact        = 0

write(*,*) "Read reaction file"
reactfile = "reactfile" ! default
call getin_p("reactfile",reactfile) ! default path
write(*,*) " reactfile = ",trim(reactfile)

call count_react(reactfile,nlines,nreact,nb_phot,nb_reaction_4,nb_reaction_3,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)

!!!!!!!!!!! Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!!!!!
nb_phot       = nb_phot       + nphot_hard_coding
nb_reaction_4 = nb_reaction_4 + n4_hard_coding
nb_reaction_3 = nb_reaction_3 + n3_hard_coding
!!!!!!!!!!! END Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!

write(*,*)'number of reaction in reaction files is = ',nreact
print*, 'nb_phot       = ', nb_phot
print*, 'nb_reaction_4 = ', nb_reaction_4
print*, 'nb_reaction_3 = ', nb_reaction_3
print*, '****************'
print*, 'nb_hv         = ', nb_hv
print*, 'nb_pfunc1     = ', nb_pfunc1
print*, 'nb_pfunc2     = ', nb_pfunc2
print*, 'nb_pfunc3     = ', nb_pfunc3
nb_phot_max       = nb_phot
nb_reaction_4_max = nb_reaction_4
nb_reaction_3_max = nb_reaction_3
nb_hv_max         = nb_hv
nb_pfunc1_max     = nb_pfunc1
nb_pfunc2_max     = nb_pfunc2
nb_pfunc3_max     = nb_pfunc3

! Allocate
allocate(indice_phot(nb_phot_max))
allocate(indice_4(nb_reaction_4_max))
allocate(indice_3(nb_reaction_3_max))
allocate(reactions(nreact))
allocate(jlabel(nb_phot_hv_max,2))
allocate(jparamline(nb_hv_max))
allocate(pfunc1_param(nb_pfunc1_max))
allocate(pfunc2_param(nb_pfunc2_max))
allocate(pfunc3_param(nb_pfunc3_max))

! Get reactants, products and number of species in reactfile
! inititialize variables
nbq             = 0
specheck(:)     = 0
specheckr(:)    = 0
specheckp(:)    = 0
reactions(:)    = reaction('',-1,0,.false.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
pfunc1_param(:) = rtype1(0.,0.,0.,0.,0.)
pfunc2_param(:) = rtype2(0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
pfunc3_param(:) = rtype3(0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.)
nb_phot         = 0
nb_reaction_4   = 0
nb_reaction_3   = 0
jlabel(:,:)     = ''
jparamline(:)   = ''

call get_react(reactfile,nlines,nreact,specheck,specheckr,specheckp,nbq,nb_phot,nb_reaction_4,nb_reaction_3)

! rewrite jlabel corretly for output
do ireact=1,nb_phot_hv_max
  if (reactions(ireact)%three_prod) then
    jlabel(ireact+1:nb_phot_hv_max,2) = jlabel(ireact:nb_phot_hv_max-1,2)
    jlabel(ireact+1:nb_phot_hv_max,1) = jlabel(ireact:nb_phot_hv_max-1,1)
    jlabel(ireact,1) = trim(jlabel(ireact,1))//'_a'
    jlabel(ireact+1,1) = trim(jlabel(ireact+1,1))//'_b'
  end if
end do

!!!!!!!!!!! Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!!!!!
call indice_HC(nb_phot,nb_reaction_4,nb_reaction_3)
!!!!!!!!!!! END Hard coding reaction !!!!!!!!!!!!!!!!!!!!!!!

write(*,*)'number of species in reaction files is = ',nbq
write(*,*)'species in reaction files are:'

do iq=1,nesp
  if (specheck(iq)==1) print*, chemnoms(iq)
end do

!===========================================================
!  check species only destroy or product
!===========================================================

do iq=1,nesp
  if (specheckr(iq)/=specheckp(iq)) then
    if (specheckr(iq)==0 .and. specheckp(iq)==1) then
      print*, 'WARNING: ', chemnoms(iq),' is only product'
    else if (specheckr(iq)==1 .and. specheckp(iq)==0) then
      print*, 'WARNING: ', chemnoms(iq),' is only destroy'
    else
      print*, 'Error in indice: bad value in specheckr or specheckp'
      call abort
    end if
  end if
end do 

!===========================================================
!  check stochiometry
!===========================================================

! If you found a way

!===========================================================
!  check dimensions 
!===========================================================

if (jonline) then
  nd = nb_hv_max
else if (nb_phot_hv_max /= 0) then
  print*,'calchim: Read photolysis lookup table'
  call read_phototable
end if

if ((nb_phot /= nb_phot_max)             .or.  &
    (nb_reaction_3 /= nb_reaction_3_max) .or.  &
    (nb_reaction_4 /= nb_reaction_4_max) .or.  &
    (nd /= nb_hv_max)) then
   print*, 'nb_phot       = ', nb_phot
   print*, 'nb_reaction_4 = ', nb_reaction_4
   print*, 'nb_reaction_3 = ', nb_reaction_3
   print*, 'nd            = ', nd
   print*, 'wrong dimensions in indice' 
   print*, 'nb_phot_max       = ', nb_phot_max
   print*, 'nb_reaction_4_max = ', nb_reaction_4_max
   print*, 'nb_reaction_3_max = ', nb_reaction_3_max
   print*, 'nb_phot_hv_max    = ', nb_phot_hv_max
   print*, 'nb_hv_max         = ', nb_hv_max
   call abort
end if  

end subroutine indice

!*****************************************************************

      subroutine gcmtochim(nlayer, zycol, lswitch, nesp, dens, rm, c) 

!*****************************************************************

      use callkeys_mod

      implicit none


!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     input:
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      
      integer, intent(in) :: nlayer ! number of atmospheric layers
      integer, intent(in) :: nesp   ! number of species in the chemistry
      integer :: lswitch            ! interface level between chemistries

      real :: zycol(nlayer,nesp)    ! volume mixing ratios in the gcm
      real :: dens(nlayer)          ! total number density (molecule.cm-3) 

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     output:
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      real, dimension(nlayer,nesp)            :: rm ! volume mixing ratios 
      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3) 
      
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     local:
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      integer      :: l, iesp
      
      rm(:,:) = 0.
      
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     initialise mixing ratios
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      do l = 1,lswitch-1
         rm(l,:) = zycol(l,:)
      end do 

      where (rm(:,:) < 1.e-30)
         rm(:,:) = 0.
      end where

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     initialise number densities
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      do iesp = 1,nesp
         do l = 1,lswitch-1
            c(l,iesp) = rm(l,iesp)*dens(l)
         end do
      end do

      end subroutine gcmtochim

!*****************************************************************
 
      subroutine chimtogcm(nlayer, zycol, lswitch, nesp, dens, c) 

 
!*****************************************************************

      use callkeys_mod

      implicit none


!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     inputs:
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
      integer, intent(in) :: nlayer  ! number of atmospheric layers
      integer, intent(in) :: nesp    ! number of species in the chemistry
      integer :: lswitch             ! interface level between chemistries

      real :: dens(nlayer)     ! total number density (molecule.cm-3) 
      real (kind = 8), dimension(nlayer,nesp) :: c  ! number densities (molecule.cm-3) 

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     output:
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       
      real zycol(nlayer,nesp)  ! volume mixing ratios in the gcm

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     local:
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      integer l
      
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!     save mixing ratios for the gcm
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      do l = 1,lswitch-1
         zycol(l,:)  = c(l,:)/dens(l) 
      end do 

      end subroutine chimtogcm

!*****************************************************************

    subroutine split_str(line,words,n,nmax)

!*****************************************************************

        implicit none
        character(*), intent(in)  :: line
        integer,      intent(in)  :: nmax
        character(*), intent(out) :: words(nmax)
        integer,      intent(out) :: n          ! number of words at the end
        integer :: ios
        character(100) :: buf(100)  ! large buffer!
    
        n = 0
        do
            n = n + 1
            if (n>nmax) then
              print*,'Error in split_str: to much words'
              print*,'limit = ',nmax
              print*,'change it, if you want, with increasing nmax'
              print*,'line is:',line
              call abort
            end if
            read(line,*,iostat=ios) buf(1:n)  ! use list-directed input
            if (ios==0) then
                words(1:n) = buf(1:n)   ! if success, copy to the original array
            else
                n = n-1
                exit       ! if all the words are obtained, finish
            endif
        enddo
    end subroutine split_str

!*****************************************************************

    subroutine find_vtype(reactline,vtype)

!*****************************************************************

        use tracer_h,     only: nesp
        use chimiedata_h, only: indexchim

        implicit none

        character(len = 50), intent(in)    :: reactline ! all reactants of one reaction
        character(*),        intent(inout) :: vtype     ! "v3", "v4" or "vphot"

        integer             :: nwr                      ! number of reactant for a reaction
        integer,parameter   :: nmax = 5                 ! number max of reactants or products
        character(len = 24) :: words(nmax)              ! to get words in reactants and products line
        integer             :: stochio(nesp)            ! stoichiometry of reaction
        integer             :: iword

        ! init
        stochio(:) = 0

        ! split reactant
        call split_str(reactline,words,nwr,nmax)

        ! set stochio
        do iword = 1,nwr
            if (trim(words(iword))=="M") exit
            if (trim(words(iword))/="hv"         &
                .and. trim(words(iword))/="dummy") stochio(indexchim(words(iword))) = stochio(indexchim(words(iword))) + 1
        end do

        ! set vtype
        if (sum(stochio)==1) then
            vtype = "vphot"
        else if (sum(stochio)==2) then
            if (any(stochio==2)) then
                vtype = "v3"
            else
                vtype = "v4"
            end if
        else if (sum(stochio)==3) then
            if (any(stochio==2)) then
                vtype = "v4"
            else
                print*,'Error in photochemistry_asis (find_vtype):'
                print*,'3 different reactants not implemented'
                call abort
            end if
        else
            print*,'Error in photochemistry_asis (find_vtype):'
            print*,'more than 2 different reactants not implemented'
            call abort
        end if

    end subroutine find_vtype

!*****************************************************************

    subroutine count_react(reactfile,nlines,nreact,nb_phot,nb_reaction_4,nb_reaction_3,nb_hv,nb_pfunc1,nb_pfunc2,nb_pfunc3)

!*****************************************************************

        use types_asis, only : nb_phot_hv_max

        implicit none
        character(*), intent(in)     :: reactfile     ! name of the file to read 
        integer,      intent(inout)  :: nlines        ! number of lines in the file
        integer,      intent(out)    :: nreact        ! real number of reaction
        integer,      intent(inout)  :: nb_phot       ! dimension of "vphot" reaction
        integer,      intent(inout)  :: nb_reaction_4 ! dimension of "v4" reaction
        integer,      intent(inout)  :: nb_reaction_3 ! dimension of "v3" reaction
        integer,      intent(inout)  :: nb_hv         ! number of typerate 0 reaction
        integer,      intent(inout)  :: nb_pfunc1     ! number of typerate 1 reaction
        integer,      intent(inout)  :: nb_pfunc2     ! number of typerate 2 reaction
        integer,      intent(inout)  :: nb_pfunc3     ! number of typerate 3 reaction

        ! local
        character(len = 50) :: reactline              ! all reactants of one reaction
        character(len = 50) :: prodline               ! all produts of one reaction
        integer             :: typerate               ! get the type of the rate reaction coefficient (pfunc_i)
        character(5)        :: vtype                  ! "v3", "v4" or "vphot"
        integer             :: nwp                    ! number of products for a reaction
        integer,parameter   :: nmax = 5               ! number max of reactants or products
        character(len = 24) :: words(nmax)            ! to get words in reactants and products line
        integer             :: ierr

        nreact = 0
        vtype  = ''

        open(402, form = 'formatted', status = 'old',                &
              file ='chemnetwork/'//trim(reactfile),iostat=ierr)

        if (ierr /= 0) THEN
           write(*,*)'Error : cannot open chemical reaction file : chemnetwork/'//trim(reactfile)
           write(*,*)'It should be in your simulation directory in chemnetwork directory'
           write(*,*)'   You can change the input chemical reactions file name in'
           write(*,*)'   callphys.def with:'
           write(*,*)'   monoreact=filename or bimolreact=filename or quadrareact=filename'
           write(*,*)'   depending on what chemical reaction type it is'
           call abort
        end if
    
        do
          read(402,'(A,A,I2)',iostat=ierr) reactline, prodline, typerate
          if (ierr<0) exit
          ! count the number of lines
          nlines = nlines + 1
          if (reactline(1:1)/='!' .and. reactline(1:1)/='') then
            ! count the number of reaction
            nreact = nreact + 1

            call find_vtype(reactline,vtype)
            call split_str(prodline,words,nwp,nmax)

            ! count the dimension of each rate reaction coefficient type
            if (trim(vtype)=="vphot") then
                nb_phot = nb_phot + 1
                ! check three product reaction
                if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
                          .and. trim(words(1))/=trim(words(3)) &
                          .and. trim(words(2))/=trim(words(3))) nb_phot = nb_phot + 1
            else if (trim(vtype)=="v4") then
                nb_reaction_4 = nb_reaction_4 + 1
                ! check three product reaction
                if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
                          .and. trim(words(1))/=trim(words(3)) &
                          .and. trim(words(2))/=trim(words(3))) nb_reaction_4 = nb_reaction_4 + 1
            else if (trim(vtype)=="v3") then
                nb_reaction_3 = nb_reaction_3 + 1
                ! check three product reaction
                if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
                          .and. trim(words(1))/=trim(words(3)) &
                          .and. trim(words(2))/=trim(words(3))) nb_reaction_3 = nb_reaction_3 + 1
            else
                print*,'Error in photochemistry_asis (count_react):'
                print*,'vtype not found'
                call abort
            end if

            ! count the number of each rate reaction coefficient type
            if (typerate==0) then
              nb_hv = nb_hv + 1
              nb_phot_hv_max = nb_phot_hv_max + 1
              if (nwp>2 .and. trim(words(1))/=trim(words(2)) &
                        .and. trim(words(1))/=trim(words(3)) &
                        .and. trim(words(2))/=trim(words(3))) nb_phot_hv_max = nb_phot_hv_max + 1
            else if (typerate==1) then
              nb_pfunc1 = nb_pfunc1 + 1
            else if (typerate==2) then
              nb_pfunc2 = nb_pfunc2 + 1
            else if (typerate==3) then
              nb_pfunc3 = nb_pfunc3 + 1
            else
              print*, 'Error in indice: wrong index coefficient rate line ',nlines
              print*, 'in file : chemnetwork/'//trim(reactfile)
              print*, 'It should be 0 for photolysis reations and 1 or 2 for the others'
              print*, 'And not : ', typerate
              call abort
            end if

          end if

        end do

        close(402)

    end subroutine count_react

!*****************************************************************

    subroutine get_react(reactfile,nlines,nreact,specheck,specheckr,specheckp, &
                         nbq,nb_phot,nb_reaction_4,nb_reaction_3)

!*****************************************************************

        use types_asis
        use tracer_h
        use chimiedata_h, only: indexchim
        use callkeys_mod, only: jonline

        implicit none
        character(*), intent(in)     :: reactfile          ! name of the file to read 
        integer,      intent(in)     :: nlines             ! number of lines in the file
        integer,      intent(in)     :: nreact             ! real number of reaction in the file
        integer,      intent(out)    :: nb_phot            ! number of "vphot" calculation reactions
        integer,      intent(out)    :: nb_reaction_4      ! number of "vphot" calculation reactions
        integer,      intent(out)    :: nb_reaction_3      ! number of "vphot" calculation reactions
        integer,      intent(inout)  :: specheck(nesp)     ! to count the species of traceur.def in reactions file
        integer,      intent(inout)  :: specheckr(nesp)    ! to count the reactant species of traceur.def in reactions file
        integer,      intent(inout)  :: specheckp(nesp)    ! to count the product species of traceur.def in reactions file
        integer,      intent(inout)  :: nbq                ! number of species in reactions file

        ! local
        character(len = 50)  :: reactline                  ! all reactants of one reaction
        character(len = 50)  :: prodline                   ! all produts of one reaction
        integer              :: nwr                        ! number of reactants for a reaction
        integer              :: nwp                        ! number of products for a reaction
        integer,parameter    :: nmax = 5                   ! number max of reactants or products
        character(len = 24)  :: words(nmax)                ! to get words in reactants and products line
        real                 :: nindice(2*nmax)            ! stoichiometry of species (for indice variables)
        integer              :: iindice(2*nmax)            ! indice of species (for indice variables)
        character(len = 500) :: paramline                  ! line of reactions parameters
        integer              :: stochio(nesp)              ! stoichiometry of reaction
        integer              :: ierr, ilines, ireact, i_dummy, iw, ispe, iloc
        integer              :: nb_phot_hv, nb_hv, nb_pfunc1, nb_pfunc2, nb_pfunc3

        i_dummy       = 1
        ireact        = 0
        nb_phot       = 0
        nb_phot_hv    = 0
        nb_hv         = 0
        nb_pfunc1     = 0
        nb_pfunc2     = 0
        nb_pfunc3     = 0

        open(402, form = 'formatted', status = 'old',                &
              file ='chemnetwork/'//trim(reactfile),iostat=ierr)

        if (ierr /= 0) THEN
           write(*,*)'Error : cannot open chemical reaction file : chemnetwork/'//trim(reactfile)
           write(*,*)'It should be in your simulation directory in chemnetwork directory'
           write(*,*)'   You can change the input chemical reactions file name in'
           write(*,*)'   callphys.def with:'
           write(*,*)'   monoreact=filename or bimolreact=filename or quadrareact=filename'
           write(*,*)'   depending on what chemical reaction type it is'
           call abort
        end if

        do ilines=1,nlines
          paramline = ''

          read(402,'(A,A,A)') reactline, prodline, paramline

          ! continue only if it's not a comment line
          if (reactline(1:1)/='!' .and. reactline(1:1)/='') then

            ! increment number of reactions
            ireact = ireact + 1

            ! fill reaction vtype
            call find_vtype(reactline,reactions(ireact)%vtype)
            if (trim(reactions(ireact)%vtype)=="v4") then
                nb_reaction_4 = nb_reaction_4 + 1
            else if (trim(reactions(ireact)%vtype)=="v3") then
                nb_reaction_3 = nb_reaction_3 + 1
            else if (trim(reactions(ireact)%vtype)=="vphot") then
                nb_phot = nb_phot + 1
            else
                print*,'Error in photochemistry_asis (get_react):'
                print*,'vtype not found'
                call abort
            end if

            ! fill reaction rate type and parameters
            if (trim(paramline)=='') then
                print*, 'Error in reactfile: where are the parameters - line',ilines
                call abort
            else
                read(paramline,*) reactions(ireact)%rtype
                if (reactions(ireact)%rtype==1) then
                  nb_pfunc1 = nb_pfunc1 + 1
                  read(paramline,*) reactions(ireact)%rtype, pfunc1_param(nb_pfunc1)
                else if (reactions(ireact)%rtype==0) then
                  nb_hv = nb_hv + 1
                  nb_phot_hv = nb_phot_hv + 1
                  if (jonline) then
                    read(paramline,'(I5,A,A)') reactions(ireact)%rtype, jlabel(nb_hv,1), jparamline(nb_hv)
                  else
                    read(paramline,*) reactions(ireact)%rtype, jlabel(nb_hv,1)
                  end if
                else if (reactions(ireact)%rtype==2) then
                  nb_pfunc2 = nb_pfunc2 + 1
                  read(paramline,*) reactions(ireact)%rtype, pfunc2_param(nb_pfunc2)
                else if (reactions(ireact)%rtype==3) then
                  nb_pfunc3 = nb_pfunc3 + 1
                  read(paramline,*) reactions(ireact)%rtype, pfunc3_param(nb_pfunc3)
                end if
            end if

            ! WARNING: the implementation is limited to 3 different reactants (for now only 2) and 3 different products
            nindice(:) = 0.0     ! 3 first indice for reactants, 3 following for products
            iindice(:) = i_dummy ! 3 first indice for reactants, 3 following for products

            !-----------------!
            !--- reactants ---!
            !-----------------!

            ! init for reactant
            stochio(:) = 0
            ! split reactant
            call split_str(reactline,words,nwr,nmax)
            ! set species that are photolysed
            if (reactions(ireact)%rtype==0) jlabel(nb_hv,2) = words(1)
            ! find reactant stochio
            do iw = 1,nwr
                ! fill third body index
                if (trim(words(iw))=="M") then
                    if (iw==nwr) then
                        exit
                    else if (iw==nwr-1) then
                        reactions(ireact)%third_body = indexchim(words(iw+1))
                        exit
                    else
                        print*, 'Error in reactfile: just only one specie can be after M corresponding to the third body - line',ilines
                        call abort
                    end if
                end if
                ! count stochio
                if (trim(words(iw))/="hv"         &
                    .and. trim(words(iw))/="dummy") stochio(indexchim(words(iw))) = stochio(indexchim(words(iw))) + 1
            end do
            ! get reaction reactant stochio and indice
            iloc = 0
            do ispe = 1,nesp
                if (stochio(ispe)==0) cycle
                iloc = iloc + 1
                nindice(iloc) = stochio(ispe)
                iindice(iloc) = ispe
                ! check up: to count the species used in 'reactfile'
                if (specheck(ispe)==0) then
                    specheckr(ispe) = 1
                    specheck(ispe) = 1
                    nbq = nbq + 1
                else if (specheckr(ispe)==0) then
                    specheckr(ispe) = 1
                end if
            end do
            ! fill reaction reactant stochio and indice
            reactions(ireact)%ir1 = nindice(1)
            reactions(ireact)%r1  = iindice(1)
            reactions(ireact)%ir2 = nindice(2)
            reactions(ireact)%r2  = iindice(2)
            reactions(ireact)%ir3 = nindice(3)
            reactions(ireact)%r3  = iindice(3)

            !----------------!
            !--- products ---!
            !----------------!

            ! init for products
            stochio(:) = 0
            ! split products
            call split_str(prodline,words,nwp,nmax)
            ! find products stochio
            do iw = 1,nwp
                stochio(indexchim(words(iw))) = stochio(indexchim(words(iw))) + 1
            end do
            ! get reaction product stochio and indice
            iloc = 3
            do ispe = 1,nesp
                if (stochio(ispe)==0) cycle
                iloc = iloc + 1
                nindice(iloc) = stochio(ispe)
                iindice(iloc) = ispe
                ! check up: to count the species used in 'reactfile'
                if (specheck(ispe)==0) then
                    specheckp(ispe) = 1
                    specheck(ispe) = 1
                    nbq = nbq + 1
                else if (specheckp(ispe)==0) then
                    specheckp(ispe) = 1
                end if
            end do
            ! fill reaction product stochio and indice
            reactions(ireact)%ip1 = nindice(4)
            reactions(ireact)%p1  = iindice(4)
            reactions(ireact)%ip2 = nindice(5)
            reactions(ireact)%p2  = iindice(5)
            reactions(ireact)%ip3 = nindice(6)
            reactions(ireact)%p3  = iindice(6)
            ! set reaction three prod to true with checking in position 6 if there is three different products
            if (nindice(6)/=0) reactions(ireact)%three_prod = .true.


            !-------------------------!
            !--- fill indice array ---!
            !-------------------------!

            if (trim(reactions(ireact)%vtype)=="v4") then
                if (nindice(6)==0) then ! reaction with 1 or 2 products
                    ! z4spec(ir1,r1,ir2,r2,ip1,p1,ip2,p2)
                    indice_4(nb_reaction_4) = z4spec(nindice(1), iindice(1), nindice(2), iindice(2), nindice(4), iindice(4), nindice(5), iindice(5))
                else ! reaction with 3 different products
                    indice_4(nb_reaction_4) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(4), iindice(4), nindice(5)/2., iindice(5))
                    nb_reaction_4 = nb_reaction_4 + 1
                    indice_4(nb_reaction_4) = z4spec(nindice(1)/2., iindice(1), nindice(2)/2., iindice(2), nindice(6), iindice(6), nindice(5)/2., iindice(5))
                end if
            else if (trim(reactions(ireact)%vtype)=="v3") then
                if (nindice(6)==0) then ! reaction with 1 or 2 products
                    ! z3spec(ir1,r1,ip1,p1,ip2,p2)
                    indice_3(nb_reaction_3) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
                else ! reaction with 3 different products
                    indice_3(nb_reaction_3) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
                    nb_reaction_3 = nb_reaction_3 + 1
                    indice_3(nb_reaction_3) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
                end if
            else if (trim(reactions(ireact)%vtype)=="vphot") then
                if (reactions(ireact)%rtype==0) then
                    if (nindice(6)==0) then ! reaction with 1 or 2 products
                        ! z3spec(ir1,r1,ip1,p1,ip2,p2)
                        indice_phot(nb_phot_hv) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
                    else ! reaction with 3 different products
                        indice_phot(nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
                        nb_phot = nb_phot + 1
                        nb_phot_hv = nb_phot_hv + 1
                        indice_phot(nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
                    end if
                else
                    if (nindice(6)==0) then ! reaction with 1 or 2 products
                        ! z3spec(ir1,r1,ip1,p1,ip2,p2)
                        indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1), iindice(1), nindice(4), iindice(4), nindice(5), iindice(5))
                    else ! reaction with 3 different products
                        indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(4), iindice(4), 0.0, i_dummy)
                        nb_phot = nb_phot + 1
                        indice_phot(nb_phot_hv_max+nb_phot-nb_phot_hv) = z3spec(nindice(1)/2., iindice(1), nindice(5), iindice(5), nindice(6), iindice(6))
                    end if
                end if
            else
                print*,'Error in photochemistry_asis (get_react):'
                print*,'vtype',reactions(ireact)%vtype,' not implemented'
                call abort
            end if
        
          end if ! end if comment line

        end do ! end loop on lines

        close(402)

    end subroutine get_react

end subroutine photochemistry_asis