Chapter 3#
MICM Box Model with Multiple Grid Cells in Fortran#
In this third MUSICA Fortran example, we will extend the MICM box model from Chapter 2 to handle multiple grid cells in a single call to the MICM solver. This demonstrates how MUSICA can efficiently solve chemical mechanisms for multiple independent air masses simultaneously.
Each grid cell represents an independent set of well-mixed air masses, allowing us to model different atmospheric conditions (temperature, pressure, concentrations) and observe how the same chemical mechanism evolves differently under various conditions.
We will use the same analytical configuration files from Chapter 2
(config.json, species.json, and reactions.json)
saved in the subdirectory configs/v0/analytical.
For reference, these files define:
A system of six chemical species: A, B, C, D, E, and F
Four reactions including two Arrhenius-type reactions and two user-defined rate reactions
The same mechanism will be solved simultaneously across multiple grid cells
To create a multiple grid cells box model, save the following Fortran code to a file named micm_multiple_grid_cells.F90:
program test_micm_multiple_grid_cells use, intrinsic :: iso_c_binding use, intrinsic :: ieee_arithmetic use iso_fortran_env, only : real64 use musica_util, only: assert, error_t, string_t, mapping_t use musica_micm, only: micm_t, solver_stats_t use musica_micm, only: Rosenbrock, RosenbrockStandardOrder use musica_state, only: conditions_t, state_t #define ASSERT( expr ) call assert( expr, __FILE__, __LINE__ ) #define ASSERT_EQ( a, b ) call assert( a == b, __FILE__, __LINE__ ) #define ASSERT_NE( a, b ) call assert( a /= b, __FILE__, __LINE__ ) #define ASSERT_NEAR( a, b, tol ) call assert( abs(a - b) < abs(a + b) * tol, __FILE__, __LINE__ ) #define ASSERT_GT( a, b ) call assert( a > b, __FILE__, __LINE__ ) implicit none call multiple_grid_cells_example() contains !> Runs a multiple grid cells box model using the MICM solver subroutine multiple_grid_cells_example() character(len=256) :: config_path integer :: solver_type integer :: num_grid_cells real(real64), parameter :: GAS_CONSTANT = 8.31446261815324_real64 ! J mol-1 K-1 real(real64) :: time_step type(string_t) :: solver_state type(solver_stats_t) :: solver_stats type(error_t) :: error type(micm_t), pointer :: micm type(state_t), pointer :: state integer :: i, j, cell_id config_path = "configs/v0/analytical" solver_type = RosenbrockStandardOrder num_grid_cells = 3 ! Use 3 grid cells to demonstrate multiple cells write(*,*) "Creating MICM solver with", num_grid_cells, "grid cells..." micm => micm_t(config_path, solver_type, error) ASSERT( error%is_success() ) write(*,*) "Creating State for multiple grid cells..." state => micm%get_state(num_grid_cells, error) ASSERT( error%is_success() ) time_step = 200 ! Set up conditions for each grid cell do cell_id = 1, num_grid_cells state%conditions(cell_id)%temperature = 273.0 + (cell_id - 1) * 10.0 ! Different temperatures state%conditions(cell_id)%pressure = 1.0e5 state%conditions(cell_id)%air_density = state%conditions(cell_id)%pressure / & (GAS_CONSTANT * state%conditions(cell_id)%temperature) end do ! Set initial concentrations for each grid cell associate( cell_stride => state%species_strides%grid_cell, & var_stride => state%species_strides%variable ) do cell_id = 1, num_grid_cells ! Set different initial concentrations for each cell do i = 1, 6 ! Assuming 6 species ! Cell 1: all species start at 1.0 ! Cell 2: all species start at 2.0 ! Cell 3: all species start at 0.5 if (cell_id == 1) then state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride) = 1.0 else if (cell_id == 2) then state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride) = 2.0 else state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride) = 0.5 end if end do end do end associate ! Set rate parameters for each grid cell associate( cell_stride => state%rate_parameters_strides%grid_cell, & var_stride => state%rate_parameters_strides%variable ) do cell_id = 1, num_grid_cells ! Set the same rate parameters for all cells state%rate_parameters(1 + (cell_id - 1) * cell_stride + 0 * var_stride) = 0.001 state%rate_parameters(1 + (cell_id - 1) * cell_stride + 1 * var_stride) = 0.002 end do end associate ! Print species information write(*,*) "Species in the mechanism:" do i = 1, state%species_ordering%size() write(*,*) "Species Name:", trim(state%species_ordering%name(i)), & ", Index:", state%species_ordering%index(i) end do ! Print initial concentrations for each grid cell write(*,*) "" write(*,*) "Initial concentrations by grid cell:" associate( cell_stride => state%species_strides%grid_cell, & var_stride => state%species_strides%variable ) do cell_id = 1, num_grid_cells write(*,'(A,I0,A,F6.1,A)') "Grid Cell ", cell_id, " (T=", & state%conditions(cell_id)%temperature, "K):" do i = 1, 6 write(*,'(A,F8.3)', advance='no') " ", & state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride) end do write(*,*) "" end do end associate write(*,*) "" write(*,*) "Solving for all grid cells simultaneously..." call micm%solve(time_step, state, solver_state, solver_stats, error) ASSERT( error%is_success() ) ! Print final concentrations for each grid cell write(*,*) "" write(*,*) "Final concentrations by grid cell:" associate( cell_stride => state%species_strides%grid_cell, & var_stride => state%species_strides%variable ) do cell_id = 1, num_grid_cells write(*,'(A,I0,A,F6.1,A)') "Grid Cell ", cell_id, " (T=", & state%conditions(cell_id)%temperature, "K):" do i = 1, 6 write(*,'(A,F8.3)', advance='no') " ", & state%concentrations(1 + (cell_id - 1) * cell_stride + (i - 1) * var_stride) end do write(*,*) "" end do end associate write(*,*) "" write(*,*) "Solver completed successfully for all", num_grid_cells, "grid cells!" deallocate( micm ) deallocate( state ) end subroutine multiple_grid_cells_example end program test_micm_multiple_grid_cells
Key differences from the single grid cell example in Chapter 2:
Grid Cell Configuration:
num_grid_cellsis set to 3 instead of 1Each grid cell can have different environmental conditions (temperature, pressure)
Each grid cell can start with different initial concentrations
Setting Conditions and Initial Values:
The state object now manages data for multiple grid cells using strides:
state%species_strides%variableprovides the stride for accessing species data across cellsstate%rate_parameters_strides%variableprovides the stride for rate parameter dataConcentrations are stored in a flat array with stride-based indexing
! Different temperatures for each cell
do cell_id = 1, num_grid_cells
state%conditions(cell_id)%temperature = 273.0 + (cell_id - 1) * 10.0
state%conditions(cell_id)%pressure = 1.0e5
state%conditions(cell_id)%air_density = state%conditions(cell_id)%pressure / &
(GAS_CONSTANT * state%conditions(cell_id)%temperature)
end do
Accessing Multi-Cell Data:
To access concentration data for a specific species in a specific grid cell:
cell_stride = state%species_strides%grid_cell
var_stride = state%species_strides%variable
cell_concentration = state%concentrations(1 + (cell_id - 1) * cell_stride + (species_index - 1) * var_stride)
Fortran’s 1-based indexing requires subtracting 1 from the 0-based cell_id
and species_index to adjust offsets.
Simultaneous Solving:
A call to micm%solve with a state that represents multiple grid cells is only
required once per time step; all grid cells are solved simultaneously in a multi-grid cell micm state.
call micm%solve(time_step, state, solver_state, solver_stats, error)
This is much more efficient than calling the solver separately for each grid cell, especially when dealing with large numbers of atmospheric grid cells in climate or weather models.
To build and run the example, follow the same build instructions as in Chapter 1 and 2. The compilation command would be:
gfortran -o micm_multiple_grid_cells micm_multiple_grid_cells.F90 -I<MUSICA_DIR>/include -L<MUSICA_DIR>/lib64 -lmusica-fortran -lmusica -lstdc++ -lyaml-cpp
Assuming you name the executable micm_multiple_grid_cells, you can run the program as follows:
$ ./micm_multiple_grid_cells
Creating MICM solver with 3 grid cells...
Creating State for multiple grid cells...
Species in the mechanism:
Species Name:A, Index: 1
Species Name:B, Index: 2
Species Name:C, Index: 5
Species Name:D, Index: 3
Species Name:E, Index: 4
Species Name:F, Index: 6
Initial concentrations by grid cell:
Grid Cell 1 (T= 273.0K):
1.000 1.000 1.000 1.000 1.000 1.000
Grid Cell 2 (T= 283.0K):
2.000 2.000 2.000 2.000 2.000 2.000
Grid Cell 3 (T= 293.0K):
0.500 0.500 0.500 0.500 0.500 0.500
Solving for all grid cells simultaneously...
Final concentrations by grid cell:
Grid Cell 1 (T= 273.0K):
0.382 1.468 0.670 1.116 1.150 1.214
Grid Cell 2 (T= 283.0K):
0.764 2.936 1.340 2.232 2.300 2.428
Grid Cell 3 (T= 293.0K):
0.191 0.734 0.335 0.558 0.575 0.607
Solver completed successfully for all 3 grid cells!
The results differ in each grid cell due to the varying conditions.