CAMP 1.0.0
Chemistry Across Multiple Phases
rxn_CMAQ_H2O2.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_CMAQ_H2O2 module.
7
8!> \page camp_rxn_CMAQ_H2O2 CAMP: Special CMAQ Reaction for H2O2
9!!
10!! For the reactions:
11!!
12!! \f[\ce{
13!! HO2 + HO2 ->[k] H2O2
14!! }\f]
15!! \f[\ce{
16!! HO2 + HO2 + H2O ->[k] H2O2
17!! }\f]
18!!
19!! CMAQ rate constants are calculated as:
20!!
21!! \f[
22!! k=k_1+k_2[\mbox{M}]
23!! \f]
24!!
25!! where \f$k_1\f$ and \f$k_2\f$ are \ref camp_rxn_arrhenius "Arrhenius" rate
26!! constants with \f$D=300\f$ and \f$E=0\f$, and \f$[\mbox{M}]\f$ is the
27!! concentration of air \cite Gipson.
28!!
29!! Input data for CMAQ \f$\ce{H2O2}\f$ reactions have the following
30!! format:
31!! \code{.json}
32!! {
33!! "type" : "CMAQ_H2O2",
34!! "k1_A" : 5.6E-12,
35!! "k1_B" : -1.8,
36!! "k1_C" : 180.0,
37!! "k2_A" : 3.4E-12,
38!! "k2_B" : -1.6,
39!! "k2_C" : 104.1,
40!! "time unit" : "MIN",
41!! "reactants" : {
42!! "spec1" : {},
43!! "spec2" : { "qty" : 2 },
44!! ...
45!! },
46!! "products" : {
47!! "spec3" : {},
48!! "spec4" : { "yield" : 0.65 },
49!! ...
50!! }
51!! }
52!! \endcode
53!! The key-value pairs \b reactants, and \b products are required. Reactants
54!! without a \b qty value are assumed to appear once in the reaction equation.
55!! Products without a specified \b yield are assumed to have a \b yield of
56!! 1.0.
57!!
58!! The two sets of parameters beginning with \b k1_ and \b k2_ are the
59!! \ref camp_rxn_arrhenius "Arrhenius" parameters for the \f$k_1\f$ and
60!! \f$k_2\f$ rate constants, respectively. When not present, \b _A
61!! parameters are assumed to be 1.0, \b _B to be 0.0, and \b _C to be 0.0.
62!!
63!! The unit for time is assumed to be s, but inclusion of the optional
64!! key-value pair \b time \b unit = \b MIN can be used to indicate a rate
65!! with min as the time unit.
66
67!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
68
69!> The rxn_CMAQ_H2O2_t type and associated functions.
71
75 use camp_constants, only: const
79 use camp_util, only: i_kind, dp, to_string, &
81
82 implicit none
83 private
84
85#define NUM_REACT_ this%condensed_data_int(1)
86#define NUM_PROD_ this%condensed_data_int(2)
87#define k1_A_ this%condensed_data_real(1)
88#define k1_B_ this%condensed_data_real(2)
89#define k1_C_ this%condensed_data_real(3)
90#define k2_A_ this%condensed_data_real(4)
91#define k2_B_ this%condensed_data_real(5)
92#define k2_C_ this%condensed_data_real(6)
93#define CONV_ this%condensed_data_real(7)
94#define NUM_INT_PROP_ 2
95#define NUM_REAL_PROP_ 7
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_cmaq_h2o2_t
104
105 !> Generic test reaction data type
106 type, extends(rxn_data_t) :: rxn_cmaq_h2o2_t
107 contains
108 !> Reaction initialization
109 procedure :: initialize
110 !> Finalize the reaction
112 end type rxn_cmaq_h2o2_t
113
114 !> Constructor for rxn_CMAQ_H2O2_t
116 procedure :: constructor
117 end interface rxn_cmaq_h2o2_t
118
119contains
120
121!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
122
123 !> Constructor for CMAQ H2O2 reaction
124 function constructor() result(new_obj)
125
126 !> A new reaction instance
127 type(rxn_cmaq_h2o2_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_cmaq_h2o2_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(255324828, &
161 "Missing property set needed to initialize reaction")
162 key_name = "reactants"
163 call assert_msg(940945637, &
164 this%property_set%get_property_t(key_name, reactants), &
165 "CMAQ H2O2 reaction is missing reactants")
166 key_name = "products"
167 call assert_msg(435739232, &
168 this%property_set%get_property_t(key_name, products), &
169 "CMAQ H2O2 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(952475195, 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 = "k1_A"
208 if (.not. this%property_set%get_real(key_name, k1_a_)) then
209 k1_a_ = 1.0
210 end if
211 key_name = "k1_B"
212 if (.not. this%property_set%get_real(key_name, k1_b_)) then
213 k1_b_ = 0.0
214 end if
215 key_name = "k1_C"
216 if (.not. this%property_set%get_real(key_name, k1_c_)) then
217 k1_c_ = 0.0
218 end if
219 key_name = "k2_A"
220 if (.not. this%property_set%get_real(key_name, k2_a_)) then
221 k2_a_ = 1.0
222 end if
223 key_name = "k2_B"
224 if (.not. this%property_set%get_real(key_name, k2_b_)) then
225 k2_b_ = 0.0
226 end if
227 key_name = "k2_C"
228 if (.not. this%property_set%get_real(key_name, k2_c_)) then
229 k2_c_ = 0.0
230 end if
231 key_name = "time unit"
232 if (this%property_set%get_string(key_name, string_val)) then
233 if (trim(string_val).eq."MIN") then
234 k1_a_ = k1_a_ / 60.0
235 k2_a_ = k2_a_ / 60.0
236 end if
237 endif
238
239 ! Include the multiplication of [M]*k2 into k2_A_
240 k2_a_ = k2_a_ * real(1.0d6, kind=dp)
241
242 ! Get the indices and chemical properties for the reactants
243 call reactants%iter_reset()
244 i_spec = 1
245 do while (reactants%get_key(spec_name))
246
247 ! Save the index of this species in the state variable array
248 react_(i_spec) = chem_spec_data%gas_state_id(spec_name)
249
250 ! Make sure the species exists
251 call assert_msg(345360993, react_(i_spec).gt.0, &
252 "Missing CMAQ H2O2 reactant: "//spec_name)
253
254 ! Get properties included with this reactant in the reaction data
255 call assert(796763915, reactants%get_property_t(val=spec_props))
256 key_name = "qty"
257 if (spec_props%get_int(key_name, temp_int)) then
258 do i_qty = 1, temp_int - 1
259 react_(i_spec + i_qty) = react_(i_spec)
260 end do
261 i_spec = i_spec + temp_int - 1
262 end if
263
264 call reactants%iter_next()
265 i_spec = i_spec + 1
266 end do
267
268 ! Get the indices and chemical properties for the products
269 call products%iter_reset()
270 i_spec = 1
271 do while (products%get_key(spec_name))
272
273 ! Save the index of this species in the state variable array
274 prod_(i_spec) = chem_spec_data%gas_state_id(spec_name)
275
276 ! Make sure the species exists
277 call assert_msg(234948182, prod_(i_spec).gt.0, &
278 "Missing CMAQ H2O2 product: "//spec_name)
279
280 ! Get properties included with this product in the reaction data
281 call assert(267035567, products%get_property_t(val=spec_props))
282 key_name = "yield"
283 if (spec_props%get_real(key_name, temp_real)) then
284 yield_(i_spec) = temp_real
285 else
286 yield_(i_spec) = 1.0
287 end if
288
289 call products%iter_next()
290 i_spec = i_spec + 1
291 end do
292
293 end subroutine initialize
294
295!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
296
297 !> Finalize the reaction
298 subroutine finalize(this)
299
300 !> Reaction data
301 type(rxn_cmaq_h2o2_t), intent(inout) :: this
302
303 if (associated(this%property_set)) &
304 deallocate(this%property_set)
305 if (allocated(this%condensed_data_real)) &
306 deallocate(this%condensed_data_real)
307 if (allocated(this%condensed_data_int)) &
308 deallocate(this%condensed_data_int)
309
310 end subroutine finalize
311
312!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
313
314 !> Finalize an array of reactions
315 subroutine finalize_array(this)
316
317 !> Array of reaction data
318 type(rxn_cmaq_h2o2_t), intent(inout) :: this(:)
319
320 integer(kind=i_kind) :: i_rxn
321
322 do i_rxn = 1, size(this)
323 call finalize(this(i_rxn))
324 end do
325
326 end subroutine finalize_array
327
328!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
329
330end module camp_rxn_cmaq_h2o2
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_CMAQ_H2O2_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