#include "dims.h"
      module crates_module
      implicit none 
#include "params.h"
!
! Chemical reaction rates.
! Temperature independent rates are global real parameters,
!   (module data), set at compile time.
! Temperature dependent rates are dimensioned (nlon,nlev) and are
!   in thread private common (crates_tdep.h). Dependent rates are 
!   set by sub rates_tdep at every latitude at each time step. 
!
! Temperature independent global rate parameters:
      real,parameter ::
     |  RK4   = 1.0E-10,
     |  RK5   = 4.4E-10,
     |  RK6   = 4.0E-10,
     |  RK7   = 2.0E-10,
     |  RK8   = 1.0E-12,
     |  RK9   = 6.0E-11,
     |  RK10  = 1.3E-10,
     |  RK16  = 4.8E-10,
     |  RK17  = 1.0E-10,
     |  RK18  = 4.0E-10,
     |  RK21  = 0.047,
     |  RK22  = 0.171,
     |  RK23  = 8.E-10,
     |  RK24  = 5.0E-12,
     |  RK26  = 7.E-10,
     |  RK27  = 7.7E-5,
     |  BETA2 = 5.0E-12,
!    |  BETA4 = 6.9E-13,
     |  BETA4 = 5.0E-13,
     |  BETA6 = 7.0E-11,
     |  BETA7 = 1.06E-5
      contains
!-----------------------------------------------------------------------
      subroutine chemrates_tdep
!
! Time-gcm wrapper for rates_tdep (temperature-dependent chemical rate 
!   coefficients, called at all latitudes and every time step).
! This is called from parallel latitude loop (dynamics.f)
!
#include "fcom.h"
#include "index.h"
#include "buff.h"
!     logical,save :: isfirst=.true.
!
!     if (isfirst) then
!       write(6,"('chemrates_tdep: calling rates_tdep.')")
!       isfirst = .false.
!     endif
!
! Note zkmx means top level zkmxp of rates is not defined:
      call rates_tdep(f(1,nj+nt),f(1,nj+nte),f(1,nj+nti),f(1,nno2),
     |  f(1,nnvo2),zimxp,zkmx)
      end subroutine chemrates_tdep
!-----------------------------------------------------------------------
      subroutine rates_tdep(tn,te,ti,fno2,fnvo2,nlon,nlev)
      use input_module,only: f107  ! 10.7 cm flux (from input and/or gpi)
      use init_module,only: sfeps  ! flux variation from orbital excentricity
!
! Set temperature-dependent reaction rates in thread-private common in
!   crates_tdep.h. This routine is called from dynamics.F, from inside
!   multi-tasked latitude scan.
!
#include "crates_tdep.h"
!
! Args:
      integer,intent(in) ::
     |  nlon,                ! number of longitudes (# of columns)
     |  nlev                 ! number of vertical levels in each column
      real,intent(in) ::
     |  tn(nlon,nlev),       ! neutral temperature (deg K)
     |  te(nlon,nlev),       ! electron temperature (deg K)
     |  ti(nlon,nlev),       ! ion temperature (deg K)
     |  fno2(nlon,nlev+1),   ! o2 line integral (see chapmn.F)
     |  fnvo2(nlon,nlev+1)   ! o2 column number density (see chapmn.F)
!
! Local:
      integer :: i,k
      real :: ti1(nlon,nlev), ti2(nlon,nlev), ti3(nlon,nlev),
     |        etvib(nlon,nlev)
!
! Loop over number of columns (nlon), and number of levels in each 
! column:
      do k=1,nlev
        do i=1,nlon
!
! ti1 = T1/300. (S15)    ti2 = T2/300. (S14)   ti3 = TR/300. (S13)
!
          ti1(i,k) = (0.667* ti(i,k)+0.333 *tn(i,k))/300.
          ti2(i,k) = (0.6363*ti(i,k)+0.3637*tn(i,k))/300.
          ti3(i,k) = .5*(ti(i,k)+tn(i,k))/300.
!
! Rate coefficient for O+ + O2 -> O + O2+ (HIERL)
          rk1(i,k)=1.6E-11*ti1(i,k)**(-0.52)+5.5E-11*
     |      exp(-22.85/ti1(i,k))
!
! k2 assumes t2/300 in ti2:
          rk2(i,k)=(8.6E-14*ti2(i,k)-5.92E-13)*ti2(i,k)+1.533E-12
!
! Rate coefficient for O+ + N2 -> N + NO+
! See also sub altv (altv.f). Tvib is output of sub altv, if it
! is called (see dynamics.f).
!
          tvib(i,k) = tn(i,k)                ! tvib in crates_tdep.h
          etvib(i,k) = exp(-3353./tvib(i,k)) ! etvib is local
          rk2(i,k) = (((((270.*etvib(i,k)+220.)*etvib(i,k)+85.)*
     |      etvib(i,k)+38.)*etvib(i,k)+1.)*rk2(i,k)*etvib(i,k)+
     |      rk2(i,k))*(1.-etvib(i,k))
!
          if (300.*ti3(i,k) >= 1500.) then
            rk3(i,k)=5.2E-11*ti3(i,k)**0.2
          else
            rk3(i,k)=1.4E-10*ti3(i,k)**(-0.44)
          endif
          rk19(i,k) = 4.0E-8*sqrt(300./te(i,k))
          rk20(i,k) = 1.5E-7*sqrt(300./te(i,k))
          rk25(i,k) = 6.6E-8*sqrt(300./te(i,k))
          rkm12(i,k) = 9.59E-34*exp(480./tn(i,k))
          ra1(i,k)=4.2E-7*(300./te(i,k))**0.85
          if (te(i,k) >= 1200.) then
            ra2(i,k)=1.6E-7*(300./te(i,k))**0.55
          else
            ra2(i,k)=2.7E-7*(300./te(i,k))**0.7
          endif
          ra3(i,k)=1.8E-7*(300./te(i,k))**0.39
          beta1(i,k) = 1.5E-11*exp(-3600./tn(i,k))
!         beta3(i,k) = 1.6E-10*exp(-460./tn(i,k))
          beta3(i,k) = 2.5E-10*sqrt(tn(i,k)/300.)*exp(-600./tn(i,k))
          beta5(i,k) = 3.6E-10*sqrt(te(i,k)/300.)
          beta8(i,k) = 4.5E-6*(1.+0.11*(f107-65.)/165.)*
     |      exp(-1.E-8*fno2(i,k)**0.38)*sfeps
          beta9(i,k)=2.91E11*(1.+0.2*(f107-65.)/100.)*2.E-18*
     |      exp(-8.E-21*fno2(i,k))*sfeps
          beta9n(i,k)=5.E9*(1.+0.2*(f107-65.)/100.)*2.E-18*
     |      exp(-8.E-21*fnvo2(i,k))*sfeps
          beta17(i,k)=1.E-32*sqrt(300./te(i,k))
        enddo ! i=1,nlon
      enddo ! k=1,nlev
!
! Top level (nlev+1 == zkmxp):
      do i=1,nlon
        beta8(i,nlev+1)=4.5E-6*(1.+0.11*(f107-65.)/165.)*
     |    exp(-1.E-8*fno2(i,nlev+1)**0.38)
        beta9(i,nlev+1)=2.91E11*(1.+0.2*(f107-65.)/100.)*2.E-18*
     |    exp(-8.E-21*fno2(i,nlev+1))
        beta9n(i,nlev+1)=5.E9*(1.+0.2*(f107-65.)/100.)*2.E-18*
     |    exp(-8.E-21*fnvo2(i,nlev+1))
      enddo
      end subroutine rates_tdep
      end module crates_module