From a6f1da8a6f1ff1877793970b8a6ce4ec0acb19d8 Mon Sep 17 00:00:00 2001 From: David Fillmore Date: Mon, 23 Mar 2026 20:07:34 -0600 Subject: [PATCH 1/3] MUSICA-MICM chemistry solver integration 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) --- src/core_atmosphere/Registry.xml | 18 ++ .../chemistry/mpas_atm_chemistry.F | 21 ++- .../chemistry/musica/mpas_musica.F | 169 ++++++++++++++++-- src/core_atmosphere/mpas_atm_core.F | 2 +- src/core_init_atmosphere/Registry.xml | 9 + 5 files changed, 199 insertions(+), 20 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index c0bdcdd05d..88f61285f9 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1737,6 +1737,15 @@ description="Graupel mixing ratio" packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/> + + + + + + @@ -2103,6 +2112,15 @@ description="Tendency of graupel mass per unit volume divided by d(zeta)/dz" packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/> + + + + + + diff --git a/src/core_atmosphere/chemistry/mpas_atm_chemistry.F b/src/core_atmosphere/chemistry/mpas_atm_chemistry.F index 7fda1f27a8..d0d43910a1 100644 --- a/src/core_atmosphere/chemistry/mpas_atm_chemistry.F +++ b/src/core_atmosphere/chemistry/mpas_atm_chemistry.F @@ -19,6 +19,8 @@ !------------------------------------------------------------------------- module mpas_atm_chemistry + use mpas_kind_types, only : StrKIND, RKIND + implicit none private @@ -86,16 +88,31 @@ end subroutine chemistry_init !> This routine steps the chemistry packages, updating their state !> based on the current simulation time and conditions. !------------------------------------------------------------------------ - subroutine chemistry_step() + subroutine chemistry_step(dt) #ifdef MPAS_USE_MUSICA + use iso_fortran_env, only: real64 use mpas_musica, only: musica_step #endif use mpas_log, only : mpas_log_write + use mpas_derived_types, only: MPAS_LOG_CRIT + + real (kind=RKIND), intent(in) :: dt #ifdef MPAS_USE_MUSICA + real(real64) :: time_step + integer :: error_code + character(len=:), allocatable :: error_message + + time_step = dt + call mpas_log_write('Stepping chemistry packages...') - ! call musica_step() + + call musica_step(time_step, error_code, error_message) + + if (error_code /= 0) then + call mpas_log_write(error_message, messageType=MPAS_LOG_CRIT) + end if #endif end subroutine chemistry_step diff --git a/src/core_atmosphere/chemistry/musica/mpas_musica.F b/src/core_atmosphere/chemistry/musica/mpas_musica.F index a531a2414e..43156e6c38 100644 --- a/src/core_atmosphere/chemistry/musica/mpas_musica.F +++ b/src/core_atmosphere/chemistry/musica/mpas_musica.F @@ -44,21 +44,30 @@ module mpas_musica !> setting up necessary parameters and data structures for the simulation. !> For now, we will load fixed configurations for MICM and TUV-x. !> Later, we will gradually replace the fixed configuration elements - !> with runtime updates using the MUSICA API + !> with runtime updates using the MUSICA API. !------------------------------------------------------------------------ subroutine musica_init(filename_of_micm_configuration, & number_of_grid_cells, & error_code, error_message) + use iso_fortran_env, only: real64 + use musica_micm, only : RosenbrockStandardOrder use musica_util, only : error_t, string_t use mpas_log, only : mpas_log_write - character(len=*), intent(in) :: filename_of_micm_configuration - integer, intent(in) :: number_of_grid_cells - integer, intent(out) :: error_code - character(len=:), allocatable, intent(out) :: error_message + character(len=*), intent(in) :: filename_of_micm_configuration + integer, intent(in) :: number_of_grid_cells + integer, intent(out) :: error_code + character(len=:), allocatable, intent(out) :: error_message + + ! TEMPORARY: ABBA specific initial concentrations + real(real64), allocatable :: init_ab(:) + real(real64), allocatable :: init_a(:) + real(real64), allocatable :: init_b(:) + + character(len=*), parameter :: INITIAL_CONCENTRATION_PROPERTY = "__initial concentration" type(error_t) :: error type(string_t) :: micm_version @@ -99,25 +108,48 @@ end subroutine musica_init ! subroutine musica_step ! !> \brief Steps the MUSICA chemistry package - !> \author CheMPAS-A Developers - !> \date 20 August 2025 + !> \author David Fillmore + !> \date 17 November 2025 !> \details - !> This subroutine steps the MUSICA chemistry package, updating its state - !> based on the current simulation time and conditions. - !> It first calls the TUV-x package to compute photolysis rates, - !> then calls the MICM package to solve the ODEs for chemical species - !> concentrations. + !> This subroutine steps the MICM ODE solver. !------------------------------------------------------------------------ - subroutine musica_step() + subroutine musica_step(time_step, error_code, error_message) + + use iso_fortran_env, only: real64 + + use musica_micm, only : solver_stats_t + use musica_util, only : error_t, string_t use mpas_log, only : mpas_log_write + real(real64), intent(in) :: time_step + integer, intent(out) :: error_code + character(len=:), allocatable, intent(out) :: error_message + + type(string_t) :: solver_state + type(solver_stats_t) :: solver_stats + type(error_t) :: error + if (.not. musica_is_initialized) return - call mpas_log_write('Stepping MUSICA chemistry package...') + call mpas_log_write('[MUSICA] Stepping MICM solver...') + + call micm%solve(time_step, state, solver_state, solver_stats, error) + if (has_error_occurred(error, error_message, error_code)) return + + if (solver_state%get_char_array() /= 'Converged') then + call mpas_log_write('[MUSICA Warning] MICM solver failure') + end if + + call log_solver_stats(solver_stats) - ! Here we would typically call the TUV-x and MICM packages to perform - ! the necessary computations, but for now we will just log the step. + ! TEMPORARY: ABBA specific logging + call mpas_log_write('[MUSICA] MICM AB concentrations') + call log_concentrations(state, state%number_of_grid_cells, 'AB', error) + call mpas_log_write('[MUSICA] MICM A concentrations') + call log_concentrations(state, state%number_of_grid_cells, 'A', error) + call mpas_log_write('[MUSICA] MICM B concentrations') + call log_concentrations(state, state%number_of_grid_cells, 'B', error) end subroutine musica_step @@ -143,6 +175,110 @@ subroutine musica_finalize() end subroutine musica_finalize + subroutine assign_rate_parameters(state, number_of_grid_cells) + ! Provide a default value for any rate parameters (none for ABBA by default). + use iso_fortran_env, only : real64 + + type(state_t), intent(inout) :: state + integer, intent(in) :: number_of_grid_cells + + integer :: cell_stride, var_stride, rp_index, cell, idx + + if (state%number_of_rate_parameters == 0) return + + cell_stride = state%rate_parameters_strides%grid_cell + var_stride = state%rate_parameters_strides%variable + + do rp_index = 1, state%rate_parameters_ordering%size() + do cell = 1, number_of_grid_cells + idx = 1 + (cell - 1) * cell_stride + (rp_index - 1) * var_stride + state%rate_parameters(idx) = 1.0_real64 + end do + end do + end subroutine assign_rate_parameters + + subroutine set_species_profile(state, species_name, values, err) + ! Write a 1-D profile into the contiguous concentrations array. + use iso_fortran_env, only : real64 + use musica_util, only : error_t, string_t + + type(state_t), intent(inout) :: state + character(len=*), intent(in) :: species_name + real(real64), intent(in) :: values(:) + type(error_t), intent(inout) :: err + integer :: species_index, cell_stride, var_stride, cell, idx + + species_index = state%species_ordering%index(species_name, err) + if (.not. err%is_success()) return + + cell_stride = state%species_strides%grid_cell + var_stride = state%species_strides%variable + + do cell = 1, size(values) + idx = 1 + (cell - 1) * cell_stride + (species_index - 1) * var_stride + state%concentrations(idx) = values(cell) + end do + end subroutine set_species_profile + + subroutine log_concentrations(state, number_of_grid_cells, species_name, err) + use iso_fortran_env, only : real64 + use mpas_log, only : mpas_log_write + use mpas_kind_types, only : RKIND + use musica_util, only : error_t, string_t + + type(state_t), intent(in) :: state + integer, intent(in) :: number_of_grid_cells + character(len=*), intent(in) :: species_name + type(error_t), intent(inout) :: err + + integer :: species_index, cell_stride, var_stride, cell, idx + real (kind=RKIND) :: concentration + + species_index = state%species_ordering%index(species_name, err) + if (.not. err%is_success()) return + + cell_stride = state%species_strides%grid_cell + var_stride = state%species_strides%variable + + do cell = 1, number_of_grid_cells + idx = 1 + (cell - 1) * cell_stride + (species_index - 1) * var_stride + concentration = state%concentrations(idx) + call mpas_log_write(' MICM concentration: $r', realArgs=[concentration]) + end do + end subroutine log_concentrations + + subroutine log_solver_stats(solver_stats) + use mpas_log, only : mpas_log_write + use mpas_kind_types, only : RKIND + use musica_micm, only : solver_stats_t + + type(solver_stats_t), intent(in) :: solver_stats + + integer :: function_calls, jacobian_updates, number_of_steps + integer :: accepted, rejected, decompositions, solves + + real (kind=RKIND) :: final_time + + function_calls = solver_stats%function_calls() + jacobian_updates = solver_stats%jacobian_updates() + number_of_steps = solver_stats%number_of_steps() + accepted = solver_stats%accepted() + rejected = solver_stats%rejected() + decompositions = solver_stats%decompositions() + solves = solver_stats%solves() + final_time = solver_stats%final_time() + + call mpas_log_write('[MUSICA] MICM Solver statistics ...') + call mpas_log_write(' MICM function calls: $i', intArgs=[function_calls]) + call mpas_log_write(' MICM jacobian updates: $i', intArgs=[jacobian_updates]) + call mpas_log_write(' MICM number of steps: $i', intArgs=[number_of_steps]) + call mpas_log_write(' MICM accepted: $i', intArgs=[accepted]) + call mpas_log_write(' MICM rejected: $i', intArgs=[rejected]) + call mpas_log_write(' MICM LU decompositions: $i', intArgs=[decompositions]) + call mpas_log_write(' MICM linear solves: $i', intArgs=[solves]) + call mpas_log_write(' MICM final time (s): $r', realArgs=[final_time]) + end subroutine log_solver_stats + !------------------------------------------------------------------------- ! function has_error_occurred ! @@ -175,7 +311,6 @@ logical function has_error_occurred(error, error_message, error_code) error_message = '[MUSICA Error]: ' // error%category( ) // '[' // & trim( adjustl( error_code_str ) ) // ']: ' // error%message( ) has_error_occurred = .true. - end function has_error_occurred end module mpas_musica diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index eef1951e36..effa48beeb 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -1038,7 +1038,7 @@ subroutine atm_do_timestep(domain, dt, itimestep) #endif #ifdef DO_CHEMISTRY - call chemistry_step() + call chemistry_step(dt) #endif call atm_timestep(domain, dt, currTime, itimestep, exchange_halo_group) diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index fc9be9fb65..9c84da47cd 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -1163,6 +1163,15 @@ + + + + + + From 7a15848d78d360ea1b990344617633d6801f8ec0 Mon Sep 17 00:00:00 2001 From: David Fillmore Date: Sat, 18 Apr 2026 11:56:25 -0600 Subject: [PATCH 2/3] Address PR #1425 review feedback - 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) --- src/core_atmosphere/Registry.xml | 38 +++++++++---------- .../chemistry/mpas_atm_chemistry.F | 3 +- .../chemistry/musica/mpas_musica.F | 30 ++++++++++++--- src/core_init_atmosphere/Registry.xml | 18 ++++----- 4 files changed, 53 insertions(+), 36 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 88f61285f9..0c356eaf58 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1737,15 +1737,6 @@ description="Graupel mixing ratio" packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/> - - - - - - @@ -1768,6 +1759,15 @@ + + + + + + #endif @@ -2112,15 +2112,6 @@ description="Tendency of graupel mass per unit volume divided by d(zeta)/dz" packages="mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/> - - - - - - @@ -2141,8 +2132,17 @@ description="Tendency of water-friendly aerosol number concentration multiplied by dry air density divided by d(zeta)/dz" packages="mp_thompson_aers_in"/> - + + + + + + #endif diff --git a/src/core_atmosphere/chemistry/mpas_atm_chemistry.F b/src/core_atmosphere/chemistry/mpas_atm_chemistry.F index d0d43910a1..1d707bb3ed 100644 --- a/src/core_atmosphere/chemistry/mpas_atm_chemistry.F +++ b/src/core_atmosphere/chemistry/mpas_atm_chemistry.F @@ -19,8 +19,6 @@ !------------------------------------------------------------------------- module mpas_atm_chemistry - use mpas_kind_types, only : StrKIND, RKIND - implicit none private @@ -94,6 +92,7 @@ subroutine chemistry_step(dt) use iso_fortran_env, only: real64 use mpas_musica, only: musica_step #endif + use mpas_kind_types, only : RKIND use mpas_log, only : mpas_log_write use mpas_derived_types, only: MPAS_LOG_CRIT diff --git a/src/core_atmosphere/chemistry/musica/mpas_musica.F b/src/core_atmosphere/chemistry/musica/mpas_musica.F index 43156e6c38..e8c35efcd0 100644 --- a/src/core_atmosphere/chemistry/musica/mpas_musica.F +++ b/src/core_atmosphere/chemistry/musica/mpas_musica.F @@ -55,6 +55,7 @@ subroutine musica_init(filename_of_micm_configuration, & use musica_micm, only : RosenbrockStandardOrder use musica_util, only : error_t, string_t + use mpas_kind_types, only : RKIND use mpas_log, only : mpas_log_write character(len=*), intent(in) :: filename_of_micm_configuration @@ -62,11 +63,6 @@ subroutine musica_init(filename_of_micm_configuration, & integer, intent(out) :: error_code character(len=:), allocatable, intent(out) :: error_message - ! TEMPORARY: ABBA specific initial concentrations - real(real64), allocatable :: init_ab(:) - real(real64), allocatable :: init_a(:) - real(real64), allocatable :: init_b(:) - character(len=*), parameter :: INITIAL_CONCENTRATION_PROPERTY = "__initial concentration" type(error_t) :: error @@ -74,6 +70,9 @@ subroutine musica_init(filename_of_micm_configuration, & ! TEMPORARY: Hard-coded options for the MICM solver integer :: solver_type = RosenbrockStandardOrder integer :: i_species + character(len=:), allocatable :: species_name + real(real64) :: initial_concentration + real(real64), allocatable :: initial_profile(:) ! Skip MUSICA initialization if no configuration file is provided if (len_trim(filename_of_micm_configuration) == 0) then @@ -96,10 +95,29 @@ subroutine musica_init(filename_of_micm_configuration, & state => micm%get_state(number_of_grid_cells, error) if (has_error_occurred(error, error_message, error_code)) return + allocate(initial_profile(number_of_grid_cells)) + do i_species = 1, state%species_ordering%size() - call mpas_log_write('MICM species: ' // state%species_ordering%name(i_species)) + species_name = state%species_ordering%name(i_species) + call mpas_log_write('MICM species: ' // species_name) + + initial_concentration = micm%get_species_property_double(species_name, & + INITIAL_CONCENTRATION_PROPERTY, & + error) + if (has_error_occurred(error, error_message, error_code)) return + + initial_profile(:) = initial_concentration + call set_species_profile(state, species_name, initial_profile, error) + if (has_error_occurred(error, error_message, error_code)) return + + call mpas_log_write(' initial concentration: $r', & + realArgs=[real(initial_concentration, kind=RKIND)]) end do + deallocate(initial_profile) + + call assign_rate_parameters(state, number_of_grid_cells) + musica_is_initialized = .true. end subroutine musica_init diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index 9c84da47cd..b0e3c981bc 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -1163,15 +1163,6 @@ - - - - - - @@ -1181,6 +1172,15 @@ packages="mp_thompson_aers_in"/> + + + + + + From 752f89599459038c086c0a842487714aaf2a489d Mon Sep 17 00:00:00 2001 From: David Fillmore Date: Sat, 18 Apr 2026 12:02:52 -0600 Subject: [PATCH 3/3] Propagate MICM solver failures as errors 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) --- src/core_atmosphere/chemistry/musica/mpas_musica.F | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/core_atmosphere/chemistry/musica/mpas_musica.F b/src/core_atmosphere/chemistry/musica/mpas_musica.F index e8c35efcd0..f67b80b227 100644 --- a/src/core_atmosphere/chemistry/musica/mpas_musica.F +++ b/src/core_atmosphere/chemistry/musica/mpas_musica.F @@ -155,11 +155,20 @@ subroutine musica_step(time_step, error_code, error_message) call micm%solve(time_step, state, solver_state, solver_stats, error) if (has_error_occurred(error, error_message, error_code)) return + call log_solver_stats(solver_stats) + if (solver_state%get_char_array() /= 'Converged') then - call mpas_log_write('[MUSICA Warning] MICM solver failure') + error_code = 1 + error_message = '[MUSICA Error] MICM solver did not converge (state=' // & + solver_state%get_char_array() // ')' + return end if - call log_solver_stats(solver_stats) + if (solver_stats%final_time() < time_step) then + error_code = 1 + error_message = '[MUSICA Error] MICM solver returned before reaching the requested time step' + return + end if ! TEMPORARY: ABBA specific logging call mpas_log_write('[MUSICA] MICM AB concentrations')