Skip to content

MPAS MUSICA MICM State

Jump to bottom
David Fillmore edited this page Jan 6, 2026 · 1 revision

Plan: Bidirectional MICM-MPAS State Coupling

Goal

Implement chemistry_from_MPAS and chemistry_to_MPAS routines to couple MICM chemistry state with MPAS model state, following the physics interface pattern. Focus on ABBA test species (qAB, qA, qB) with abstraction layer for future extension.

Scope: ABBA Model Only

This implementation focuses exclusively on the ABBA test chemistry model (species: A, B, AB). Full tracer mappings for production chemistry mechanisms will come later, pending the MPAS runtime variables feature (see below).

Future: MPAS Runtime Variables

A new MPAS feature, runtime variables, is under development but not yet ready. This feature will allow MPAS tracers to be defined at runtime (via configuration) rather than at compile-time (via Registry.xml).

Once available:

  • MICM will use runtime variables to dynamically register its species as MPAS tracers
  • No consistency check between YAML and Registry.xml will be needed
  • Full tracer mappings for arbitrary chemistry mechanisms will be supported

Current placeholder: The existing check_registry_tracer_consistency() routine validates that MICM species match MPAS tracers at initialization. This is a temporary solution until runtime variables are available.

Species Naming Convention (Current)

  • MICM species (runtime, YAML): A, B, AB
  • MPAS tracers (compile-time, Registry.xml): qA, qB, qAB (prepend "q")
  • check_registry_tracer_consistency() validates this mapping at init (placeholder)

Architecture

mpas_atm_chemistry.F (Chemistry-Agnostic Layer)
├── chemistry_from_MPAS()  - Extract MPAS state, call MICM-specific routine
├── chemistry_to_MPAS()    - Update MPAS state from MICM results
└── chemistry_step()       - Modified to orchestrate coupling

mpas_musica.F (MICM-Specific Layer)
├── MICM_from_chemistry()  - MPAS scalars → MICM concentrations (with unit conversion)
└── MICM_to_chemistry()    - MICM concentrations → MPAS scalars (with unit conversion)

Data Flow (Each Timestep)

chemistry_step(dt, mesh, state, diag, dimensions, time_lev)
  │
  ├─> chemistry_from_MPAS()
  │     └─> MICM_from_chemistry()   # q [kg/kg] → C [mol/m³]
  │
  ├─> musica_step()                 # MICM solves chemistry
  │
  └─> chemistry_to_MPAS()
        └─> MICM_to_chemistry()     # C [mol/m³] → q [kg/kg]

Unit Conversions

  • MPAS → MICM: C = q × rho_air / M_species
  • MICM → MPAS: q = C × M_species / rho_air

Where:

  • rho_air = zz × rho_zz × (1 + qv) [kg/m³]
  • M_species = molar mass [kg/mol] (e.g., M_AB=0.058, M_A=M_B=0.029)

Files to Modify

1. src/core_atmosphere/chemistry/musica/mpas_musica.F

Add new subroutines:

subroutine MICM_from_chemistry(scalars, rho_air, temperature, pressure, &
                                nCells, nVertLevels, &
                                index_qAB, index_qA, index_qB, &
                                error_code, error_message)
    ! Convert MPAS scalars to MICM state
    ! - Set state%conditions (T, P, air_density)
    ! - Convert mixing ratios to concentrations using q->C formula
    ! - Store in state%concentrations using strided access
    ! - Uses existing species mapping validated by check_registry_tracer_consistency
subroutine MICM_to_chemistry(scalars, rho_air, nCells, nVertLevels, &
                              index_qAB, index_qA, index_qB, &
                              error_code, error_message)
    ! Convert MICM state to MPAS scalars
    ! - Read state%concentrations using strided access
    ! - Convert concentrations to mixing ratios using C->q formula
    ! - Update scalars array (qA, qB, qAB)

Add module constants:

real(real64), parameter :: M_AB = 0.058_real64  ! kg/mol
real(real64), parameter :: M_A  = 0.029_real64  ! kg/mol
real(real64), parameter :: M_B  = 0.029_real64  ! kg/mol
real(real64), parameter :: M_AIR = 0.0289644_real64  ! kg/mol

2. src/core_atmosphere/chemistry/mpas_atm_chemistry.F

Modify chemistry_step signature:

subroutine chemistry_step(dt, mesh, state, diag, dimensions, time_lev)

Add new subroutines:

subroutine chemistry_from_MPAS(mesh, state, diag, dimensions, time_lev, &
                                nCells, nVertLevels)
    ! 1. Get scalars, rho_zz, zz, theta_m, exner from pools
    ! 2. Calculate rho_air and temperature
    ! 3. Call MICM_from_chemistry()
subroutine chemistry_to_MPAS(mesh, state, diag, dimensions, time_lev, &
                              nCells, nVertLevels)
    ! 1. Get scalars and rho_air
    ! 2. Call MICM_to_chemistry()

Modify chemistry_step implementation:

call chemistry_from_MPAS(...)     ! Before solver - calls MICM_from_chemistry
call musica_step(dt, ...)         ! Solve
call chemistry_to_MPAS(...)       ! After solver - calls MICM_to_chemistry

3. src/core_atmosphere/mpas_atm_core.F

Modify call to chemistry_step (around line 1045-1058):

block_chem => domain % blocklist
do while (associated(block_chem))
   call mpas_pool_get_subpool(block_chem % structs, 'mesh', mesh_chem)
   call mpas_pool_get_subpool(block_chem % structs, 'state', state_chem)
   call mpas_pool_get_subpool(block_chem % structs, 'diag', diag_chem)

   time_lev_chem = 1
   call chemistry_step(dt, mesh_chem, state_chem, diag_chem, &
                       block_chem % dimensions, time_lev_chem)

   block_chem => block_chem % next
end do

Key Implementation Patterns

Strided MICM Array Access

cell_stride = state%species_strides%grid_cell
var_stride  = state%species_strides%variable
spec_idx = state%species_ordering%index("AB", err)

do k = 1, nVertLevels
    idx = 1 + (k-1)*cell_stride + (spec_idx-1)*var_stride
    state%concentrations(idx) = conc_value
end do

MPAS Pool Access (following physics pattern)

call mpas_pool_get_array(state, 'scalars', scalars, time_lev)
call mpas_pool_get_array(state, 'rho_zz', rho_zz, time_lev)
call mpas_pool_get_array(mesh, 'zz', zz)
call mpas_pool_get_array(diag, 'exner', exner)
call mpas_pool_get_dimension(state, 'index_qAB', index_qAB)

Air Density and Temperature

! Moist air density [kg/m³]
rho_air(k, iCell) = zz(k, iCell) * rho_zz(k, iCell) * (1.0 + qv(k, iCell))

! Temperature [K]
theta = theta_m(k, iCell) / (1.0 + rvord * qv(k, iCell))
T = theta * exner(k, iCell)

Processing Strategy

Column-by-column approach:

  • Loop over cells
  • For each cell, copy column data to/from MICM state
  • MICM state sized for nVertLevels (single column)

This matches current MICM initialization and is simpler to implement correctly.

Implementation Order

  1. Add MICM_from_chemistry and MICM_to_chemistry to mpas_musica.F
  2. Add chemistry_from_MPAS and chemistry_to_MPAS to mpas_atm_chemistry.F
  3. Modify chemistry_step signature and implementation
  4. Update call site in mpas_atm_core.F
  5. Build and test with supercell case

Testing

  1. Build: make-mpas
  2. Run: $HOME/bin/run-mpas
  3. Verify in log: ABBA concentrations evolve over timesteps
  4. Inspect output: ncks -v qAB,qA,qB output.nc

Clone this wiki locally