!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: USR_RATES                                              !
!                                                                      !
!  Purpose:                                            !
!                                                                      !
!  Author: M.Denison                                  Date: 26-Oct-18  !
!                                                                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE USR_RATES(IJK, RATES)

      USE compar
      USE constant
      USE energy
      USE fldvar
      USE fun_avg
      USE functions
      USE funits
      USE geometry
      USE indices
      USE parallel
      USE param
      USE param1
      USE physprop
      USE run
      USE rxns
      USE sendrecv
      USE toleranc


!      USE usr

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: IJK
	  INTEGER ir

      DOUBLE PRECISION, DIMENSION(NO_OF_RXNS), INTENT(OUT) :: RATES

! Reaction specific variables:
!`````````````````````````````````````````````````````````````````````//
! Bounded phase temperatures (K)
      DOUBLE PRECISION xTg   ! gas phase
      DOUBLE PRECISION xTs   ! solids phase
! Gas phase pressure in atm
!      DOUBLE PRECISION Pg_KPa
      DOUBLE PRECISION Pg_atm
! Gas phase concentrations (mol/m^3)
      DOUBLE PRECISION c_CO, c_O2, c_CO2, c_H2, c_H2O, c_CH4
! Partial pressures (atm)
      DOUBLE PRECISION p_O2, p_H2O, p_H2, p_CO2, p_CO, p_CH4
! Total solids surface area per volume (1/m)
      DOUBLE PRECISION Sa
! Diffusion coefficient for SiH2 (m^2/sec)
!     DOUBLE PRECISION DIFF_SiH2
! SiH2 gas constant (KPa.m^3)/(kg-SiH2.K)
!      DOUBLE PRECISION R_SiH2
! Calculate film diffusion resistance. (kg-SiH2/(KPa.sec.m^2)
!      DOUBLE PRECISION K_f
! Apparent first order heterogeneous reation rate constant  (m/sec)
!      DOUBLE PRECISION k_so
! Constant of hydroggen inhibition of heterogeneous reaction (1/KPa)
!      DOUBLE PRECISION K_H
! Rate constants
      DOUBLE PRECISION k13, k14, k15, k16, kox, kdry, kdevol, k1
	  DOUBLE PRECISION k6, k7, k8, k9, k10, k11
! fraction of Carbon converted to CO2
      DOUBLE PRECISION psi
! Temporary values for reversible reactions
      DOUBLE PRECISION FWD, RVS, NET
!	  INTEGER M


! Prevent reactions when mass fraction falls below c_Limiter.
      DOUBLE PRECISION, parameter :: c_Limiter = 1.0d-7

!`````````````````````````````````````````````````````````````````````//
      INCLUDE 'species.inc'
!`````````````````````````````````````````````````````````````````````//

! Calculate the bounded temperatures.
      xTg = max(min(T_g(IJK),   TMAX), 300.0d0)
      xTs = max(min(T_s(IJK,1), TMAX), 300.0d0)

! Convert gas pressure to atm
!      Pg_KPa = 1.0d-3 * P_g(IJK)  !  KPa
      Pg_atm = P_g(IJK) / 101325.0  !  atm

! Initialize partial pressures and concentrations.
	  c_CO = ZERO; c_O2 = ZERO; c_CO2 = ZERO; c_H2 = ZERO;
	  c_H2O = ZERO; c_CH4=ZERO;

! Calculate the partial pressure and molar concentration of CO
      IF(X_g(IJK,CO) > c_Limiter) then
         p_CO = Pg_atm * X_g(IJK,CO) * (MW_MIX_g(IJK) / MW_g(CO))
         c_CO = RO_g(IJK) * X_g(IJK,CO) / Mw_g(CO)
      ENDIF

! Calculate the partial pressure and molar concentration of O2
      IF(X_g(IJK,O2) > c_Limiter) then
         p_O2 = Pg_atm * X_g(IJK,O2) * (MW_MIX_g(IJK) / MW_g(O2))
         c_O2 = RO_g(IJK) * X_g(IJK,O2) / Mw_g(O2)
      ENDIF

! Calculate the partial pressure and molar concentration of CO2
      IF(X_g(IJK,CO2) > c_Limiter) then
         p_CO2 = Pg_atm * X_g(IJK,CO2) * (MW_MIX_g(IJK) / MW_g(CO2))
         c_CO2 = RO_g(IJK) * X_g(IJK,CO2) / Mw_g(CO2)
      ENDIF

! Calculate the partial pressure and molar concentration of H2
      IF(X_g(IJK,H2_REF_ELEMENT) > c_Limiter) then
         p_H2 = Pg_atm * X_g(IJK,H2_REF_ELEMENT) * (MW_MIX_g(IJK) / MW_g(H2_REF_ELEMENT))
         c_H2 = RO_g(IJK)*X_g(IJK,H2_REF_ELEMENT)/Mw_g(H2_REF_ELEMENT)
      ENDIF

! Calculate the partial pressure and molar concentration of H2O
      IF(X_g(IJK,H2O) > c_Limiter) then
         p_H2O = Pg_atm * X_g(IJK,H2O) * (MW_MIX_g(IJK) / MW_g(H2O))
         c_H2O = RO_g(IJK) * X_g(IJK,H2O) / Mw_g(H2O)
      ENDIF

! Calculate the partial pressure and molar concentration of CH4
      IF(X_g(IJK,CH4) > c_Limiter) then
         p_CH4 = Pg_atm * X_g(IJK,CH4) * (MW_MIX_g(IJK) / MW_g(CH4))
         c_CH4 = RO_g(IJK) * X_g(IJK,CH4) / Mw_g(CH4)
      ENDIF

!**********************************************************************!
!                         Homogeneous Reactions                        !
!**********************************************************************!

! CO + 1/2 O2 --> CO2  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
	  IF(xTg < 1150.0) then
	     FWD = 2.61d12 * exp(-2.2929d4/xTg) * EP_g(IJK) * c_CO * c_H2O**0.5 * c_O2**0.25
	  ELSE
	     FWD = 6.52d6 * exp(-8.0352d3/xTg) * EP_g(IJK) * c_CO * c_H2O**0.5 * c_O2**0.25	  
	  ENDIF

      RATES(CO_oxidation) =  FWD

! CO + H2O <--> CO2 + H2 (kg-mole/m^3.s)
!---------------------------------------------------------------------//
      FWD = 7.68d10 * exp(-3.664d4/xTg) * EP_g(IJK) * sqrt(c_CO)*c_H2O
      RVS = 6.4d9 * exp(-3.926d4/xTg) * EP_g(IJK) * sqrt(c_H2) * c_CO2

      NET = FWD - RVS
      IF(NET >= ZERO) THEN
         RATES(WGS_forward) =  NET
         RATES(WGS_reverse) =  ZERO
      ELSE
         RATES(WGS_forward) =  ZERO
         RATES(WGS_reverse) = -NET
      ENDIF

! H2 + 1/2 O2 --> H2O  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
      FWD = 1.8d13 * exp(-1.7614d4/xTg) * EP_g(IJK) * c_H2 * sqrt(c_O2)

      RATES(H2_oxidation) =  FWD

! CH4 + 2 O2 --> CO2 + 2 H2O  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!      FWD = 1.3d8 * exp(-2.4367d4/xTg) * EP_g(IJK) * c_O2**1.3 * c_CH4**(-0.3)
      FWD = 1.3d8 * exp(-2.4367d4/xTg) * EP_g(IJK) * c_O2**1.3 * c_CH4**(0.3)
!	  print*,EP_g(IJK),c_O2,c_CH4

      RATES(CH4_oxidation) =  FWD

! C_g + 1/2 O2 --> CO2  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
      FWD = 1.0d10 * exp(-1.7614d4/xTg) * EP_g(IJK) * c_O2

      RATES(C_oxidation) =  FWD


!**********************************************************************!
!                        Heterogeneous Reactions                       !
!**********************************************************************!
      IF(EP_s(IJK,1) > ZERO_EP_s) THEN

! Cu2O + 1/2 O2 --> 2 CuO  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(Cu2Ooxidation)%Species(CU2O)%pmap
         k13 = 6.3145d13 * exp(-4.1497d4/xTg) * EP_g(IJK) * ROP_s(IJK,1) * X_s(IJK,1,CU2O)
	     k14 = 5.6d3 * exp(-8179.0/xTg) *EP_g(IJK) *ROP_s(IJK,1) * X_s(IJK,1,CU2O)
! Reaction rate:
         RATES(Cu2Ooxidation) = (k13 - k14 * p_O2**1.3)/MW_s(1,CU2O)

! CuO --> 1/2 Cu2O+ 1/4 O2  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(CuOdecomp)%Species(CUO)%pmap
         k15 = 3.0 * exp(-2350.0/xTg) * EP_g(IJK) * ROP_s(IJK,1) * X_s(IJK,1,CUO)
	     k16 = 1.62d8 * exp(-2.795d4/xTg) *EP_g(IJK) *ROP_s(IJK,1) * X_s(IJK,1,CUO)
! Reaction rate:
         RATES(CuOdecomp) = (k15 * p_O2 - k16)/MW_s(1,CUO)
		 
	   ENDIF

	   IF(EP_s(IJK,2) > ZERO_EP_s) THEN
! Drying H2O_L_1 --> H2O  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(Drying)%Species(H2O_L_1)%pmap
         kdry = 1.3d7 * exp(-20133.0/xTg) * ROP_s(IJK,2) * X_s(IJK,2,H2O_L_1)
! Reaction rate:
         RATES(Drying) = kdry/MW_s(2,H2O_L_1)

! Coal --> Volatiles  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(VolatileRelease)%Species(Coal)%pmap
         kdevol = 4.28d14 * exp(-27477.3/xTg) * ROP_s(IJK,2) * X_s(IJK,2,Coal)
! Reaction rate:
         RATES(VolatileRelease) = 0.536068 * kdevol / MW_s(2,Coal)

! Coal --> Char  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
! Reaction rate:
         RATES(char_production) = 0.463932 * kdevol / MW_s(2,Coal)

! CGR_REF_ELEMENT + O2 --> CO2  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(Char_C_oxidation)%Species(CGR_REF_ELEMENT)%pmap

! Total surface area of solids in computational cell (m^2/m^3)
         Sa = 6.0d0 * EP_s(IJK,2) / D_p0(2)   ! (1/m)
!         Sa = 6.0d0 * EP_s(IJK,M) / 0.0001   ! (1/m)

         kox = 1.45d6 * exp(-10049.3/xTg) * EP_g(IJK) * Sa / MW_s(2,Coal)
		 
! fraction of carbon converted to CO2
         k1 = 200.0 * exp(-4529.4/xTg)
		 
		 psi = 1.0 / (1.0 + k1)

! Reaction rate:
         RATES(Char_C_oxidation) = 0.51926 * psi * kox * p_O2
!		 print*,RATES(Char_C_oxidation),psi,kox,p_O2,Sa,EP_s(IJK,2),MW_s(2,Coal)

! CGR_REF_ELEMENT + 1/2 O2 --> CO  (kg-mole/m^3.s)
!---------------------------------------------------------------------//

! Reaction rate:
         RATES(Char_to_CO_oxidation) = 0.51926 * (1.0 - psi) * kox * p_O2
		 
! Other char oxidation
         RATES(Char_other_oxidation) = 0.48074 * kox * p_O2

! CGR_REF_ELEMENT + H2O --> CO + H2_REF_ELEMENT  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(steam_gasification)%Species(CGR_REF_ELEMENT)%pmap
		 
         k6 = 7.0d1 * exp(-15096.6/xTg) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT)
! Reaction rate:
         RATES(steam_gasification) = k6 * p_H2O /MW_s(2,CGR_REF_ELEMENT)

! CO + H2_REF_ELEMENT --> CGR_REF_ELEMENT + H2O  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(steam_gasification_reverse)%Species(CGR_REF_ELEMENT)%pmap
		 
         k7 = 7.0d1 * exp(1233.4/xTg - 17.29) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT)
! Reaction rate:
         RATES(steam_gasification_reverse) = k7 * p_H2 * p_H2O /MW_s(2,CGR_REF_ELEMENT)

! Reaction rate:
         RATES(Char_to_CO_oxidation) = 0.51926 * (1.0 - psi) * kox * p_O2/MW_s(2,CGR_REF_ELEMENT)

! CGR_REF_ELEMENT + CO2 --> 2 CO (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(CO2_gasification)%Species(CGR_REF_ELEMENT)%pmap
		 
         k8 = 7.0d1 * exp(-15096.6/xTg) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT)
! Reaction rate:
         RATES(CO2_gasification) = k8 * p_CO2 /MW_s(2,CGR_REF_ELEMENT)

! 2 CO --> CGR_REF_ELEMENT + CO2  (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!        M = Reaction(CO2_gasification_reverse)%Species(CGR_REF_ELEMENT)%pmap
		 
         k9 = 7.0d1 * exp(5185.4/xTg - 20.92) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT)
! Reaction rate:
         RATES(CO2_gasification_reverse) = k9 * p_CO * p_CO /MW_s(2,CGR_REF_ELEMENT)

! CGR_REF_ELEMENT + 2 H2 --> CH4 (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(Char_methanation)%Species(CGR_REF_ELEMENT)%pmap
		 
         k10 = exp(-8078.0/xTg-7.087) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT)
! Reaction rate:
         RATES(Char_methanation) = k10 * p_CO2 /MW_s(2,CGR_REF_ELEMENT)

! CH4 --> CGR_REF_ELEMENT + 2 H2 (kg-mole/m^3.s)
!---------------------------------------------------------------------//
!         M = Reaction(Char_methanation_reverse)%Species(CGR_REF_ELEMENT)%pmap
		 
         k11 = exp(-13578.0/xTg-0.372) * EP_g(IJK) * ROP_s(IJK,2) * X_s(IJK,2,CGR_REF_ELEMENT)
! Reaction rate:
         RATES(Char_methanation_reverse) = k11 * sqrt(p_CH4) / MW_s(2,CGR_REF_ELEMENT)

      END IF
  
!         do ir=1,20
!            print*,'RATES ',ir,RATES(ir)
!         enddo
! Save reaction rates for visualization
!      IF(nrr>=RX1F) ReactionRates(IJK,RX1F) = RATES(RX1F)
!      IF(nrr>=RX1R) ReactionRates(IJK,RX1R) = RATES(RX1R)
!      IF(nrr>=RX2F) ReactionRates(IJK,RX2F) = RATES(RX2F)
!      IF(nrr>=RX2R) ReactionRates(IJK,RX2R) = RATES(RX2R)
!      IF(nrr>=RX3)  ReactionRates(IJK,RX3)  = RATES(RX3)
!      IF(nrr>=RX4)  ReactionRates(IJK,RX4)  = RATES(RX4)

      RETURN

      END SUBROUTINE USR_RATES
