CAMP 1.0.0
Chemistry Across Multiple Phases
rxn_photolysis.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_photolysis module.
7
8!> \page camp_rxn_photolysis CAMP: Photolysis
9!!
10!! Photolysis reactions take the form:
11!!
12!! \f[\ce{
13!! X + h $\nu$ -> Y_1 ( + Y_2 \dots )
14!! }\f]
15!!
16!! where \f$\ce{X}\f$ is the species being photolyzed, and
17!! \f$\ce{Y_n}\f$ are the photolysis products.
18!!
19!! Photolysis rate constants (including the \f$\ce{h $\nu$}\f$ term)
20!! can be constant or set from an external photolysis module using the
21!! \c camp_rxn_photolysis::rxn_update_data_photolysis_t object.
22!! External modules can use the
23!! \c camp_rxn_photolysis::rxn_photolysis_t::get_property_set() function during
24!! initilialization to access any needed reaction parameters to identify
25!! certain photolysis reactions.
26!! An \c camp_rxn_photolysis::update_data_photolysis_t object should be
27!! initialized for each photolysis reaction. These objects can then be used
28!! during solving to update the photolysis rate from an external module.
29!!
30!! Input data for photolysis reactions have the following format :
31!! \code{.json}
32!! {
33!! "type" : "PHOTOLYSIS",
34!! "reactants" : {
35!! "spec1" : {}
36!! },
37!! "products" : {
38!! "spec2" : {},
39!! "spec3" : { "yield" : 0.65 },
40!! ...
41!! },
42!! "scaling factor" : 1.2,
43!! ...
44!! }
45!! \endcode
46!! The key-value pairs \b reactants, and \b products are required. There must
47!! be exactly one key-value pair in the \b reactants object whose name is the
48!! species being photolyzed and whose value is an empty \c json object. Any
49!! number of products may be present. Products without a specified \b yield
50!! are assumed to have a \b yield of 1.0. The \b scaling \b factor is
51!! optional, and can be used to set a constant scaling factor for the rate
52!! constant. When the \b scaling \b factor is not provided, it is assumed to
53!! be 1.0. All other data is optional and will be available to external
54!! photolysis modules during initialization. Rate constants are in units of
55!! \f$s^{-1}\f$, and must be set using a \c rxn_photolysis_update_data_t
56!! object.
57
58!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
59
60!> The rxn_photolysis_t type and associated functions.
62
66 use camp_constants, only: const
68 use camp_mpi
71 use camp_util, only: i_kind, dp, to_string, &
73
74 use iso_c_binding
75
76 implicit none
77 private
78
79#define NUM_REACT_ this%condensed_data_int(1)
80#define NUM_PROD_ this%condensed_data_int(2)
81#define RXN_ID_ this%condensed_data_int(3)
82#define SCALING_ this%condensed_data_real(1)
83#define NUM_INT_PROP_ 3
84#define NUM_REAL_PROP_ 1
85#define NUM_ENV_PARAM_ 2
86#define REACT_(x) this%condensed_data_int(NUM_INT_PROP_ + x)
87#define PROD_(x) this%condensed_data_int(NUM_INT_PROP_ + NUM_REACT_ + x)
88#define DERIV_ID_(x) this%condensed_data_int(NUM_INT_PROP_ + NUM_REACT_ + NUM_PROD_ + x)
89#define JAC_ID_(x) this%condensed_data_int(NUM_INT_PROP_ + 2*(NUM_REACT_+NUM_PROD_) + x)
90#define YIELD_(x) this%condensed_data_real(NUM_REAL_PROP_ + x)
91
93
94 !> Generic test reaction data type
95 type, extends(rxn_data_t) :: rxn_photolysis_t
96 contains
97 !> Reaction initialization
98 procedure :: initialize
99 !> Get the reaction property set
100 procedure :: get_property_set
101 !> Initialize update data
103 !> Finalize the reaction
105 end type rxn_photolysis_t
106
107 !> Constructor for rxn_photolysis_t
109 procedure :: constructor
110 end interface rxn_photolysis_t
111
112 !> Photolysis rate update object
114 private
115 !> Flag indicating whether the update data as been allocated
116 logical :: is_malloced = .false.
117 !> Unique id for finding reactions during model initialization
118 integer(kind=i_kind) :: rxn_unique_id = 0
119 contains
120 !> Update the rate data
121 procedure :: set_rate => update_data_rate_set
122 !> Determine the pack size of the local update data
124 !> Pack the local update data to a binary
125 procedure :: internal_bin_pack
126 !> Unpack the local update data from a binary
128 !> Finalize the rate update data
131
132 !> Interface to c reaction functions
133 interface
134
135 !> Allocate space for a rate update
136 function rxn_photolysis_create_rate_update_data() result (update_data) &
137 bind (c)
138 use iso_c_binding
139 !> Allocated update_data object
140 type(c_ptr) :: update_data
142
143 !> Set a new photolysis rate
144 subroutine rxn_photolysis_set_rate_update_data(update_data, photo_id, &
145 base_rate) bind (c)
146 use iso_c_binding
147 !> Update data
148 type(c_ptr), value :: update_data
149 !> Photo id
150 integer(kind=c_int), value :: photo_id
151 !> New pre-scaling base photolysis rate
152 real(kind=c_double), value :: base_rate
154
155 !> Free an update rate data object
156 pure subroutine rxn_free_update_data(update_data) bind (c)
157 use iso_c_binding
158 !> Update data
159 type(c_ptr), value, intent(in) :: update_data
160 end subroutine rxn_free_update_data
161
162 end interface
163
164contains
165
166!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
167
168 !> Constructor for Photolysis reaction
169 function constructor() result(new_obj)
170
171 !> A new reaction instance
172 type(rxn_photolysis_t), pointer :: new_obj
173
174 allocate(new_obj)
175 new_obj%rxn_phase = gas_rxn
176
177 end function constructor
178
179!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
180
181 !> Initialize the reaction data, validating component data and loading
182 !! any required information into the condensed data arrays for use during
183 !! solving
184 subroutine initialize(this, chem_spec_data, aero_phase, aero_rep, n_cells)
185
186 !> Reaction data
187 class(rxn_photolysis_t), intent(inout) :: this
188 !> Chemical species data
189 type(chem_spec_data_t), intent(in) :: chem_spec_data
190 !> Aerosol phase data
191 type(aero_phase_data_ptr), intent(in) :: aero_phase(:)
192 !> Aerosol representations
193 type(aero_rep_data_ptr), pointer, intent(in) :: aero_rep(:)
194 !> Number of grid cells to solve simultaneously
195 integer(kind=i_kind), intent(in) :: n_cells
196
197 type(property_t), pointer :: spec_props, reactants, products
198 character(len=:), allocatable :: key_name, spec_name
199 integer(kind=i_kind) :: i_spec, i_qty
200
201 integer(kind=i_kind) :: temp_int
202 real(kind=dp) :: temp_real
203
204 ! Get the species involved
205 if (.not. associated(this%property_set)) call die_msg(408416753, &
206 "Missing property set needed to initialize reaction")
207 key_name = "reactants"
208 call assert_msg(415586594, &
209 this%property_set%get_property_t(key_name, reactants), &
210 "Photolysis reaction is missing reactants")
211 key_name = "products"
212 call assert_msg(180421290, &
213 this%property_set%get_property_t(key_name, products), &
214 "Photolysis reaction is missing products")
215
216 ! Count the number of reactants (including those with a qty specified)
217 call reactants%iter_reset()
218 i_spec = 0
219 do while (reactants%get_key(spec_name))
220 ! Get properties included with this reactant in the reaction data
221 call assert(240165383, reactants%get_property_t(val=spec_props))
222 key_name = "qty"
223 if (spec_props%get_int(key_name, temp_int)) i_spec = i_spec+temp_int-1
224 call reactants%iter_next()
225 i_spec = i_spec + 1
226 end do
227
228 ! Allocate space in the condensed data arrays
229 ! Space in this example is allocated for two sets of inidices for the
230 ! reactants and products, one molecular property for each reactant,
231 ! yields for the products and three reaction parameters.
232 allocate(this%condensed_data_int(num_int_prop_ + &
233 (i_spec + 2) * (i_spec + products%size())))
234 allocate(this%condensed_data_real(num_real_prop_ + products%size()))
235 this%condensed_data_int(:) = int(0, kind=i_kind)
236 this%condensed_data_real(:) = real(0.0, kind=dp)
237
238 ! Save space for the environment-dependent parameters
239 this%num_env_params = num_env_param_
240
241 ! Save the size of the reactant and product arrays (for reactions where
242 ! these can vary)
243 num_react_ = i_spec
244 num_prod_ = products%size()
245
246 ! Get reaction parameters (it might be easiest to keep these at the
247 ! beginning of the condensed data array, so they can be accessed using
248 ! compliler flags)
249 key_name = "scaling factor"
250 if (.not. this%property_set%get_real(key_name, scaling_)) then
251 scaling_ = real(1.0, kind=dp)
252 end if
253
254 ! Get the indices and chemical properties for the reactants
255 call reactants%iter_reset()
256 i_spec = 1
257 do while (reactants%get_key(spec_name))
258
259 ! Save the index of this species in the state variable array
260 react_(i_spec) = chem_spec_data%gas_state_id(spec_name)
261
262 ! Make sure the species exists
263 call assert_msg(747277322, react_(i_spec).gt.0, &
264 "Missing photolysis reactant: "//spec_name)
265
266 ! Get properties included with this reactant in the reaction data
267 call assert(231542303, reactants%get_property_t(val=spec_props))
268 key_name = "qty"
269 if (spec_props%get_int(key_name, temp_int)) then
270 do i_qty = 1, temp_int - 1
271 react_(i_spec + i_qty) = react_(i_spec)
272 end do
273 i_spec = i_spec + temp_int - 1
274 end if
275
276 call reactants%iter_next()
277 i_spec = i_spec + 1
278 end do
279
280 ! Make sure exactly one reactant is present
281 call assert_msg(908486656, i_spec.eq.2, "Incorrect number of reactants"//&
282 " for Photolysis reaction: "//to_string(i_spec-1))
283
284 ! Get the indices and chemical properties for the products
285 call products%iter_reset()
286 i_spec = 1
287 do while (products%get_key(spec_name))
288
289 ! Save the index of this species in the state variable array
290 prod_(i_spec) = chem_spec_data%gas_state_id(spec_name)
291
292 ! Make sure the species exists
293 call assert_msg(360988742, prod_(i_spec).gt.0, &
294 "Missing photolysis product: "//spec_name)
295
296 ! Get properties included with this product in the reaction data
297 call assert(173703744, products%get_property_t(val=spec_props))
298 key_name = "yield"
299 if (spec_props%get_real(key_name, temp_real)) then
300 yield_(i_spec) = temp_real
301 else
302 yield_(i_spec) = 1.0
303 end if
304
305 call products%iter_next()
306 i_spec = i_spec + 1
307 end do
308
309 ! Initialize the reaction id
310 rxn_id_ = -1
311
312 end subroutine initialize
313
314!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
315
316 !> Get the reaction properties. (For use by external photolysis modules.)
317 function get_property_set(this) result(prop_set)
318
319 !> Reaction properties
320 type(property_t), pointer :: prop_set
321 !> Reaction data
322 class(rxn_photolysis_t), intent(in) :: this
323
324 prop_set => this%property_set
325
326 end function get_property_set
327
328!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
329
330 !> Finalize the reaction
331 subroutine finalize(this)
332
333 !> Reaction data
334 type(rxn_photolysis_t), intent(inout) :: this
335
336 if (associated(this%property_set)) &
337 deallocate(this%property_set)
338 if (allocated(this%condensed_data_real)) &
339 deallocate(this%condensed_data_real)
340 if (allocated(this%condensed_data_int)) &
341 deallocate(this%condensed_data_int)
342
343 end subroutine finalize
344
345!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
346
347 !> Finalize an array of reactions
348 subroutine finalize_array(this)
349
350 !> Array of reaction data
351 type(rxn_photolysis_t), intent(inout) :: this(:)
352
353 integer(kind=i_kind) :: i
354
355 do i = 1, size(this)
356 call finalize(this(i))
357 end do
358
359 end subroutine finalize_array
360
361!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
362
363 !> Set packed update data for photolysis rate constants
364 subroutine update_data_rate_set(this, base_rate)
365
366 !> Update data
367 class(rxn_update_data_photolysis_t), intent(inout) :: this
368 !> Updated pre-scaling photolysis rate
369 real(kind=dp), intent(in) :: base_rate
370
371 call rxn_photolysis_set_rate_update_data(this%get_data(), &
372 this%rxn_unique_id, base_rate)
373
374 end subroutine update_data_rate_set
375
376!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
377
378 !> Initialize update data
379 subroutine update_data_initialize(this, update_data, rxn_type)
380
381 use camp_rand, only : generate_int_id
382
383 !> The reaction to update
384 class(rxn_photolysis_t), intent(inout) :: this
385 !> Update data object
386 class(rxn_update_data_photolysis_t), intent(out) :: update_data
387 !> Reaction type id
388 integer(kind=i_kind), intent(in) :: rxn_type
389
390 ! If a reaction id has not yet been generated, do it now
391 if (rxn_id_.eq.-1) then
392 rxn_id_ = generate_int_id()
393 endif
394
395 update_data%rxn_unique_id = rxn_id_
396 update_data%rxn_type = int(rxn_type, kind=c_int)
397 update_data%update_data = rxn_photolysis_create_rate_update_data()
398 update_data%is_malloced = .true.
399
400 end subroutine update_data_initialize
401
402!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
403
404 !> Determine the size of a binary required to pack the reaction data
405 integer(kind=i_kind) function internal_pack_size(this, comm) &
406 result(pack_size)
407
408 !> Reaction update data
409 class(rxn_update_data_photolysis_t), intent(in) :: this
410 !> MPI communicator
411 integer, intent(in) :: comm
412
413 pack_size = &
414 camp_mpi_pack_size_logical(this%is_malloced, comm) + &
415 camp_mpi_pack_size_integer(this%rxn_unique_id, comm)
416
417 end function internal_pack_size
418
419!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
420
421 !> Pack the given value to the buffer, advancing position
422 subroutine internal_bin_pack(this, buffer, pos, comm)
423
424 !> Reaction update data
425 class(rxn_update_data_photolysis_t), intent(in) :: this
426 !> Memory buffer
427 character, intent(inout) :: buffer(:)
428 !> Current buffer position
429 integer, intent(inout) :: pos
430 !> MPI communicator
431 integer, intent(in) :: comm
432
433#ifdef CAMP_USE_MPI
434 integer :: prev_position
435
436 prev_position = pos
437 call camp_mpi_pack_logical(buffer, pos, this%is_malloced, comm)
438 call camp_mpi_pack_integer(buffer, pos, this%rxn_unique_id, comm)
439 call assert(649543400, &
440 pos - prev_position <= this%pack_size(comm))
441#endif
442
443 end subroutine internal_bin_pack
444
445!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
446
447 !> Unpack the given value from the buffer, advancing position
448 subroutine internal_bin_unpack(this, buffer, pos, comm)
449
450 !> Reaction update data
451 class(rxn_update_data_photolysis_t), intent(inout) :: this
452 !> Memory buffer
453 character, intent(inout) :: buffer(:)
454 !> Current buffer position
455 integer, intent(inout) :: pos
456 !> MPI communicator
457 integer, intent(in) :: comm
458
459#ifdef CAMP_USE_MPI
460 integer :: prev_position
461
462 prev_position = pos
463 call camp_mpi_unpack_logical(buffer, pos, this%is_malloced, comm)
464 call camp_mpi_unpack_integer(buffer, pos, this%rxn_unique_id, comm)
465 call assert(254749806, &
466 pos - prev_position <= this%pack_size(comm))
467 this%update_data = rxn_photolysis_create_rate_update_data()
468#endif
469
470 end subroutine internal_bin_unpack
471
472!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
473
474 !> Finalize an update data object
475 subroutine update_data_finalize(this)
476
477 !> Update data object to free
478 type(rxn_update_data_photolysis_t), intent(inout) :: this
479
480 if (this%is_malloced) call rxn_free_update_data(this%update_data)
481
482 end subroutine update_data_finalize
483
484!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
485
486 !> Finalize an array of update data objects
488
489 !> Array of update data objects to free
490 type(rxn_update_data_photolysis_t), intent(inout) :: this(:)
491
492 integer(kind=i_kind) :: i
493
494 do i = 1, size(this)
495 call update_data_finalize(this(i))
496 end do
497
498 end subroutine update_data_finalize_array
499
500!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
501
502end module camp_rxn_photolysis
Initialize the aerosol representation data, validating component data and loading any required inform...
Extending-type binary pack function (Internal use only)
Extending-type binary unpack function (Internal use only)
Extending-type binary pack size (internal use only)
Get the non-unique name of a chemical species by its unique name.
Free an update rate data object.
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.
class(property_t) function, pointer get_property_set(this)
Get the aerosol phase property set.
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.
integer(kind=i_kind) function pack_size(this, comm)
Determine the size of a binary required to pack the aerosol representation 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
Wrapper functions for MPI.
Definition mpi.F90:13
subroutine camp_mpi_pack_logical(buffer, position, val, comm)
Packs the given value into the buffer, advancing position.
Definition mpi.F90:792
subroutine camp_mpi_unpack_integer(buffer, position, val, comm)
Unpacks the given value from the buffer, advancing position.
Definition mpi.F90:1023
subroutine camp_mpi_unpack_logical(buffer, position, val, comm)
Unpacks the given value from the buffer, advancing position.
Definition mpi.F90:1131
integer function camp_mpi_pack_size_logical(val, comm)
Determines the number of bytes required to pack the given value.
Definition mpi.F90:484
subroutine camp_mpi_pack_integer(buffer, position, val, comm)
Packs the given value into the buffer, advancing position.
Definition mpi.F90:691
integer function camp_mpi_pack_size_integer(val, comm)
Determines the number of bytes required to pack the given value.
Definition mpi.F90:398
The property_t structure and associated subroutines.
Definition property.F90:9
Random number generators.
Definition rand.F90:9
integer(kind=i_kind) function generate_int_id()
Generate an integer id Ids will be sequential, and can only be generated by the primary process.
Definition rand.F90:435
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_photolysis_t type and associated functions.
subroutine update_data_initialize(this, update_data, rxn_type)
Initialize update data.
subroutine update_data_rate_set(this, base_rate)
Set packed update data for photolysis rate constants.
subroutine update_data_finalize(this)
Finalize an update data object.
subroutine update_data_finalize_array(this)
Finalize an array of update data objects.
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.
Abstract reaction data type.
Definition rxn_data.F90:99
Generic test reaction data type.