Chapter 2

An MICM Box Model Fortran Example

In this next MUSICA Fortran example, we will setup a MICM solver, starting with a set of MICM configuration files, and run the solver for a single integration time step.

The MICM configuration is specified in a top-level config.json file, which simply lists the chemical species configuration file followed by the reactions configuration file.

{"camp-files": ["species.json", "reactions.json"]}

For this example, we will have a system of three chemical species A, B, and C, defined in the JSON file species.json as follows:

{
  "camp-data":
  [
    {
      "name": "A", "type": "CHEM_SPEC"
    },
    {
      "name": "B", "type": "CHEM_SPEC"
    },
    {
      "name": "C", "type": "CHEM_SPEC"
    },
    {
      "name": "D", "type": "CHEM_SPEC"
    },
    {
      "name": "E", "type": "CHEM_SPEC"
    },
    {
      "name": "F", "type": "CHEM_SPEC"
    }
  ]
}

The reactions.json specifies a mechanism, or a set of reactions for the system. Here, we will introduce two Arrhenius type reactions, the first with B evolving to C, and specifying all five reaction parameters, and the second reaction with A evolving to B and using only two reaction parameters. The mechanism configuration might then be set up as:

{
  "camp-data":
  [
    {
      "type": "MECHANISM",
      "name": "music box interactive configuration",
      "reactions":
      [
        {
          "type": "USER_DEFINED",
          "MUSICA name" : "reaction 2",
          "reactants": {"E": {"qty": 1}},
          "products": {"F": {"yield": 1}}
        },
        {
          "type": "USER_DEFINED",
          "MUSICA name" : "reaction 1",
          "reactants": {"D": {"qty": 1}},
          "products": {"E": {"yield": 1}}
        },
        {
          "type": "ARRHENIUS", "A": 0.004, "C" : 50,
          "reactants": {"A": {"qty": 1}},
          "products": {"B": {"yield": 1}}
        },
        {
          "type": "ARRHENIUS", "A": 0.012, "B": -2, "C" : 75, "D": 50, "E": 1.0e-6,
          "reactants": {"B": {"qty": 1}},
          "products": {"C": {"yield": 1}}
        }
      ]
    }
  ]
}

More information on MICM configurations and reactions can be found in the MICM documentation at `https://ncar.github.io/micm/user_guide/`_

The Fortran example code is shown below in full:

program test_micm_box_model

  use, intrinsic :: iso_c_binding
  use, intrinsic :: ieee_arithmetic

  use musica_util, only: error_t, string_t, mapping_t
  use musica_micm, only: micm_t, solver_stats_t
  use musica_micm, only: Rosenbrock, RosenbrockStandardOrder

  implicit none

  call box_model()

contains

  subroutine box_model()

    character(len=256) :: config_path
    integer(c_int)     :: solver_type
    integer(c_int)     :: num_grid_cells

    real(c_double), parameter :: GAS_CONSTANT = 8.31446261815324_c_double ! J mol-1 K-1

    real(c_double) :: time_step
    real(c_double), target :: temperature(1)
    real(c_double), target :: pressure(1)
    real(c_double), target :: air_density(1)
    real(c_double), target :: concentrations(6)
    real(c_double), target :: user_defined_reaction_rates(2)

    type(string_t)       :: solver_state
    type(solver_stats_t) :: solver_stats
    type(error_t)        :: error

    type(micm_t), pointer :: micm

    integer :: i

    config_path = "configs/analytical"
    solver_type = RosenbrockStandardOrder
    num_grid_cells = 1

    time_step = 200
    temperature(1) = 273.0
    pressure(1) = 1.0e5
    air_density(:) = pressure(:) / (GAS_CONSTANT * temperature(:))

    concentrations = (/ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 /)
    user_defined_reaction_rates = (/ 0.001, 0.002 /)

    write(*,*) "Creating MICM solver..."
    micm => micm_t(config_path, solver_type, num_grid_cells, error)

    do i = 1, micm%species_ordering%size()
      print *, "Species Name:", micm%species_ordering%name(i), &
               ", Index:", micm%species_ordering%index(i)
    end do

    write(*,*) "Solving starts..."
    call micm%solve(time_step, temperature, pressure, air_density, concentrations, &
        user_defined_reaction_rates, solver_state, solver_stats, error)
    write(*,*) "After solving, concentrations", concentrations

    deallocate( micm )

  end subroutine box_model

end program test_micm_box_model

From the musica_util module we need the Fortran types error_t, string_t, and mapping_t. A pointer to a musica_micm::micm_t will serve as the interface to the MICM solver (in the example the pointer name is micm). Note that the config_path in the code sample has been set to configs/analytical, so that subdir should be created relative to the main program and contain the MICM JSON configuration files, or otherwise the config_path should be modified appropriately. The initial species concentrations are initialized in the concentrations array, which is an argument to the MICM solver.

Finally, a single time step solution is obtained through a call to micm%solve, after which the updated concentrations may be displayed.

$ ./test_micm_box_model
  Creating MICM solver...
  Species Name:A, Index:           1
  Species Name:B, Index:           2
  Species Name:C, Index:           3
  Solving starts...
  After solving, concentrations  0.38  1.61E-009  2.62
$