CAMP 1.0.0
Chemistry Across Multiple Phases
rxn_troe.F90
Go to the documentation of this file.
1! Copyright (C) 2021 Barcelona Supercomputing Center and University of
2! Illinois at Urbana-Champaign
3! SPDX-License-Identifier: MIT
4
5!> \file
6!> The camp_rxn_troe module.
7
8!> \page camp_rxn_troe CAMP: Troe Reaction
9!!
10!! Troe (fall-off) reaction rate constant equations take the form:
11!!
12!! \f[
13!! \frac{k_0[\mbox{M}]}{1+k_0[\mbox{M}]/k_{\inf}}F_C^{(1+1/N[log_{10}(k_0[\mbox{M}]/k_{\inf})]^2)^{-1}}
14!! \f]
15!!
16!! where \f$k_0\f$ is the low-pressure limiting rate constant, \f$k_{\inf}\f$
17!! is the high-pressure limiting rate constant, \f$[\mbox{M}]\f$ is the
18!! density of air, and \f$F_C\f$ and \f$N\f$ are parameters
19!! that determine the shape of the fall-off curve, and are typically 0.6 and
20!! 1.0, respectively \cite Finlayson-Pitts2000 \cite Gipson. \f$k_0\f$ and
21!! \f$k_{\inf}\f$ are calculated as \ref camp_rxn_arrhenius "Arrhenius" rate
22!! constants with \f$D=300\f$ and \f$E=0\f$.
23!!
24!! Input data for Troe reactions have the following format :
25!! \code{.json}
26!! {
27!! "type" : "TROE",
28!! "k0_A" : 5.6E-12,
29!! "k0_B" : -1.8,
30!! "k0_C" : 180.0,
31!! "kinf_A" : 3.4E-12,
32!! "kinf_B" : -1.6,
33!! "kinf_C" : 104.1,
34!! "Fc" : 0.7,
35!! "N" : 0.9,
36!! "time unit" : "MIN",
37!! "reactants" : {
38!! "spec1" : {},
39!! "spec2" : { "qty" : 2 },
40!! ...
41!! },
42!! "products" : {
43!! "spec3" : {},
44!! "spec4" : { "yield" : 0.65 },
45!! ...
46!! }
47!! }
48!! \endcode
49!! The key-value pairs \b reactants, and \b products are required. Reactants
50!! without a \b qty value are assumed to appear once in the reaction equation.
51!! Products without a specified \b yield are assumed to have a \b yield of
52!! 1.0.
53!!
54!! The two sets of parameters beginning with \b k0_ and \b kinf_ are the
55!! \ref camp_rxn_arrhenius "Arrhenius" parameters for the \f$k_0\f$ and
56!! \f$k_{\inf}\f$ rate constants, respectively. When not present, \b _A
57!! parameters are assumed to be 1.0, \b _B to be 0.0, \b _C to be 0.0, \b Fc
58!! to be 0.6 and \b N to be 1.0.
59!!
60!! The unit for time is assumed to be s, but inclusion of the optional
61!! key-value pair \b time \b unit = \b MIN can be used to indicate a rate
62!! with min as the time unit.
63
64!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
65
66!> The rxn_troe_t type and associated functions.
68
72 use camp_constants, only: const
76 use camp_util, only: i_kind, dp, to_string, &
78
79 implicit none
80 private
81
82#define NUM_REACT_ this%condensed_data_int(1)
83#define NUM_PROD_ this%condensed_data_int(2)
84#define K0_A_ this%condensed_data_real(1)
85#define K0_B_ this%condensed_data_real(2)
86#define K0_C_ this%condensed_data_real(3)
87#define KINF_A_ this%condensed_data_real(4)
88#define KINF_B_ this%condensed_data_real(5)
89#define KINF_C_ this%condensed_data_real(6)
90#define FC_ this%condensed_data_real(7)
91#define N_ this%condensed_data_real(8)
92#define SCALING_ this%condensed_data_real(9)
93#define CONV_ this%condensed_data_real(10)
94#define NUM_INT_PROP_ 2
95#define NUM_REAL_PROP_ 10
96#define NUM_ENV_PARAM_ 1
97#define REACT_(x) this%condensed_data_int(NUM_INT_PROP_ + x)
98#define PROD_(x) this%condensed_data_int(NUM_INT_PROP_ + NUM_REACT_ + x)
99#define DERIV_ID_(x) this%condensed_data_int(NUM_INT_PROP_ + NUM_REACT_ + NUM_PROD_ + x)
100#define JAC_ID_(x) this%condensed_data_int(NUM_INT_PROP_ + 2*(NUM_REACT_+NUM_PROD_) + x)
101#define YIELD_(x) this%condensed_data_real(NUM_REAL_PROP_ + x)
102
103public :: rxn_troe_t
104
105 !> Generic test reaction data type
106 type, extends(rxn_data_t) :: rxn_troe_t
107 contains
108 !> Reaction initialization
109 procedure :: initialize
110 !> Finalize the reaction
112 end type rxn_troe_t
113
114 !> Constructor for rxn_troe_t
115 interface rxn_troe_t
116 procedure :: constructor
117 end interface rxn_troe_t
118
119contains
120
121!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
122
123 !> Constructor for Troe reaction
124 function constructor() result(new_obj)
125
126 !> A new reaction instance
127 type(rxn_troe_t), pointer :: new_obj
128
129 allocate(new_obj)
130 new_obj%rxn_phase = gas_rxn
131
132 end function constructor
133
134!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
135
136 !> Initialize the reaction data, validating component data and loading
137 !! any required information into the condensed data arrays for use during
138 !! solving
139 subroutine initialize(this, chem_spec_data, aero_phase, aero_rep, n_cells)
140
141 !> Reaction data
142 class(rxn_troe_t), intent(inout) :: this
143 !> Chemical species data
144 type(chem_spec_data_t), intent(in) :: chem_spec_data
145 !> Aerosol phase data
146 type(aero_phase_data_ptr), intent(in) :: aero_phase(:)
147 !> Aerosol representations
148 type(aero_rep_data_ptr), pointer, intent(in) :: aero_rep(:)
149 !> Number of grid cells to solve simultaneously
150 integer(kind=i_kind), intent(in) :: n_cells
151
152 type(property_t), pointer :: spec_props, reactants, products
153 character(len=:), allocatable :: key_name, spec_name, string_val
154 integer(kind=i_kind) :: i_spec, i_qty
155
156 integer(kind=i_kind) :: temp_int
157 real(kind=dp) :: temp_real
158
159 ! Get the species involved
160 if (.not. associated(this%property_set)) call die_msg(510658779, &
161 "Missing property set needed to initialize reaction")
162 key_name = "reactants"
163 call assert_msg(852878121, &
164 this%property_set%get_property_t(key_name, reactants), &
165 "Troe reaction is missing reactants")
166 key_name = "products"
167 call assert_msg(965196466, &
168 this%property_set%get_property_t(key_name, products), &
169 "Troe reaction is missing products")
170
171 ! Count the number of reactants (including those with a qty specified)
172 call reactants%iter_reset()
173 i_spec = 0
174 do while (reactants%get_key(spec_name))
175 ! Get properties included with this reactant in the reaction data
176 call assert(844081716, reactants%get_property_t(val=spec_props))
177 key_name = "qty"
178 if (spec_props%get_int(key_name, temp_int)) i_spec = i_spec+temp_int-1
179 call reactants%iter_next()
180 i_spec = i_spec + 1
181 end do
182
183 ! Allocate space in the condensed data arrays
184 ! Space in this example is allocated for two sets of inidices for the
185 ! reactants and products, one molecular property for each reactant,
186 ! yields for the products and three reaction parameters.
187 allocate(this%condensed_data_int(num_int_prop_ + &
188 (i_spec + 2) * (i_spec + products%size())))
189 allocate(this%condensed_data_real(num_real_prop_ + products%size()))
190 this%condensed_data_int(:) = int(0, kind=i_kind)
191 this%condensed_data_real(:) = real(0.0, kind=dp)
192
193 ! Save space for the environment-dependent parameters
194 this%num_env_params = num_env_param_
195
196 ! Save the size of the reactant and product arrays (for reactions where
197 ! these can vary)
198 num_react_ = i_spec
199 num_prod_ = products%size()
200
201 ! Set the #/cc -> ppm conversion prefactor
202 conv_ = const%avagadro / const%univ_gas_const * 10.0d0**(-12.0d0)
203
204 ! Get reaction parameters (it might be easiest to keep these at the
205 ! beginning of the condensed data array, so they can be accessed using
206 ! compliler flags)
207 key_name = "k0_A"
208 if (.not. this%property_set%get_real(key_name, k0_a_)) then
209 k0_a_ = 1.0
210 end if
211 key_name = "k0_B"
212 if (.not. this%property_set%get_real(key_name, k0_b_)) then
213 k0_b_ = 0.0
214 end if
215 key_name = "k0_C"
216 if (.not. this%property_set%get_real(key_name, k0_c_)) then
217 k0_c_ = 0.0
218 end if
219 key_name = "kinf_A"
220 if (.not. this%property_set%get_real(key_name, kinf_a_)) then
221 kinf_a_ = 1.0
222 end if
223 key_name = "kinf_B"
224 if (.not. this%property_set%get_real(key_name, kinf_b_)) then
225 kinf_b_ = 0.0
226 end if
227 key_name = "kinf_C"
228 if (.not. this%property_set%get_real(key_name, kinf_c_)) then
229 kinf_c_ = 0.0
230 end if
231 key_name = "Fc"
232 if (.not. this%property_set%get_real(key_name, fc_)) then
233 fc_ = 0.6
234 end if
235 key_name = "N"
236 if (.not. this%property_set%get_real(key_name, n_)) then
237 n_ = 1.0
238 end if
239 key_name = "time unit"
240 scaling_ = real(1.0, kind=dp)
241 if (this%property_set%get_string(key_name, string_val)) then
242 if (trim(string_val).eq."MIN") then
243 scaling_ = real(1.0d0/60.0d0, kind=dp)
244 end if
245 endif
246
247 ! Include [M] in K0_A_
248 k0_a_ = k0_a_ * real(1.0d6, kind=dp)
249
250 ! Get the indices and chemical properties for the reactants
251 call reactants%iter_reset()
252 i_spec = 1
253 do while (reactants%get_key(spec_name))
254
255 ! Save the index of this species in the state variable array
256 react_(i_spec) = chem_spec_data%gas_state_id(spec_name)
257
258 ! Make sure the species exists
259 call assert_msg(595701751, react_(i_spec).gt.0, &
260 "Missing Troe reactant: "//spec_name)
261
262 ! Get properties included with this reactant in the reaction data
263 call assert(221292659, reactants%get_property_t(val=spec_props))
264 key_name = "qty"
265 if (spec_props%get_int(key_name, temp_int)) then
266 do i_qty = 1, temp_int - 1
267 react_(i_spec + i_qty) = react_(i_spec)
268 end do
269 i_spec = i_spec + temp_int - 1
270 end if
271
272 call reactants%iter_next()
273 i_spec = i_spec + 1
274 end do
275
276 ! Get the indices and chemical properties for the products
277 call products%iter_reset()
278 i_spec = 1
279 do while (products%get_key(spec_name))
280
281 ! Save the index of this species in the state variable array
282 prod_(i_spec) = chem_spec_data%gas_state_id(spec_name)
283
284 ! Make sure the species exists
285 call assert_msg(480024633, prod_(i_spec).gt.0, &
286 "Missing Troe product: "//spec_name)
287
288 ! Get properties included with this product in the reaction data
289 call assert(393355097, products%get_property_t(val=spec_props))
290 key_name = "yield"
291 if (spec_props%get_real(key_name, temp_real)) then
292 yield_(i_spec) = temp_real
293 else
294 yield_(i_spec) = 1.0
295 end if
296
297 call products%iter_next()
298 i_spec = i_spec + 1
299 end do
300
301 end subroutine initialize
302
303!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
304
305 !> Finalize the reaction
306 subroutine finalize(this)
307
308 !> Reaction data
309 type(rxn_troe_t), intent(inout) :: this
310
311 if (associated(this%property_set)) &
312 deallocate(this%property_set)
313 if (allocated(this%condensed_data_real)) &
314 deallocate(this%condensed_data_real)
315 if (allocated(this%condensed_data_int)) &
316 deallocate(this%condensed_data_int)
317
318 end subroutine finalize
319
320!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
321
322 !> Finalize an array of reactions
323 subroutine finalize_array(this)
324
325 !> Array of reaction data
326 type(rxn_troe_t), intent(inout) :: this(:)
327
328 integer(kind=i_kind) :: i_rxn
329
330 do i_rxn = 1, size(this)
331 call finalize(this(i_rxn))
332 end do
333
334 end subroutine finalize_array
335
336!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
337
338end module camp_rxn_troe
Initialize the aerosol representation data, validating component data and loading any required inform...
Get the non-unique name of a chemical species by its unique name.
Interface for to_string functions.
Definition util.F90:32
The abstract aero_phase_data_t structure and associated subroutines.
subroutine finalize_array(this)
Finalize the aerosol phase data.
type(aero_phase_data_t) function, pointer constructor(phase_name, init_size)
Constructor for aero_phase_data_t.
subroutine finalize(this)
Finalize the aerosol phase data.
The abstract aero_rep_data_t structure and associated subroutines.
The camp_state_t structure and associated subroutines.
Definition camp_state.F90:9
The chem_spec_data_t structure and associated subroutines.
Physical constants.
Definition constants.F90:9
integer, parameter dp
Kind of a double precision real number.
Definition constants.F90:16
type(const_t), save const
Fixed variable for accessing the constant's values.
Definition constants.F90:77
integer, parameter i_kind
Kind of an integer.
Definition constants.F90:21
The property_t structure and associated subroutines.
Definition property.F90:9
The rxn_data_t structure and associated subroutines.
Definition rxn_data.F90:60
integer(kind=i_kind), parameter, public gas_rxn
Gas-phase reaction.
Definition rxn_data.F90:85
The rxn_troe_t type and associated functions.
Definition rxn_troe.F90:67
Common utility subroutines.
Definition util.F90:9
subroutine assert(code, condition_ok)
Errors unless condition_ok is true.
Definition util.F90:165
subroutine die_msg(code, error_msg)
Error immediately.
Definition util.F90:196
subroutine assert_msg(code, condition_ok, error_msg)
Errors unless condition_ok is true.
Definition util.F90:130
#define n_
Pointer type for building arrays.
Pointer to aero_rep_data_t extending types.
Abstract reaction data type.
Definition rxn_data.F90:99
Generic test reaction data type.
Definition rxn_troe.F90:115