ESCOMP / ESCOMP/atmospheric_physics
Provide tropospheric surface area density calculations to musica chemistry
- Dominant language
- Fortran
- Stars
- 12
- Forks
- 38
- Avg merge
- 11h 16m
- Merged PRs (30d)
- 5
Description
CAM-SIMA, through the abstract aerosol interface, now supports some aerosol properties. Those aerosol properties can be used to drive some gas-phase aerosol reactions in musica. This issue details a plan to make those connections at runtime for tropospheric surface area density calculations.
- CAM-SIMA can provide the [surface area density](https://github.com/ESCOMP/CAM-SIMA/blob/3db176f5a4613e1606ab09f4a9d2fb46c45312b9/src/aerosol/aerosol_state_mod.F90#L350) and [effective radius](https://github.com/ESCOMP/CAM-SIMA/blob/3db176f5a4613e1606ab09f4a9d2fb46c45312b9/src/aerosol/aerosol_state_mod.F90#L351) for a particular list of aerosol species types
- It provides them for a list of species types, aggregated over modes, with the surface of each internally mixed mode partitioned by the mass fraction of the listed types ([modal_aerosol_state_mod.F90:848-1010](https://github.com/ESCOMP/CAM-SIMA/blob/3db176f5a4613e1606ab09f4a9d2fb46c45312b9/src/aerosol/modal_aerosol_state_mod.F90#L848-L1010)).
- MUSICA needs the number concentration and effective particle radius for a particular list of aerosol species types
- MUSICA needs these values to be set in our [rate parameters](https://github.com/NCAR/musica/blob/8316648b8a335352e770644031941a3a412c3b4a/fortran/micm/state.F90#L131) using the **reaction name** for our [surface reaction type](https://mechanismconfiguration.readthedocs.io/en/latest/v1/reactions/surface.html)
- the `N` we back out is **not** the mode's particle number concentration. It is the number of particles of radius `d_m/2` that would carry the mode's reactive surface, so it already includes the lognormal surface enhancement `exp(2 ln^2 sigma_g)` and the mass fraction of the mapped species types. Call it a "surface-equivalent number" so nobody confuses it with a physical number concentration
- Therefore we need
- a mapping from a list of CAM-SIMA aerosol species types to a MUSICA surface reaction
- Or, we map the CAM-SIMA aerosol species to species in MUSICA. From there, we can map the musica species to the particular reactions, which will then allow us to calculate and set the aerosol information we need. Either one works, as long as the aerosol side stays a list of species types and not a single species
**Acceptance criteria**
- Configured CAM-SIMA aerosol species can be mapped to configured MUSICA surface reactions
- Properly mapped species have their aerosol properties for surface area density and effective radius either used in the musica ccpp wrapper to set the number concentration and effective radius for each musica surface reaction at each chemistry timestep
- Tests are added which have valid gas-phase surface reactions that match CAM-SIMA aerosols for MAM
- CARMA is not yet implemented. BAM is not the target of this work. Focus on MAM
- **mechanism-parameter check.** Diff the mechanism configuration against the constants in CAM's `mo_usrrxt`: `reaction probability` per reaction equals CAM's `gamma_*` constant; `diffusion coeff` is `1e-5 m2/s` for every reacting gas (CAM's `dg = 0.1 cm2/s`), except effectively infinite for glyoxal, whose CAM rate has no diffusion term; `molecular weight` matches the value behind CAM's hard-coded mean-speed coefficient for that gas
- **summed-rate check.** Given a MAM4 state, the sum over modes of MICM's rates should be close enough to CAM's `hetrxtrate` to **1e-3 relative** ([mo_usrrxt.F90:3215-3231](https://github.com/ESCOMP/CAM/blob/06cf015bb555fdcd2d2fc34e53bebe977230db86/src/chemistry/mozart/mo_usrrxt.F90#L3215-L3231)). Exact agreement is not possible: CAM rounds its mean-speed coefficients (for example `1.40e3` for N2O5) to three digits, which limits agreement to about 1e-4 on its own
- Consider saving the current rates calculated in CAM to a file (describe **exactly** how they were obtained, which commit, etc.) and add a regression test against these saved values.
- Tests are added which catch missing/invalid/incorrect mappings of CAM-SIMA aerosols and musica mechanisms
**Ideas**
- Create one MICM `SURFACE` reaction per mode per gas. Create one configuration file for each mam configuration
- naming convention: `.mode` for mode `l` of `nbins` (for example `N2O5_sulfate.mode1` ... `.mode4`), all sharing the same gas-phase reactant/products, reaction probability, and species-type list — only `(N_l, r_l)` differs per mode
- we can use the [mechanism configuration file list format](https://mechanismconfiguration.readthedocs.io/en/latest/v1/overview.html#file-list-format) to reduce the amount of duplication. Place each surface reaction in a file for each mam configuration and share the same core parts of TS1
- no special handling needed on the MICM side to sum these back together, multiple reactions consuming the same gas-phase species already add in the solver, same as CAM's sum over modes does
- Set the gas-phase diffusion coefficient in the mechanism/rate-constant config to `1e-5 m2 s-1`
(matching CAM's constant `D_g = 0.1 cm2/s`, used for every species) for parity with CAM, unless a
species-specific, physically-based value is deliberately chosen instead
- Calculate the number concentration from the surface area density $SAD=4\pi r_{eff}^2 N$, where $N$ is the number concentration
- Do this per mode, not on the aggregate sad/reff outputs. CAM's mo_usrrxt/hetrxtrate never uses the aggregate sad/reff. Those are diagnostic-only in CAM (SAD_TROP, REFF_TROP history fields). The real rate sums per mode, using the optional per-mode `sfc(:,:,l)` (cm2/cm3) and `dm_aer(:,:,l)` (cm) outputs of surf_area_dens. Per mode: `r_l = dm_l / 2` and `N_l = sfc_l / (pi * dm_l^2)` reproduces CAM's per-mode term exactly.
- note surf_area_dens skips the primary_carbon mode entirely, so black-c only shows up through the accumulation mode
- make sure to get the units right. CAM-SIMA provides surface area density in [cm2 cm-3] and effective radius in [cm]. Musica needs the number concentration in [particles m-3], and the effective radius in [m]
- Create a mapping for CAM-SIMA aerosols and musica reactions
- consider implementing a similar mapping to what we do for tuvx
- we create a [mapping file stored alongside our mechanism and tuvx data](https://github.com/NCAR/cam-sima-chemistry-data/blob/main/mechanisms/chapman/tuvx_micm_mapping.json) and then [read that file to create a mapping from tuvx to musica](https://github.com/ESCOMP/atmospheric_physics/blob/cb0992cbb5fcd2a3b22418478ebef222e7208a4e/schemes/musica/tuvx/musica_ccpp_tuvx.F90#L482-L493)
- we could do something similar, either in the namelist for cam or just another json file since we already have APIs to read these, that map CAM-SIMA MAM aerosol species names to the corresponding surface reaction name in a particular mechanism they should be mapped to
- default the mapping to the two groupings CAM already uses today, for parity: `sad_chem_spec_types` (default sulfate, s-organic, p-organic, black-c, ammonium) for the general tropospheric heterogeneous reactions, and `sad_seasalt_spec_types` (default seasalt) for the sea-salt-specific ones (dehalogenation, ClNO2 yield). both are CAM `aerosol_nl` namelist variables today and are runtime configurable, we should keep that
- **each mapping entry should map to a *list* of CAM-SIMA aerosol species types, not a single species name.** `aero_state%surf_area_dens`'s `types_list` argument ([aerosol_state_mod.F90:340](https://github.com/ESCOMP/CAM-SIMA/blob/3db176f5a4613e1606ab09f4a9d2fb46c45312b9/src/aerosol/aerosol_state_mod.F90#L340)) is itself an array. A MAM mode is internally mixed, and [`modal_aerosol_state_mod.F90`](https://github.com/ESCOMP/CAM-SIMA/blob/development/src/aerosol/modal_aerosol_state_mod.F90#L848-L1010) computes surface area as the mass-fraction-weighted share attributable to whichever species types are listed. Some surface reactions may need more than one species type combined into one aggregate surface (for example sulfate and nitrate together), so the mapping schema should be `reaction_name -> [species_type, ...]` from the start, for example `{"N2O5_sulfate": ["sulfate"]}`, rather than `{"N2O5_sulfate": "sulfate"}`, which cannot express the combined case later without a schema change
- Add a new file inside the musica tree, for example `atmospheric_physics/schemes/musica/aerosol/musica_ccpp_surface_area_density.F90`.
- Use the aerosol interface to fetch the MAM `aero_state` / `aero_props`, similar to what [ndrop does](https://github.com/ESCOMP/atmospheric_physics/blob/cb0992cbb5fcd2a3b22418478ebef222e7208a4e/schemes/microp_aero/ccpp/ndrop_bam_ccpp.F90#L102-L135).
- Alternatively, select and cache the MAM aerosol-model index once, at init, the way [`nucleate_ice_ccpp_init` does](https://github.com/ESCOMP/atmospheric_physics/blob/13d634758a26127ac574972006b4b87656355322/schemes/microp_aero/ccpp/nucleate_ice_ccpp.F90#L113-L133): then loop `aerosol_instances_get_num_models()` once and keep the matching index as module state.
- At each chemistry timestep, use the cached index for a direct `aerosol_instances_get_props`/`aerosol_instances_get_state` lookup (no search), the way `nucleate_ice_ccpp_run` does at [lines 377-379](https://github.com/ESCOMP/atmospheric_physics/blob/13d634758a26127ac574972006b4b87656355322/schemes/microp_aero/ccpp/nucleate_ice_ccpp.F90#L377-L379)
- **call `surf_area_dens` once per species-type list, not once per label.** Each call loops every column, level, mode, and species, and allocates work arrays. CAM makes only two calls per timestep total (one per list). TS1 has about 40 uptake reactions times four modes, so 160 labels. Calling once per label costs about 80x more than grouping the labels by their shared species-type list and filling every label in a group from one call
- **make this a separate CCPP scheme, not an internal call from `musica_ccpp_run`.** `musica_ccpp` is now split into `musica_ccpp_photolysis` and `musica_ccpp_chemistry` ([atmospheric_physics#443](https://github.com/ESCOMP/atmospheric_physics/pull/443)). Place the new scheme between them, taking `micm_rate_parameters` as `inout`. The slot is already earmarked as `Slot for a SAD -> MICM rate-parameter provider:` in `test/test_suites/suite_modal_aero_wet_musica.xml`, right above `mam_surfarea_diagnostics`
- validate the mode count at init. For each base reaction, check that the number of `.mode` labels equals `aero_props%nbins()`, and abort if not. A silent mode-count mismatch would otherwise only show up as a wrong answer at run time
- do not special-case the `primary_carbon` mode out of the mechanism config. `surf_area_dens` returns zero surface for that mode, so its `.mode` reaction just receives `N = 0`. Generating it and feeding zero is correct, and simpler than excluding it
- mapping the sea-salt species type to the sea-salt surface is necessary, but it is not sufficient for the sea-salt-specific reactions. The ClNO2 yield depends on total sea-salt mass, with an HCl-limited branch. The dehalogenation rates carry land masks and a dehalogenation factor. Defer this extra physics, but say so explicitly, do not imply the mapping alone finishes these reactions
- A MAM sandbox for this work already exists: `jimmielin/CAM-SIMA` branch `hplin/modal_aero_sad_sandbox` (clone, then run `bin/git-fleximod update`). Its test suite `modal_aero_wet_musica` (`atmospheric_physics/test/test_suites/`) runs the full MAM wet-process block plus [`mam_surfarea_diagnostics`](https://github.com/ESCOMP/atmospheric_physics/tree/main/schemes/sima_diagnostics), which already calls `aero_state%surf_area_dens` the way CAM does, writes `SAD_TROP`, `REFF_TROP`, and per-mode `sad_a01..04`/`dmaer_a01..04` to the history tape, and runs the identity check from the acceptance criteria below. It has been verified to run to completion on Derecho with wet-process bit-for-bit checks intact
- Below is a list of relevant PRs/files that could be helpful for this work or provide more context
Relevant PRs
- https://github.com/ESCOMP/CAM/pull/1628
- https://github.com/ESCOMP/atmospheric_physics/pull/443 (splits `musica_ccpp` into `musica_ccpp_photolysis` and `musica_ccpp_chemistry`, the slot for the new SAD scheme)
Relevant code locations
In [cam-sima](https://github.com/ESCOMP/CAM-SIMA/tree/development)
- [cime_config/namelist_definition_cam.xml](https://github.com/ESCOMP/CAM-SIMA/blob/3db176f5a4613e1606ab09f4a9d2fb46c45312b9/cime_config/namelist_definition_cam.xml#L444-L488)
- this namelist has the definition of species that will be included for different aerosols
- this particular section seems to show that we are defining these BAM species
- `sulf`, `dust1`, `dust2`, `dust3`, `dust4`, `bcar1`, `bcar2`, `ocar1`, `ocar2`, `sslt1`, `sslt2`, `sslt3`, `sslt4`
- when a configuration like this is present, we will need to connect them to species in musica and get the properties we need for our [surface reactions](https://mechanismconfiguration.readthedocs.io/en/latest/v1/reactions/surface.html)
- this means we will need a way to connect names configured in musica configurations with names configured in the nameilst. We might do this similar to how we connect rates in tuvx with rates in micm
- [cime_config/namelist_definition_cam.xml#L430-L441](https://github.com/ESCOMP/CAM-SIMA/blob/3db176f5a4613e1606ab09f4a9d2fb46c45312b9/cime_config/namelist_definition_cam.xml#L430-L441)
- `mode_defs` is the MAM equivalent of the BAM species list above, but it defaults to empty (`''`) in cam-sima right now. there's no shipped MAM mode/species example in this repo yet, we'll need to write our own for testing (something like the mam4 `mode_defs`/`rad_climate` example from classic CAM would work as a starting point)
- [test/unit/fortran/src/aerosol/modal_test_helpers.F90](https://github.com/ESCOMP/CAM-SIMA/blob/development/test/unit/fortran/src/aerosol/modal_test_helpers.F90)
- a small hand-built 2-mode mock (accum: sulfate + black-c, coarse: dust) that our own modal aerosol unit tests already use. good starting point for testing this too
- [src/aerosol/radiative_aerosol_definitions.F90](https://github.com/ESCOMP/CAM-SIMA/blob/development/src/aerosol/radiative_aerosol_definitions.F90#L15-L34)
- This shows that there are only a couple of allowed species types in cam-sima.
- [src/aerosol/aerosol_state_mod.F90](https://github.com/ESCOMP/CAM-SIMA/blob/development/src/aerosol/aerosol_state_mod.F90)
- specifically the [surface area density and effective radius attributes](https://github.com/ESCOMP/CAM-SIMA/blob/3db176f5a4613e1606ab09f4a9d2fb46c45312b9/src/aerosol/aerosol_state_mod.F90#L350-L351)
- [src/aerosol/modal_aerosol_state_mod.F90](https://github.com/ESCOMP/CAM-SIMA/blob/development/src/aerosol/modal_aerosol_state_mod.F90)
- mam implementaiton
- [src/aerosol/aerosol_instances_mod.F90](https://github.com/ESCOMP/CAM-SIMA/blob/development/src/aerosol/aerosol_instances_mod.F90)
- allows us to get the aerosol state
- [src/aerosol/aerosol_spec_utils.F90](https://github.com/ESCOMP/CAM-SIMA/blob/development/src/aerosol/aerosol_spec_utils.F90)
- potentially useful for finding a specific species maybe?
In [atmospheric physics](https://github.com/ESCOMP/atmospheric_physics)
- [schemes/musica/musica_ccpp.F90](https://github.com/ESCOMP/atmospheric_physics/blob/main/schemes/musica/musica_ccpp.F90)
- [schemes/musica/micm/musica_ccpp_micm.F90](https://github.com/ESCOMP/atmospheric_physics/blob/main/schemes/musica/micm/musica_ccpp_micm.F90)
In [musica](https://github.com/NCAR/musica)
- [micm state]()
- [rate parameters](https://github.com/NCAR/musica/blob/8316648b8a335352e770644031941a3a412c3b4a/fortran/micm/state.F90#L131)
- [rate parameter ordering](https://github.com/NCAR/musica/blob/8316648b8a335352e770644031941a3a412c3b4a/fortran/micm/state.F90#L133)
- You can use these to set a particular rate paramter
- [getting an index](https://github.com/NCAR/musica/blob/8316648b8a335352e770644031941a3a412c3b4a/fortran/test/integration/test_micm_api.F90#L115-L120)
- [setting a rate for an index](https://github.com/NCAR/musica/blob/8316648b8a335352e770644031941a3a412c3b4a/fortran/test/integration/test_micm_api.F90#L133-L135)
- If you happen to make an array that contains every rate parameter correctly ordered, you can pass it directly to musica, `state%rate_parameters(:) = my_rate_values(:)`
- because we can solve multiple grid cells simultaneously, this would need to contain all the different rate parameters across all grid cells, corectly ordered
Contributor guide
No contributing guide indexed for this repository
Research direction
Start with test/test_suites/suite_modal_aero_wet_musica.xml and the existing mam_surfarea_diagnostics scheme, then inspect schemes/microp_aero/ccpp/nucleate_ice_ccpp.F90 and the MUSICA chemistry/photolysis slot. Done means MAM mappings use species-type lists, per-mode rates meet the CAM checks, invalid mode or mechanism mappings are rejected, and sea-salt-specific extra physics is explicitly deferred.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- fortran
- Domain
- backend, testing-qa
- Issue type
- Feature
- Difficulty
- 5/5
- Estimated time
- Over a week
- Activity status
- Active
- Clarity
- Mostly clear
- Newbie friendliness
- 48/100