MUSICA-MICM chemistry solver by dwfncar · Pull Request #1425 · MPAS-Dev/MPAS-Model

dwfncar · 2026-03-24T02:14:42Z

Added the MUSICA-MICM ODE chemistry solver.
MICM initialization, time-stepping, and finalization are called through
the existing chemistry hooks in mpas_atm_chemistry.F.
ABBA mechanism is hard coded for testing.

ABBA MICM configuration

The ABBA mechanism used for testing is two USER_DEFINED reactions (a
forward dissociation AB → A + B and a reverse recombination
A + B → AB) with prescribed scaling factors. abba.yaml:

version: "1.0.0"
name: ABBA
species:
  - name: A
    __absolute tolerance: 1e-4
    __long name: atom_A
    __do advect: true
    __initial concentration: 0.0
    __molar mass: 0.029
  - name: B
    __absolute tolerance: 1e-4
    __long name: atom_B
    __do advect: true
    __initial concentration: 0.0
    __molar mass: 0.029
  - name: AB
    __absolute tolerance: 1e-4
    __long name: molecule_AB
    __do advect: true
    __initial concentration: 1.0
    __molar mass: 0.058

phases:
  - name: gas
    species:
      - A
      - B
      - AB

reactions:
  - type: USER_DEFINED
    gas phase: gas
    reactants:
      - species name: AB
        coefficient: 1
    products:
      - species name: A
        coefficient: 1
      - species name: B
        coefficient: 1
    scaling factor: 2.0e-3
    name: forward_AB_to_A_B

  - type: USER_DEFINED
    gas phase: gas
    reactants:
      - species name: A
        coefficient: 1
      - species name: B
        coefficient: 1
    products:
      - species name: AB
        coefficient: 1
    scaling factor: 1.0e-3
    name: reverse_A_B_to_AB

Smoke test

Built on Ubuntu with the conda-forge mpas env (gfortran / openmpi / netcdf / pnetcdf) and MUSICA=true, linked against MUSICA-Fortran 0.14.5 / MICM 3.11.0. Ran the supercell idealized case with config_micm_file='abba.yaml', 30 s simulated time, 8 MPI ranks. Results:

Species	t=0 (YAML init)	t=3 s	t=6 s
AB	1.0	0.994018	0.988072
A	0.0	0.005982	0.011921
B	0.0	0.005982	0.011921

Atom conservation [X] + [AB] = 1.0 holds exactly for X ∈ {A, B}. Solver: 9 steps, 9 accepted / 0 rejected, final_time = 3.0 s matches config_dt. Zero model errors.

Integrate the MUSICA-MICM ODE chemistry solver with the MPAS atmosphere core. MICM initialization, time-stepping, and finalization are called through the existing chemistry hooks in mpas_atm_chemistry.F. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

mgduda · 2026-04-08T21:05:24Z

+
+                        <var name="qB" array_group="passive" units="kg kg^{-1}"
+                             description="Atomic B mixing ratio"/>
+


These variables are declared here in the Registry.xml file for the init_atmosphere core, but it doesn't look like there's any corresponding code changes to actually provide values for them. Is there some offline script that's used to initialize the qAB, qA, and qB fields? Or are they just expected to be zero in the initial conditions?

I use an offline python script to initialize the test tracers, e.g. to a sinusoidal pattern.

For symmetry with core_atmosphere/Registry.xml, these have also been moved to the end of the scalars var_array with array_group="chem_test" in 7a15848.

mgduda · 2026-04-08T21:07:17Z

+
+                        <var name="qB" array_group="passive" units="kg kg^{-1}"
+                             description="Atomic B mixing ratio"/>
+


In principle, the order shouldn't matter too much, but it might be nicer to add these at the end of the scalars var_array. We could also give them a more descriptive group name, perhaps something like chem_test?

Moved to the end of scalars with array_group="chem_test" in 7a15848.

mgduda · 2026-04-08T21:08:01Z

+
+                        <var name="tend_qB" name_in_code="qB" array_group="passive" units="kg m^{-3} s^{-1}"
+                             description="Tendency of atomic B mass per unit volume divided by d(zeta)/dz"/>
+


Similar to the comment for the scalars var_array, perhaps we can add these at the end of the var_array using the chem_test array group?

Moved to the end of scalars_tend with array_group="chem_test" in 7a15848.

mgduda · 2026-04-08T21:13:15Z

 !-------------------------------------------------------------------------
 module mpas_atm_chemistry

+    use mpas_kind_types, only : StrKIND, RKIND


It doesn't look like StrKIND is only used in chemistry_init, where there's already an use statement to import StrKIND.
Also, it looks like RKIND is only used in the chemistry_step routine. Would it be worth moving the use statement to import RKIND to chemistry_step for consistency with the import of real64?

Good catch. The module-level use mpas_kind_types was redundant (chemistry_init already imports StrKIND locally, and RKIND is only needed in chemistry_step). In 7a15848 the module-level import is removed and use mpas_kind_types, only : RKIND is moved into chemistry_step.

mgduda · 2026-04-08T21:16:43Z

+        ! TEMPORARY: ABBA specific initial concentrations
+        real(real64), allocatable :: init_ab(:)
+        real(real64), allocatable :: init_a(:)
+        real(real64), allocatable :: init_b(:)


It doesn't look like these variables are actually used. Is there otherwise a reason for declaring them here?

These were placeholders for initial-concentration seeding that had not been wired up yet. In 7a15848 they are replaced with a single reused initial_profile(:) populated per species in a loop over state%species_ordering, using micm%get_species_property_double() with the INITIAL_CONCENTRATION_PROPERTY constant. MICM species concentrations are now initialized from the __initial concentration field in the MICM YAML, so the hard-coded ABBA scaffolding is gone.

mgduda · 2026-04-08T21:17:20Z

+        real(real64), allocatable :: init_a(:)
+        real(real64), allocatable :: init_b(:)
+
+        character(len=*), parameter :: INITIAL_CONCENTRATION_PROPERTY = "__initial concentration"


It doesn't look like this parameter is used anywhere.

Now used in 7a15848 as the property key passed to micm%get_species_property_double() to read initial concentrations from the MICM YAML. See the reply on the init_ab/init_a/init_b thread for the full pattern.

mgduda · 2026-04-08T23:20:19Z

@dwfncar Is there a MICM configuration file that's needed in order to run with the ABBA mechanism? If so (and assuming it's short), could you include the contents of that configuration file in the PR description?

- Registry (core_atmosphere, core_init_atmosphere): move qAB/qA/qB and tend_qAB/tend_qA/tend_qB to the end of their scalars/scalars_tend var_arrays with array_group="chem_test" (per mgduda review comments). - mpas_atm_chemistry.F: move RKIND import from module level into chemistry_step, where it is the only consumer. - mpas_musica.F: wire up INITIAL_CONCENTRATION_PROPERTY through micm%get_species_property_double() so initial species concentrations from the MICM YAML (__initial concentration) are loaded into state%concentrations at init. Replace unused init_ab/init_a/init_b placeholders with a single reused initial_profile(:). Also call assign_rate_parameters after the species loop so USER_DEFINED reactions in the ABBA mechanism have a defined rate. Smoke test on the supercell case now shows physically correct AB -> A + B decay with atom conservation. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

musica_step now treats non-convergence and early termination as fatal errors that set error_code and return with a descriptive error_message, instead of just logging a warning. The caller (mpas_atm_chemistry:chemistry_step) already routes non-zero error_code through MPAS_LOG_CRIT, so the model now aborts on a silent MICM failure rather than advancing with stale state. Two checks added: - solver_state must be 'Converged' (previously only log-warned) - solver_stats%final_time() must reach the requested time_step (catches partial/interrupted solves that would otherwise be invisible) log_solver_stats is now called before the checks so the stats leading to a failure are always in the log. Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

dwfncar · 2026-04-18T18:09:28Z

Done — the full abba.yaml is now in the PR description.

dwfncar · 2026-04-18T18:09:39Z

@mgduda Thanks for the review. All six review threads addressed in 7a15848; individual replies left on each.

One extra, unsolicited change in a follow-up commit (752f895) that I want to flag because it changes semantics: musica_step now propagates non-convergence and early termination as errors (error_code = 1 + descriptive error_message) instead of just logging a warning. Since chemistry_step already routes non-zero error_code through MPAS_LOG_CRIT, the model will now abort on silent MICM solver failures rather than advancing with stale state. Happy to split that into a separate PR if you would prefer to review it in isolation.

dwfncar marked this pull request as draft March 24, 2026 02:14

K20shores mentioned this pull request Mar 24, 2026

MICM Chemistry Solver Integration NCAR/MPAS-Model-ACOM-dev#8

Open

dwfncar marked this pull request as ready for review March 30, 2026 17:48

mgduda added Atmosphere Chemistry Changes specific to chemistry in MPAS labels Apr 7, 2026

mgduda self-requested a review April 7, 2026 22:11

mgduda requested changes Apr 8, 2026

View reviewed changes

dwfncar and others added 2 commits April 18, 2026 11:56

dwfncar requested a review from mgduda April 18, 2026 18:13


		<var name="qB" array_group="passive" units="kg kg^{-1}"
		description="Atomic B mixing ratio"/>


		<var name="tend_qB" name_in_code="qB" array_group="passive" units="kg m^{-3} s^{-1}"
		description="Tendency of atomic B mass per unit volume divided by d(zeta)/dz"/>

Conversation

dwfncar commented Mar 24, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

ABBA MICM configuration

Smoke test

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mgduda commented Apr 8, 2026

Uh oh!

dwfncar commented Apr 18, 2026

Uh oh!

dwfncar commented Apr 18, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

dwfncar commented Mar 24, 2026 •

edited

Loading