diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index f1a55bc3a3..2257a1621f 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -1742,6 +1742,15 @@
+
+
+
+
+
+
#endif
@@ -2106,9 +2115,18 @@
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 7fda1f27a8..1d707bb3ed 100644
--- a/src/core_atmosphere/chemistry/mpas_atm_chemistry.F
+++ b/src/core_atmosphere/chemistry/mpas_atm_chemistry.F
@@ -86,16 +86,32 @@ 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_kind_types, only : RKIND
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..f67b80b227 100644
--- a/src/core_atmosphere/chemistry/musica/mpas_musica.F
+++ b/src/core_atmosphere/chemistry/musica/mpas_musica.F
@@ -44,27 +44,35 @@ 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_kind_types, only : RKIND
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
+
+ character(len=*), parameter :: INITIAL_CONCENTRATION_PROPERTY = "__initial concentration"
type(error_t) :: error
type(string_t) :: micm_version
! 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
@@ -87,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
@@ -99,25 +126,57 @@ 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...')
- ! 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.
+ 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
+ error_code = 1
+ error_message = '[MUSICA Error] MICM solver did not converge (state=' // &
+ solver_state%get_char_array() // ')'
+ return
+ end if
+
+ 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')
+ 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 +202,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 +338,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 d45a115e33..0ee2e73049 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -1173,6 +1173,15 @@
packages="mp_thompson_aers_in"/>
+
+
+
+
+
+