!> \file kinetics_module.F90
!!  Contains type definitions for MICM kinetics

module kinetics_module

use machine, only: r8 => kind_phys

!*****************************************************************************************
! This module is completely auto-generated by Configurator
!*****************************************************************************************

use ODE_defs_module

implicit none

private
public :: kinetics_type

! k_rateConst are computed at the beginning of the 
!   chemistry_box_solver time step.
!   They are not computed for internal steps of the
!   box-model time-step advancer
! k_rateConst will be thread-safe memory provided elsewhere.
! rate_constant_store will be an accessor to memory
! For now, it is allocated here. It is not thread safe

! Filter with CPP for PGI compiler
#ifndef __PGI
!> \section arg_table_kinetics_type
!! | local_name | standard_name                                    | long_name                               | units   | rank | type          |    kind   | intent | optional |
!! |------------|--------------------------------------------------|-----------------------------------------|---------|------|---------------|-----------|--------|----------|
!! | theKinetics | kinetics_data                                   | chemistry kinetics                      | DDT     |    0 | kinetics_type |           | none   | F        |
!! | nkRxt       | no_gas_phase_reactions                          | ngas_phase_reactions                    | count   |    0 | integer       |           | none   | F        |
!!

#endif

type, extends(baseOdeDefs) ::  kinetics_type
! type  kinetics_type
  integer               :: nkReact
  real(r8), allocatable :: k_rateConst(:)
  real(r8), allocatable :: j_rateConst(:)
contains
  private
  procedure, public :: k_rateConst_init 
  procedure, public :: k_rateConst_run
  procedure, public :: k_rateConst_print
  procedure, public :: forcing => force
  procedure, public :: jacobian => jac
  final             :: DasEnder
end type kinetics_type

contains

  !------------------------------------------------------
  ! allocate k_rateConst member array
  !------------------------------------------------------
  subroutine k_rateConst_init( this, nkRxt )
      
    class(kinetics_type) :: this
    integer, intent(in)  :: nkRxt          ! total no. gas phase rxtions

    this%nkReact = nkRxt
    if( .not. allocated(this%k_rateConst) ) then
      allocate( this%k_rateConst(nkRxt) )
    else
      write(*,*) 'k_rateConst_init: k_rateConst already allocated'
    endif

  end subroutine k_rateConst_init

  !------------------------------------------------------
  ! Compute k_rateConst
  ! Execute once for the chemistry-time-step advance
  ! Not called from the solver
  !------------------------------------------------------
  subroutine k_rateConst_run( this )
  
  class(kinetics_type) :: this

  ! Rate Constants
  ! Y0_a
  this%k_rateConst(1) = 0.04_r8
  ! Y1_Y2_M_b
  this%k_rateConst(2) = 1.e4_r8
  ! Y1_Y1_a
  this%k_rateConst(3) = 1.5e7_r8
  
  end subroutine k_rateConst_run

  subroutine k_rateConst_print( this )

    class(kinetics_type) :: this

    write(*,*) 'rate constants'
    write(*,'(1p,5(1x,g0))') this%k_rateConst(:)

  end subroutine k_rateConst_print

  !---------------------------
  !  cleanup when k_rateConst type is removed
  !---------------------------
  subroutine DasEnder( this )

  type(kinetics_type) :: this

  if( allocated( this%k_rateConst ) ) then
    deallocate( this%k_rateConst )
  endif

  end subroutine DasEnder
  
  !---------------------------
  ! Compute time rate of change of each molecule (vmr) given reaction rates
  !---------------------------
  function force( this, data )
  
    class(kinetics_type) :: this
    real(r8), intent(in)::  data(:)              ! volume mixing ratios of each component in order
    real(r8)            ::  force(size(data))    ! rate of change of each molecule

    real(r8) ::  k_rates(this%nkReact)          ! rates of each reaction
  
    call compute_rates(this, data, k_rates)
 
   ! set the forcing array
    force(1) = k_rates(2) - k_rates(1)
    force(2) = k_rates(1) - k_rates(2) - 2._r8 * k_rates(3)
    force(3) = 2._r8 * k_rates(3)

  end function force
  
  !---------------------------
  ! Compute sensitivity of molecular forcing to each vmr (derivative of force w.r.t. each vmr)
  !---------------------------
  function jac( this, data )
  
    class(kinetics_type) :: this
    real(r8), intent(in)::  data(:)         ! volume mixing ratios of each component in order
    real(r8) :: jac(size(data),size(data))   ! sensitivity of forcing to changes in each vmr
  
    jac(:,:) = 0._r8
  
    ! Jacobian  
    jac(1,1)  = jac(1,1) - this%k_rateConst(1)
    jac(1,2)  = jac(1,2) + this%k_rateConst(2) * data(3)
    jac(1,3)  = jac(1,3) + this%k_rateConst(2) * data(2)

    jac(2,1)  = jac(2,1) + this%k_rateConst(1)
    jac(2,2)  = jac(2,2) - this%k_rateConst(2) * data(3)
    jac(2,3)  = jac(2,3) - this%k_rateConst(2) * data(2)
    jac(2,2)  = jac(2,2) - 4._r8 * this%k_rateConst(3) * data(2)

    jac(3,2)  = jac(3,2) + 4._r8 * this%k_rateConst(3) * data(2)

  end function jac

  !---------------------------
  ! Compute reaction rates, given vmr of each species and rate constants
  !---------------------------
  subroutine compute_rates(this, vmr, k_rates)
  
    class(kinetics_type)  :: this
    real(r8), intent(in)  ::  vmr(:)           ! volume mixing ratios of each component in order
    real(r8), intent(out) ::  k_rates(:)   ! rates for each reaction (sometimes called velocity of reaction)
  
! Rates
   ! Y0_a
   k_rates(1) = this%k_rateConst(1) * vmr(1)
   ! Y1_Y2_M_b
   k_rates(2) = this%k_rateConst(2) * vmr(2) * vmr(3)
   ! Y1_Y1_a
   k_rates(3) = this%k_rateConst(3) * vmr(2) * vmr(2)

  end subroutine compute_rates
  
end module kinetics_module