diff --git a/.gitignore b/.gitignore
new file mode 100644
index 00000000..83c5d6be
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,7 @@
+**/*.DS_Store
+**/*.pyc
+**/build/
+**/CMakeCache.txt
+**/CMakeFiles/
+.vscode
+xcode/
\ No newline at end of file
diff --git a/schemes/musica/micm/musica_ccpp_micm.F90 b/schemes/musica/micm/musica_ccpp_micm.F90
index 1ef50d5a..7b9fbae9 100644
--- a/schemes/musica/micm/musica_ccpp_micm.F90
+++ b/schemes/musica/micm/musica_ccpp_micm.F90
@@ -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
@@ -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)
@@ -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
@@ -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
diff --git a/schemes/musica/micm/musica_ccpp_micm_util.F90 b/schemes/musica/micm/musica_ccpp_micm_util.F90
index 47b92c81..2e12de18 100644
--- a/schemes/musica/micm/musica_ccpp_micm_util.F90
+++ b/schemes/musica/micm/musica_ccpp_micm_util.F90
@@ -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
diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90
index 1ef2ea63..d3512f4a 100644
--- a/schemes/musica/musica_ccpp.F90
+++ b/schemes/musica/musica_ccpp.F90
@@ -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
@@ -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
@@ -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()
@@ -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)
diff --git a/schemes/musica/musica_ccpp_namelist.xml b/schemes/musica/musica_ccpp_namelist.xml
index a8417876..5015f2d3 100644
--- a/schemes/musica/musica_ccpp_namelist.xml
+++ b/schemes/musica/musica_ccpp_namelist.xml
@@ -103,5 +103,17 @@
UNSET_PATH
+
+ char*512
+ musica_ccpp
+ musica_ccpp
+ filename_of_tuvx_micm_mapping_configuration
+ none
+
+ A configuration file for the mapping of TUV-x photolysis rates to MICM custom rate parameters
+
+
+ UNSET_PATH
+
diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90
index bcc7be09..5a12a3ee 100644
--- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90
+++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90
@@ -2,10 +2,11 @@ module musica_ccpp_tuvx
! Note: "tuvx_t" is included in an external pre-built TUV-x library that the host
! model is responsible for linking to during compilation
- use musica_tuvx, only: tuvx_t, grid_t, profile_t
- use musica_ccpp_util, only: has_error_occurred
use ccpp_kinds, only: kind_phys
use musica_ccpp_namelist, only: filename_of_tuvx_configuration
+ use musica_ccpp_util, only: has_error_occurred
+ use musica_tuvx, only: tuvx_t, grid_t, profile_t
+ use musica_util, only: mappings_t, index_mappings_t
implicit none
private
@@ -18,13 +19,18 @@ module musica_ccpp_tuvx
type(profile_t), pointer :: temperature_profile => null()
type(profile_t), pointer :: surface_albedo_profile => null()
+ type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( )
+ integer :: number_of_photolysis_rate_constants = 0
+
contains
!> Initializes TUV-x
subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
- wavelength_grid_interfaces, errmsg, errcode)
+ wavelength_grid_interfaces, micm_rate_parameter_ordering, &
+ errmsg, errcode)
use musica_tuvx, only: grid_map_t, profile_map_t, radiator_map_t
- use musica_util, only: error_t
+ use musica_util, only: error_t, configuration_t
+ use musica_ccpp_namelist, only: filename_of_tuvx_micm_mapping_configuration
use musica_ccpp_tuvx_util, only: tuvx_deallocate
use musica_ccpp_tuvx_height_grid, &
only: create_height_grid, height_grid_label, height_grid_unit
@@ -38,14 +44,17 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
integer, intent(in) :: vertical_layer_dimension ! (count)
integer, intent(in) :: vertical_interface_dimension ! (count)
real(kind_phys), intent(in) :: wavelength_grid_interfaces(:) ! m
+ type(mappings_t), intent(in) :: micm_rate_parameter_ordering ! index mappings for MICM rate parameters
character(len=512), intent(out) :: errmsg
integer, intent(out) :: errcode
! local variables
- type(grid_map_t), pointer :: grids
- type(profile_map_t), pointer :: profiles
- type(radiator_map_t), pointer :: radiators
- type(error_t) :: error
+ type(grid_map_t), pointer :: grids
+ type(profile_map_t), pointer :: profiles
+ type(radiator_map_t), pointer :: radiators
+ type(configuration_t) :: config
+ type(mappings_t), pointer :: photolysis_rate_constants_ordering
+ type(error_t) :: error
grids => grid_map_t( error )
if (has_error_occurred( error, errmsg, errcode )) return
@@ -122,7 +131,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
return
end if
- tuvx => tuvx_t( filename_of_tuvx_configuration, grids, profiles, &
+ tuvx => tuvx_t( trim(filename_of_tuvx_configuration), grids, profiles, &
radiators, error )
if (has_error_occurred( error, errmsg, errcode )) then
call tuvx_deallocate( grids, profiles, radiators, null(), height_grid, &
@@ -134,11 +143,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
wavelength_grid, temperature_profile, surface_albedo_profile )
grids => tuvx%get_grids( error )
- if (has_error_occurred( error, errmsg, errcode )) then
- deallocate( tuvx )
- tuvx => null()
- return
- end if
+ if (has_error_occurred( error, errmsg, errcode )) return
height_grid => grids%get( height_grid_label, height_grid_unit, error )
if (has_error_occurred( error, errmsg, errcode )) then
@@ -179,6 +184,25 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, &
call tuvx_deallocate( grids, profiles, null(), null(), null(), null(), &
null(), null() )
+ photolysis_rate_constants_ordering => &
+ tuvx%get_photolysis_rate_constants_ordering( error )
+ if (has_error_occurred( error, errmsg, errcode )) return
+ number_of_photolysis_rate_constants = photolysis_rate_constants_ordering%size()
+
+ call config%load_from_file( trim(filename_of_tuvx_micm_mapping_configuration), error )
+ if (has_error_occurred( error, errmsg, errcode )) then
+ deallocate( photolysis_rate_constants_ordering )
+ photolysis_rate_constants_ordering => null()
+ return
+ end if
+
+ photolysis_rate_constants_mapping => &
+ index_mappings_t( config, photolysis_rate_constants_ordering, &
+ micm_rate_parameter_ordering, error )
+ deallocate( photolysis_rate_constants_ordering )
+ photolysis_rate_constants_ordering => null()
+ if (has_error_occurred( error, errmsg, errcode )) return
+
end subroutine tuvx_init
!> Calculates photolysis rate constants for the current model conditions
@@ -188,34 +212,43 @@ subroutine tuvx_run(temperature, dry_air_density, &
surface_temperature, surface_geopotential, &
surface_albedo, &
standard_gravitational_acceleration, &
- photolysis_rate_constants, errmsg, errcode)
- use musica_util, only: error_t
- use musica_ccpp_tuvx_height_grid, only: set_height_grid_values, calculate_heights
- use musica_ccpp_tuvx_temperature, only: set_temperature_values
+ rate_parameters, errmsg, errcode)
+ use musica_util, only: error_t
+ use musica_ccpp_tuvx_height_grid, only: set_height_grid_values, calculate_heights
+ use musica_ccpp_tuvx_temperature, only: set_temperature_values
+ use musica_ccpp_util, only: has_error_occurred
use musica_ccpp_tuvx_surface_albedo, only: set_surface_albedo_values
- real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer)
- real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer)
- 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
- ! temporarily set to Chapman mechanism and 1 dimension
- ! until mapping between MICM and TUV-x is implemented
- real(kind_phys), intent(out) :: photolysis_rate_constants(:) ! s-1 (column, reaction)
- character(len=512), intent(out) :: errmsg
- integer, intent(out) :: errcode
+ real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer)
+ real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer)
+ 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
+ real(kind_phys), intent(inout) :: rate_parameters(:,:,:) ! various units (column, layer, reaction)
+ character(len=512), intent(out) :: errmsg
+ integer, intent(out) :: errcode
! local variables
real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_midpoint, dim = 2)) :: height_midpoints
real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_interface, dim = 2)) :: height_interfaces
real(kind_phys) :: reciprocal_of_gravitational_acceleration ! s2 m-1
- integer :: i_col
+ real(kind_phys), dimension(size(rate_parameters, dim=2)+2, &
+ number_of_photolysis_rate_constants) :: photolysis_rate_constants, & ! s-1
+ heating_rates ! K s-1 (TODO: check units)
+ real(kind_phys) :: solar_zenith_angle ! degrees
+ real(kind_phys) :: earth_sun_distance ! AU
+ type(error_t) :: error
+ integer :: i_col, i_level
reciprocal_of_gravitational_acceleration = 1.0_kind_phys / standard_gravitational_acceleration
+ ! surface albedo with respect to direct UV/visible radiation
+ call set_surface_albedo_values( surface_albedo_profile, surface_albedo, errmsg, errcode )
+ if (errcode /= 0) return
+
do i_col = 1, size(temperature, dim=1)
call calculate_heights( geopotential_height_wrt_surface_at_midpoint(i_col,:), &
geopotential_height_wrt_surface_at_interface(i_col,:), &
@@ -229,14 +262,28 @@ subroutine tuvx_run(temperature, dry_air_density, &
call set_temperature_values( temperature_profile, temperature(i_col,:), &
surface_temperature(i_col), errmsg, errcode )
if (errcode /= 0) return
- end do
-
- ! surface albedo with respect to direct UV/visible radiation
- call set_surface_albedo_values( surface_albedo_profile, surface_albedo, errmsg, errcode )
- if (errcode /= 0) return
- ! stand-in until actual photolysis rate constants are calculated
- photolysis_rate_constants(:) = 1.0e-6_kind_phys
+ ! temporary values until these are available from the host model
+ solar_zenith_angle = 0.0_kind_phys
+ earth_sun_distance = 1.0_kind_phys
+
+ ! calculate photolysis rate constants and heating rates
+ call tuvx%run( solar_zenith_angle, earth_sun_distance, &
+ photolysis_rate_constants(:,:), heating_rates(:,:), &
+ error )
+ if (has_error_occurred( error, errmsg, errcode )) return
+
+ ! filter out negative photolysis rate constants
+ photolysis_rate_constants(:,:) = &
+ max( photolysis_rate_constants(:,:), 0.0_kind_phys )
+
+ ! map photolysis rate constants to the host model's rate parameters and vertical grid
+ do i_level = 1, size(rate_parameters, dim=2)
+ call photolysis_rate_constants_mapping%copy_data( &
+ photolysis_rate_constants(size(rate_parameters, dim=2)-i_level+2,:), &
+ rate_parameters(i_col,i_level,:) )
+ end do
+ end do
end subroutine tuvx_run
@@ -273,6 +320,11 @@ subroutine tuvx_final(errmsg, errcode)
tuvx => null()
end if
+ if (associated( photolysis_rate_constants_mapping )) then
+ deallocate( photolysis_rate_constants_mapping )
+ photolysis_rate_constants_mapping => null()
+ end if
+
end subroutine tuvx_final
end module musica_ccpp_tuvx
\ No newline at end of file
diff --git a/test/docker/Dockerfile.musica b/test/docker/Dockerfile.musica
index 65940f1a..83dea485 100644
--- a/test/docker/Dockerfile.musica
+++ b/test/docker/Dockerfile.musica
@@ -5,7 +5,7 @@
FROM ubuntu:22.04
-ARG MUSICA_GIT_TAG=dbbdb22f5f2807e27c2695db85291951ba178634
+ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8
RUN apt update \
&& apt install -y sudo \
@@ -68,6 +68,12 @@ RUN cd musica \
COPY . atmospheric_physics
RUN sudo chown -R test_user:test_user atmospheric_physics
+# Clone the MUSICA chemistry data set repository
+RUN git clone https://github.com/NCAR/cam-sima-chemistry-data.git \
+ && cd cam-sima-chemistry-data \
+ && git checkout 144d982715f5bd6fa61e76f255a52db1843c55f2 \
+ && mv mechanisms /home/test_user/atmospheric_physics/schemes/musica/configurations
+
# Must make ccpp-framework available before building test
RUN cd atmospheric_physics/test \
&& mkdir lib \
@@ -84,10 +90,4 @@ RUN cd atmospheric_physics/test \
-D CCPP_ENABLE_MEMCHECK=ON \
&& cmake --build ./build
-RUN cd atmospheric_physics/test \
- && mkdir third_party \
- && cd third_party \
- && git clone --depth 1 https://github.com/NCAR/tuv-x.git \
- && cp -r tuv-x/data ../build/data
-
WORKDIR /home/test_user/atmospheric_physics/test/build
\ No newline at end of file
diff --git a/test/docker/Dockerfile.musica.no_install b/test/docker/Dockerfile.musica.no_install
index b59489af..e8f70fb0 100644
--- a/test/docker/Dockerfile.musica.no_install
+++ b/test/docker/Dockerfile.musica.no_install
@@ -8,7 +8,7 @@
FROM ubuntu:22.04
-ARG MUSICA_GIT_TAG=dbbdb22f5f2807e27c2695db85291951ba178634
+ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8
RUN apt update \
&& apt install -y sudo \
@@ -55,6 +55,12 @@ ENV MUSICA_GIT_TAG=${MUSICA_GIT_TAG}
COPY . atmospheric_physics
RUN sudo chown -R test_user:test_user atmospheric_physics
+# Clone the MUSICA chemistry data set repository
+RUN git clone https://github.com/NCAR/cam-sima-chemistry-data.git \
+ && cd cam-sima-chemistry-data \
+ && git checkout 144d982715f5bd6fa61e76f255a52db1843c55f2 \
+ && mv mechanisms /home/test_user/atmospheric_physics/schemes/musica/configurations
+
# Must make ccpp-framework available before building test
RUN cd atmospheric_physics/test \
&& mkdir lib \
@@ -71,10 +77,4 @@ RUN cd atmospheric_physics/test \
-D CCPP_ENABLE_MEMCHECK=ON \
&& cmake --build ./build
-RUN cd atmospheric_physics/test \
- && mkdir third_party \
- && cd third_party \
- && git clone --depth 1 https://github.com/NCAR/tuv-x.git \
- && cp -r tuv-x/data ../build/data
-
WORKDIR /home/test_user/atmospheric_physics/test/build
\ No newline at end of file
diff --git a/test/musica/CMakeLists.txt b/test/musica/CMakeLists.txt
index e41405c1..74b29693 100644
--- a/test/musica/CMakeLists.txt
+++ b/test/musica/CMakeLists.txt
@@ -56,14 +56,10 @@ add_test(
add_memory_check_test(test_musica_api $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY})
+# copy MUSICA configuration data into the build folder
add_custom_target(
copy_micm_configs ALL ${CMAKE_COMMAND} -E copy_directory
- ${CMAKE_CURRENT_SOURCE_DIR}/micm/configs/chapman ${CMAKE_BINARY_DIR}/chapman
-)
-
-add_custom_target(
- copy_tuvx_configs ALL ${CMAKE_COMMAND} -E copy_directory
- ${CMAKE_CURRENT_SOURCE_DIR}/tuvx/configs ${CMAKE_BINARY_DIR}/configs
+ ${CMAKE_CURRENT_SOURCE_DIR}/../../schemes/musica/configurations ${CMAKE_BINARY_DIR}/musica_configurations
)
add_subdirectory(micm)
diff --git a/test/musica/micm/configs/chapman/config.json b/test/musica/micm/configs/chapman/config.json
deleted file mode 100644
index 04d0ef28..00000000
--- a/test/musica/micm/configs/chapman/config.json
+++ /dev/null
@@ -1,6 +0,0 @@
-{
- "camp-files": [
- "species.json",
- "reactions.json"
- ]
-}
diff --git a/test/musica/micm/configs/chapman/reactions.json b/test/musica/micm/configs/chapman/reactions.json
deleted file mode 100644
index 32af8982..00000000
--- a/test/musica/micm/configs/chapman/reactions.json
+++ /dev/null
@@ -1,97 +0,0 @@
-{
- "camp-data": [
- {
- "name": "Chapman",
- "type": "MECHANISM",
- "reactions": [
- {
- "type": "PHOTOLYSIS",
- "reactants": {
- "O2": {}
- },
- "products": {
- "O": {
- "yield": 2.0
- }
- },
- "MUSICA name": "jO2"
- },
- {
- "type": "ARRHENIUS",
- "A": 8.018e-17,
- "reactants": {
- "O": {},
- "O2": {}
- },
- "products": {
- "O3": {}
- },
- "MUSICA name": "R2"
- },
- {
- "type": "PHOTOLYSIS",
- "reactants": {
- "O3": {}
- },
- "products": {
- "O": {},
- "O2": {}
- },
- "MUSICA name": "jO3->O"
- },
- {
- "type": "ARRHENIUS",
- "A": 1.576e-15,
- "reactants": {
- "O": {},
- "O3": {}
- },
- "products": {
- "O2": {
- "yield": 2.0
- }
- },
- "MUSICA name": "R4"
- },
- {
- "type": "PHOTOLYSIS",
- "reactants": {
- "O3": {}
- },
- "products": {
- "O1D": {},
- "O2": {}
- },
- "MUSICA name": "jO3->O1D"
- },
- {
- "type": "ARRHENIUS",
- "A": 7.11e-11,
- "reactants": {
- "O1D": {},
- "M": {}
- },
- "products": {
- "O": {},
- "M": {}
- },
- "MUSICA name": "R6"
- },
- {
- "type": "ARRHENIUS",
- "A": 1.2e-10,
- "reactants": {
- "O1D": {},
- "O3": {}
- },
- "products": {
- "O2": {
- "yield": 2.0
- }
- },
- "MUSICA name": "R7"
- }
- ]
- }
- ]
-}
\ No newline at end of file
diff --git a/test/musica/micm/configs/chapman/species.json b/test/musica/micm/configs/chapman/species.json
deleted file mode 100644
index aa2c2dcf..00000000
--- a/test/musica/micm/configs/chapman/species.json
+++ /dev/null
@@ -1,39 +0,0 @@
-{
- "camp-data": [
- {
- "name": "M",
- "type": "CHEM_SPEC",
- "tracer type": "THIRD_BODY",
- "molecular weight [kg mol-1]": 0.029,
- "__is advected": false
- },
- {
- "name": "O2",
- "type": "CHEM_SPEC",
- "tracer type": "CONSTANT",
- "molecular weight [kg mol-1]": 0.032,
- "__is advected": false
- },
- {
- "name": "O",
- "type": "CHEM_SPEC",
- "absolute tolerance": 1e-12,
- "molecular weight [kg mol-1]": 0.016,
- "__is advected": false
- },
- {
- "name": "O1D",
- "type": "CHEM_SPEC",
- "absolute tolerance": 1e-12,
- "molecular weight [kg mol-1]": 0.016,
- "__is advected": false
- },
- {
- "name": "O3",
- "type": "CHEM_SPEC",
- "absolute tolerance": 1e-12,
- "molecular weight [kg mol-1]": 0.048,
- "__is advected": true
- }
- ]
-}
diff --git a/test/musica/micm/test_micm_util.F90 b/test/musica/micm/test_micm_util.F90
index 798f1408..3b5cf4fa 100644
--- a/test/musica/micm/test_micm_util.F90
+++ b/test/musica/micm/test_micm_util.F90
@@ -7,71 +7,10 @@ program test_micm_util
#define ASSERT(x) if (.not.(x)) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: x"; stop 1; endif
#define ASSERT_NEAR( a, b, abs_error ) if( (abs(a - b) >= abs_error) .and. (abs(a - b) /= 0.0) ) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: a, b"; stop 1; endif
- call test_reshape()
call test_unit_conversion()
contains
- subroutine test_reshape()
- use iso_c_binding, only: c_double
- use ccpp_kinds, only: kind_phys
-
- integer, parameter :: NUM_SPECIES = 4
- integer, parameter :: NUM_COLUMNS = 2
- integer, parameter :: NUM_LAYERS = 2
- integer, parameter :: NUM_GRID_CELLS = NUM_COLUMNS * NUM_LAYERS
- real, parameter :: ABS_ERROR = 1e-7
- real(c_double) :: expected_conditions(NUM_GRID_CELLS) = [100.0, 200.0, 300.0, 400.0]
- real(c_double) :: expected_constituents(NUM_GRID_CELLS*NUM_SPECIES) = &
- [0.1, 0.2, 0.3, 0.4, 0.21, 0.22, 0.23, 0.24, 0.41, 0.42, 0.43, 0.44, 0.31, 0.32, 0.33, 0.34]
- real(kind_phys) :: temperature(NUM_COLUMNS,NUM_LAYERS)
- real(kind_phys) :: pressure(NUM_COLUMNS,NUM_LAYERS)
- real(kind_phys) :: dry_air_density(NUM_COLUMNS,NUM_LAYERS)
- real(kind_phys) :: constituents(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES)
- real(c_double) :: micm_temperature(NUM_GRID_CELLS)
- real(c_double) :: micm_pressure(NUM_GRID_CELLS)
- real(c_double) :: micm_dry_air_density(NUM_GRID_CELLS)
- real(c_double) :: micm_constituents(NUM_GRID_CELLS*NUM_SPECIES)
- integer :: i_column, i_layer, i_elem, i_arr
-
- temperature(:,1) = (/ 100._kind_phys, 200._kind_phys /)
- temperature(:,2) = (/ 300._kind_phys, 400._kind_phys /)
- pressure(:,1) = (/ 100._kind_phys, 200._kind_phys /)
- pressure(:,2) = (/ 300._kind_phys, 400._kind_phys /)
- dry_air_density(:,1) = (/ 100._kind_phys, 200._kind_phys /)
- dry_air_density(:,2) = (/ 300._kind_phys, 400._kind_phys /)
- constituents(1,1,:) = (/ 0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys, 0.4_kind_phys /)
- constituents(1,2,:) = (/ 0.41_kind_phys, 0.42_kind_phys, 0.43_kind_phys, 0.44_kind_phys /)
- constituents(2,1,:) = (/ 0.21_kind_phys, 0.22_kind_phys, 0.23_kind_phys, 0.24_kind_phys /)
- constituents(2,2,:) = (/ 0.31_kind_phys, 0.32_kind_phys, 0.33_kind_phys, 0.34_kind_phys /)
-
- call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, &
- micm_temperature, micm_pressure, micm_dry_air_density, micm_constituents)
-
- do i_elem = 1, NUM_GRID_CELLS
- ASSERT(micm_temperature(i_elem) == expected_conditions(i_elem))
- ASSERT(micm_pressure(i_elem) == expected_conditions(i_elem))
- ASSERT(micm_dry_air_density(i_elem) == expected_conditions(i_elem))
- end do
-
- do i_elem = 1, size(micm_constituents)
- ASSERT_NEAR(micm_constituents(i_elem), expected_constituents(i_elem), ABS_ERROR)
- end do
-
- call reshape_into_ccpp_arr(micm_constituents, constituents)
-
- i_arr = 1
- do i_layer = 1, NUM_LAYERS
- do i_column = 1, NUM_COLUMNS
- do i_elem = 1, NUM_SPECIES
- ASSERT_NEAR(constituents(i_column, i_layer, i_elem), expected_constituents(i_arr), ABS_ERROR)
- i_arr = i_arr + 1
- end do
- end do
- end do
-
- end subroutine test_reshape
-
subroutine test_unit_conversion()
use ccpp_kinds, only: kind_phys
diff --git a/test/musica/musica_ccpp_namelist.F90 b/test/musica/musica_ccpp_namelist.F90
index 69c803b3..d27cc754 100644
--- a/test/musica/musica_ccpp_namelist.F90
+++ b/test/musica/musica_ccpp_namelist.F90
@@ -4,9 +4,9 @@ module musica_ccpp_namelist
implicit none
private
- public :: filename_of_micm_configuration, filename_of_tuvx_configuration
-
- character(len=*), parameter :: filename_of_micm_configuration = 'chapman'
- character(len=*), parameter :: filename_of_tuvx_configuration = 'configs/ts1_tsmlt.json'
+
+ character(len=250), public :: filename_of_micm_configuration = 'musica_configurations/chapman/micm/config.json'
+ character(len=250), public :: filename_of_tuvx_configuration = 'musica_configurations/chapman/tuvx/config.json'
+ character(len=250), public :: filename_of_tuvx_micm_mapping_configuration = 'musica_configurations/chapman/tuvx_micm_mapping.json'
end module musica_ccpp_namelist
diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90
index 6b4d4154..61cca2f7 100644
--- a/test/musica/test_musica_api.F90
+++ b/test/musica/test_musica_api.F90
@@ -5,51 +5,25 @@ program run_test_musica_ccpp
implicit none
#define ASSERT(x) if (.not.(x)) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: x"; stop 1; endif
+#define ASSERT_NEAR( a, b, abs_error ) if( (abs(a - b) >= abs_error) .and. (abs(a - b) /= 0.0) ) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: a, b"; stop 1; endif
- call test_musica_ccpp_api()
+ call test_chapman()
+ call test_terminator()
contains
- subroutine test_musica_ccpp_api()
- use musica_micm, only: Rosenbrock, RosenbrockStandardOrder
- use ccpp_kinds, only: kind_phys
- use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
- use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t
+ !> Returns the wavelength edges used in the CAM-Chem photolysis rate constant lookup table
+ !! These are the values that will be used in CAM-SIMA and correspond to the wavelength
+ !! bins used in the CAM-Chem photolysis rate constant lookup table.
+ !!
+ !! We're using the actual values here because several of the TS1/TSMLT photolysis
+ !! rate constant configurations are sensitive to the wavelength grid.
+ subroutine get_wavelength_edges(edges)
+ use ccpp_kinds, only: kind_phys
- integer, parameter :: NUM_SPECIES = 4
- integer, parameter :: NUM_COLUMNS = 2
- integer, parameter :: NUM_LAYERS = 2
- integer, parameter :: NUM_WAVELENGTH_BINS = 102
- integer :: num_grid_cells = NUM_COLUMNS * NUM_LAYERS
- integer :: solver_type = Rosenbrock
- integer :: errcode
- character(len=512) :: errmsg
- real(kind_phys) :: time_step = 60._kind_phys ! s
- real(kind_phys), dimension(NUM_WAVELENGTH_BINS+1) :: photolysis_wavelength_grid_interfaces ! m
- real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: geopotential_height_wrt_surface_at_midpoint ! m
- real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m
- real(kind_phys), dimension(NUM_COLUMNS) :: surface_temperature ! K
- real(kind_phys), dimension(NUM_COLUMNS) :: surface_geopotential ! m2 s-2
- real(kind_phys) :: surface_albedo ! unitless
- real(kind_phys) :: standard_gravitational_acceleration ! s2 m-1
- real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: temperature ! K
- real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: pressure ! Pa
- real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: dry_air_density ! kg m-3
- real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: constituents ! kg kg-1
- type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:)
- type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:)
- type(ccpp_constituent_properties_t), pointer :: const_prop
- real(kind_phys) :: molar_mass
- character(len=512) :: species_name, units
- logical :: tmp_bool, is_advected
- integer :: i
-
- ! These are the values that will be used in CAM-SIMA and correspond to the wavelength
- ! bins used in the CAM-Chem photolysis rate constant lookup table.
- !
- ! We're using the actual values here because several of the TS1/TSMLT photolysis
- ! rate constant configurations are sensitive to the wavelength grid.
- photolysis_wavelength_grid_interfaces = (/ &
+ real(kind_phys), dimension(:), intent(out) :: edges
+
+ edges = (/ &
120.0e-9_kind_phys, &
121.4e-9_kind_phys, &
121.9e-9_kind_phys, &
@@ -154,6 +128,60 @@ subroutine test_musica_ccpp_api()
700.0e-9_kind_phys, &
750.0e-9_kind_phys &
/)
+
+ end subroutine get_wavelength_edges
+
+ !> Tests the Chapman chemistry scheme
+ subroutine test_chapman()
+ use musica_micm, only: Rosenbrock, RosenbrockStandardOrder
+ use ccpp_kinds, only: kind_phys
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t
+ use musica_ccpp_micm, only: micm
+ use musica_ccpp_namelist, only: filename_of_micm_configuration, &
+ filename_of_tuvx_configuration, &
+ filename_of_tuvx_micm_mapping_configuration
+
+ implicit none
+
+ integer, parameter :: NUM_SPECIES = 5
+ ! This test requires that the number of grid cells = 4, which is the default
+ ! vector dimension for MICM. This restriction will be removed once
+ ! https://github.com/NCAR/musica/issues/217 is finished.
+ integer, parameter :: NUM_COLUMNS = 2
+ integer, parameter :: NUM_LAYERS = 2
+ integer, parameter :: NUM_WAVELENGTH_BINS = 102
+ integer :: NUM_GRID_CELLS = NUM_COLUMNS * NUM_LAYERS
+ integer :: solver_type = Rosenbrock
+ integer :: errcode
+ character(len=512) :: errmsg
+ real(kind_phys) :: time_step = 60._kind_phys ! s
+ real(kind_phys), dimension(NUM_WAVELENGTH_BINS+1) :: photolysis_wavelength_grid_interfaces ! m
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: geopotential_height_wrt_surface_at_midpoint ! m
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m
+ real(kind_phys), dimension(NUM_COLUMNS) :: surface_temperature ! K
+ real(kind_phys), dimension(NUM_COLUMNS) :: surface_geopotential ! m2 s-2
+ real(kind_phys) :: surface_albedo ! unitless
+ real(kind_phys) :: standard_gravitational_acceleration ! s2 m-1
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: temperature ! K
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: pressure ! Pa
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: dry_air_density ! kg m-3
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: constituents ! kg kg-1
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: initial_constituents ! kg kg-1
+ type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:)
+ type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:)
+ type(ccpp_constituent_properties_t), pointer :: const_prop
+ real(kind_phys) :: molar_mass, base_conc
+ character(len=512) :: species_name, units
+ character(len=:), allocatable :: micm_species_name
+ logical :: tmp_bool, is_advected
+ integer :: i, j, k
+ integer :: N2_index, O2_index, O_index, O1D_index, O3_index
+ real(kind_phys) :: total_O, total_O_init
+
+ call get_wavelength_edges(photolysis_wavelength_grid_interfaces)
+ solver_type = Rosenbrock
+ time_step = 60._kind_phys
geopotential_height_wrt_surface_at_midpoint(1,:) = (/ 2000.0_kind_phys, 500.0_kind_phys /)
geopotential_height_wrt_surface_at_midpoint(2,:) = (/ 2000.0_kind_phys, -500.0_kind_phys /)
geopotential_height_wrt_surface_at_interface(1,:) = (/ 3000.0_kind_phys, 1000.0_kind_phys, 0.0_kind_phys /)
@@ -168,12 +196,16 @@ subroutine test_musica_ccpp_api()
pressure(:,2) = (/ 8000.04_kind_phys, 9000.04_kind_phys /)
dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /)
dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /)
- constituents(1,1,:) = (/ 0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys, 0.4_kind_phys /)
- constituents(1,2,:) = (/ 0.41_kind_phys, 0.42_kind_phys, 0.43_kind_phys, 0.44_kind_phys /)
- constituents(2,1,:) = (/ 0.21_kind_phys, 0.22_kind_phys, 0.23_kind_phys, 0.24_kind_phys /)
- constituents(2,2,:) = (/ 0.31_kind_phys, 0.32_kind_phys, 0.33_kind_phys, 0.34_kind_phys /)
- call musica_ccpp_register(solver_type, num_grid_cells, constituent_props, errmsg, errcode)
+ filename_of_micm_configuration = 'musica_configurations/chapman/micm/config.json'
+ filename_of_tuvx_configuration = 'musica_configurations/chapman/tuvx/config.json'
+ filename_of_tuvx_micm_mapping_configuration = 'musica_configurations/chapman/tuvx_micm_mapping.json'
+
+ call musica_ccpp_register(solver_type, NUM_GRID_CELLS, constituent_props, errmsg, errcode)
+ if (errcode /= 0) then
+ write(*,*) trim(errmsg)
+ stop 3
+ endif
ASSERT(allocated(constituent_props))
ASSERT(size(constituent_props) == NUM_SPECIES)
do i = 1, size(constituent_props)
@@ -185,13 +217,213 @@ subroutine test_musica_ccpp_api()
ASSERT(errcode == 0)
call constituent_props(i)%is_advected(is_advected, errcode, errmsg)
ASSERT(errcode == 0)
- tmp_bool = (trim(species_name) == "O2" .and. molar_mass == 0.032_kind_phys .and. .not. is_advected) .or. &
- (trim(species_name) == "O" .and. molar_mass == 0.016_kind_phys .and. .not. is_advected) .or. &
- (trim(species_name) == "O1D" .and. molar_mass == 0.016_kind_phys .and. .not. is_advected) .or. &
- (trim(species_name) == "O3" .and. molar_mass == 0.048_kind_phys .and. is_advected)
+ tmp_bool = (trim(species_name) == "O2" .and. molar_mass == 0.0319988_kind_phys .and. is_advected) .or. &
+ (trim(species_name) == "O" .and. molar_mass == 0.0159994_kind_phys .and. .not. is_advected) .or. &
+ (trim(species_name) == "O1D" .and. molar_mass == 0.0159994_kind_phys .and. .not. is_advected) .or. &
+ (trim(species_name) == "O3" .and. molar_mass == 0.0479982_kind_phys .and. is_advected) .or. &
+ (trim(species_name) == "N2" .and. molar_mass == 0.0280134_kind_phys .and. is_advected)
ASSERT(tmp_bool)
call constituent_props(i)%units(units, errcode, errmsg)
+ if (errcode /= 0) then
+ write(*,*) errcode, trim(errmsg)
+ stop 3
+ endif
+ ASSERT(trim(units) == 'kg kg-1')
+ end do
+ if (errcode /= 0) then
+ write(*,*) errcode, trim(errmsg)
+ stop 3
+ end if
+
+ allocate(constituent_props_ptr(size(constituent_props)))
+ do i = 1, size(constituent_props)
+ const_prop => constituent_props(i)
+ call constituent_props_ptr(i)%set(const_prop, errcode, errmsg)
+ end do
+
+ call musica_ccpp_init(NUM_LAYERS, NUM_LAYERS+1, photolysis_wavelength_grid_interfaces, &
+ errmsg, errcode)
+ if (errcode /= 0) then
+ write(*,*) trim(errmsg)
+ stop 3
+ endif
+
+ do i = 1, micm%species_ordering%size()
+ micm_species_name = micm%species_ordering%name(i)
+ if (micm_species_name == "O2") then
+ O2_index = i
+ base_conc = 0.21_kind_phys
+ else if (micm_species_name == "O") then
+ O_index = i
+ base_conc = 1.0e-9_kind_phys
+ else if (micm_species_name == "O1D") then
+ O1D_index = i
+ base_conc = 1.0e-9_kind_phys
+ else if (micm_species_name == "O3") then
+ O3_index = i
+ base_conc = 1.0e-4_kind_phys
+ else if (micm_species_name == "N2") then
+ N2_index = i
+ base_conc = 0.79_kind_phys
+ else
+ write(*,*) "Unknown species: ", micm_species_name
+ stop 3
+ endif
+ do j = 1, NUM_COLUMNS
+ do k = 1, NUM_LAYERS
+ constituents(j,k,i) = base_conc * (1.0 + 0.1 * (j-1) + 0.01 * (k-1))
+ end do
+ end do
+ end do
+ initial_constituents(:,:,:) = constituents(:,:,:)
+
+ write(*,*) "[MUSICA INFO] Initial Time Step"
+ write(*,fmt="(1x,f10.2)") time_step
+ write(*,*) "[MUSICA INFO] Initial Temperature"
+ write(*,fmt="(4(1x,f10.4))") temperature
+ write(*,*) "[MUSICA INFO] Initial Pressure"
+ write(*,fmt="(4(1x,f10.4))") pressure
+ write(*,*) "[MUSICA INFO] Initial Concentrations"
+ write(*,fmt="(4(3x,e13.6))") constituents
+
+ call musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props_ptr, &
+ constituents, geopotential_height_wrt_surface_at_midpoint, &
+ geopotential_height_wrt_surface_at_interface, surface_temperature, &
+ surface_geopotential, surface_albedo, standard_gravitational_acceleration, &
+ errmsg, errcode)
+ if (errcode /= 0) then
+ write(*,*) trim(errmsg)
+ stop 3
+ endif
+
+ write(*,*) "[MUSICA INFO] Solved Concentrations"
+ write(*,fmt="(4(3x,e13.6))") constituents
+
+ call musica_ccpp_final(errmsg, errcode)
+
+ if (errcode /= 0) then
+ write(*,*) trim(errmsg)
+ stop 3
+ endif
+
+ ! Check the results
+ do i = 1, NUM_COLUMNS
+ do j = 1, NUM_LAYERS
+ ! N2 should be unchanged
+ ASSERT_NEAR(constituents(i,j,N2_index), initial_constituents(i,j,N2_index), 1.0e-13)
+ ! O3 and O2 should be relatively unchanged
+ ASSERT_NEAR(constituents(i,j,O3_index), initial_constituents(i,j,O3_index), 1.0e-4)
+ ASSERT_NEAR(constituents(i,j,O2_index), initial_constituents(i,j,O2_index), 1.0e-4)
+ ! O and O1D should be < 1e-10 kg kg-1
+ ASSERT(constituents(i,j,O_index) < 1.0e-10)
+ ASSERT(constituents(i,j,O1D_index) < 1.0e-10)
+ ! total O mass should be conserved
+ total_O = constituents(i,j,O_index) + constituents(i,j,O1D_index) + &
+ constituents(i,j,O2_index) + constituents(i,j,O3_index)
+ total_O_init = initial_constituents(i,j,O_index) + initial_constituents(i,j,O1D_index) + &
+ initial_constituents(i,j,O2_index) + initial_constituents(i,j,O3_index)
+ ASSERT_NEAR(total_O, total_O_init, 1.0e-13)
+ end do
+ end do
+
+ deallocate(constituent_props_ptr)
+
+ end subroutine test_chapman
+
+ !> Tests the simple Terminator chemistry scheme
+ subroutine test_terminator()
+ use musica_micm, only: Rosenbrock, RosenbrockStandardOrder
+ use ccpp_kinds, only: kind_phys
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
+ use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t
+ use musica_ccpp_micm, only: micm
+ use musica_ccpp_namelist, only: filename_of_micm_configuration, &
+ filename_of_tuvx_configuration, &
+ filename_of_tuvx_micm_mapping_configuration
+
+ implicit none
+
+ integer, parameter :: NUM_SPECIES = 2
+ ! This test requires that the number of grid cells = 4, which is the default
+ ! vector dimension for MICM. This restriction will be removed once
+ ! https://github.com/NCAR/musica/issues/217 is finished.
+ integer, parameter :: NUM_COLUMNS = 2
+ integer, parameter :: NUM_LAYERS = 2
+ integer, parameter :: NUM_WAVELENGTH_BINS = 102
+ integer :: NUM_GRID_CELLS = NUM_COLUMNS * NUM_LAYERS
+ integer :: solver_type = Rosenbrock
+ integer :: errcode
+ character(len=512) :: errmsg
+ real(kind_phys) :: time_step = 60._kind_phys ! s
+ real(kind_phys), dimension(NUM_WAVELENGTH_BINS+1) :: photolysis_wavelength_grid_interfaces ! m
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: geopotential_height_wrt_surface_at_midpoint ! m
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m
+ real(kind_phys), dimension(NUM_COLUMNS) :: surface_temperature ! K
+ real(kind_phys), dimension(NUM_COLUMNS) :: surface_geopotential ! m2 s-2
+ real(kind_phys) :: surface_albedo ! unitless
+ real(kind_phys) :: standard_gravitational_acceleration ! s2 m-1
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: temperature ! K
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: pressure ! Pa
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: dry_air_density ! kg m-3
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: constituents ! kg kg-1
+ real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: initial_constituents ! kg kg-1
+ type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:)
+ type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:)
+ type(ccpp_constituent_properties_t), pointer :: const_prop
+ real(kind_phys) :: molar_mass, base_conc
+ character(len=512) :: species_name, units
+ character(len=:), allocatable :: micm_species_name
+ logical :: tmp_bool, is_advected
+ integer :: i, j, k
+ integer :: Cl_index, Cl2_index
+ real(kind_phys) :: total_Cl, total_Cl_init
+
+ call get_wavelength_edges(photolysis_wavelength_grid_interfaces)
+ solver_type = Rosenbrock
+ time_step = 60._kind_phys
+ geopotential_height_wrt_surface_at_midpoint(1,:) = (/ 2000.0_kind_phys, 500.0_kind_phys /)
+ geopotential_height_wrt_surface_at_midpoint(2,:) = (/ 2000.0_kind_phys, -500.0_kind_phys /)
+ geopotential_height_wrt_surface_at_interface(1,:) = (/ 3000.0_kind_phys, 1000.0_kind_phys, 0.0_kind_phys /)
+ geopotential_height_wrt_surface_at_interface(2,:) = (/ 3000.0_kind_phys, 500.0_kind_phys, -1500.0_kind_phys /)
+ surface_temperature = (/ 300.0_kind_phys, 300.0_kind_phys /)
+ surface_geopotential = (/ 100.0_kind_phys, 200.0_kind_phys /)
+ surface_albedo = 0.10_kind_phys
+ standard_gravitational_acceleration = 10.0_kind_phys
+ temperature(:,1) = (/ 100._kind_phys, 200._kind_phys /)
+ temperature(:,2) = (/ 300._kind_phys, 400._kind_phys /)
+ pressure(:,1) = (/ 6000.04_kind_phys, 7000.04_kind_phys /)
+ pressure(:,2) = (/ 8000.04_kind_phys, 9000.04_kind_phys /)
+ dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /)
+ dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /)
+
+ filename_of_micm_configuration = 'musica_configurations/terminator/micm/config.json'
+ filename_of_tuvx_configuration = 'musica_configurations/terminator/tuvx/config.json'
+ filename_of_tuvx_micm_mapping_configuration = 'musica_configurations/terminator/tuvx_micm_mapping.json'
+
+ call musica_ccpp_register(solver_type, NUM_GRID_CELLS, constituent_props, errmsg, errcode)
+ if (errcode /= 0) then
+ write(*,*) trim(errmsg)
+ stop 3
+ endif
+ ASSERT(allocated(constituent_props))
+ ASSERT(size(constituent_props) == NUM_SPECIES)
+ do i = 1, size(constituent_props)
+ ASSERT(constituent_props(i)%is_instantiated(errcode, errmsg))
+ ASSERT(errcode == 0)
+ call constituent_props(i)%standard_name(species_name, errcode, errmsg)
+ ASSERT(errcode == 0)
+ call constituent_props(i)%molar_mass(molar_mass, errcode, errmsg)
ASSERT(errcode == 0)
+ call constituent_props(i)%is_advected(is_advected, errcode, errmsg)
+ ASSERT(errcode == 0)
+ tmp_bool = (trim(species_name) == "Cl" .and. molar_mass == 0.035453_kind_phys .and. is_advected) .or. &
+ (trim(species_name) == "Cl2" .and. molar_mass == 0.070906_kind_phys .and. is_advected)
+ ASSERT(tmp_bool)
+ call constituent_props(i)%units(units, errcode, errmsg)
+ if (errcode /= 0) then
+ write(*,*) errcode, trim(errmsg)
+ stop 3
+ endif
ASSERT(trim(units) == 'kg kg-1')
end do
if (errcode /= 0) then
@@ -212,6 +444,26 @@ subroutine test_musica_ccpp_api()
stop 3
endif
+ do i = 1, micm%species_ordering%size()
+ micm_species_name = micm%species_ordering%name(i)
+ if (micm_species_name == "Cl") then
+ Cl_index = i
+ base_conc = 1.0e-10_kind_phys
+ else if (micm_species_name == "Cl2") then
+ Cl2_index = i
+ base_conc = 1.0e-6_kind_phys
+ else
+ write(*,*) "Unknown species: ", micm_species_name
+ stop 3
+ endif
+ do j = 1, NUM_COLUMNS
+ do k = 1, NUM_LAYERS
+ constituents(j,k,i) = base_conc * (1.0 + 0.1 * (j-1) + 0.01 * (k-1))
+ end do
+ end do
+ end do
+ initial_constituents(:,:,:) = constituents(:,:,:)
+
write(*,*) "[MUSICA INFO] Initial Time Step"
write(*,fmt="(1x,f10.2)") time_step
write(*,*) "[MUSICA INFO] Initial Temperature"
@@ -241,8 +493,22 @@ subroutine test_musica_ccpp_api()
stop 3
endif
+ ! Check the results
+ do i = 1, NUM_COLUMNS
+ do j = 1, NUM_LAYERS
+ ! Cl2 should not be unchanged
+ ASSERT(constituents(i,j,Cl2_index) .ne. initial_constituents(i,j,Cl2_index))
+ ! Cl should not be unchanged
+ ASSERT(constituents(i,j,Cl_index) .ne. initial_constituents(i,j,Cl_index))
+ ! total Cl mass should be conserved
+ total_Cl = constituents(i,j,Cl_index) + constituents(i,j,Cl2_index)
+ total_Cl_init = initial_constituents(i,j,Cl_index) + initial_constituents(i,j,Cl2_index)
+ ASSERT_NEAR(total_Cl, total_Cl_init, 1.0e-13)
+ end do
+ end do
+
deallocate(constituent_props_ptr)
- end subroutine test_musica_ccpp_api
+ end subroutine test_terminator
end program run_test_musica_ccpp
\ No newline at end of file