!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: USR_RATES                                              !
!                                                                      !
!  Purpose: Hook for user defined reaction rates.                      !
!                                                                      !
!  Author: J.Musser                                   Date: 10-Oct-12  !
!                                                                      !
!  Comments: Write reaction rates in amount of chemical substance      !
!  per time. For CGS this has units mol/s. For _MFIX_ SI this has      !
!  units kmol/s.                                                       !
!                                                                      !
!  WARNING: Only discrete phase reactions should be specifed here.     !
!  Homogeneous gas phase reactions in DEM simulations must be given    !
!  in the continuum reaction hook (usr_rates.f).                       !
!                                                                      !
!  The call to usr_rates_des is made from inside a particle loop which !
!  is nested inside an IJK loop. Fluid grid calculations independent   !
!  of particle properties can be carried out in des/usr4_des.f to      !
!  reduce redundant calculations.                                      !
!                                                                      !
!  Example: Evaporation                                                !
!                                                                      !
!  mfix.dat input:                                                     !
!``````````````````````````````````````````````````````````````````````!
!    @(DES_RXNS)                                                       !
!      Evap { chem_eq = "Liquid --> Vapor" }                           !
!    @(DES_END)                                                        !
!``````````````````````````````````````````````````````````````````````!
!                                                                      !
!  des/usr_rates_des.f input:                                          !
!``````````````````````````````````````````````````````````````````````!
!    Sa = Pi*(2.0d0*DES_RADIUS(NP))**2    ! Particle surface area      !
!    H2O_xfr = ...  ! An expression for mass transfer coefficient      !
!    Cmg_H2O = ...  ! Molar concentration grad of water vapor          !
!    DES_RATES(EVAP) = Sa * H2O_xfr * Cmg_H2O                          !
!``````````````````````````````````````````````````````````````````````!
!  * Species alias and reaction names given in the data file can be    !
!    used in reference to the reaction index in DES_RATES and a        !
!    species index in gas/solids phase variables.                      !
!                                                                      !
!  * Additional information is provided in section 4.11 of the code    !
!    Readme.                                                           !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE USR_RATES_DES(NP, pM, IJK, DES_RATES)

      USE compar
      USE constant
      USE des_thermo
      USE des_rxns
      USE discretelement
      USE energy
      USE fldvar
      USE funits
      USE geometry
      USE indices
      USE param
      USE param1
      USE physprop
      USE rxns
      USE run
      USE usr
      USE fun_avg
      USE functions

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NP  ! Global index of particle
      INTEGER, INTENT(IN) :: pM  ! Solid phase index of particle NP
      INTEGER, INTENT(IN) :: IJK ! Fluid cell index containing NP

! Calculated reaction rates. (reacted moles per sec)
      DOUBLE PRECISION, INTENT(OUT) :: DES_RATES(NO_OF_DES_RXNS)


! Reaction specific variables:
!`````````````````````````````````````````````````````````````````````//

!     Uncomment the following line when customizing this usr_rates_des.f to include generated file species.inc:
!
      INCLUDE 'species.inc'


! Reaction rates:
!`````````````````````````````````````````````````````````````````````//
! Include reaction rates here. Reaction rates should be stored in the
! variable DES_RATES. The reaction name given in the data file can be
! used to store the rate in the appropriate array location. Additional
! input format parameters are given in Section 4.11 of the code Readme.
! Calculation variables:
!`````````````````````````````````````````````````````````````````````//
! Bounded phase temperatures (K)
      DOUBLE PRECISION :: TgX   ! gas phase
      DOUBLE PRECISION :: TsX   ! solids phase 1
      DOUBLE PRECISION :: TgsX  ! Average of gas/solids 1     

      DOUBLE PRECISION :: X_gTar
      DOUBLE PRECISION :: X_sWood

      ! Prevent reactions when mass fraction falls below c_Limiter.
      DOUBLE PRECISION, parameter :: c_Limiter = 1.0d-6
! Reaction specific variables:
!`````````````````````````````````````````````````````````````````````//
!A = pre-factor (1/s) and E = activation energy (J/mol)
      DOUBLE PRECISION, PARAMETER :: A1 = 1.08e11,   E1 = 148.0e3
      DOUBLE PRECISION, PARAMETER :: A2 = 4.38e9,    E2 = 152.7e3
      DOUBLE PRECISION, PARAMETER :: A3 = 3.75e6,    E3 = 111.7e3
      DOUBLE PRECISION, PARAMETER :: A4 = 1.00e5,    E4 = 108.0e3
      DOUBLE PRECISION, PARAMETER :: A5 = 4.28e6,    E5 = 108.0e3
      DOUBLE PRECISION, PARAMETER :: Adrying = 5.13e6,    Edrying =  87.9e3
      DOUBLE PRECISION :: K1, K2, K3, K4, K5, Kdrying
      DOUBLE PRECISION, PARAMETER :: RGAS_LOCAL = 8.314472 !J/(K.mol)
      DOUBLE PRECISION, PARAMETER :: TMAX = 2000.d0
      DOUBLE PRECISION :: woodMole, H2OMole, scale_oil, scale_gas, scale_char
! calculate variables for reactions:
!`````````````````````````````````````````````````````````````````````//

! skip sands
      if(pM==2) then
         DES_RATES(:) = 0.d0
         des_usr_var(1:4,NP) = 0.d0
	 !print *, 'llu skip sand' 
         return     
      endif
      
      
! Bound the gas and solids phase temperatures.
      TgX  = min(TMAX, T_g(IJK))
      TsX  = min(TMAX, DES_T_s(NP))


      !TsX = -0.327538*time**4 + 7.846755*time**3 - 68.227705*time**2 + 263.972943*time + 341.367649
      DES_T_s(NP) = TsX
      TgX = TsX
      TgsX = HALF * (TgX + TsX)
      ! Alias the wood/biooil phase volume fractions
 !-------------- drying   -----------
      H2OMole  = DES_X_s(NP,H2O_L)*pmass(NP)/MW_s(1,H2O_L) 
      Kdrying = Adrying * Exp(-Edrying/(RGAS_LOCAL*TsX))
      DES_RATES(W2Char) = Kdrying * H2OMole
!------------------------------------------------------------------------------   

!-------------- Biomass to biooil,gas, and char (g) conversion   -----------
      X_sWood  = DES_X_s(NP,wood)
      woodMole  = DES_X_s(NP,wood)*pmass(NP)/MW_s(1,wood) 
      IF( woodMole .LE. ZERO) THEN
        DES_RATES(W2BIOOIL)   = ZERO
        DES_RATES(W2Gas)   = ZERO
        DES_RATES(W2Char) = ZERO
      ELSE
        scale_char  = 1.05
        scale_oil  = 0.082
        scale_gas  = 2.5
        !K1 = A1 * Exp(-E1/(RGAS_LOCAL*T_s(IJK,1)))
        K1 = A1 * Exp(-E1/(RGAS_LOCAL*TsX))
        DES_RATES(W2BIOOIL) = K1 * woodMole *scale_oil

        K2 = A2 * Exp(-E2/(RGAS_LOCAL*TsX)) *scale_gas
        DES_RATES(W2Gas) = K2 * woodMole

        K3 = A3 * Exp(-E3/(RGAS_LOCAL*TsX)) *scale_char
        DES_RATES(W2Char) = K3 * woodMole
      END IF
!------------------------------------------------------------------------------
      X_gTar = X_g(IJK,biooil)
      IF( X_g(IJK,biooil) .LE. c_Limiter) THEN
        DES_RATES(BIOOIL2CHAR) = ZERO
      ELSE
        ! Tar to char conversion
        K4 = A4 * Exp(-E4/(RGAS_LOCAL*TgX))
        DES_RATES(BIOOIL2CHAR) = K4 * Rop_g(IJK)*X_g(IJK,biooil)/MW_g(biooil)*PVOL(NP)/EP_S(IJK,1)
      END IF
 
      if(EP_S(IJK,1)>0.d0) then
 
      else
         DES_RATES(:) = 0.d0
          print *, 'lludebug, ',__FILE__,ep_s(ijk,1)      
      endif

! print *, 'llurates', DES_RATES(:),EP_S(IJK,1),Rop_s(IJK,1),X_s(IJK,1,wood)
!  Write the reaction rate data into .SPA file and for visualizing.
      des_usr_var(1,NP) = DES_RATES(W2BIOOIL)
      des_usr_var(2,NP) = DES_RATES(W2Gas)
      des_usr_var(3,NP) = DES_RATES(W2Char)
      des_usr_var(4,NP) = DES_RATES(BIOOIL2CHAR)            
      RETURN

      END SUBROUTINE USR_RATES_DES
