CAMP 1.0.0
Chemistry Across Multiple Phases
camp_box_model_data.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_camp_box_model_data_t type and related functions
7
8!> A simple box model for \ref index "CAMP" mechanisms
9!! \todo{ Modify this to run full PartMC scenarios with CAMP enabled
10!! so that aerosol physical processes can be included. }
12
17 use camp_util
18
19 implicit none
20 private
21
22 public :: camp_box_model_data_t
23
24 ! New-line character
25 character(len=*), parameter :: new_line = char(10)
26 ! File unit to use for scipts output
27 integer(kind=i_kind), parameter :: scripts_file_unit = 8
28
29 !> Property time profile
30 type :: profile_t
31 !> Name of the profile
32 character(len=:), allocatable :: name
33 !> Current time step in the profile
34 integer(kind=i_kind) :: curr_time_step
35 !> Index for the current position in the transition arrays
36 integer(kind=i_kind) :: curr_transition
37 !> Time step for the next transition, starting with 0 (initial conditions)
38 integer(kind=i_kind), allocatable :: transition_time_step(:)
39 !> New values for each transition, starting with the initial value
40 real(kind=dp), allocatable :: transition_value(:)
41 contains
42 !> Reset the profile
43 procedure :: reset
44 !> Advance the profile by one time step
45 procedure :: advance
46 !> Get the current profile value
47 procedure :: current_value
48 !> Print the profile configuration
49 procedure :: print => profile_do_print
50 end type profile_t
51
52#ifdef CAMP_USE_JSON
53 !> Constructor for profile_t
54 interface profile_t
55 procedure :: profile_constructor
56 end interface profile_t
57#endif
58
59 !> Reaction rate time profile
60 type, extends(profile_t) :: rxn_profile_t
61 !> Reaction update object
62 class(rxn_update_data_t), pointer :: update_data
63 contains
64 !> Update the reaction rate
65 procedure :: update_rxn
66 !> Finalize the reaction profile
68 end type rxn_profile_t
69
70#ifdef CAMP_USE_JSON
71 !> Constructor for rxn_profile_t
72 interface rxn_profile_t
73 procedure :: rxn_profile_constructor
74 end interface rxn_profile_t
75#endif
76
77 !> Aerosol representation time profile
78 type, extends(profile_t) :: aero_rep_profile_t
79 !> Index for update type
80 integer(kind=i_kind) :: update_type
81 !> Section id (modal/binned)
82 integer(kind=i_kind) :: section_id = 0
83 !> Aerosol representation update object
84 class(aero_rep_update_data_t), pointer :: update_data
85 contains
86 !> Update the aerosol representation property of interest
87 procedure :: update_aero_rep
88 !> Finalize the profile
90 end type aero_rep_profile_t
91
92#ifdef CAMP_USE_JSON
93 !> Constructor for aero_rep_profile_t
95 procedure :: aero_rep_profile_constructor
96 end interface aero_rep_profile_t
97#endif
98
99 !> CAMP Box model
101 private
102 !> CAMP core
103 type(camp_core_t), pointer :: camp_core
104 !> CAMP state
105 type(camp_state_t), pointer :: camp_state
106 !> Initial CAMP state
107 type(camp_state_t), pointer :: initial_camp_state
108 !> State species names
109 type(string_t), allocatable :: spec_names(:)
110 !> Number of time steps
111 integer(kind=i_kind) :: num_steps
112 !> Time step length [s]
113 real(kind=dp) :: time_step__s
114 !> Total integration time [s]
115 real(kind=dp) :: total_time__s
116 !> Temperature [K] profile
117 type(profile_t) :: temperature__k
118 !> Pressure [Pa] profile
119 type(profile_t) :: pressure__pa
120 !> Reaction profiles
121 type(rxn_profile_t), allocatable :: rxn_profiles(:)
122 !> Aerosol representation profiles
123 type(aero_rep_profile_t), allocatable :: aero_rep_profiles(:)
124 contains
125 !> Run the box model
126 procedure :: run
127 !> Create a plotting script for the results
129 !> Print the box model configuration
130 procedure :: print => do_print
131 !> Finalize the box model
132 final :: finalize
133 end type camp_box_model_data_t
134
135 !> Constructor for camp_box_model_data_t
137 procedure :: constructor
138 end interface camp_box_model_data_t
139
140contains
141
142!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
143
144 !> Constructor for the CAMP box model
145 function constructor( config_file ) result( new_obj )
146
147#ifdef CAMP_USE_JSON
148 use json_module
149#endif
150
151 !> Pointer to a new box model object
152 type(camp_box_model_data_t), pointer :: new_obj
153 !> Box model configuration file
154 character(len=*), intent(in) :: config_file
155#ifdef CAMP_USE_JSON
156 type(json_core), target :: json
157 type(json_file) :: j_file
158 type(json_value), pointer :: j_obj, j_next, j_box_config
159
160 logical(kind=json_lk) :: found, valid
161 character(kind=json_ck, len=:), allocatable :: unicode_str_val
162 character(kind=json_ck, len=:), allocatable :: json_err_msg
163 character(kind=json_ck, len=:), allocatable :: spec_name
164 integer(kind=json_ik) :: int_value
165 real(kind=json_rk) :: real_value
166
167 logical :: file_exists
168 integer(kind=i_kind) :: num_rates, i_rate
169 integer(kind=i_kind) :: num_aero_reps, i_aero_rep
170 integer(kind=i_kind) :: spec_id
171
172 ! Create the new box model
173 allocate( new_obj )
174
175 ! Create the CAMP core
176 new_obj%camp_core => camp_core_t( config_file )
177 call new_obj%camp_core%initialize( )
178
179 ! Get the state species names
180 new_obj%spec_names = new_obj%camp_core%unique_names( )
181
182 ! Load the configuration data from the json file
183 call j_file%initialize( )
184 call j_file%get_core( json )
185 json_err_msg = ""
186 call assert_msg( 135902099, trim( config_file ).ne."", &
187 "Received empty string for file path" )
188 inquire( file = trim( config_file ), exist = file_exists )
189 call assert_msg( 181758805, file_exists, "Cannot find file: "// &
190 trim( config_file ) )
191 call j_file%load_file( filename = trim( config_file ) )
192
193 ! Get the box model configuration data
194 call j_file%get( "camp-box-model", j_box_config, found )
195 call assert_msg( 198537480, found, &
196 "Missing 'camp-box-model' configuration data "// &
197 "in configuration file.")
198 call json%validate( j_box_config, valid, json_err_msg )
199 call assert_msg( 992222859, valid, &
200 "Bad JSON format in 'camp-box-model' configuration "// &
201 "data: "//trim( json_err_msg ) )
202
203 ! Get the number of time steps
204 call json%get_child( j_box_config, "output time step [s]", j_obj, found )
205 call assert_msg( 912341318, found, &
206 "Missing 'output time step [s]' "// &
207 "in box model data" )
208 call json%get( j_obj, real_value )
209 new_obj%time_step__s = real_value
210
211 ! Get the integration time
212 call json%get_child( j_box_config, "total integration time [s]", j_obj, &
213 found )
214 call assert_msg( 442175882, found, &
215 "Missing 'total integration time [s]' "// &
216 "in box model data" )
217 call json%get( j_obj, real_value )
218 new_obj%num_steps = ceiling( real_value / new_obj%time_step__s )
219 new_obj%total_time__s = real_value
220
221 ! Get the temperature profile
222 call json%get_child( j_box_config, "temperature [K]", j_obj, found )
223 call assert_msg( 212013359, found, &
224 "Missing 'temperature [K]' "// &
225 "in box model data" )
226 new_obj%temperature__K = profile_t( json, j_obj )
227
228 ! Get the pressure profile
229 call json%get_child( j_box_config, "pressure [Pa]", j_obj, found )
230 call assert_msg( 115322763, found, &
231 "Missing 'pressure [Pa]' "// &
232 "in box model data" )
233 new_obj%pressure__Pa = profile_t( json, j_obj )
234
235 ! Get the reaction rate profiles
236 call json%get_child( j_box_config, "rates", j_obj, found )
237 if( found ) then
238 call json%info( j_obj, n_children = num_rates )
239 allocate( new_obj%rxn_profiles( num_rates ) )
240 j_next => null( )
241 call json%get( j_box_config, "rates(1)", j_obj, found )
242 call assert( 693229278, found )
243 i_rate = 0
244 do while( associated( j_obj ) )
245 i_rate = i_rate + 1
246 new_obj%rxn_profiles( i_rate ) = &
247 rxn_profile_t( new_obj%camp_core, json, j_obj )
248 call json%get_next( j_obj, j_next )
249 j_obj => j_next
250 end do
251 call assert( 763815888, i_rate .eq. num_rates)
252 end if
253
254 ! Get the aerosol representation profiles
255 call json%get_child( j_box_config, "aerosol representations", j_obj, &
256 found )
257 if( found ) then
258 call json%info( j_obj, n_children = num_aero_reps )
259 allocate( new_obj%aero_rep_profiles( num_aero_reps ) )
260 j_next => null( )
261 call json%get( j_box_config, "aerosol representations(1)", j_obj, &
262 found )
263 call assert( 877171198, found )
264 i_aero_rep = 0
265 do while( associated( j_obj ) )
266 i_aero_rep = i_aero_rep + 1
267 new_obj%aero_rep_profiles( i_aero_rep ) = &
268 aero_rep_profile_t( new_obj%camp_core, json, j_obj )
269 call json%get_next( j_obj, j_next )
270 j_obj => j_next
271 end do
272 call assert( 472301285, i_aero_rep .eq. num_aero_reps )
273 end if
274
275 ! Initialize the solver now that all the update data objects are attached
276 ! to reactions
277 call new_obj%camp_core%solver_initialize( )
278
279 ! Create a CAMP state for the box model runs and initial conditions
280 new_obj%camp_state => new_obj%camp_core%new_state( )
281 new_obj%initial_camp_state => new_obj%camp_core%new_state( )
282
283 ! Get the initial species concentrations
284 call json%get( j_box_config, "initial state(1)", j_obj, found )
285 if( found ) then
286 j_next => null( )
287 do while( associated( j_obj ) )
288 call json%info( j_obj, name = spec_name )
289 call json%get( j_obj, real_value )
290 call assert_msg( 301872431, &
291 new_obj%camp_core%spec_state_id( spec_name, &
292 spec_id ), &
293 "Cannot find species '"//trim( spec_name )// &
294 "' in the chemical mechanism" )
295 new_obj%initial_camp_state%state_var( spec_id ) = real_value
296 call json%get_next( j_obj, j_next )
297 j_obj => j_next
298 end do
299 end if
300
301 ! Set the initial environtmental conditions
302 call new_obj%initial_camp_state%env_states(1)%set_temperature_K( &
303 new_obj%temperature__K%current_value( ) )
304 call new_obj%initial_camp_state%env_states(1)%set_pressure_Pa( &
305 new_obj%pressure__Pa%current_value( ) )
306 call new_obj%initial_camp_state%update_env_state( )
307
308 ! free the json file
309 call j_file%destroy( )
310 call json%destroy( )
311
312#else
313 new_obj => null( )
314 call die_msg( 444795237, "JSON must be enabled for the CAMP box model" )
315#endif
316
317 end function constructor
318
319!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
320
321 !> Run the camp-chem box model
322 subroutine run( this, output_file_unit )
323
325
326 !> CAMP box model
327 class(camp_box_model_data_t), intent(inout) :: this
328 !> Output file unit
329 integer(kind=i_kind), intent(in), optional :: output_file_unit
330
331 type(solver_stats_t) :: solver_stats
332
333 integer(kind=i_kind) :: f_unit
334 integer(kind=i_kind) :: i_time, i_rate
335 real(kind=dp) :: step_size__s, curr_time__s
336
337 f_unit = 6
338 if( present( output_file_unit ) ) f_unit = output_file_unit
339
340 ! Initialize the model state
341 this%camp_state%state_var( : ) = this%initial_camp_state%state_var( : )
342 this%camp_state%env_var( : ) = this%initial_camp_state%env_var( : )
343
344 ! Output file header and intitial conditions
345 write(f_unit,*) 0.0, this%camp_state%env_var( : ), &
346 this%camp_state%state_var( : )
347
348 ! Reset the profiles
349 call this%temperature__K%reset( )
350 call this%pressure__Pa%reset( )
351 do i_rate = 1, size( this%rxn_profiles )
352 call this%rxn_profiles( i_rate )%reset( )
353 end do
354
355 ! Solve the box model
356 do i_time = 1, this%num_steps
357
358 ! Set the current time and step size
359 curr_time__s = i_time * this%time_step__s
360 if( curr_time__s .gt. this%total_time__s ) then
361 step_size__s = this%time_step__s - &
362 (curr_time__s - this%total_time__s )
363 else
364 step_size__s = this%time_step__s
365 end if
366
367 ! Update the model conditions
368 call this%temperature__K%advance( )
369 call this%camp_state%env_states(1)%set_temperature_K( &
370 this%temperature__K%current_value( ) )
371 call this%pressure__Pa%advance( )
372 call this%camp_state%env_states(1)%set_pressure_Pa( &
373 this%pressure__Pa%current_value( ) )
374 do i_rate = 1, size( this%rxn_profiles )
375 call this%rxn_profiles( i_rate )%advance( )
376 call this%rxn_profiles( i_rate )%update_rxn( this%camp_core )
377 end do
378
379 ! Solve the chemistry
380 call this%camp_core%solve( this%camp_state, step_size__s, &
381 solver_stats = solver_stats )
382
383 ! Output the model state
384 write(f_unit,*) curr_time__s, this%camp_state%env_var( : ), &
385 this%camp_state%state_var( : )
386
387 end do
388
389 end subroutine run
390
391!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
392
393 !> Create a gnuplot configuration file for plotting box model results
394 subroutine create_gnuplot_config_file( this, file_prefix )
395
396 !> Box model data
397 class(camp_box_model_data_t), intent(in) :: this
398 !> Prefix for the results file (assumes a results file name of
399 !> `'file_prefix'_results.txt`)
400 character(len=*), intent(in) :: file_prefix
401
402 character(len=:), allocatable :: spec_name, file_name
403 integer(kind=i_kind) :: f_unit, i_char, i_spec
404
405 f_unit = scripts_file_unit
406
407 ! Create the gnuplot script
408 file_name = file_prefix//".conf"
409 open(unit=f_unit, file=file_name, status="replace", action="write")
410 write(f_unit,*) "# GNU Plot configuration for CAMP box model output"
411 write(f_unit,*) "# Run as: gnuplot "//file_name
412 write(f_unit,*) "set terminal png truecolor"
413 write(f_unit,*) "set autoscale"
414 write(f_unit,*) "set xrange [ 0.0:", this%total_time__s, "]"
415
416 ! output environmental conditions
417 write(f_unit,*) "set output '"//file_prefix//"_temperature.png"
418 write(f_unit,*) "plot\"
419 write(f_unit,*) " '"//file_prefix//"_results.txt'\"
420 write(f_unit,*) " using 1:2 title 'Temperature [K]'"
421 write(f_unit,*) "set output '"//file_prefix//"_pressure.png"
422 write(f_unit,*) "plot\"
423 write(f_unit,*) " '"//file_prefix//"_results.txt'\"
424 write(f_unit,*) " using 1:3 title 'Pressure [Pa]'"
425
426 ! output species concentrations and parameter values
427 do i_spec = 1, size( this%spec_names )
428 spec_name = this%spec_names( i_spec )%string
429 forall( i_char = 1 : len( spec_name ), &
430 spec_name( i_char : i_char ) .eq. '/') &
431 spec_name( i_char:i_char ) = '_'
432 write(f_unit,*) "set output '"//file_prefix//"_"//spec_name//".png'"
433 write(f_unit,*) "plot\"
434 write(f_unit,*) " '"//file_prefix//"_results.txt'\"
435 write(f_unit,*) " using 1:"//trim( to_string( i_spec + 3 ) )// &
436 " title '"//spec_name//"'"
437 end do
438 close(f_unit)
439
440 deallocate( spec_name )
441 deallocate( file_name )
442
443 end subroutine create_gnuplot_config_file
444
445!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
446
447 !> Print out the configuration of the box model
448 subroutine do_print( this, file_unit )
449
450 !> Box model data
451 class(camp_box_model_data_t), intent(in) :: this
452 !> Output file unit
453 integer(kind=i_kind), intent(in), optional :: file_unit
454
455 character(len=*), parameter :: fmt_blank_line = "('')"
456 character(len=*), parameter :: fmt_state_hdr = &
457 "(' | ',A50,' | ',A13,' |')"
458 character(len=*), parameter :: fmt_state_data = &
459 "(' | ',A50,' | ',ES13.3,' |')"
460 integer(kind=i_kind) :: f_unit, i_rate, i_spec
461
462 f_unit = 6
463 if( present( file_unit ) ) f_unit = file_unit
464
465 write(f_unit,*) "**********************"
466 write(f_unit,*) "*** CAMP Box Model ***"
467 write(f_unit,*) "**********************"
468 write(f_unit,fmt_blank_line)
469 write(f_unit,*) "total integration time", this%total_time__s, "s"
470 write(f_unit,*) "time step size ", this%time_step__s, "s"
471 write(f_unit,*) "number of time steps ", this%num_steps
472 write(f_unit,fmt_blank_line)
473 write(f_unit,*) "** temperature profile **"
474 call this%temperature__K%print( f_unit )
475 write(f_unit,*) "** end temperature profile **"
476 write(f_unit,fmt_blank_line)
477 write(f_unit,*) "** pressure [Pa] profile **"
478 call this%pressure__Pa%print( f_unit )
479 write(f_unit,*) "** end pressure [Pa] profile **"
480 write(f_unit,fmt_blank_line)
481 write(f_unit,*) "****************************"
482 write(f_unit,*) "** Reaction Rate Profiles **"
483 write(f_unit,*) "****************************"
484 write(f_unit,fmt_blank_line)
485 do i_rate = 1, size( this%rxn_profiles )
486 call this%rxn_profiles( i_rate )%print( f_unit )
487 write(f_unit,fmt_blank_line)
488 end do
489 write(f_unit,*) "****************************"
490 write(f_unit,*) "** Reaction Rate Profiles **"
491 write(f_unit,*) "****************************"
492 write(f_unit,fmt_blank_line)
493 write(f_unit,*) "*******************"
494 write(f_unit,*) "** Initial State **"
495 write(f_unit,*) "*******************"
496 write(f_unit,fmt_blank_line)
497 write(f_unit,fmt_state_hdr) "Species", "Value"
498 do i_spec = 1, size( this%spec_names )
499 write(f_unit,fmt_state_data) this%spec_names( i_spec )%string, &
500 this%initial_camp_state%state_var( i_spec )
501 end do
502 write(f_unit,fmt_blank_line)
503 call this%camp_core%print( f_unit )
504 write(f_unit,fmt_blank_line)
505 write(f_unit,*) "***********************"
506 write(f_unit,*) "** End Initial State **"
507 write(f_unit,*) "***********************"
508 write(f_unit,fmt_blank_line)
509 write(f_unit,*) "**************************"
510 write(f_unit,*) "*** End CAMP Box Model ***"
511 write(f_unit,*) "**************************"
512
513 end subroutine do_print
514
515!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
516
517 !> Finalize the box model
518 subroutine finalize( this )
519
520 !> CAMP box model
521 type(camp_box_model_data_t), intent(inout) :: this
522
523 if( associated( this%camp_core ) ) &
524 deallocate( this%camp_core )
525 if( associated( this%camp_state ) ) &
526 deallocate( this%camp_state )
527 if( associated( this%initial_camp_state) ) &
528 deallocate( this%initial_camp_state )
529 if( allocated( this%spec_names ) ) &
530 deallocate( this%spec_names )
531 if( allocated( this%rxn_profiles ) ) &
532 deallocate( this%rxn_profiles )
533 if( allocated( this%aero_rep_profiles ) ) &
534 deallocate( this%aero_rep_profiles )
535
536 end subroutine finalize
537
538!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
539
540#ifdef CAMP_USE_JSON
541 !> Constructor for profile_t
542 function profile_constructor( json, j_obj ) result( new_obj )
543
544 use json_module
545
546 !> New profile_t object
547 type(profile_t) :: new_obj
548 !> JSON core
549 type(json_core), intent(inout) :: json
550 !> JSON object to build profile from
551 type(json_value), pointer, intent(inout) :: j_obj
552
553 type(json_value), pointer :: j_child, j_next, j_val
554 integer(kind=json_ik) :: num_trans
555 logical(kind=json_lk) :: found
556 character(kind=json_ck, len=:), allocatable :: profile_name
557
558 integer(kind=i_kind) :: i_trans
559
560 ! Get the name of the profile
561 call json%info( j_obj, name = profile_name )
562 new_obj%name = profile_name
563
564 ! Get the transitions
565 call json%get_child( j_obj, "transitions", j_child, found )
566 if( found ) then
567 call json%info( j_obj, n_children = num_trans )
568 allocate( new_obj%transition_time_step( 0:num_trans ) )
569 allocate( new_obj%transition_value( 0:num_trans ) )
570 j_next => null( )
571 call json%get( j_obj, "transitions(1)", j_child, found )
572 i_trans = 0
573 do while( associated( j_child ) )
574 i_trans = i_trans + 1
575
576 ! get the time step
577 call json%get_child( j_child, "time step", j_val, found )
578 call assert_msg( 325008338, found, &
579 "Missing 'time step' from element "// &
580 trim( to_string( i_trans ) )//" in transition '"// &
581 trim( new_obj%name )//"'" )
582 call json%get( j_val, new_obj%transition_time_step( i_trans ) )
583
584 ! Get the value
585 call json%get_child( j_child, "value", j_val, found )
586 call assert_msg( 483183389, found, &
587 "Missing 'value' from element "// &
588 trim( to_string( i_trans ) )//" in transition '"// &
589 trim( new_obj%name )//"'" )
590 call json%get( j_val, new_obj%transition_value( i_trans ) )
591 call json%get_next( j_child, j_next )
592 j_child => j_next
593 end do
594 call assert( 249923619, i_trans .eq. num_trans )
595 else
596 allocate( new_obj%transition_time_step( 0:0 ) )
597 allocate( new_obj%transition_value( 0:0 ) )
598 end if
599
600 ! Get the initial state
601 call json%get_child( j_obj, "initial", j_val, found )
602 call assert_msg( 382107479, found, &
603 "Missing 'initial' value in profile '"// &
604 trim( new_obj%name )//"'" )
605 call json%get( j_val, new_obj%transition_value( 0 ) )
606 new_obj%transition_time_step( 0 ) = 1
607
608 call new_obj%reset( )
609
610 end function profile_constructor
611#endif
612!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
613
614 !> Reset the profile to the initial state
615 subroutine reset( this )
616
617 !> Property profile
618 class(profile_t), intent(inout) :: this
619
620 this%curr_transition = 0
621 this%curr_time_step = 1
622
623 end subroutine reset
624
625!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
626
627 !> Advance the profile by one time step
628 subroutine advance( this )
629
630 !> Property profile
631 class(profile_t), intent(inout) :: this
632
633 this%curr_time_step = this%curr_time_step + 1
634
635 if( this%curr_transition .eq. &
636 size( this%transition_time_step ) - 1 ) return
637
638 if( this%transition_time_step( this%curr_transition + 1 ) .eq. &
639 this%curr_time_step ) &
640 this%curr_transition = this%curr_transition + 1
641
642 end subroutine advance
643
644!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
645
646 !> Get the current value of the profile
647 function current_value( this )
648
649 !> Current value of the profile
650 real(kind=dp) :: current_value
651 !> Property profile
652 class(profile_t), intent(inout) :: this
653
654 current_value = this%transition_value( this%curr_transition )
655
656 end function current_value
657
658!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
659
660 !> Print the profile configuration
661 subroutine profile_do_print( this, file_unit )
662
663 !> Property profile
664 class(profile_t), intent(in) :: this
665 !> Output file unit
666 integer(kind=i_kind), intent(in), optional :: file_unit
667
668 character(len=*), parameter :: fmt_trans_hdr = &
669 "(' | ',A13,' | ',A13,' |')"
670 character(len=*), parameter :: fmt_trans_data = &
671 "(' | ',I13,' | ',ES13.3,' |')"
672 integer(kind=i_kind) :: f_unit, i_trans
673
674 f_unit = 6
675 if( present( file_unit ) ) f_unit = file_unit
676
677 write(f_unit,*) "Profile name: "//trim(this%name)
678 write(f_unit,*) " ** Transitions **"
679 write(f_unit,fmt_trans_hdr) "Time Step", "Value"
680 do i_trans = 0, size( this%transition_time_step ) - 1
681 write(f_unit,fmt_trans_data) this%transition_time_step( i_trans ), &
682 this%transition_value( i_trans )
683 end do
684 write(f_unit,*) " ** End Transitions **"
685
686 end subroutine profile_do_print
687
688!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
689#ifdef CAMP_USE_JSON
690 !> Constructor for rxn_profile_t
691 function rxn_profile_constructor( camp_core, json, j_obj ) result( new_obj )
692
698 use json_module
699
700 !> New reaction profile
701 type(rxn_profile_t) :: new_obj
702 !> CAMP core
703 type(camp_core_t), intent(inout) :: camp_core
704 !> JSON core
705 type(json_core), intent(inout) :: json
706 !> JSON object to build the profile from
707 type(json_value), pointer, intent(inout) :: j_obj
708
709 type(mechanism_data_t), pointer :: mech
710 class(rxn_data_t), pointer :: rxn
711 type(profile_t) :: base_profile
712
713 integer(kind=i_kind) :: i_mech, i_rxn
714 character(len=:), allocatable :: rxn_label
715 logical :: found
716
717 found = .false.
718
719 ! Load the standard profile data
720 base_profile = profile_t( json, j_obj )
721 new_obj%name = base_profile%name
722 new_obj%transition_time_step = base_profile%transition_time_step
723 new_obj%transition_value = base_profile%transition_value
724
725 ! Find the reaction and initialize the update data object
726 do i_mech = 1, size( camp_core%mechanism )
727 mech => camp_core%mechanism( i_mech )%val
728 do i_rxn = 1, mech%size( )
729 rxn => mech%get_rxn( i_rxn )
730 if( .not. rxn%property_set%get_string( "camp-box-model-id", &
731 rxn_label ) ) cycle
732 if( .not. trim( rxn_label ) .eq. new_obj%name ) cycle
733 select type( rxn )
734 class is( rxn_emission_t)
735 allocate( rxn_update_data_emission_t::new_obj%update_data )
736 class is( rxn_first_order_loss_t)
737 allocate( rxn_update_data_first_order_loss_t::new_obj%update_data )
738 class is( rxn_photolysis_t)
739 allocate( rxn_update_data_photolysis_t::new_obj%update_data )
740 class is( rxn_wet_deposition_t)
741 allocate( rxn_update_data_wet_deposition_t::new_obj%update_data )
742 class default
743 call die_msg( 592455681, "Invalid reaction for rate updates" )
744 end select
745 call camp_core%initialize_update_object( rxn, new_obj%update_data )
746 found = .true.
747 exit
748 end do
749 if( found ) exit
750 end do
751
752 call assert_msg( 800298506, found, "Could not find reaction label '"// &
753 trim( new_obj%name ) )
754
755 end function rxn_profile_constructor
756#endif
757!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
758
759 !> Update a reaction with the current rate from the profile
760 subroutine update_rxn( this, camp_core )
761
766
767 !> Reaction rate profile
768 class(rxn_profile_t) :: this
769 !> CAMP core
770 type(camp_core_t), intent(inout) :: camp_core
771
772 class(rxn_update_data_t), pointer :: ud
773
774 select type( ud => this%update_data )
776 call ud%set_rate( this%current_value( ) )
778 call ud%set_rate( this%current_value( ) )
780 call ud%set_rate( this%current_value( ) )
782 call ud%set_rate( this%current_value( ) )
783 class default
784 call die( 497670619 )
785 end select
786
787 call camp_core%update_data( this%update_data )
788
789 end subroutine update_rxn
790
791!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
792
793 !> Finalize the reaction profile
794 subroutine rxn_profile_finalize( this )
795
796 !> Reaction profile
797 type(rxn_profile_t), intent(inout) :: this
798
799 if( associated( this%update_data ) ) deallocate( this%update_data )
800
801 end subroutine rxn_profile_finalize
802
803!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
804#ifdef CAMP_USE_JSON
805 !> Constructor for aero_rep_profile_t
806 function aero_rep_profile_constructor( camp_core, json, j_obj ) &
807 result( new_obj )
808
812 use json_module
813
814 !> New aerosol representation profile
815 type(aero_rep_profile_t) :: new_obj
816 !> CAMP core
817 type(camp_core_t), intent(inout) :: camp_core
818 !> JSON core
819 type(json_core), intent(inout) :: json
820 !> JSON object to build the profile from
821 type(json_value), pointer, intent(inout) :: j_obj
822
823 class(aero_rep_data_t), pointer :: aero_rep
824 type(profile_t) :: base_profile
825
826 character(kind=json_ck, len=:), allocatable :: update_type
827 character(kind=json_ck, len=:), allocatable :: section_name
828 logical(kind=json_lk) :: found
829
830 integer(kind=i_kind) :: i_aero_rep
831 character(len=:), allocatable :: aero_rep_label
832
833 found = .false.
834
835 ! Load the standard profile data
836 base_profile = profile_t( json, j_obj )
837 new_obj%name = base_profile%name
838 new_obj%transition_time_step = base_profile%transition_time_step
839 new_obj%transition_value = base_profile%transition_value
840
841 ! Get the type of value being updated
842 call json%get( j_obj, "update type", update_type, found )
843 call assert_msg( 591311579, found, &
844 "Missing update type for aerosol representation "// &
845 "profile '"//new_obj%name//"'" )
846
847 ! Find the reaction and initialize the update data object
848 do i_aero_rep = 1, size( camp_core%aero_rep )
849 aero_rep => camp_core%aero_rep( i_aero_rep )%val
850 if( .not. aero_rep%property_set%get_string( "camp-box-model-id", &
851 aero_rep_label ) ) cycle
852 if( .not. trim( aero_rep_label ) .eq. new_obj%name ) cycle
853 select type( aero_rep )
855
856 ! Modal/binned updates must include a section id
857 call json%get( j_obj, "section name", section_name, found )
858 call assert_msg( 454366901, found, &
859 "Missing section name for modal/binned "// &
860 "aerosol representation '"//new_obj%name//"'" )
861 call assert_msg( 112599854, aero_rep%get_section_id( section_name, &
862 new_obj%section_id ), &
863 "Cannot find aerosol section '"//section_name// &
864 "' in aerosol representation '"// &
865 new_obj%name//"'" )
866 select case( update_type )
867 case( "UPDATE_GMD" )
868 allocate( aero_rep_update_data_modal_binned_mass_gmd_t::new_obj%update_data )
869 case( "UPDATE_GSD" )
870 allocate( aero_rep_update_data_modal_binned_mass_gsd_t::new_obj%update_data )
871 case default
872 call die_msg( 513320382, "Invalid update type for "// &
873 "modal/binned aerosol rep: "//update_type )
874 end select
876 select case( update_type )
877 case( "UPDATE_NUMBER" )
878 allocate( aero_rep_update_data_single_particle_number_t::new_obj%update_data )
879 case default
880 call die_msg( 379396160, "Invalid update type for "// &
881 "single particle aerosol rep: "//update_type )
882 end select
883 class default
884 call die_msg( 203974949, "Invalid aerosol representation " )
885 end select
886 call camp_core%initialize_update_object( aero_rep, new_obj%update_data )
887 found = .true.
888 exit
889 end do
890
891 call assert_msg( 712992422, found, "Could not find aerosol "// &
892 "representation label '"// trim( new_obj%name ) )
893
895#endif
896!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
897
898 !> Update a reaction with the current rate from the profile
899 subroutine update_aero_rep( this, camp_core )
900
903
904 !> Aerosol representation profile
905 class(aero_rep_profile_t) :: this
906 !> CAMP core
907 type(camp_core_t), intent(inout) :: camp_core
908
909 class(aero_rep_update_data_t), pointer :: ud
910
911 select type( ud => this%update_data )
913 call assert( 832366630, this%section_id.gt.0 )
914 call ud%set_GMD( this%section_id, this%current_value( ) )
916 call assert( 264057359, this%section_id.gt.0 )
917 call ud%set_GSD( this%section_id, this%current_value( ) )
919 call ud%set_number__n_m3( 1, this%current_value( ) )
920 class default
921 call die( 623019372 )
922 end select
923
924 call camp_core%update_data( this%update_data )
925
926 end subroutine update_aero_rep
927
928!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
929
930 !> Finalize the aerosol representation profile
931 subroutine aero_rep_profile_finalize( this )
932
933 !> Aerosol representation profile
934 type(aero_rep_profile_t), intent(inout) :: this
935
936 if( associated( this%update_data ) ) deallocate( this%update_data )
937
938 end subroutine aero_rep_profile_finalize
939
940!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
941
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_rep_data_t structure and associated subroutines.
subroutine do_print(this, file_unit)
Print the aerosol representation data.
The abstract aero_rep_modal_binned_mass_t structure and associated subroutines.
The aero_rep_single_particle_t type and associated subroutines.
A simple box model for CAMP mechanisms.
subroutine rxn_profile_finalize(this)
Finalize the reaction profile.
subroutine profile_do_print(this, file_unit)
Print the profile configuration.
type(camp_box_model_data_t) function, pointer constructor(config_file)
Constructor for the CAMP box model.
character(len= *), parameter new_line
integer(kind=i_kind), parameter scripts_file_unit
subroutine finalize(this)
Finalize the box model.
subroutine update_rxn(this, camp_core)
Update a reaction with the current rate from the profile.
subroutine advance(this)
Advance the profile by one time step.
type(aero_rep_profile_t) function aero_rep_profile_constructor(camp_core, json, j_obj)
Constructor for aero_rep_profile_t.
type(rxn_profile_t) function rxn_profile_constructor(camp_core, json, j_obj)
Constructor for rxn_profile_t.
subroutine reset(this)
Reset the profile to the initial state.
real(kind=dp) function current_value(this)
Get the current value of the profile.
subroutine create_gnuplot_config_file(this, file_prefix)
Create a gnuplot configuration file for plotting box model results.
subroutine run(this, output_file_unit)
Run the camp-chem box model.
subroutine aero_rep_profile_finalize(this)
Finalize the aerosol representation profile.
type(profile_t) function profile_constructor(json, j_obj)
Constructor for profile_t.
subroutine update_aero_rep(this, camp_core)
Update a reaction with the current rate from the profile.
The camp_core_t structure and associated subroutines.
Definition camp_core.F90:87
The camp_state_t structure and associated subroutines.
Definition camp_state.F90:9
The mechanism_data_t structure and associated subroutines.
The rxn_data_t structure and associated subroutines.
Definition rxn_data.F90:60
The rxn_emission_t type and associated functions.
The rxn_first_order_loss_t type and associated functions.
The rxn_photolysis_t type and associated functions.
The rxn_wet_deposition_t type and associated functions.
The solver_stats_t type and associated subroutines.
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 die(code)
Error immediately.
Definition util.F90:184
subroutine assert_msg(code, condition_ok, error_msg)
Errors unless condition_ok is true.
Definition util.F90:130
Abstract aerosol representation data type.
Aerosol representation time profile.
Part-MC model data.
Abstract reaction data type.
Definition rxn_data.F90:98
Generic test reaction data type.
Generic test reaction data type.
String type for building arrays of string of various size.
Definition util.F90:53