CAMP 1.0.0
Chemistry Across Multiple Phases
rxn_arrhenius.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_arrhenius module.
7
8!> \page camp_rxn_arrhenius CAMP: Arrhenius Reaction
9!!
10!! Arrhenius-like reaction rate constant equations are calculated as follows:
11!!
12!! \f[
13!! Ae^{(\frac{-E_a}{k_bT})}(\frac{T}{D})^B(1.0+E*P)
14!! \f]
15!!
16!! where \f$A\f$ is the pre-exponential factor
17!! (# \f$(\mbox{cm}^{-3})^{-(n-1)}\mbox{s}^{-1}\f$),
18!! \f$n\f$ is the number of reactants, \f$E_a\f$ is the activation energy (J),
19!! \f$k_b\f$ is the Boltzmann constant (J/K), \f$D\f$ (K), \f$B\f$ (unitless)
20!! and \f$E\f$ (\f$Pa^{-1}\f$) are reaction parameters, \f$T\f$ is the
21!! temperature (K), and \f$P\f$ is the pressure (Pa). The first two terms are
22!! described in Finlayson-Pitts and Pitts (2000) \cite Finlayson-Pitts2000 .
23!! The final term is included to accomodate CMAQ EBI solver type 7 rate
24!! constants.
25!!
26!! Input data for Arrhenius equations has the following format:
27!! \code{.json}
28!! {
29!! "type" : "ARRHENIUS",
30!! "A" : 123.45,
31!! "Ea" : 123.45,
32!! "B" : 1.3,
33!! "D" : 300.0,
34!! "E" : 0.6E-5,
35!! "time unit" : "MIN",
36!! "reactants" : {
37!! "spec1" : {},
38!! "spec2" : { "qty" : 2 },
39!! ...
40!! },
41!! "products" : {
42!! "spec3" : {},
43!! "spec4" : { "yield" : 0.65 },
44!! ...
45!! }
46!! }
47!! \endcode
48!! The key-value pairs \b reactants, and \b products are required. Reactants
49!! without a \b qty value are assumed to appear once in the reaction equation.
50!! Products without a specified \b yield are assumed to have a \b yield of
51!! 1.0.
52!!
53!! Optionally, a parameter \b C may be included, and is taken to equal
54!! \f$\frac{-E_a}{k_b}\f$. Note that either \b Ea or \b C may be included, but
55!! not both. When neither \b Ea or \b C are included, they are assumed to be
56!! 0.0. When \b A is not included, it is assumed to be 1.0, when \b D is not
57!! included, it is assumed to be 300.0 K, when \b B is not included, it is
58!! assumed to be 0.0, and when \b E is not included, it is assumed to be 0.0.
59!! The unit for time is assumed to be s, but inclusion of the optional
60!! key-value pair \b time \b unit = \b MIN can be used to indicate a rate
61!! with min as the time unit.
62
63!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
64
65!> The rxn_arrhenius_t type and associated functions.
67
71 use camp_constants, only: const
75 use camp_util, only: i_kind, dp, to_string, &
77
78 implicit none
79 private
80
81#define NUM_REACT_ this%condensed_data_int(1)
82#define NUM_PROD_ this%condensed_data_int(2)
83#define A_ this%condensed_data_real(1)
84#define B_ this%condensed_data_real(2)
85#define C_ this%condensed_data_real(3)
86#define D_ this%condensed_data_real(4)
87#define E_ this%condensed_data_real(5)
88#define CONV_ this%condensed_data_real(6)
89#define NUM_INT_PROP_ 2
90#define NUM_REAL_PROP_ 6
91#define NUM_ENV_PARAM_ 1
92#define REACT_(x) this%condensed_data_int(NUM_INT_PROP_ + x)
93#define PROD_(x) this%condensed_data_int(NUM_INT_PROP_ + NUM_REACT_ + x)
94#define DERIV_ID_(x) this%condensed_data_int(NUM_INT_PROP_ + NUM_REACT_ + NUM_PROD_ + x)
95#define JAC_ID_(x) this%condensed_data_int(NUM_INT_PROP_ + 2*(NUM_REACT_+NUM_PROD_) + x)
96#define YIELD_(x) this%condensed_data_real(NUM_REAL_PROP_ + x)
97
98 public :: rxn_arrhenius_t
99
100 !> Generic test reaction data type
101 type, extends(rxn_data_t) :: rxn_arrhenius_t
102 contains
103 !> Reaction initialization
104 procedure :: initialize
105 !> Finalize the reaction
107 end type rxn_arrhenius_t
108
109 !> Constructor for rxn_arrhenius_t
111 procedure :: constructor
112 end interface rxn_arrhenius_t
113
114contains
115
116!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
117
118 !> Constructor for Arrhenius reaction
119 function constructor() result(new_obj)
120
121 !> A new reaction instance
122 type(rxn_arrhenius_t), pointer :: new_obj
123
124 allocate(new_obj)
125 new_obj%rxn_phase = gas_rxn
126
127 end function constructor
128
129!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
130
131 !> Initialize the reaction data, validating component data and loading
132 !! any required information into the condensed data arrays for use during
133 !! solving
134 subroutine initialize(this, chem_spec_data, aero_phase, aero_rep, n_cells)
135
136 !> Reaction data
137 class(rxn_arrhenius_t), intent(inout) :: this
138 !> Chemical species data
139 type(chem_spec_data_t), intent(in) :: chem_spec_data
140 !> Aerosol phase data
141 type(aero_phase_data_ptr), intent(in) :: aero_phase(:)
142 !> Aerosol representations
143 type(aero_rep_data_ptr), pointer, intent(in) :: aero_rep(:)
144 !> Number of grid cells being solved simultaneously
145 integer(kind=i_kind), intent(in) :: n_cells
146
147 type(property_t), pointer :: spec_props, reactants, products
148 character(len=:), allocatable :: key_name, spec_name, string_val
149 integer(kind=i_kind) :: i_spec, i_qty
150
151 integer(kind=i_kind) :: temp_int
152 real(kind=dp) :: temp_real
153
154 ! Get the species involved
155 if (.not. associated(this%property_set)) call die_msg(255324828, &
156 "Missing property set needed to initialize reaction")
157 key_name = "reactants"
158 call assert_msg(250060521, &
159 this%property_set%get_property_t(key_name, reactants), &
160 "Arrhenius reaction is missing reactants")
161 key_name = "products"
162 call assert_msg(304540307, &
163 this%property_set%get_property_t(key_name, products), &
164 "Arrhenius reaction is missing products")
165
166 ! Count the number of reactants (including those with a qty specified)
167 call reactants%iter_reset()
168 i_spec = 0
169 do while (reactants%get_key(spec_name))
170 ! Get properties included with this reactant in the reaction data
171 call assert(243342975, reactants%get_property_t(val=spec_props))
172 key_name = "qty"
173 if (spec_props%get_int(key_name, temp_int)) i_spec = i_spec+temp_int-1
174 call reactants%iter_next()
175 i_spec = i_spec + 1
176 end do
177
178 ! Allocate space in the condensed data arrays
179 allocate(this%condensed_data_int(num_int_prop_ + &
180 (i_spec + 2) * (i_spec + products%size())))
181 allocate(this%condensed_data_real(num_real_prop_ + products%size()))
182 this%condensed_data_int(:) = int(0, kind=i_kind)
183 this%condensed_data_real(:) = real(0.0, kind=dp)
184
185 ! Save space for the environment dependent parameters
186 this%num_env_params = num_env_param_
187
188 ! Save the size of the reactant and product arrays (for reactions where
189 ! these can vary)
190 num_react_ = i_spec
191 num_prod_ = products%size()
192
193 ! Set the #/cc -> ppm conversion prefactor
194 conv_ = const%avagadro / const%univ_gas_const * 10.0d0**(-12.0d0)
195
196 ! Get reaction parameters (it might be easiest to keep these at the
197 ! beginning of the condensed data array, so they can be accessed using
198 ! compliler flags)
199 key_name = "A"
200 if (.not. this%property_set%get_real(key_name, a_)) then
201 a_ = 1.0
202 end if
203 key_name = "time unit"
204 if (this%property_set%get_string(key_name, string_val)) then
205 if (trim(string_val).eq."MIN") then
206 a_ = a_ / 60.0
207 end if
208 endif
209 key_name = "Ea"
210 if (this%property_set%get_real(key_name, temp_real)) then
211 c_ = -temp_real/const%boltzmann
212 key_name = "C"
213 call assert_msg(297370315, &
214 .not.this%property_set%get_real(key_name, temp_real), &
215 "Received both Ea and C parameter for Arrhenius equation")
216 else
217 key_name = "C"
218 if (.not. this%property_set%get_real(key_name, c_)) then
219 c_ = 0.0
220 end if
221 end if
222 key_name = "D"
223 if (.not. this%property_set%get_real(key_name, d_)) then
224 d_ = 300.0
225 end if
226 key_name = "B"
227 if (.not. this%property_set%get_real(key_name, b_)) then
228 b_ = 0.0
229 end if
230 key_name = "E"
231 if (.not. this%property_set%get_real(key_name, e_)) then
232 e_ = 0.0
233 end if
234
235 call assert_msg(344705857, .not. ((b_.ne.real(0.0, kind=dp)) &
236 .and.(d_.eq.real(0.0, kind=dp))), &
237 "D cannot be zero if B is non-zero in Arrhenius reaction.")
238
239 ! Get the indices and chemical properties for the reactants
240 call reactants%iter_reset()
241 i_spec = 1
242 do while (reactants%get_key(spec_name))
243
244 ! Save the index of this species in the state variable array
245 react_(i_spec) = chem_spec_data%gas_state_id(spec_name)
246
247 ! Make sure the species exists
248 call assert_msg(751684145, react_(i_spec).gt.0, &
249 "Missing Arrhenius reactant: "//spec_name)
250
251 ! Get properties included with this reactant in the reaction data
252 call assert(796763915, reactants%get_property_t(val=spec_props))
253 key_name = "qty"
254 if (spec_props%get_int(key_name, temp_int)) then
255 do i_qty = 1, temp_int - 1
256 react_(i_spec + i_qty) = react_(i_spec)
257 end do
258 i_spec = i_spec + temp_int - 1
259 end if
260
261 call reactants%iter_next()
262 i_spec = i_spec + 1
263 end do
264
265 ! Get the indices and chemical properties for the products
266 call products%iter_reset()
267 i_spec = 1
268 do while (products%get_key(spec_name))
269
270 ! Save the index of this species in the state variable array
271 prod_(i_spec) = chem_spec_data%gas_state_id(spec_name)
272
273 ! Make sure the species exists
274 call assert_msg(234495887, prod_(i_spec).gt.0, &
275 "Missing Arrhenius product: "//spec_name)
276
277 ! Get properties included with this product in the reaction data
278 call assert(451185800, products%get_property_t(val=spec_props))
279 key_name = "yield"
280 if (spec_props%get_real(key_name, temp_real)) then
281 yield_(i_spec) = temp_real
282 else
283 yield_(i_spec) = 1.0
284 end if
285
286 call products%iter_next()
287 i_spec = i_spec + 1
288 end do
289
290 end subroutine initialize
291
292!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
293
294 !> Finalize the reaction
295 subroutine finalize(this)
296
297 !> Reaction data
298 type(rxn_arrhenius_t), intent(inout) :: this
299
300 if (associated(this%property_set)) &
301 deallocate(this%property_set)
302 if (allocated(this%condensed_data_real)) &
303 deallocate(this%condensed_data_real)
304 if (allocated(this%condensed_data_int)) &
305 deallocate(this%condensed_data_int)
306
307 end subroutine finalize
308
309!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
310
311 !> Finalize an array of reactions
312 subroutine finalize_array(this)
313
314 !> Array of reaction data
315 type(rxn_arrhenius_t), intent(inout) :: this(:)
316
317 integer(kind=i_kind) :: i
318
319 do i = 1, size(this)
320 call finalize(this(i))
321 end do
322
323 end subroutine finalize_array
324
325!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
326
327end module camp_rxn_arrhenius
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_arrhenius_t type and associated functions.
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
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
Pointer type for building arrays.
Pointer to aero_rep_data_t extending types.
Generic test reaction data type.
Abstract reaction data type.
Definition rxn_data.F90:99