added support for pressure-dependent-Arrhenius reactions

photochem/cython/PhotochemWrk.pyx

@@ -138,6 +138,17 @@ cdef class PhotochemWrk:
       cdef ndarray arr = np.empty((dim1, dim2), np.double, order="F")
       wrk_pxd.photochemwrk_densities_get(&self._ptr, &dim1, &dim2, <double *>arr.data)
       return arr
+
+  property rx_rates:
+    """ndarray[double,dim=2], shape (nz,nrT). Reaction rate constants in various units
+    involving molecules cm^3 and s. These rates include 3rd body contributions.
+    """
+    def __get__(self):
+      cdef int dim1, dim2
+      wrk_pxd.photochemwrk_rx_rates_get_size(&self._ptr, &dim1, &dim2)
+      cdef ndarray arr = np.empty((dim1, dim2), np.double, order="F")
+      wrk_pxd.photochemwrk_rx_rates_get(&self._ptr, &dim1, &dim2, <double *>arr.data)
+      return arr
 
   property mubar:
     """ndarray[double,dim=1], shape (nz). The mean molar mass of each atmospheric layer

photochem/cython/PhotochemWrk_pxd.pxd

@@ -33,6 +33,9 @@ cdef extern void photochemwrk_density_get(void *ptr, int *dim1, double *arr)
 cdef extern void photochemwrk_densities_get_size(void *ptr, int *dim1, int *dim2)
 cdef extern void photochemwrk_densities_get(void *ptr, int *dim1, int *dim2, double *arr)
 
+cdef extern void photochemwrk_rx_rates_get_size(void *ptr, int *dim1, int *dim2)
+cdef extern void photochemwrk_rx_rates_get(void *ptr, int *dim1, int *dim2, double *arr)
+
 cdef extern void photochemwrk_mubar_get_size(void *ptr, int *dim1)
 cdef extern void photochemwrk_mubar_get(void *ptr, int *dim1, double *arr)
 

photochem/fortran/PhotochemWrk_wrapper.f90

@@ -176,6 +176,24 @@ subroutine photochemwrk_densities_get(ptr, dim1, dim2, arr) bind(c)
     arr = wrk%densities
   end subroutine
 
+  subroutine photochemwrk_rx_rates_get_size(ptr, dim1, dim2) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(out) :: dim1, dim2
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    dim1 = size(wrk%rx_rates,1)
+    dim2 = size(wrk%rx_rates,2)
+  end subroutine
+
+  subroutine photochemwrk_rx_rates_get(ptr, dim1, dim2, arr) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(in) :: dim1, dim2
+    real(c_double), intent(out) :: arr(dim1, dim2)
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    arr = wrk%rx_rates
+  end subroutine
+
   subroutine photochemwrk_mubar_get_size(ptr, dim1) bind(c)
     type(c_ptr), intent(in) :: ptr
     integer(c_int), intent(out) :: dim1

photochem/utils/_convert_cantera.py

@@ -34,6 +34,11 @@ def photochem2cantera_main(data):
             rxt = "elementary"
         else:
             rxt = data['reactions'][i]['type']
+
+        if rxt == 'falloff':
+            if ' + M ' in data['reactions'][i]['equation']:
+                data['reactions'][i]['equation'] = data['reactions'][i]['equation'].replace('+ M','(+ M)')
+
         if rxt != 'photolysis':
             reactions_new.append(data['reactions'][i])
 
@@ -63,7 +68,7 @@ def photochem2cantera_main(data):
             data['species'][i]['thermo']['temperature-ranges'] = data['species'][i]['thermo']['temperature-ranges'][:3]
             data['species'][i]['thermo']['data'] = data['species'][i]['thermo']['data'][:2] 
 
-    units = flowmap({"length": "cm", "quantity": "molec", "activation-energy": "K"})
+    units = flowmap({"length": "cm", "quantity": "molec", "activation-energy": "K", "pressure": "dyn/cm^2"})
     phases = [{}]
     phases[0]["name"] = "gas name"
     phases[0]["thermo"] = 'ideal-gas'

photochem/utils/_format.py

@@ -96,7 +96,8 @@ def FormatReactions_main(data):
     # Reactions
     if 'reactions' in data:
         for i in range(len(data['reactions'])):
-            order = ['equation','type','rate-constant','low-P-rate-constant','high-P-rate-constant','duplicate','efficiencies','JPL','citation']
+            order = ['equation','type','rate-constant','rate-constants','low-P-rate-constant',
+                     'high-P-rate-constant','duplicate','efficiencies','JPL','citation']
             copy = data['reactions'][i].copy()
             data['reactions'][i].clear()
             for key in order:
@@ -107,6 +108,10 @@ def FormatReactions_main(data):
             for key in flowstyle:
                 if key in data['reactions'][i].keys():
                     data['reactions'][i][key] = flowmap(data['reactions'][i][key])
+
+            if 'rate-constants' in data['reactions'][i]:
+                for j in range(len(data['reactions'][i]['rate-constants'])):
+                    data['reactions'][i]['rate-constants'][j] = flowmap(data['reactions'][i]['rate-constants'][j])
 
     return data
 

src/atmosphere/photochem_atmosphere_rhs.f90

@@ -540,7 +540,7 @@ module subroutine prep_all_background_gas(self, usol_in, err)
     enddo
 
     !!! reaction rates
-    call reaction_rates(self%dat, self%var, wrk%density, wrk%densities, wrk%rx_rates)
+    call reaction_rates(self%dat, self%var, wrk%pressure, wrk%density, wrk%densities, wrk%rx_rates)
 
     ! Update the photon_flux if the function is associated.
     ! we use time wrk%tn, which MUST be updated.

src/evoatmosphere/photochem_evoatmosphere_rhs.f90

@@ -612,7 +612,7 @@ module subroutine prep_all_evo_gas(self, usol_in, err)
     enddo
 
     !!! reaction rates
-    call reaction_rates(self%dat, self%var, wrk%density, wrk%densities, wrk%rx_rates)
+    call reaction_rates(self%dat, self%var, wrk%pressure, wrk%density, wrk%densities, wrk%rx_rates)
 
     ! Update the photon_flux if the function is associated.
     ! we use time wrk%tn, which MUST be updated.

src/input/photochem_input_read.f90

@@ -1208,7 +1208,7 @@ subroutine compare_rxtype_string(tmp, eqr, eqp, reverse, rxtype_int, err)
     logical, intent(in) :: reverse
     integer, intent(in) :: rxtype_int
     character(:), allocatable, intent(out) :: err
-    character(len=s_str_len) :: rxtype
+    character(:), allocatable :: rxtype
     integer i
     logical k, j, m, l, kk, jj
     l = .false.
@@ -1225,6 +1225,11 @@ subroutine compare_rxtype_string(tmp, eqr, eqp, reverse, rxtype_int, err)
       rxtype = 'three-body'
     elseif (rxtype_int == 3) then
       rxtype = 'falloff'
+    elseif (rxtype_int == 4) then
+      rxtype = 'pressure-dependent-Arrhenius'
+    else
+      err = 'Internal error in Photochem involving "compare_rxtype_string". Report this bug!'
+      return
     endif
 
     if ((trim(rxtype) == 'three-body') .or. (trim(rxtype) == 'falloff')) then
@@ -1265,7 +1270,7 @@ subroutine compare_rxtype_string(tmp, eqr, eqp, reverse, rxtype_int, err)
                 ' can not contain "hv". Only photolysis reactions can.'
         return
       endif
-    elseif (trim(rxtype) == 'elementary') then
+    elseif (trim(rxtype) == 'elementary' .or. trim(rxtype) == 'pressure-dependent-Arrhenius') then
       do i = 1,size(eqr)
         if (trim(eqr(i)) == 'M') j = .true.
         if (trim(eqr(i)) == 'hv') m = .true.
@@ -1530,7 +1535,7 @@ function is_rx_reverse(rx_string, err) result(reverse)
 
   subroutine get_rateparams(dat, reaction_d, infile, rx, err)
     use photochem_enum, only: NoFalloff, TroeWithoutT2Falloff, TroeWithT2Falloff, JPLFalloff
-    use photochem_types, only: Reaction, BaseRate, ElementaryRate, ThreeBodyRate, FalloffRate, PhotolysisRate
+    use photochem_types, only: Reaction, BaseRate, ElementaryRate, ThreeBodyRate, FalloffRate, PhotolysisRate, PressDependentRate
     type(PhotochemData), intent(in) :: dat
     class(type_dictionary), intent(in) :: reaction_d
     character(len=*), intent(in) :: infile
@@ -1541,10 +1546,11 @@ subroutine get_rateparams(dat, reaction_d, infile, rx, err)
     type (type_error), allocatable :: io_err
     character(len=str_len) :: rxtype_str
     type(type_dictionary), pointer :: dict
+    type(type_list), pointer :: list
+    class(type_list_item), pointer :: item
     logical :: use_jpl
     real(dp) :: T2
-
-
+
     rxtype_str = reaction_d%get_string("type", default="elementary", error = io_err) 
 
     if (rxtype_str == 'photolysis') then
@@ -1559,6 +1565,9 @@ subroutine get_rateparams(dat, reaction_d, infile, rx, err)
     elseif (rxtype_str == 'falloff') then
       allocate(FalloffRate::rx%rp)
       rx%rp%rxtype = 3
+    elseif (rxtype_str == 'pressure-dependent-Arrhenius') then
+      allocate(PressDependentRate::rx%rp)
+      rx%rp%rxtype = 4
     else
       err = 'IOError: reaction type '//trim(rxtype_str)//' is not a valid reaction type.'
       return
@@ -1642,9 +1651,88 @@ subroutine get_rateparams(dat, reaction_d, infile, rx, err)
         rp%T3 = dict%get_real('T3',error = io_err)
         if (allocated(io_err)) then; err = trim(infile)//trim(io_err%message); return; endif
       endif
+
+    class is (PressDependentRate); block
+      use futils, only: argsort
+      use photochem_types, only: MultiArrheniusRate
+      type(MultiArrheniusRate) :: rate
+      real(dp), allocatable :: P(:), A(:), b(:), Ea(:)
+      integer, allocatable :: inds(:)
+      integer :: i, j
+
+      list => reaction_d%get_list('rate-constants',.true.,error=io_err)
+      if (allocated(io_err)) then; err = trim(infile)//trim(io_err%message); return; endif
+
+      i = list%size()
+      allocate(P(i), A(i), b(i), Ea(i), inds(i))
+
+      i = 1
+      item => list%first
+      do while (associated(item))
+        select type (listitem => item%node)
+        class is (type_dictionary)
+
+          P(i) = listitem%get_real('P',error = io_err)
+          if (allocated(io_err)) then; err = trim(io_err%message); return; endif
+
+          if (P(i) <= 0.0_dp) then
+            err = 'Pressures must be positive in pressure-dependent-Arrhenius reaction '// &
+            trim(reaction_d%get_string("equation",error=io_err))
+            return
+          endif
+
+          call get_arrhenius(listitem, A(i), b(i), Ea(i), err)
+          if (allocated(err)) return
+
+        class default
+          err = '"rate-constants" must be a list of dictionaries for reaction '// &
+          trim(reaction_d%get_string("equation",error=io_err))
+          return
+        end select
+        i = i + 1
+        item => item%next
+      enddo
+
+      ! sort
+      inds = argsort(P)
+      P = P(inds)
+      A = A(inds)
+      b = b(inds)
+      Ea = Ea(inds)
+
+      allocate(rp%logP(0))
+      allocate(rp%rate(0))
+      i = 1
+      do while (i <= size(P))
+        rp%logP = [rp%logP, log(P(i))]
+        rate = MultiArrheniusRate([A(i)], [b(i)], [Ea(i)])  
+        j = i + 1
+        do
+          if (j > size(P)) then
+            i = j
+            exit
+          endif
+          if (P(j) == P(i)) then
+            rate%A = [rate%A, A(j)]
+            rate%b = [rate%b, b(j)]
+            rate%Ea = [rate%Ea, Ea(j)]
+          else
+            i = j
+            exit
+          endif
+          j = j + 1
+        enddo
+        rp%rate = [rp%rate, rate]
+      enddo
+
+      if (size(rp%logP) < 2) then
+        err = 'pressure-dependent-Arrhenius reaction '// &
+        trim(reaction_d%get_string("equation",error=io_err))// &
+        ' must have > 1 unique pressures.'
+        return
+      endif
 
-
-    end select
+    endblock; end select
 
   end subroutine
 

src/photochem_common.f90

@@ -17,25 +17,26 @@ module photochem_common
 
 contains
 
-  pure subroutine reaction_rates(dat, var, density, densities, rx_rates)
+  subroutine reaction_rates(dat, var, pressure, density, densities, rx_rates)
     use photochem_enum, only: NoFalloff, TroeWithoutT2Falloff, TroeWithT2Falloff, JPLFalloff
-    use photochem_types, only: ElementaryRate, ThreeBodyRate, FalloffRate
-    use photochem_eqns, only: arrhenius_rate, Troe_noT2, Troe_withT2, falloff_rate
+    use photochem_types, only: ElementaryRate, ThreeBodyRate, FalloffRate, PressDependentRate
+    use photochem_eqns, only: arrhenius_rate, Troe_noT2, Troe_withT2, falloff_rate, searchsorted
     use photochem_const, only: Rgas, k_boltz, smallest_real ! constants
 
     type(PhotochemData), intent(in) :: dat
     type(PhotochemVars), intent(in) :: var
+    real(dp), intent(in) :: pressure(:) ! (nz)
     real(dp), intent(in) :: density(:) ! (nz)
     real(dp), intent(in) :: densities(:,:) ! (nsp+1,nz)
     real(dp), intent(out) :: rx_rates(:,:) ! (nz,nrT) 
 
     integer :: i, j, k, n, l, m
     real(dp) :: eff_den(var%nz), F(var%nz)
     real(dp) :: k0, kinf(var%nz), Pr(var%nz)
+    real(dp) :: logP, k1, k2, logk
     real(dp) :: gibbR_forward, gibbP_forward
     real(dp) :: Dg_forward
 
-
     do i = 1,dat%nrF
       select type(rp => dat%rx(i)%rp)
       class is (ElementaryRate)
@@ -101,6 +102,57 @@ pure subroutine reaction_rates(dat, var, density, densities, rx_rates)
           rx_rates(j,i) = falloff_rate(kinf(j), Pr(j), F(j))
         enddo
 
+      class is (PressDependentRate)
+        n = size(rp%logP) ! number of pressures in grid
+
+        do j = 1,var%nz
+          logP = log(pressure(j))
+
+          if (logP <= rp%logP(1)) then
+            ! If P below grid, then constant extrapolation below.
+            k1 = 0.0_dp
+            do k = 1,size(rp%rate(1)%A)
+              k1 = k1 + &
+                   arrhenius_rate(rp%rate(1)%A(k), rp%rate(1)%b(k), &
+                                  rp%rate(1)%Ea(k), var%temperature(j))
+            enddo
+            rx_rates(j,i) = k1
+          elseif (logP >= rp%logP(n)) then
+            ! If P above grid, then constant extrapolation above.
+            k1 = 0.0_dp
+            do k = 1,size(rp%rate(n)%A)
+              k1 = k1 + &
+                   arrhenius_rate(rp%rate(n)%A(k), rp%rate(n)%b(k), &
+                                  rp%rate(n)%Ea(k), var%temperature(j))
+            enddo
+            rx_rates(j,i) = k1
+          else
+            ! P lies within grid
+            l = searchsorted(rp%logP, logP) - 1
+
+            ! Compute rate at lower pressure
+            k1 = 0.0_dp
+            do k = 1,size(rp%rate(l)%A)
+              k1 = k1 + &
+                   arrhenius_rate(rp%rate(l)%A(k), rp%rate(l)%b(k), &
+                                  rp%rate(l)%Ea(k), var%temperature(j))
+            enddo
+            k1 = log(k1)
+
+            ! Compute rate at higher pressure
+            k2 = 0.0_dp
+            do k = 1,size(rp%rate(l+1)%A)
+              k2 = k2 + &
+                   arrhenius_rate(rp%rate(l+1)%A(k), rp%rate(l+1)%b(k), &
+                                  rp%rate(l+1)%Ea(k), var%temperature(j))
+            enddo
+            k2 = log(k2)
+
+            ! log-interpolate to get rate at intermediate pressure.
+            logk = k1 + (k2 - k1)*((logP - rp%logP(l))/(rp%logP(l+1) - rp%logP(l)))
+            rx_rates(j,i) = exp(logk)
+          endif
+        enddo
       end select
 
     enddo ! end loop over forward reactions

src/photochem_eqns.f90

@@ -364,4 +364,38 @@ pure function zahnle_Hescape_coeff(S1) result(coeff)
 
   end function
 
+  !> Mimics numpy.searchsorted
+  pure function searchsorted(arr, val) result(ind)
+    real(dp), intent(in) :: arr(:) !! Input sorted array
+    real(dp), intent(in) :: val !! Value to compare to arr
+    integer :: ind !! Index that satisfies arr(i-1) < val <= arr(i)
+
+    integer :: low, high, mid
+
+    if (val <= arr(1)) then
+      ind = 1
+      return
+    endif
+
+    if (val > arr(size(arr))) then
+      ind = size(arr) + 1
+      return
+    endif
+
+    low = 1
+    high = size(arr)
+    do
+      mid = (low + high)/2
+      if (val > arr(mid)) then
+        low = mid
+      else
+        high = mid
+      endif
+      if (high-1 == low) exit
+    enddo
+
+    ind = high
+
+  end function
+
 end module