Skip to content

Commit

Permalink
Add Chapman and Terminator configurations for chemistry (#151)
Browse files Browse the repository at this point in the history 
Adds configurations for Terminator and Chapman chemistry and basic
integration tests.

This includes some data files needed for TUV-x to calculate photolysis
rate constants for these mechanisms, but I wasn't sure if this is the
best place for them. I'm happy to remove the configuration data from
this PR if there is a better place to keep them.

closes #150

Also:
- Adds calculation of photolysis rate constants from TUV-x
- Maps TUV-x rate constants to MICM rate parameters
- Removes some array copying using updated MUSICA functions that accept
pointers to Fortran arrays
- Sets up a gitignore file for the repo
  • Loading branch information
@mattldawson
mattldawson authored Nov 14, 2024
1 parent 026c4c5 commit 0ddf4d2
Show file tree
Hide file tree
Showing 15 changed files with 504 additions and 459 deletions.
7 changes: 7 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
**/*.DS_Store
**/*.pyc
**/build/
**/CMakeCache.txt
**/CMakeFiles/
.vscode
xcode/
47 changes: 23 additions & 24 deletions schemes/musica/micm/musica_ccpp_micm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,17 +2,18 @@ module musica_ccpp_micm

! Note: "micm_t" is included in an external pre-built MICM library that the host
! model is responsible for linking to during compilation
use musica_micm, only: micm_t
use ccpp_kinds, only: kind_phys
use musica_ccpp_util, only: has_error_occurred
use musica_ccpp_namelist, only: filename_of_micm_configuration
use ccpp_kinds, only: kind_phys
use musica_micm, only: micm_t

implicit none
private

public :: micm_register, micm_init, micm_run, micm_final
public :: micm_register, micm_init, micm_run, micm_final, micm, number_of_rate_parameters

type(micm_t), pointer :: micm => null( )
integer :: number_of_rate_parameters = 0

contains

Expand All @@ -36,7 +37,7 @@ subroutine micm_register(solver_type, num_grid_cells, constituent_props, errmsg,
logical :: is_advected
integer :: i, species_index

micm => micm_t(filename_of_micm_configuration, solver_type, num_grid_cells, error)
micm => micm_t(trim(filename_of_micm_configuration), solver_type, num_grid_cells, error)
if (has_error_occurred(error, errmsg, errcode)) return

allocate(constituent_props(micm%species_ordering%size()), stat=errcode)
Expand Down Expand Up @@ -73,6 +74,7 @@ subroutine micm_register(solver_type, num_grid_cells, constituent_props, errmsg,
if (errcode /= 0) return
end associate ! map
end do
number_of_rate_parameters = micm%user_defined_reaction_rates%size()

end subroutine micm_register

Expand All @@ -91,34 +93,31 @@ subroutine micm_run(time_step, temperature, pressure, dry_air_density, &
user_defined_rate_parameters, constituents, errmsg, errcode)
use musica_micm, only: solver_stats_t
use musica_util, only: string_t, error_t
use iso_c_binding, only: c_double
use iso_c_binding, only: c_double, c_loc

real(kind_phys), intent(in) :: time_step ! s
real(c_double), intent(in) :: temperature(:) ! K
real(c_double), intent(in) :: pressure(:) ! Pa
real(c_double), intent(in) :: dry_air_density(:) ! kg m-3
real(c_double), intent(in) :: user_defined_rate_parameters(:) ! various units
real(c_double), intent(inout) :: constituents(:) ! mol m-3
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode
real(kind_phys), intent(in) :: time_step ! s
real(kind_phys), target, intent(in) :: temperature(:,:) ! K
real(kind_phys), target, intent(in) :: pressure(:,:) ! Pa
real(kind_phys), target, intent(in) :: dry_air_density(:,:) ! kg m-3
real(kind_phys), target, intent(in) :: user_defined_rate_parameters(:,:,:) ! various units
real(kind_phys), target, intent(inout) :: constituents(:,:,:) ! mol m-3
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode

! local variables
type(string_t) :: solver_state
type(solver_stats_t) :: solver_stats
type(error_t) :: error
real(c_double) :: c_time_step
integer :: i_elem

c_time_step = real(time_step, c_double)

call micm%solve(c_time_step, &
temperature, &
pressure, &
dry_air_density, &
constituents, &
user_defined_rate_parameters, &
solver_state, &
solver_stats, &
call micm%solve(real(time_step, kind=c_double), &
c_loc(temperature), &
c_loc(pressure), &
c_loc(dry_air_density), &
c_loc(constituents), &
c_loc(user_defined_rate_parameters), &
solver_state, &
solver_stats, &
error)
if (has_error_occurred(error, errmsg, errcode)) return

Expand Down
72 changes: 1 addition & 71 deletions schemes/musica/micm/musica_ccpp_micm_util.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,80 +2,10 @@ module musica_ccpp_micm_util
implicit none

private
public :: reshape_into_micm_arr, reshape_into_ccpp_arr, convert_to_mol_per_cubic_meter, &
convert_to_mass_mixing_ratio
public :: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio

contains

!> Reshape array (2D/3D -> 1D) and convert type (kind_phys -> c_double)
subroutine reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, &
micm_temperature, micm_pressure, micm_dry_air_density, &
micm_constituents)
use iso_c_binding, only: c_double
use ccpp_kinds, only: kind_phys

real(kind_phys), intent(in) :: temperature(:,:) ! K
real(kind_phys), intent(in) :: pressure(:,:) ! Pa
real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3
real(kind_phys), intent(in) :: constituents(:,:,:) ! kg kg-1
real(c_double), intent(out) :: micm_temperature(:) ! K
real(c_double), intent(out) :: micm_pressure(:) ! Pa
real(c_double), intent(out) :: micm_dry_air_density(:) ! kg m-3
real(c_double), intent(out) :: micm_constituents(:) ! kg kg-1

! local variables
integer :: num_columns, num_layers, num_constituents
integer :: i_column, i_layer, i_elem, i_constituents

num_columns = size(constituents, dim=1)
num_layers = size(constituents, dim=2)
num_constituents = size(constituents, dim=3)

! Reshape into 1-D arry in species-column first order, referring to
! state.variables_[i_cell][i_species] = concentrations[i_species_elem++]
i_elem = 1
i_constituents = 1
do i_layer = 1, num_layers
do i_column = 1, num_columns
micm_temperature(i_elem) = real(temperature(i_column, i_layer), c_double)
micm_pressure(i_elem) = real(pressure(i_column, i_layer), c_double)
micm_dry_air_density(i_elem) = real(dry_air_density(i_column, i_layer), c_double)
micm_constituents(i_constituents : i_constituents + num_constituents - 1) &
= real(constituents(i_column, i_layer, :), c_double)
i_elem = i_elem + 1
i_constituents = i_constituents + num_constituents
end do
end do

end subroutine reshape_into_micm_arr

!> Reshape array (1D -> 3D) and convert type (c_double -> kind_phys)
subroutine reshape_into_ccpp_arr(micm_constituents, constituents)
use iso_c_binding, only: c_double
use ccpp_kinds, only: kind_phys

real(c_double), intent(in) :: micm_constituents(:) ! kg kg-1
real(kind_phys), intent(inout) :: constituents(:,:,:) ! kg kg-1

! local variables
integer :: num_columns, num_layers, num_constituents
integer :: i_column, i_layer, i_constituents

num_columns = size(constituents, dim=1)
num_layers = size(constituents, dim=2)
num_constituents = size(constituents, dim=3)

i_constituents = 1
do i_layer = 1, num_layers
do i_column = 1, num_columns
constituents(i_column, i_layer, :) &
= real(micm_constituents(i_constituents : i_constituents + num_constituents - 1), kind_phys)
i_constituents = i_constituents + num_constituents
end do
end do

end subroutine reshape_into_ccpp_arr

! Convert CAM-SIMA unit to MICM unit (kg kg-1 -> mol m-3)
subroutine convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents)
use ccpp_kinds, only: kind_phys
Expand Down
82 changes: 34 additions & 48 deletions schemes/musica/musica_ccpp.F90
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,9 @@
!> Top-level wrapper for MUSICA chemistry components
module musica_ccpp
use musica_ccpp_micm, only: micm_register, micm_init, micm_run, micm_final
use musica_ccpp_tuvx, only: tuvx_init, tuvx_run, tuvx_final
use musica_ccpp_micm, only: micm_register, micm_init, micm_run, micm_final
use musica_ccpp_namelist, only: filename_of_tuvx_micm_mapping_configuration
use musica_ccpp_tuvx, only: tuvx_init, tuvx_run, tuvx_final
use musica_util, only: index_mappings_t

implicit none
private
Expand Down Expand Up @@ -30,18 +32,20 @@ end subroutine musica_ccpp_register
subroutine musica_ccpp_init(vertical_layer_dimension, vertical_interface_dimension, &
photolysis_wavelength_grid_interfaces, errmsg, errcode)
use ccpp_kinds, only : kind_phys

use musica_ccpp_micm, only: micm
use musica_ccpp_util, only: has_error_occurred
integer, intent(in) :: vertical_layer_dimension ! (count)
integer, intent(in) :: vertical_interface_dimension ! (count)
real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! m
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode

call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
photolysis_wavelength_grid_interfaces, errmsg, errcode)
if (errcode /= 0) return
call micm_init(errmsg, errcode)
if (errcode /= 0) return
call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
photolysis_wavelength_grid_interfaces, &
micm%user_defined_reaction_rates, errmsg, errcode)
if (errcode /= 0) return

end subroutine musica_ccpp_init

Expand All @@ -56,53 +60,42 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
geopotential_height_wrt_surface_at_interface, surface_temperature, &
surface_geopotential, surface_albedo, &
standard_gravitational_acceleration, errmsg, errcode)
use musica_ccpp_micm_util, only: reshape_into_micm_arr, reshape_into_ccpp_arr
use musica_ccpp_micm_util, only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio
use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
use ccpp_kinds, only: kind_phys
use iso_c_binding, only: c_double
use musica_ccpp_micm, only: number_of_rate_parameters
use musica_ccpp_micm_util, only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio

real(kind_phys), intent(in) :: time_step ! s
real(kind_phys), intent(in) :: temperature(:,:) ! K
real(kind_phys), intent(in) :: pressure(:,:) ! Pa
real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3
real(kind_phys), intent(in) :: time_step ! s
real(kind_phys), target, intent(in) :: temperature(:,:) ! K
real(kind_phys), target, intent(in) :: pressure(:,:) ! Pa
real(kind_phys), target, intent(in) :: dry_air_density(:,:) ! kg m-3
type(ccpp_constituent_prop_ptr_t), &
intent(in) :: constituent_props(:)
real(kind_phys), intent(inout) :: constituents(:,:,:) ! kg kg-1
real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_midpoint(:,:) ! m (column, layer)
real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface)
real(kind_phys), intent(in) :: surface_temperature(:) ! K
real(kind_phys), intent(in) :: surface_geopotential(:) ! m2 s-2
real(kind_phys), intent(in) :: surface_albedo ! unitless
real(kind_phys), intent(in) :: standard_gravitational_acceleration ! m s-2
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode
intent(in) :: constituent_props(:)
real(kind_phys), target, intent(inout) :: constituents(:,:,:) ! kg kg-1
real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_midpoint(:,:) ! m (column, layer)
real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface)
real(kind_phys), intent(in) :: surface_temperature(:) ! K
real(kind_phys), intent(in) :: surface_geopotential(:) ! m2 s-2
real(kind_phys), intent(in) :: surface_albedo ! unitless
real(kind_phys), intent(in) :: standard_gravitational_acceleration ! m s-2
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode

! local variables
real(c_double), dimension(size(temperature, dim=1) &
* size(temperature, dim=2)) :: micm_temperature
real(c_double), dimension(size(pressure, dim=1) &
* size(pressure, dim=2)) :: micm_pressure
real(c_double), dimension(size(dry_air_density, dim=1) &
* size(dry_air_density, dim=2)) :: micm_dry_air_density
real(c_double), dimension(size(constituents, dim=1) &
* size(constituents, dim=2) &
* size(constituents, dim=3)) :: micm_constituents ! mol m-3
real(kind_phys), dimension(size(constituents, dim=3)) :: molar_mass_arr ! kg mol-1

! temporarily dimensioned to Chapman mechanism until mapping between MICM and TUV-x is implemented
real(c_double), dimension(size(constituents, dim=1) &
* size(constituents, dim=2) &
* 3) :: photolysis_rate_constants ! s-1
real(kind_phys), dimension(size(constituents, dim=1), &
size(constituents, dim=2), &
number_of_rate_parameters) :: rate_parameters ! various units
integer :: i_elem

! Calculate photolysis rate constants using TUV-x
call tuvx_run(temperature, dry_air_density, &
geopotential_height_wrt_surface_at_midpoint, &
geopotential_height_wrt_surface_at_interface, &
surface_temperature, surface_geopotential, &
surface_albedo, &
standard_gravitational_acceleration, &
photolysis_rate_constants, &
rate_parameters, &
errmsg, errcode)

! Get the molar mass that is set in the call to instantiate()
Expand All @@ -127,16 +120,9 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
! Convert CAM-SIMA unit to MICM unit (kg kg-1 -> mol m-3)
call convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents)

! Reshape array (2D/3D -> 1D) and convert type (kind_phys -> c_double)
call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, &
micm_temperature, micm_pressure, micm_dry_air_density, micm_constituents)

! temporarily pass in unmapped photolysis rate constants until mapping between MICM and TUV-x is implemented
call micm_run(time_step, micm_temperature, micm_pressure, micm_dry_air_density, &
photolysis_rate_constants, micm_constituents, errmsg, errcode)

! Reshape array (1D -> 3D) and convert type (c_double -> kind_phys)
call reshape_into_ccpp_arr(micm_constituents, constituents)
! Solve chemistry at the current time step
call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, &
constituents, errmsg, errcode)

! Convert MICM unit back to CAM-SIMA unit (mol m-3 -> kg kg-1)
call convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constituents)
Expand Down
12 changes: 12 additions & 0 deletions schemes/musica/musica_ccpp_namelist.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,5 +103,17 @@
<value>UNSET_PATH</value>
</values>
</entry>
<entry id="filename_of_tuvx_micm_mapping_configuration">
<type>char*512</type>
<category>musica_ccpp</category>
<group>musica_ccpp</group>
<standard_name>filename_of_tuvx_micm_mapping_configuration</standard_name>
<units>none</units>
<desc>
A configuration file for the mapping of TUV-x photolysis rates to MICM custom rate parameters
</desc>
<values>
<value>UNSET_PATH</value>
</values>

</entry_id_pg>
Loading

0 comments on commit 0ddf4d2

Please sign in to comment.