From e8a29b37fb86af36035d585c6c87dd815e007927 Mon Sep 17 00:00:00 2001 From: Jiwon Gim <55209567+boulderdaze@users.noreply.github.com> Date: Tue, 14 Jan 2025 10:01:02 -0700 Subject: [PATCH] Set gas-species profiles in TUV-x and map indices between constituents and MICM (#184) Originator(s): @boulderdaze Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): - Closes #98 - Closes #165 Describe any changes made to the namelist: N/A List all files eliminated and why: N/A List all files added and what they do: ``` A schemes/musica/musica_ccpp_species.F90 A schemes/musica/tuvx/musica_ccpp_tuvx_gas_species.F90 A schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 A test/musica/test_musica_species.F90 A test/musica/tuvx/test_tuvx_gas_species.F90 A test/musica/tuvx/test_tuvx_load_species.F90 ``` List all existing files that have been modified, and describe the changes: ``` M schemes/musica/micm/musica_ccpp_micm.F90 M schemes/musica/musica_ccpp.F90 M schemes/musica/musica_ccpp.meta M schemes/musica/tuvx/musica_ccpp_tuvx.F90 M schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 M test/docker/Dockerfile.musica M test/docker/Dockerfile.musica.no_install M test/musica/CMakeLists.txt M test/musica/test_musica_api.F90 M test/musica/tuvx/CMakeLists.txt M test/musica/tuvx/test_tuvx_height_grid.F90 M to_be_ccppized/ccpp_tuvx_utils.F90 ``` List any test failures: N/A Is this a science-changing update? New physics package, algorithm change, tuning changes, etc? No --- schemes/musica/micm/musica_ccpp_micm.F90 | 24 +- schemes/musica/musica_ccpp.F90 | 142 +++--- schemes/musica/musica_ccpp.meta | 2 +- schemes/musica/musica_ccpp_species.F90 | 323 ++++++++++++ schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 251 +++++---- .../tuvx/musica_ccpp_tuvx_gas_species.F90 | 145 ++++++ .../tuvx/musica_ccpp_tuvx_height_grid.F90 | 12 +- .../tuvx/musica_ccpp_tuvx_load_species.F90 | 219 ++++++++ test/docker/Dockerfile.musica | 2 +- test/docker/Dockerfile.musica.no_install | 2 +- test/musica/CMakeLists.txt | 39 +- test/musica/test_musica_api.F90 | 77 ++- test/musica/test_musica_species.F90 | 217 ++++++++ test/musica/tuvx/CMakeLists.txt | 67 +++ test/musica/tuvx/test_tuvx_gas_species.F90 | 476 ++++++++++++++++++ test/musica/tuvx/test_tuvx_height_grid.F90 | 11 +- test/musica/tuvx/test_tuvx_load_species.F90 | 362 +++++++++++++ to_be_ccppized/ccpp_tuvx_utils.F90 | 2 +- 18 files changed, 2187 insertions(+), 186 deletions(-) create mode 100644 schemes/musica/musica_ccpp_species.F90 create mode 100644 schemes/musica/tuvx/musica_ccpp_tuvx_gas_species.F90 create mode 100644 schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 create mode 100644 test/musica/test_musica_species.F90 create mode 100644 test/musica/tuvx/test_tuvx_gas_species.F90 create mode 100644 test/musica/tuvx/test_tuvx_load_species.F90 diff --git a/schemes/musica/micm/musica_ccpp_micm.F90 b/schemes/musica/micm/musica_ccpp_micm.F90 index c2e40d95..d51b51cc 100644 --- a/schemes/musica/micm/musica_ccpp_micm.F90 +++ b/schemes/musica/micm/musica_ccpp_micm.F90 @@ -19,15 +19,16 @@ module musica_ccpp_micm !> Registers MICM constituent properties with the CCPP subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, & - errmsg, errcode) + micm_species, errmsg, errcode) use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t - use musica_micm, only: Rosenbrock, RosenbrockStandardOrder + use musica_ccpp_species, only: musica_species_t use musica_util, only: error_t use iso_c_binding, only: c_int integer(c_int), intent(in) :: solver_type integer(c_int), intent(in) :: number_of_grid_cells type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:) + type(musica_species_t), allocatable, intent(out) :: micm_species(:) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -36,6 +37,7 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, & real(kind=kind_phys) :: molar_mass character(len=:), allocatable :: species_name logical :: is_advected + integer :: number_of_species integer :: i, species_index if (associated( micm )) then @@ -46,13 +48,20 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, & number_of_grid_cells, error) if (has_error_occurred(error, errmsg, errcode)) return - allocate(constituent_props(micm%species_ordering%size()), stat=errcode) + number_of_species = micm%species_ordering%size() + allocate(constituent_props(number_of_species), stat=errcode) if (errcode /= 0) then errmsg = "[MUSICA Error] Failed to allocate memory for constituent properties." return end if - do i = 1, micm%species_ordering%size() + allocate(micm_species(number_of_species), stat=errcode) + if (errcode /= 0) then + errmsg = "[MUSICA Error] Failed to allocate memory for micm species." + return + end if + + do i = 1, number_of_species associate( map => micm%species_ordering ) species_name = map%name(i) species_index = map%index(i) @@ -78,6 +87,13 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, & errcode = errcode, & errmsg = errmsg) if (errcode /= 0) return + + ! Species are ordered to match the sequence of the MICM state array + micm_species(species_index) = musica_species_t( & + name = species_name, & + unit = 'kg kg-1', & + molar_mass = molar_mass, & + index_musica_species = species_index ) end associate ! map end do number_of_rate_parameters = micm%user_defined_reaction_rates%size() diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index 4c79511f..70d9791d 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -14,16 +14,20 @@ module musica_ccpp !> \section arg_table_musica_ccpp_register Argument Table !! \htmlinclude musica_ccpp_register.html - subroutine musica_ccpp_register(constituent_props, errmsg, & - errcode) - use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t - use musica_ccpp_namelist, only: micm_solver_type + subroutine musica_ccpp_register(constituent_props, errmsg, errcode) + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_namelist, only: micm_solver_type + use musica_ccpp_species, only: musica_species_t, register_musica_species + use musica_ccpp_tuvx_load_species, only: check_tuvx_species_initialization type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode + ! local variables type(ccpp_constituent_properties_t), allocatable :: constituent_props_subset(:) + type(musica_species_t), allocatable :: micm_species(:) + type(musica_species_t), allocatable :: tuvx_species(:) integer :: number_of_grid_cells ! Temporary fix until the number of grid cells is only needed to create a MICM state @@ -32,15 +36,20 @@ subroutine musica_ccpp_register(constituent_props, errmsg, & ! the solver when the number of grid cells is known at the init stage. number_of_grid_cells = 1 call micm_register(micm_solver_type, number_of_grid_cells, constituent_props_subset, & - errmsg, errcode) + micm_species, errmsg, errcode) if (errcode /= 0) return constituent_props = constituent_props_subset deallocate(constituent_props_subset) - call tuvx_register(constituent_props_subset, errmsg, errcode) + call tuvx_register(micm_species, tuvx_species, constituent_props_subset, & + errmsg, errcode) if (errcode /= 0) return constituent_props = [ constituent_props, constituent_props_subset ] + call register_musica_species(micm_species, tuvx_species) + call check_tuvx_species_initialization(errmsg, errcode) + if (errcode /= 0) return + end subroutine musica_ccpp_register !> \section arg_table_musica_ccpp_init Argument Table @@ -48,34 +57,46 @@ end subroutine musica_ccpp_register subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, & vertical_interface_dimension, & photolysis_wavelength_grid_interfaces, & - constituent_props, errmsg, errcode) - use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t, ccpp_constituent_prop_ptr_t - use ccpp_kinds, only : kind_phys - use musica_ccpp_micm, only: micm - use musica_ccpp_namelist, only: micm_solver_type - use musica_ccpp_util, only: has_error_occurred + constituent_props_ptr, errmsg, errcode) + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t, ccpp_constituent_properties_t + use ccpp_kinds, only: kind_phys + use musica_ccpp_micm, only: micm + use musica_ccpp_namelist, only: micm_solver_type + use musica_ccpp_util, only: has_error_occurred + use musica_ccpp_species, only: initialize_musica_species_indices, initialize_molar_mass_array, & + check_initialization, musica_species_t + integer, intent(in) :: horizontal_dimension ! (count) integer, intent(in) :: vertical_layer_dimension ! (count) integer, intent(in) :: vertical_interface_dimension ! (count) real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! m - type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) + type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props_ptr(:) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode - integer :: number_of_grid_cells - type(ccpp_constituent_properties_t), allocatable :: micm_species_props(:) + ! local variables + type(ccpp_constituent_properties_t), allocatable :: constituent_props(:) + type(musica_species_t), allocatable :: micm_species(:) + integer :: number_of_grid_cells ! Temporary fix until the number of grid cells is only needed to create a MICM state ! instead of when the solver is created. ! Re-create the MICM solver with the correct number of grid cells number_of_grid_cells = horizontal_dimension * vertical_layer_dimension - call micm_register(micm_solver_type, number_of_grid_cells, micm_species_props, errmsg, errcode) + call micm_register(micm_solver_type, number_of_grid_cells, constituent_props, & + micm_species, errmsg, errcode) 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, & - constituent_props, errmsg, errcode) + photolysis_wavelength_grid_interfaces, & + micm%user_defined_reaction_rates, errmsg, errcode) + if (errcode /= 0) return + + call initialize_musica_species_indices(constituent_props_ptr, errmsg, errcode) + if (errcode /= 0) return + call initialize_molar_mass_array(constituent_props_ptr, errmsg, errcode) + if (errcode /= 0) return + call check_initialization(errmsg, errcode) if (errcode /= 0) return end subroutine musica_ccpp_init @@ -98,6 +119,9 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co use ccpp_kinds, only: kind_phys 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 + use musica_ccpp_species, only: number_of_micm_species, number_of_tuvx_species, & + micm_indices_constituent_props, tuvx_indices_constituent_props, micm_molar_mass_array, & + extract_subset_constituents, update_constituents real(kind_phys), intent(in) :: time_step ! s real(kind_phys), target, intent(in) :: temperature(:,:) ! K @@ -122,69 +146,71 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co integer, intent(out) :: errcode ! local variables - real(kind_phys), dimension(size(constituents, dim=3)) :: molar_mass_arr ! kg mol-1 real(kind_phys), dimension(size(constituents, dim=1), & size(constituents, dim=2), & - number_of_rate_parameters) :: rate_parameters ! various units - integer :: i_elem + number_of_rate_parameters) :: rate_parameters ! various units + real(kind_phys), dimension(size(constituents, dim=1), & + size(constituents, dim=2), & + number_of_micm_species) :: constituents_micm_species ! kg kg-1 + real(kind_phys), dimension(size(constituents, dim=1), & + size(constituents, dim=2), & + number_of_tuvx_species) :: constituents_tuvx_species ! kg kg-1 + + call extract_subset_constituents(tuvx_indices_constituent_props, constituents, & + constituents_tuvx_species, errmsg, errcode) + if (errcode /= 0) return ! Calculate photolysis rate constants using TUV-x - call tuvx_run(temperature, dry_air_density, & - constituents, & - geopotential_height_wrt_surface_at_midpoint, & - geopotential_height_wrt_surface_at_interface, & - surface_geopotential, surface_temperature, & - surface_albedo, & - photolysis_wavelength_grid_interfaces, & - extraterrestrial_flux, & - standard_gravitational_acceleration, & - cloud_area_fraction, & - air_pressure_thickness, & - solar_zenith_angle, & - earth_sun_distance, & - rate_parameters, & + call tuvx_run(temperature, dry_air_density, & + constituents_tuvx_species, & + geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, & + surface_geopotential, surface_temperature, & + surface_albedo, & + photolysis_wavelength_grid_interfaces, & + extraterrestrial_flux, & + standard_gravitational_acceleration, & + cloud_area_fraction, & + air_pressure_thickness, & + solar_zenith_angle, & + earth_sun_distance, & + rate_parameters, & errmsg, errcode) - ! Get the molar mass that is set in the call to instantiate() - do i_elem = 1, size(molar_mass_arr) - call constituent_props(i_elem)%molar_mass(molar_mass_arr(i_elem), errcode, errmsg) - if (errcode /= 0) then - errmsg = "[MUSICA Error] Unable to get molar mass." - return - end if - end do - - ! TODO(jiwon) Check molar mass is non zero as it becomes a denominator for unit converison - ! this code will be deleted when the framework does the check - do i_elem = 1, size(molar_mass_arr) - if (molar_mass_arr(i_elem) <= 0) then - errcode = 1 - errmsg = "[MUSICA Error] Molar mass must be greater than zero." - return - end if - end do + call update_constituents(tuvx_indices_constituent_props, constituents_tuvx_species, & + constituents, errmsg, errcode) + if (errcode /= 0) return + call extract_subset_constituents(micm_indices_constituent_props, constituents, & + constituents_micm_species, errmsg, errcode) + if (errcode /= 0) return ! 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) + call convert_to_mol_per_cubic_meter(dry_air_density, micm_molar_mass_array, constituents_micm_species) ! Solve chemistry at the current time step call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, & - constituents, errmsg, errcode) + constituents_micm_species, 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) + call convert_to_mass_mixing_ratio(dry_air_density, micm_molar_mass_array, constituents_micm_species) + call update_constituents(micm_indices_constituent_props, constituents_micm_species, & + constituents, errmsg, errcode) + if (errcode /= 0) return + end subroutine musica_ccpp_run !> \section arg_table_musica_ccpp_final Argument Table !! \htmlinclude musica_ccpp_final.html subroutine musica_ccpp_final(errmsg, errcode) + use musica_ccpp_species, only: cleanup_musica_species character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode + call cleanup_musica_species() call tuvx_final(errmsg, errcode) call micm_final(errmsg, errcode) end subroutine musica_ccpp_final -end module musica_ccpp +end module musica_ccpp \ No newline at end of file diff --git a/schemes/musica/musica_ccpp.meta b/schemes/musica/musica_ccpp.meta index 801c3c79..bad6f271 100644 --- a/schemes/musica/musica_ccpp.meta +++ b/schemes/musica/musica_ccpp.meta @@ -53,7 +53,7 @@ type = real | kind = kind_phys dimensions = (photolysis_wavelength_grid_interface_dimension) intent = in -[ constituent_props ] +[ constituent_props_ptr ] standard_name = ccpp_constituent_properties units = None type = ccpp_constituent_prop_ptr_t diff --git a/schemes/musica/musica_ccpp_species.F90 b/schemes/musica/musica_ccpp_species.F90 new file mode 100644 index 00000000..3f746a2f --- /dev/null +++ b/schemes/musica/musica_ccpp_species.F90 @@ -0,0 +1,323 @@ +module musica_ccpp_species + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: cleanup_musica_species, register_musica_species, initialize_musica_species_indices, & + initialize_molar_mass_array, extract_subset_constituents, update_constituents, & + check_initialization + + integer, parameter, public :: MUSICA_INT_UNASSIGNED = -99999 + + !> Definition of musica species object + type, public :: musica_species_t + character(len=:), allocatable :: name + character(len=:), allocatable :: unit + real(kind_phys) :: molar_mass = 0.0_kind_phys ! kg mol-1 + integer :: index_musica_species = MUSICA_INT_UNASSIGNED + integer :: index_constituent_props = MUSICA_INT_UNASSIGNED + logical :: profiled = .false. ! TUV-x gas species optional + real(kind_phys) :: scale_height = 0.0_kind_phys ! km, TUV-x gas species optional + contains + ! Deallocates memory associated with this musica species object + procedure :: deallocate => musica_species_t_deallocate + end type musica_species_t + + interface musica_species_t + procedure species_t_constructor + end interface musica_species_t + + ! Species are ordered to match the sequence of the MICM state array + type(musica_species_t), allocatable, protected, public :: micm_species_set(:) + type(musica_species_t), allocatable, protected, public :: tuvx_species_set(:) + integer, allocatable, protected, public :: micm_indices_constituent_props(:) + integer, allocatable, protected, public :: tuvx_indices_constituent_props(:) + real(kind_phys), allocatable, protected, public :: micm_molar_mass_array(:) ! kg mol-1 + integer, protected, public :: number_of_micm_species = MUSICA_INT_UNASSIGNED + integer, protected, public :: number_of_tuvx_species = MUSICA_INT_UNASSIGNED + +contains + + !> Constructor for musica species object + function species_t_constructor(name, unit, molar_mass, scale_height, & + index_musica_species, index_constituent_props) result( this ) + character(len=*), intent(in) :: name + character(len=*), intent(in) :: unit + real(kind_phys), intent(in) :: molar_mass ! kg mol-1 + real(kind_phys), intent(in) :: scale_height ! km + integer, intent(in) :: index_musica_species + integer, intent(in) :: index_constituent_props + type(musica_species_t) :: this + + this%name = name + this%unit = unit + this%molar_mass = molar_mass + this%scale_height = scale_height + this%index_musica_species = index_musica_species + this%index_constituent_props = index_constituent_props + + end function species_t_constructor + + !> Deallocates memory associated with this musica species object + subroutine musica_species_t_deallocate(this) + class(musica_species_t), intent(inout) :: this + + if (allocated(this%name)) deallocate(this%name) + if (allocated(this%unit)) deallocate(this%unit) + this%molar_mass = 0.0_kind_phys + this%index_musica_species = MUSICA_INT_UNASSIGNED + this%index_constituent_props = MUSICA_INT_UNASSIGNED + this%profiled = .false. + this%scale_height = 0.0_kind_phys + + end subroutine musica_species_t_deallocate + + !> Deallocates the memory associated with the MUSICA species object, its indices array, + ! and its molar mass array + subroutine cleanup_musica_species() + integer :: i + + if (allocated( micm_species_set )) then + do i = 1, size(micm_species_set) + call micm_species_set(i)%deallocate() + end do + deallocate( micm_species_set ) + end if + + if (allocated( tuvx_species_set )) then + do i = 1, size(tuvx_species_set) + call tuvx_species_set(i)%deallocate() + end do + deallocate( tuvx_species_set ) + end if + + if (allocated( micm_indices_constituent_props )) deallocate( micm_indices_constituent_props ) + if (allocated( tuvx_indices_constituent_props )) deallocate( tuvx_indices_constituent_props ) + if (allocated( micm_molar_mass_array )) deallocate( micm_molar_mass_array ) + + end subroutine cleanup_musica_species + + !> Allocates memory and initializes the species array for MICM and TUV-x + subroutine register_musica_species(micm_species, tuvx_species) + type(musica_species_t), intent(in) :: micm_species(:) + type(musica_species_t), intent(in) :: tuvx_species(:) + + number_of_micm_species = size(micm_species) + allocate( micm_species_set( number_of_micm_species ) ) + micm_species_set = micm_species + + number_of_tuvx_species = size(tuvx_species) + allocate( tuvx_species_set( number_of_tuvx_species ) ) + tuvx_species_set = tuvx_species + + end subroutine register_musica_species + + !> Retrieves the species indices from the constituents array and store them + subroutine find_musica_species_indices(constituent_props, musica_species_set, & + indices_constituent_props, errmsg, errcode) + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use ccpp_const_utils, only: ccpp_const_get_idx + + type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) + type(musica_species_t), intent(inout) :: musica_species_set(:) + integer, intent(inout) :: indices_constituent_props(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variables + integer :: i_elem, index_species + + do i_elem = 1, size(musica_species_set) + call ccpp_const_get_idx( constituent_props, musica_species_set(i_elem)%name, & + musica_species_set(i_elem)%index_constituent_props, errmsg, errcode ) + if (errcode /= 0) return + + index_species = musica_species_set(i_elem)%index_constituent_props + if (index_species == MUSICA_INT_UNASSIGNED) then + errmsg = "[MUSICA Error] Unable to find index for " // musica_species_set(i_elem)%name + errcode = 1 + return + end if + indices_constituent_props(i_elem) = index_species + end do + + end subroutine find_musica_species_indices + + !> Initializes arrays to store the species indices of the CCPP constituents + subroutine initialize_musica_species_indices(constituent_props, errmsg, errcode) + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + + type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + if (.not. allocated( micm_species_set ) .or. & + .not. allocated( tuvx_species_set )) then + errmsg = "[MUSICA Error] The MUSICA species set(s) are not allocated." + errcode = 1 + return + end if + + allocate( micm_indices_constituent_props( size(micm_species_set) ) ) + call find_musica_species_indices( constituent_props, micm_species_set, & + micm_indices_constituent_props, errmsg, errcode ) + if (errcode /= 0) return + + allocate( tuvx_indices_constituent_props( size(tuvx_species_set) ) ) + call find_musica_species_indices( constituent_props, tuvx_species_set, & + tuvx_indices_constituent_props, errmsg, errcode ) + if (errcode /= 0) return + + end subroutine initialize_musica_species_indices + + !> Iterates through the constituent property pointer array to populate the molar mass array, + ! storing molar mass values in a sequence that matches the order of the MICM state array + subroutine initialize_molar_mass_array(constituent_props, errmsg, errcode) + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + + type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variables + integer :: i_elem, index_species + + if (.not. allocated( micm_species_set )) then + errmsg = "[MUSICA Error] The MICM species set is not allocated." + errcode = 1 + return + end if + + allocate( micm_molar_mass_array( size(micm_species_set) ) ) + do i_elem = 1, size(micm_species_set) + call constituent_props( micm_species_set(i_elem)%index_constituent_props ) & + %molar_mass( micm_molar_mass_array(i_elem), errcode, errmsg ) + if (errcode /= 0) then + errmsg = "[MUSICA Error] Unable to get molar mass for " // micm_species_set(i_elem)%name + return + end if + end do + + ! TODO(jiwon) - This code block can be removed once the CCPP framework handles + ! the check for non-zero molar mass + do i_elem = 1, size(micm_molar_mass_array) + if (micm_molar_mass_array(i_elem) <= 0) then + errcode = 1 + errmsg = "[MUSICA Error] Molar mass must be greater than zero for " & + // micm_species_set(i_elem)%name + return + end if + end do + + end subroutine initialize_molar_mass_array + + !> Extracts sub-constituents array using the indices from constituents array + subroutine extract_subset_constituents(indices_constituent_props, constituents, & + subset_constituents, errmsg, errcode) + integer, intent(in) :: indices_constituent_props(:) + real(kind_phys), intent(in) :: constituents(:,:,:) ! kg kg-1 + real(kind_phys), intent(inout) :: subset_constituents(:,:,:) ! kg kg-1 + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variables + integer :: i + + errmsg = '' + errcode = 0 + + if (size(subset_constituents, dim=3) /= size(indices_constituent_props)) then + errmsg = "[MUSICA Error] The given dimension for the constituents & + doesn't match the size of species indices array." + errcode = 1 + return + end if + + do i = 1, size(indices_constituent_props) + subset_constituents(:,:,i) = constituents(:,:,indices_constituent_props(i)) + end do + + end subroutine extract_subset_constituents + + !> The updated subset of constituents is added back to the constituents array + ! using the indices array to specify where the updates should occur + subroutine update_constituents(indices_constituent_props, subset_constituents, & + constituents, errmsg, errcode) + integer, intent(in) :: indices_constituent_props(:) + real(kind_phys), intent(in) :: subset_constituents(:,:,:) ! kg kg-1 + real(kind_phys), intent(inout) :: constituents(:,:,:) ! kg kg-1 + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variables + integer :: i + + errmsg = '' + errcode = 0 + + if (size(subset_constituents, dim=3) /= size(indices_constituent_props)) then + errmsg = "[MUSICA Error] The given dimension for the constituents & + doesn't match the size of species indices array." + errcode = 1 + return + end if + + do i= 1, size(indices_constituent_props) + constituents(:,:,indices_constituent_props(i)) = subset_constituents(:,:,i) + end do + + end subroutine update_constituents + + !> Checks that the musica species-related objects are initialized, + ! including the MICM species set, TUV-x species set, their constituent property indices, + ! and the molar mass values for the MICM species. + ! This function is specifically designed to ensure that the musica species are properly + ! initialized, so they don't need to be checked during the run phase. + subroutine check_initialization(errmsg, errcode) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + errmsg = '' + errcode = 0 + + if (.not. allocated( micm_species_set )) then + errmsg = "[MUSICA Error] MICM species set has not been allocated." + errcode = 1 + return + end if + if (.not. allocated( micm_indices_constituent_props )) then + errmsg = "[MUSICA Error] MICM species indices array has not been allocated." + errcode = 1 + return + end if + if (.not. allocated( tuvx_species_set )) then + errmsg = "[MUSICA Error] TUV-X species set has not been allocated." + errcode = 1 + return + end if + + if (.not. allocated( tuvx_indices_constituent_props )) then + errmsg = "[MUSICA Error] TUV-X species indices array has not been allocated." + errcode = 1 + return + end if + if (.not. allocated( micm_molar_mass_array )) then + errmsg = "[MUSICA Error] MICM molar mass array has not been allocated." + errcode = 1 + return + end if + if (number_of_micm_species == MUSICA_INT_UNASSIGNED) then + errmsg = "[MUSICA Error] The 'number_of_micm_species' variable has not been initialized." + errcode = 1 + return + end if + if (number_of_tuvx_species == MUSICA_INT_UNASSIGNED) then + errmsg = "[MUSICA Error] The 'number_of_tuvx_species' variable has not been initialized." + errcode = 1 + return + end if + + end subroutine check_initialization + +end module musica_ccpp_species \ No newline at end of file diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index 681a1952..1732ffab 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -22,19 +22,14 @@ module musica_ccpp_tuvx type(profile_t), pointer :: temperature_profile => null() type(profile_t), pointer :: surface_albedo_profile => null() type(profile_t), pointer :: extraterrestrial_flux_profile => null() + type(profile_t), pointer :: dry_air_profile => null() + type(profile_t), pointer :: O2_profile => null() + type(profile_t), pointer :: O3_profile => null() type(radiator_t), pointer :: cloud_optics => null() type(radiator_t), pointer :: aerosol_optics => null() type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( ) integer, parameter :: DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS = 0 integer :: number_of_photolysis_rate_constants = DEFAULT_NUM_PHOTOLYSIS_RATE_CONSTANTS - integer, parameter :: DEFAULT_INDEX_NOT_FOUND = -1 - character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_LABEL = & - 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' - character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_LONG_NAME = & - 'Cloud water mass mixing ratio with respect to moist air plus all airborne condensates' - character(len=*), parameter :: CLOUD_LIQUID_WATER_CONTENT_UNITS = 'kg kg-1' - real(kind_phys), parameter :: CLOUD_LIQUID_WATER_CONTENT_MOLAR_MASS = 0.018_kind_phys ! kg mol-1 - integer :: index_cloud_liquid_water_content = DEFAULT_INDEX_NOT_FOUND contains @@ -52,7 +47,7 @@ subroutine reset_tuvx_map_state( grids, profiles, radiators ) end subroutine reset_tuvx_map_state - !> This is a helper subroutine created to deallocate objects associated with TUV-x + !> Deallocates objects associated with TUV-x subroutine cleanup_tuvx_resources() if (associated( height_grid )) then @@ -80,6 +75,21 @@ subroutine cleanup_tuvx_resources() extraterrestrial_flux_profile => null() end if + if (associated( dry_air_profile )) then + deallocate( dry_air_profile ) + dry_air_profile => null() + end if + + if (associated( O2_profile )) then + deallocate( O2_profile ) + O2_profile => null() + end if + + if (associated( O3_profile )) then + deallocate( O3_profile ) + O3_profile => null() + end if + if (associated( cloud_optics )) then deallocate( cloud_optics ) cloud_optics => null() @@ -98,33 +108,21 @@ subroutine cleanup_tuvx_resources() end subroutine cleanup_tuvx_resources !> Registers constituent properties with the CCPP needed by TUV-x - subroutine tuvx_register(constituent_props, errmsg, errcode) - use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t - use musica_util, only: error_t - - type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:) - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode - - allocate(constituent_props(1), stat=errcode) - if (errcode /= 0) then - errmsg = "[MUSICA Error] Failed to allocate memory for constituent properties." - return - end if - - ! Register cloud liquid water content needed for cloud optics calculations - call constituent_props(1)%instantiate( & - std_name = CLOUD_LIQUID_WATER_CONTENT_LABEL, & - long_name = CLOUD_LIQUID_WATER_CONTENT_LONG_NAME, & - units = CLOUD_LIQUID_WATER_CONTENT_UNITS, & - vertical_dim = "vertical_layer_dimension", & - default_value = 0.0_kind_phys, & - min_value = 0.0_kind_phys, & - molar_mass = CLOUD_LIQUID_WATER_CONTENT_MOLAR_MASS, & - advected = .true., & - errcode = errcode, & - errmsg = errmsg & - ) + subroutine tuvx_register(micm_species, tuvx_species, constituent_props, & + errmsg, errcode) + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_species, only: musica_species_t + use musica_ccpp_tuvx_load_species, only: configure_tuvx_species + use musica_util, only: error_t + + type(musica_species_t), intent(inout) :: micm_species(:) + type(musica_species_t), allocatable, intent(out) :: tuvx_species(:) + type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + call configure_tuvx_species(micm_species, tuvx_species, constituent_props, & + errmsg, errcode) if (errcode /= 0) return end subroutine tuvx_register @@ -132,8 +130,7 @@ end subroutine tuvx_register !> Initializes TUV-x subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & wavelength_grid_interfaces, micm_rate_parameter_ordering, & - constituent_props, errmsg, errcode) - use ccpp_const_utils, only: ccpp_const_get_idx + errmsg, errcode) use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t use musica_tuvx, only: grid_map_t, profile_map_t, radiator_map_t use musica_util, only: error_t, configuration_t @@ -150,18 +147,21 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & use musica_ccpp_tuvx_extraterrestrial_flux, & only: create_extraterrestrial_flux_profile, extraterrestrial_flux_label, & extraterrestrial_flux_unit + use musica_ccpp_tuvx_gas_species, & + only: create_dry_air_profile, create_O2_profile, create_O3_profile use musica_ccpp_tuvx_cloud_optics, & only: create_cloud_optics_radiator, cloud_optics_label use musica_ccpp_tuvx_aerosol_optics, & only: create_aerosol_optics_radiator, aerosol_optics_label + use musica_ccpp_tuvx_load_species, & + only: DRY_AIR_LABEL, O2_LABEL, O3_LABEL, TUVX_GAS_SPECIES_UNITS - 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 - type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode + 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 @@ -171,16 +171,6 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & type(mappings_t), pointer :: photolysis_rate_constants_ordering type(error_t) :: error - ! Get needed indices in constituents array - call ccpp_const_get_idx(constituent_props, CLOUD_LIQUID_WATER_CONTENT_LABEL, & - index_cloud_liquid_water_content, errmsg, errcode) - if (errcode /= 0) return - if (index_cloud_liquid_water_content == DEFAULT_INDEX_NOT_FOUND) then - errmsg = "[MUSICA Error] Unable to find index for cloud liquid water content." - errcode = 1 - return - end if - grids => grid_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) return @@ -264,6 +254,48 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & return end if + dry_air_profile => create_dry_air_profile( height_grid, errmsg, errcode ) + if (errcode /= 0) then + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + endif + + call profiles%add( dry_air_profile, error ) + if (has_error_occurred( error, errmsg, errcode )) then + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + end if + + O2_profile => create_O2_profile( height_grid, errmsg, errcode ) + if (errcode /= 0) then + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + endif + + call profiles%add( O2_profile, error ) + if (has_error_occurred( error, errmsg, errcode )) then + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + end if + + O3_profile => create_O3_profile( height_grid, errmsg, errcode ) + if (errcode /= 0) then + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + endif + + call profiles%add( O3_profile, error ) + if (has_error_occurred( error, errmsg, errcode )) then + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + end if + radiators => radiator_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) then call reset_tuvx_map_state( grids, profiles, null() ) @@ -377,6 +409,33 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & return end if + dry_air_profile => profiles%get( DRY_AIR_LABEL, TUVX_GAS_SPECIES_UNITS, error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( tuvx ) + tuvx => null() + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + end if + + O2_profile => profiles%get( O2_LABEL, TUVX_GAS_SPECIES_UNITS, error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( tuvx ) + tuvx => null() + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + end if + + O3_profile => profiles%get( O3_LABEL, TUVX_GAS_SPECIES_UNITS, error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( tuvx ) + tuvx => null() + call reset_tuvx_map_state( grids, profiles, null() ) + call cleanup_tuvx_resources() + return + end if + radiators => tuvx%get_radiators( error ) if (has_error_occurred( error, errmsg, errcode )) then deallocate( tuvx ) @@ -442,33 +501,36 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & end subroutine tuvx_init !> Calculates photolysis rate constants for the current model conditions - subroutine tuvx_run(temperature, dry_air_density, & - constituents, & - geopotential_height_wrt_surface_at_midpoint, & - geopotential_height_wrt_surface_at_interface, & - surface_geopotential, surface_temperature, & - surface_albedo, & - photolysis_wavelength_grid_interfaces, & - extraterrestrial_flux, & - standard_gravitational_acceleration, & - cloud_area_fraction, & - air_pressure_thickness, & - solar_zenith_angle, & - earth_sun_distance, & - rate_parameters, & + subroutine tuvx_run(temperature, dry_air_density, & + constituents, & + geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, & + surface_geopotential, surface_temperature, & + surface_albedo, & + photolysis_wavelength_grid_interfaces, & + extraterrestrial_flux, & + standard_gravitational_acceleration, & + cloud_area_fraction, & + air_pressure_thickness, & + solar_zenith_angle, & + earth_sun_distance, & + 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, PI - use musica_ccpp_tuvx_surface_albedo, only: set_surface_albedo_values - use musica_ccpp_tuvx_extraterrestrial_flux, only: set_extraterrestrial_flux_values - use musica_ccpp_tuvx_cloud_optics, only: set_cloud_optics_values - use musica_ccpp_tuvx_aerosol_optics, only: set_aerosol_optics_values + 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, PI + use musica_ccpp_tuvx_surface_albedo, only: set_surface_albedo_values + use musica_ccpp_tuvx_extraterrestrial_flux, only: set_extraterrestrial_flux_values + use musica_ccpp_tuvx_cloud_optics, only: set_cloud_optics_values + use musica_ccpp_tuvx_aerosol_optics, only: set_aerosol_optics_values + use musica_ccpp_tuvx_load_species, only: index_cloud_liquid_water_content, & + index_dry_air, index_O2, index_O3 + use musica_ccpp_tuvx_gas_species, only: set_gas_species_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) :: constituents(:,:,:) ! various (column, layer, constituent) + real(kind_phys), intent(in) :: constituents(:,:,:) ! kg kg-1 (column, layer, constituent) 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_geopotential(:) ! m2 s-2 @@ -488,6 +550,7 @@ subroutine tuvx_run(temperature, dry_air_density, & ! 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), dimension(size(height_interfaces)) :: height_deltas ! km 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) @@ -502,8 +565,8 @@ subroutine tuvx_run(temperature, dry_air_density, & call set_surface_albedo_values( surface_albedo_profile, surface_albedo, errmsg, errcode ) if (errcode /= 0) return - call set_extraterrestrial_flux_values( extraterrestrial_flux_profile, & - photolysis_wavelength_grid_interfaces, & + call set_extraterrestrial_flux_values( extraterrestrial_flux_profile, & + photolysis_wavelength_grid_interfaces, & extraterrestrial_flux, errmsg, errcode ) if (errcode /= 0) return @@ -519,9 +582,9 @@ subroutine tuvx_run(temperature, dry_air_density, & geopotential_height_wrt_surface_at_interface(i_col,:), & surface_geopotential(i_col), & reciprocal_of_gravitational_acceleration, & - height_midpoints, height_interfaces ) + height_midpoints, height_interfaces ) call set_height_grid_values( height_grid, height_midpoints, height_interfaces, & - errmsg, errcode ) + height_deltas, errmsg, errcode ) if (errcode /= 0) return call set_temperature_values( temperature_profile, temperature(i_col,:), & @@ -529,10 +592,24 @@ subroutine tuvx_run(temperature, dry_air_density, & if (errcode /= 0) return call set_cloud_optics_values( cloud_optics, cloud_area_fraction(i_col,:), & - air_pressure_thickness(i_col,:), & - constituents(i_col,:,index_cloud_liquid_water_content), & - reciprocal_of_gravitational_acceleration, & - errmsg, errcode ) + air_pressure_thickness(i_col,:), & + constituents(i_col,:,index_cloud_liquid_water_content), & + reciprocal_of_gravitational_acceleration, errmsg, errcode ) + if (errcode /= 0) return + + call set_gas_species_values( dry_air_profile, dry_air_density(i_col,:), & + constituents(i_col,:,index_dry_air), height_deltas, index_dry_air, & + errmsg, errcode ) + if (errcode /= 0) return + + call set_gas_species_values( O2_profile, dry_air_density(i_col,:), & + constituents(i_col,:,index_O2), height_deltas, index_O2, & + errmsg, errcode ) + if (errcode /= 0) return + + call set_gas_species_values( O3_profile, dry_air_density(i_col,:), & + constituents(i_col,:,index_O3), height_deltas, index_O3, & + errmsg, errcode ) if (errcode /= 0) return call set_aerosol_optics_values( aerosol_optics, errmsg, errcode ) @@ -546,7 +623,7 @@ subroutine tuvx_run(temperature, dry_air_density, & ! filter out negative photolysis rate constants photolysis_rate_constants(:,:) = & - max( photolysis_rate_constants(:,:), 0.0_kind_phys ) + max( photolysis_rate_constants(:,:), 0.0_kind_phys ) end if ! solar zenith angle check ! map photolysis rate constants to the host model's rate parameters and vertical grid diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species.F90 new file mode 100644 index 00000000..74b73b17 --- /dev/null +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_gas_species.F90 @@ -0,0 +1,145 @@ +module musica_ccpp_tuvx_gas_species + use ccpp_kinds, only: kind_phys + + implicit none + + private + public :: create_dry_air_profile, create_O2_profile, create_O3_profile, & + set_gas_species_values + + !> Conversion factor from km to cm + real(kind_phys), parameter, public :: km_to_cm = 1.0e5 + !> Conversion factor from m3 to cm3 + real(kind_phys), parameter, public :: m_3_to_cm_3 = 1.0e6 + +contains + + !> Creates TUV-x dry air profile + function create_dry_air_profile(height_grid, errmsg, errcode) & + result(profile) + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_tuvx_profile, only: profile_t + use musica_util, only: error_t + use musica_ccpp_species, only: tuvx_species_set + use musica_ccpp_tuvx_load_species, only: index_dry_air + + type(grid_t), intent(in) :: height_grid + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + type(profile_t), pointer :: profile + + ! local variables + type(error_t) :: error + + profile => profile_t( tuvx_species_set(index_dry_air)%name, & + tuvx_species_set(index_dry_air)%unit, & + height_grid, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end function create_dry_air_profile + + !> Creates TUV-x O2 profile + function create_O2_profile(height_grid, errmsg, errcode) & + result(profile) + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_tuvx_profile, only: profile_t + use musica_util, only: error_t + use musica_ccpp_species, only: tuvx_species_set + use musica_ccpp_tuvx_load_species, only: index_O2 + + type(grid_t), intent(in) :: height_grid + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + type(profile_t), pointer :: profile + + ! local variables + type(error_t) :: error + + profile => profile_t( tuvx_species_set(index_O2)%name, & + tuvx_species_set(index_O2)%unit, & + height_grid, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end function create_O2_profile + + !> Creates TUV-x O3 profile + function create_O3_profile(height_grid, errmsg, errcode) & + result(profile) + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_tuvx_profile, only: profile_t + use musica_util, only: error_t + use musica_ccpp_species, only: tuvx_species_set + use musica_ccpp_tuvx_load_species, only: index_O3 + + type(grid_t), intent(in) :: height_grid + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + type(profile_t), pointer :: profile + + ! local variables + type(error_t) :: error + + profile => profile_t( tuvx_species_set(index_O3)%name, & + tuvx_species_set(index_O3)%unit, & + height_grid, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end function create_O3_profile + + !> Sets the species constituents in the vertical layer + subroutine set_gas_species_values(profile, dry_air_density, constituents, & + height_deltas, index_species, errmsg, errcode) + use musica_ccpp_util, only: has_error_occurred + use musica_ccpp_species, only: tuvx_species_set + use musica_ccpp_tuvx_load_species, only: O3_LABEL + use musica_tuvx_profile, only: profile_t + use musica_util, only: error_t + + type(profile_t), intent(inout) :: profile + real(kind_phys), intent(in) :: dry_air_density(:) + real(kind_phys), intent(in) :: constituents(:) ! kg kg-1 + real(kind_phys), intent(in) :: height_deltas(:) ! km, change in height in each vertical layer + integer, intent(in) :: index_species + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variables + type(error_t) :: error + integer :: num_vertical_levels + real(kind_phys) :: constituent_mol_per_cm_3(size(constituents)) ! mol cm-3 + real(kind_phys) :: interfaces(size(constituents) + 2) + real(kind_phys) :: densities(size(constituents) + 1) + real(kind_phys) :: molar_mass + + molar_mass = tuvx_species_set(index_species)%molar_mass + constituent_mol_per_cm_3(:) = constituents(:) * dry_air_density(:) / molar_mass / m_3_to_cm_3 + + num_vertical_levels = size(constituents) + interfaces(1) = constituent_mol_per_cm_3(num_vertical_levels) + interfaces(2:num_vertical_levels+1) = constituent_mol_per_cm_3(num_vertical_levels:1:-1) + interfaces(num_vertical_levels+2) = constituent_mol_per_cm_3(1) + + if (tuvx_species_set(index_species)%name == O3_LABEL) then + densities(:) = height_deltas(:) * km_to_cm * 0.5_kind_phys & + * ( interfaces(1:num_vertical_levels+1) + interfaces(2:num_vertical_levels+2) ) + else + densities(:) = height_deltas(:) * km_to_cm & + * sqrt(interfaces(1:num_vertical_levels+1)) * sqrt(interfaces(2:num_vertical_levels+2)) + end if + + call profile%set_edge_values( interfaces, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + call profile%set_layer_densities( densities, error) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + call profile%calculate_exo_layer_density( & + tuvx_species_set(index_species)%scale_height, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end subroutine set_gas_species_values + +end module musica_ccpp_tuvx_gas_species \ No newline at end of file diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 index cff0d2a2..2fc38b97 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90 @@ -83,8 +83,9 @@ function create_height_grid(vertical_layer_dimension, vertical_interface_dimensi end function create_height_grid !> Sets TUV-x height grid values from the host-model height grid + ! and calculates the change in height for each vertical layer subroutine set_height_grid_values(height_grid, host_midpoints, & - host_interfaces, errmsg, errcode) + host_interfaces, height_deltas, errmsg, errcode) use ccpp_kinds, only: kind_phys use musica_ccpp_util, only: has_error_occurred use musica_tuvx_grid, only: grid_t @@ -93,6 +94,7 @@ subroutine set_height_grid_values(height_grid, host_midpoints, & type(grid_t), intent(inout) :: height_grid real(kind_phys), intent(in) :: host_midpoints(:) ! km real(kind_phys), intent(in) :: host_interfaces(:) ! km + real(kind_phys), intent(inout) :: height_deltas(:) ! km character(len=*), intent(out) :: errmsg integer, intent(out) :: errcode @@ -117,6 +119,12 @@ subroutine set_height_grid_values(height_grid, host_midpoints, & n_host_midpoints = size(host_midpoints) n_host_interfaces = size(host_interfaces) + if ( size(height_deltas) /= n_host_interfaces ) then + errmsg = "[MUSICA Error] Invalid size of TUV-x height deltas." + errcode = 1 + return + end if + interfaces(1) = host_interfaces(n_host_interfaces) interfaces(2:n_host_interfaces) = host_midpoints(n_host_midpoints:1:-1) interfaces(n_host_interfaces+1) = host_interfaces(1) @@ -127,6 +135,8 @@ subroutine set_height_grid_values(height_grid, host_midpoints, & midpoints(n_host_midpoints+1) = 0.5_kind_phys * & ( interfaces(n_host_interfaces) + interfaces(n_host_interfaces+1) ) + height_deltas(:) = interfaces(2:size(interfaces)) - interfaces(1:size(interfaces)-1) + call height_grid%set_edges( interfaces, error ) if ( has_error_occurred( error, errmsg, errcode ) ) return diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 new file mode 100644 index 00000000..916e5e97 --- /dev/null +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 @@ -0,0 +1,219 @@ +module musica_ccpp_tuvx_load_species + use ccpp_kinds, only: kind_phys + use musica_ccpp_species, only: MUSICA_INT_UNASSIGNED + + implicit none + private + + public :: configure_tuvx_species, check_tuvx_species_initialization + + integer, protected, public :: index_cloud_liquid_water_content = MUSICA_INT_UNASSIGNED + integer, protected, public :: index_dry_air = MUSICA_INT_UNASSIGNED + integer, protected, public :: index_O2 = MUSICA_INT_UNASSIGNED + integer, protected, public :: index_O3 = MUSICA_INT_UNASSIGNED + + ! Constants + ! Cloud liquid water + character(len=*), parameter, public :: CLOUD_LIQUID_WATER_CONTENT_LABEL = & + 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' + character(len=*), parameter, public :: CLOUD_LIQUID_WATER_CONTENT_LONG_NAME = & + 'cloud water mass mixing ratio with respect to moist air plus all airborne condensates' + character(len=*), parameter, public :: CLOUD_LIQUID_WATER_CONTENT_UNITS = 'kg kg-1' + real(kind_phys), parameter, public :: CLOUD_LIQUID_WATER_CONTENT_MOLAR_MASS = 0.018_kind_phys ! kg mol-1 + ! Gas species - dry air, O2, O3 + character(len=*), parameter, public :: DRY_AIR_LABEL = 'air' + character(len=*), parameter, public :: O2_LABEL = 'O2' + character(len=*), parameter, public :: O3_LABEL = 'O3' + character(len=*), parameter, public :: TUVX_GAS_SPECIES_UNITS = 'molecule cm-3' + real(kind_phys), parameter, public :: SCALE_HEIGHT_DRY_AIR = 8.01_kind_phys ! km + real(kind_phys), parameter, public :: SCALE_HEIGHT_O2 = 7.0_kind_phys ! km + real(kind_phys), parameter, public :: SCALE_HEIGHT_O3 = 7.0_kind_phys ! km + !> Molar mass value of dry air is obtained from 'CAM-SIMA/src/utils/std_atm_profile.F90' + real(kind_phys), parameter, public :: MOLAR_MASS_DRY_AIR = 0.0289644_kind_phys ! kg mol-1 + real(kind_phys), parameter, public :: MOLAR_MASS_O2 = 0.0319988_kind_phys ! kg mol-1 + real(kind_phys), parameter, public :: MOLAR_MASS_O3 = 0.0479982_kind_phys ! kg mol-1 + +contains + + !> Configures the TUV-x species and their constituent properties. + ! If the MICM configuration includes any TUV-x gas species, constituent properties + ! are not created; otherwise, new constituent properties are generated for each species. + subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props, & + errmsg, errcode) + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_species, only: musica_species_t + use musica_util, only: error_t + + type(musica_species_t), intent(inout) :: micm_species(:) + type(musica_species_t), allocatable, intent(out) :: tuvx_species(:) + type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variables + integer :: num_new_species = 4 + integer :: num_micm_species + ! temp_constituents_props is used to store TUVx-specific constituents and gas species + ! that are not registered by MICM. Its fixed array size represents the maximum number + ! of possible constituents. + type(ccpp_constituent_properties_t) :: temp_constituent_props(4) + logical :: is_dry_air_registered = .false. + logical :: is_O2_registered = .false. + logical :: is_O3_registered = .false. + integer :: i_new, i_species, i_tuvx_species + + num_micm_species = size(micm_species) + is_dry_air_registered = .false. + is_O2_registered = .false. + is_O3_registered = .false. + + ! Register cloud liquid water content needed for cloud optics calculations + i_new = 1 + call temp_constituent_props(i_new)%instantiate( & + std_name = CLOUD_LIQUID_WATER_CONTENT_LABEL, & + long_name = CLOUD_LIQUID_WATER_CONTENT_LONG_NAME, & + units = CLOUD_LIQUID_WATER_CONTENT_UNITS, & + vertical_dim = "vertical_layer_dimension", & + default_value = 0.0_kind_phys, & + min_value = 0.0_kind_phys, & + molar_mass = CLOUD_LIQUID_WATER_CONTENT_MOLAR_MASS, & + advected = .true., & + errcode = errcode, & + errmsg = errmsg ) + if (errcode /= 0) return + + ! Iterate through the MICM species to check if any TUV-x gas + ! species are included; if present, updates the scale height and profiled status. + do i_species = 1, num_micm_species + if (is_dry_air_registered .and. is_O2_registered .and. is_O3_registered) exit + + if ( micm_species(i_species)%name == DRY_AIR_LABEL ) then + is_dry_air_registered = .true. + micm_species(i_species)%profiled = .true. + micm_species(i_species)%scale_height = SCALE_HEIGHT_DRY_AIR + else if ( micm_species(i_species)%name == O2_LABEL ) then + is_O2_registered = .true. + micm_species(i_species)%profiled = .true. + micm_species(i_species)%scale_height = SCALE_HEIGHT_O2 + else if ( micm_species(i_species)%name == O3_LABEL ) then + is_O3_registered = .true. + micm_species(i_species)%profiled = .true. + micm_species(i_species)%scale_height = SCALE_HEIGHT_O3 + end if + end do + + if (.not. is_dry_air_registered) then + i_new = i_new + 1 + call temp_constituent_props(i_new)%instantiate( & + std_name = DRY_AIR_LABEL, & + long_name = DRY_AIR_LABEL, & + units = 'kg kg-1', & + vertical_dim = "vertical_layer_dimension", & + default_value = 0.0_kind_phys, & + min_value = 0.0_kind_phys, & + molar_mass = MOLAR_MASS_DRY_AIR, & + advected = .false., & + errcode = errcode, & + errmsg = errmsg ) + if (errcode /= 0) return + end if + + if (.not. is_O2_registered) then + i_new = i_new + 1 + call temp_constituent_props(i_new)%instantiate( & + std_name = O2_LABEL, & + long_name = O2_LABEL, & + units = 'kg kg-1', & + vertical_dim = "vertical_layer_dimension", & + default_value = 0.0_kind_phys, & + min_value = 0.0_kind_phys, & + molar_mass = MOLAR_MASS_O2, & + advected = .false., & + errcode = errcode, & + errmsg = errmsg ) + if (errcode /= 0) return + end if + + if (.not. is_O3_registered) then + i_new = i_new + 1 + call temp_constituent_props(i_new)%instantiate( & + std_name = O3_LABEL, & + long_name = O3_LABEL, & + units = 'kg kg-1', & + vertical_dim = "vertical_layer_dimension", & + default_value = 0.0_kind_phys, & + min_value = 0.0_kind_phys, & + molar_mass = MOLAR_MASS_O3, & + advected = .false., & + errcode = errcode, & + errmsg = errmsg ) + if (errcode /= 0) return + end if + + allocate( constituent_props(i_new) ) + constituent_props(:) = temp_constituent_props(1:i_new) + + allocate( tuvx_species(num_new_species) ) + i_tuvx_species = 1 + index_cloud_liquid_water_content = i_tuvx_species + tuvx_species(i_tuvx_species) = musica_species_t( & + name = CLOUD_LIQUID_WATER_CONTENT_LABEL, & + unit = CLOUD_LIQUID_WATER_CONTENT_UNITS, & + molar_mass = CLOUD_LIQUID_WATER_CONTENT_MOLAR_MASS, & + index_musica_species = i_tuvx_species ) + + i_tuvx_species = i_tuvx_species + 1 + index_dry_air = i_tuvx_species + tuvx_species(i_tuvx_species) = musica_species_t( & + name = DRY_AIR_LABEL, & + unit = TUVX_GAS_SPECIES_UNITS, & ! TUV-x profile unit, different from molar mass unit + molar_mass = MOLAR_MASS_DRY_AIR, & ! kg mol-1 + index_musica_species = i_tuvx_species, & + profiled = .true., & + scale_height = SCALE_HEIGHT_DRY_AIR ) + + i_tuvx_species = i_tuvx_species + 1 + index_O2 = i_tuvx_species + tuvx_species(i_tuvx_species) = musica_species_t( & + name = O2_LABEL, & + unit = TUVX_GAS_SPECIES_UNITS, & ! TUV-x profile unit, different from molar mass unit + molar_mass = MOLAR_MASS_O2, & ! kg mol-1 + index_musica_species = i_tuvx_species, & + profiled = .true., & + scale_height = SCALE_HEIGHT_O2 ) + + i_tuvx_species = i_tuvx_species + 1 + index_O3 = i_tuvx_species + tuvx_species(i_tuvx_species) = musica_species_t( & + name = O3_LABEL, & + unit = TUVX_GAS_SPECIES_UNITS, & ! TUV-x profile unit, different from molar mass unit + molar_mass = MOLAR_MASS_O3, & ! kg mol-1 + index_musica_species = i_tuvx_species, & + profiled = .true., & + scale_height = SCALE_HEIGHT_O3 ) + + end subroutine configure_tuvx_species + + !> Ensures that the indices of all TUV-x species are initialized. + ! This function is typically called during the initialization phase, + ! so that the indices can be used during the run phase without the need + ! for additional checks. + subroutine check_tuvx_species_initialization(errmsg, errcode) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + errmsg = '' + errcode = 0 + + if ((index_cloud_liquid_water_content == MUSICA_INT_UNASSIGNED) .or. & + (index_dry_air == MUSICA_INT_UNASSIGNED) .or. & + (index_O2 == MUSICA_INT_UNASSIGNED) .or. & + (index_O3== MUSICA_INT_UNASSIGNED)) then + errmsg = "[MUSICA Error] TUV-x species index (or indices) has not been initialized." + errcode = 1 + return + end if + + end subroutine check_tuvx_species_initialization + +end module musica_ccpp_tuvx_load_species \ No newline at end of file diff --git a/test/docker/Dockerfile.musica b/test/docker/Dockerfile.musica index 00ad7482..dbab331b 100644 --- a/test/docker/Dockerfile.musica +++ b/test/docker/Dockerfile.musica @@ -6,7 +6,7 @@ FROM ubuntu:22.04 ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8 -ARG CAM_SIMA_CHEMISTRY_DATA_TAG=be87bc14822aa50b1afda0059ab6f5b5bd7397e6 +ARG CAM_SIMA_CHEMISTRY_DATA_TAG=ea3539f1d7b71162e8a78d900ecbe265ba870e3d ARG BUILD_TYPE=Debug RUN apt update \ diff --git a/test/docker/Dockerfile.musica.no_install b/test/docker/Dockerfile.musica.no_install index f6440ac1..1a4acd79 100644 --- a/test/docker/Dockerfile.musica.no_install +++ b/test/docker/Dockerfile.musica.no_install @@ -9,7 +9,7 @@ FROM ubuntu:22.04 ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8 -ARG CAM_SIMA_CHEMISTRY_DATA_TAG=be87bc14822aa50b1afda0059ab6f5b5bd7397e6 +ARG CAM_SIMA_CHEMISTRY_DATA_TAG=ea3539f1d7b71162e8a78d900ecbe265ba870e3d ARG BUILD_TYPE=Debug RUN apt update \ diff --git a/test/musica/CMakeLists.txt b/test/musica/CMakeLists.txt index fbedfba3..f060ac59 100644 --- a/test/musica/CMakeLists.txt +++ b/test/musica/CMakeLists.txt @@ -13,10 +13,7 @@ set(MUSICA_ENABLE_INSTALL OFF) FetchContent_MakeAvailable(musica) -# --------------------------------------------------------- -# Create a test for MUSICA CCPP wrapper -# --------------------------------------------------------- - +# MUSICA CCPP API test add_executable(test_musica_api test_musica_api.F90 musica_ccpp_namelist.F90) file(GLOB MUSICA_CCPP_SOURCES @@ -87,4 +84,36 @@ add_test( COMMAND ${Python_EXECUTABLE} ${CMAKE_BINARY_DIR}/../$ENV{CCPP_STD_NAMES_PATH}/tools/meta_stdname_check.py --metafile-loc ${CMAKE_BINARY_DIR}/metadata_test/musica_ccpp.meta --stdname-dict ${CMAKE_BINARY_DIR}/../$ENV{CCPP_STD_NAMES_PATH}/standard_names.xml -) \ No newline at end of file +) + +# MUSICA species test +add_executable(test_musica_species test_musica_species.F90 musica_ccpp_namelist.F90) + +target_sources(test_musica_species + PUBLIC + ${MUSICA_CCPP_SOURCES} + ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_const_utils.F90 + ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_tuvx_utils.F90 + ${CCPP_SRC_PATH}/ccpp_constituent_prop_mod.F90 + ${CCPP_SRC_PATH}/ccpp_hash_table.F90 + ${CCPP_SRC_PATH}/ccpp_hashable.F90 + ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 +) + +target_link_libraries(test_musica_species + PRIVATE + musica::musica-fortran +) + +set_target_properties(test_musica_species + PROPERTIES + LINKER_LANGUAGE Fortran +) + +add_test( + NAME test_musica_species + COMMAND $ + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} +) + +add_memory_check_test(test_musica_species $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) \ No newline at end of file diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index 6ed8b0e8..76af8099 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -10,8 +10,13 @@ program run_test_musica_ccpp real(kind_phys), parameter :: DEGREE_TO_RADIAN = 3.14159265358979323846_kind_phys / 180.0_kind_phys + write(*,*) "[MUSICA Test] Running the Chapman test" call test_chapman() + write(*,*) "[MUSICA Test] Ends the Chapman test" + + write(*,*) "[MUSICA Test] Running the Terminator test" call test_terminator() + write(*,*) "[MUSICA Test] Ends the Terminator test" contains @@ -136,17 +141,19 @@ end subroutine get_wavelength_edges !> Tests the Chapman chemistry scheme subroutine test_chapman() - 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 + 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 + use musica_ccpp_tuvx_load_species, only: index_dry_air, index_O2, index_O3 implicit none - integer, parameter :: NUM_SPECIES = 5 - integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 + integer, parameter :: NUM_SPECIES = 5 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 + integer, parameter :: NUM_TUVX_ONLY_GAS_SPECIES = 1 ! 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. @@ -173,9 +180,9 @@ subroutine test_chapman() real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: cloud_area_fraction ! unitless real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: air_pressure_thickness ! Pa real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & - NUM_SPECIES+NUM_TUVX_CONSTITUENTS) :: constituents ! kg kg-1 + NUM_SPECIES+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) :: constituents ! kg kg-1 real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & - NUM_SPECIES+NUM_TUVX_CONSTITUENTS) :: initial_constituents ! kg kg-1 + NUM_SPECIES+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) :: initial_constituents ! kg kg-1 real(kind_phys), dimension(NUM_COLUMNS) :: solar_zenith_angle ! radians real(kind_phys) :: earth_sun_distance ! AU type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:) @@ -228,7 +235,7 @@ subroutine test_chapman() stop 3 endif ASSERT(allocated(constituent_props)) - ASSERT(size(constituent_props) == NUM_SPECIES+NUM_TUVX_CONSTITUENTS) + ASSERT(size(constituent_props) == NUM_SPECIES+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) do i = 1, size(constituent_props) ASSERT(constituent_props(i)%is_instantiated(errcode, errmsg)) ASSERT(errcode == 0) @@ -244,7 +251,8 @@ subroutine test_chapman() (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) .or. & (trim(species_name) == "cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water" .and. & - molar_mass == 0.018_kind_phys .and. is_advected) + molar_mass == 0.018_kind_phys .and. is_advected) .or. & + (trim(species_name) == "air" .and. molar_mass == 0.0289644_kind_phys .and. .not. is_advected) ASSERT(tmp_bool) call constituent_props(i)%units(units, errcode, errmsg) if (errcode /= 0) then @@ -304,6 +312,11 @@ subroutine test_chapman() constituents(j,k,NUM_SPECIES+1) = 1.0e-3_kind_phys * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) end do end do + do j = 1, NUM_COLUMNS + do k = 1, NUM_LAYERS + constituents(j,k,NUM_SPECIES+index_dry_air) = 1.0e-3_kind_phys * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) + end do + end do initial_constituents(:,:,:) = constituents(:,:,:) write(*,*) "[MUSICA INFO] Initial Time Step" @@ -371,17 +384,18 @@ end subroutine test_chapman !> Tests the simple Terminator chemistry scheme subroutine test_terminator() - 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 - + 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 + use musica_ccpp_tuvx_load_species, only: index_dry_air, index_O2, index_O3 implicit none - integer, parameter :: NUM_SPECIES = 2 - integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 + integer, parameter :: NUM_SPECIES = 2 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 + integer, parameter :: NUM_TUVX_ONLY_GAS_SPECIES = 3 ! 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. @@ -408,9 +422,9 @@ subroutine test_terminator() real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: cloud_area_fraction ! unitless real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: air_pressure_thickness ! Pa real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & - NUM_SPECIES+NUM_TUVX_CONSTITUENTS) :: constituents ! kg kg-1 + NUM_SPECIES+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) :: constituents ! kg kg-1 real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & - NUM_SPECIES+NUM_TUVX_CONSTITUENTS) :: initial_constituents ! kg kg-1 + NUM_SPECIES+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) :: initial_constituents ! kg kg-1 real(kind_phys), dimension(NUM_COLUMNS) :: solar_zenith_angle ! radians real(kind_phys) :: earth_sun_distance ! AU type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:) @@ -463,7 +477,7 @@ subroutine test_terminator() stop 3 endif ASSERT(allocated(constituent_props)) - ASSERT(size(constituent_props) == NUM_SPECIES+NUM_TUVX_CONSTITUENTS) + ASSERT(size(constituent_props) == NUM_SPECIES+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) do i = 1, size(constituent_props) ASSERT(constituent_props(i)%is_instantiated(errcode, errmsg)) ASSERT(errcode == 0) @@ -475,8 +489,11 @@ subroutine test_terminator() 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) .or. & - (trim(species_name) == "cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water" & - .and. molar_mass == 0.018_kind_phys .and. is_advected) + (trim(species_name) == "cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water" & + .and. molar_mass == 0.018_kind_phys .and. is_advected) .or. & + (trim(species_name) == "air" .and. molar_mass == 0.0289644_kind_phys .and. .not. is_advected) .or. & + (trim(species_name) == "O2" .and. molar_mass == 0.0319988_kind_phys .and. .not. is_advected) .or. & + (trim(species_name) == "O3" .and. molar_mass == 0.0479982_kind_phys .and. .not. is_advected) ASSERT(tmp_bool) call constituent_props(i)%units(units, errcode, errmsg) if (errcode /= 0) then @@ -527,6 +544,14 @@ subroutine test_terminator() constituents(j,k,NUM_SPECIES+1) = 1.0e-3_kind_phys * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) end do end do + ! set concentrations for TUV-x gas species + do j = 1, NUM_COLUMNS + do k = 1, NUM_LAYERS + constituents(j,k,NUM_SPECIES+index_dry_air) = 1.0e-3_kind_phys * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) + constituents(j,k,NUM_SPECIES+index_O2) = 1.0e-3_kind_phys * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) + constituents(j,k,NUM_SPECIES+index_O3) = 1.0e-3_kind_phys * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) + end do + end do initial_constituents(:,:,:) = constituents(:,:,:) write(*,*) "[MUSICA INFO] Initial Time Step" diff --git a/test/musica/test_musica_species.F90 b/test/musica/test_musica_species.F90 new file mode 100644 index 00000000..579031de --- /dev/null +++ b/test/musica/test_musica_species.F90 @@ -0,0 +1,217 @@ +program test_musica_species + + use musica_ccpp_species + + 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_register_musica_species() + call test_initialize_musica_species_indices_and_molar_mass() + call test_extract_and_update_subset_constituents() + +contains + + subroutine test_register_musica_species() + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_species, only: musica_species_t + use musica_ccpp_tuvx_load_species, only: configure_tuvx_species, & + check_tuvx_species_initialization + + integer, parameter :: NUM_MICM_SPECIES = 6 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 4 + type(musica_species_t) :: micm_species(NUM_MICM_SPECIES) + type(musica_species_t), allocatable :: tuvx_species(:) + type(ccpp_constituent_properties_t) :: micm_constituent_props(NUM_MICM_SPECIES) + type(ccpp_constituent_properties_t), allocatable :: tuvx_constituent_props(:) + character(len=512) :: errmsg + integer :: errcode + real(kind_phys) :: molar_mass_group(NUM_MICM_SPECIES) = & + [0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys, 0.4_kind_phys, 0.5_kind_phys, 0.6_kind_phys] + integer :: i_species + character(len=512) :: species_names(NUM_MICM_SPECIES) + + species_names(1) = 'N2' + species_names(2) = 'O2' ! shared species + species_names(3) = 'FOO' + species_names(4) = 'O1D' + species_names(5) = 'BAZ' + species_names(6) = 'O3' ! shared species + + do i_species = 1, NUM_MICM_SPECIES + call micm_constituent_props(i_species)%instantiate( & + std_name = trim(species_names(i_species)), & + long_name = trim(species_names(i_species)), & + units = 'kg kg-1', & + vertical_dim = 'vertical_layer_dimension', & + default_value = 0.0_kind_phys, & + min_value = 0.0_kind_phys, & + molar_mass = molar_mass_group(i_species), & + advected = .true., & + errcode = errcode, & + errmsg = errmsg) + + micm_species(i_species) = musica_species_t( & + name = species_names(i_species), & + unit = 'kg kg-1', & + molar_mass = molar_mass_group(i_species), & + index_musica_species = i_species ) + end do + + call configure_tuvx_species( micm_species, tuvx_species, tuvx_constituent_props, & + errmsg, errcode ) + ASSERT(errcode == 0) + + call register_musica_species( micm_species, tuvx_species ) + ASSERT(allocated(micm_species_set)) + ASSERT(allocated(tuvx_species_set)) + ASSERT(number_of_micm_species == NUM_MICM_SPECIES) + ASSERT(number_of_tuvx_species == NUM_TUVX_CONSTITUENTS) + + call check_tuvx_species_initialization( errmsg, errcode ) + ASSERT(errcode == 0) + + do i_species = 1, size(micm_species) + call micm_species(i_species)%deallocate() + end do + do i_species = 1, size(tuvx_species) + call tuvx_species(i_species)%deallocate() + end do + deallocate(tuvx_species) + deallocate(tuvx_constituent_props) + + call cleanup_musica_species() + + end subroutine test_register_musica_species + + subroutine test_initialize_musica_species_indices_and_molar_mass() + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t, ccpp_constituent_prop_ptr_t + use musica_ccpp, only: musica_ccpp_register + use musica_ccpp_namelist, only: filename_of_micm_configuration, & + filename_of_tuvx_configuration, & + filename_of_tuvx_micm_mapping_configuration + + type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:) + type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:) + type(ccpp_constituent_properties_t), pointer :: const_prop + character(len=512) :: errmsg + integer :: errcode + integer :: i + + 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( constituent_props, errmsg, errcode ) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + ASSERT(allocated(constituent_props)) + + 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 initialize_musica_species_indices( constituent_props_ptr, errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(allocated(micm_indices_constituent_props)) + ASSERT(allocated(tuvx_indices_constituent_props)) + + call initialize_molar_mass_array( constituent_props_ptr, errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(allocated(micm_molar_mass_array)) + + do i = 1, size(micm_species_set) + ASSERT(micm_species_set(i)%index_musica_species == i) + ASSERT(micm_species_set(i)%index_constituent_props == micm_indices_constituent_props(i)) + ASSERT(micm_species_set(i)%molar_mass == micm_molar_mass_array(i)) + end do + + do i = 1, size(tuvx_species_set) + ASSERT(tuvx_species_set(i)%index_musica_species == i) + ASSERT(tuvx_species_set(i)%index_constituent_props == tuvx_indices_constituent_props(i)) + end do + + call check_initialization( errmsg, errcode ) + ASSERT(errcode == 0) + + call cleanup_musica_species() + + end subroutine test_initialize_musica_species_indices_and_molar_mass + + + subroutine test_extract_and_update_subset_constituents() + use ccpp_kinds, only: kind_phys + + integer, parameter :: NUM_COLUMNS = 3 + integer, parameter :: NUM_LAYERS = 3 + integer, parameter :: NUM_SPECIES = 5 + integer, parameter :: NUM_SUBSET_SPECIES = 3 + integer :: indices_array(NUM_SUBSET_SPECIES) = [2, 4, 5] + real(kind_phys) :: constituents(NUM_COLUMNS, NUM_LAYERS, NUM_SPECIES) + real(kind_phys) :: subset_constituents(NUM_COLUMNS, NUM_LAYERS, NUM_SUBSET_SPECIES) + real(kind_phys) :: expected_constituents(NUM_COLUMNS, NUM_LAYERS, NUM_SPECIES) + real(kind_phys) :: expected_subset_constituents(NUM_COLUMNS, NUM_LAYERS, NUM_SUBSET_SPECIES) + character(len=512) :: errmsg + integer :: errcode + integer :: i, j, k + + constituents = reshape([1.0, 2.0, 3.0, 4.0, 5.0, & + 6.0, 7.0, 8.0, 9.0, 10.0, & + 11.0, 12.0, 13.0, 14.0, 15.0, & + 16.0, 17.0, 18.0, 19.0, 20.0, & + 21.0, 22.0, 23.0, 24.0, 25.0, & + 26.0, 27.0, 28.0, 29.0, 30.0, & + 31.0, 32.0, 33.0, 34.0, 35.0, & + 36.0, 37.0, 38.0, 39.0, 40.0, & + 41.0, 42.0, 43.0, 44.0, 45.0, & + 46.0, 47.0, 48.0, 49.0, 50.0], shape(constituents)) + expected_constituents(:,:,:) = constituents(:,:,:) + do k = 1, NUM_SUBSET_SPECIES + do j = 1, NUM_LAYERS + do i = 1, NUM_COLUMNS + expected_subset_constituents(i,j,k) = constituents(i,j,indices_array(k)) + end do + end do + end do + + call extract_subset_constituents( indices_array, constituents, & + subset_constituents, errmsg, errcode ) + ASSERT(errcode == 0) + do k = 1, NUM_SUBSET_SPECIES + do j = 1, NUM_LAYERS + do i = 1, NUM_COLUMNS + ASSERT(expected_subset_constituents(i,j,k) == subset_constituents(i,j,k)) + end do + end do + end do + + do k = 1, NUM_SUBSET_SPECIES + do j = 1, NUM_LAYERS + do i = 1, NUM_COLUMNS + subset_constituents(i,j,k) = subset_constituents(i,j,k) + 100.0_kind_phys + expected_constituents(i,j,indices_array(k)) = & + expected_constituents(i,j,indices_array(k)) + 100.0_kind_phys + end do + end do + end do + + call update_constituents( indices_array, subset_constituents, & + constituents, errmsg, errcode ) + ASSERT(errcode == 0) + do k = 1, NUM_SPECIES + do j = 1, NUM_LAYERS + do i = 1, NUM_COLUMNS + ASSERT(expected_constituents(i,j,k) == constituents(i,j,k)) + end do + end do + end do + + end subroutine test_extract_and_update_subset_constituents + +end program test_musica_species \ No newline at end of file diff --git a/test/musica/tuvx/CMakeLists.txt b/test/musica/tuvx/CMakeLists.txt index a636ec8f..c12d0d57 100644 --- a/test/musica/tuvx/CMakeLists.txt +++ b/test/musica/tuvx/CMakeLists.txt @@ -202,3 +202,70 @@ add_test( ) add_memory_check_test(test_tuvx_aerosol_optics $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) + +# TUV-x gas species profiles +add_executable(test_tuvx_gas_species test_tuvx_gas_species.F90) + +target_sources(test_tuvx_gas_species + PUBLIC + ${MUSICA_CCPP_SOURCES} + ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_const_utils.F90 + ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_tuvx_utils.F90 + ${CCPP_SRC_PATH}/ccpp_constituent_prop_mod.F90 + ${CCPP_SRC_PATH}/ccpp_hash_table.F90 + ${CCPP_SRC_PATH}/ccpp_hashable.F90 + ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 + ${CMAKE_CURRENT_SOURCE_DIR}/../musica_ccpp_namelist.F90 +) + +target_link_libraries(test_tuvx_gas_species + PRIVATE + musica::musica-fortran +) + +set_target_properties(test_tuvx_gas_species + PROPERTIES + LINKER_LANGUAGE Fortran +) + +add_test( + NAME test_tuvx_gas_species + COMMAND $ + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} +) + +add_memory_check_test(test_tuvx_gas_species $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) + +# Configure TUV-x species +add_executable(test_tuvx_load_species test_tuvx_load_species.F90) + +target_sources(test_tuvx_load_species + PUBLIC + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_load_species.F90 + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_gas_species.F90 + ${MUSICA_SRC_PATH}/musica_ccpp_species.F90 + ${MUSICA_SRC_PATH}/musica_ccpp_util.F90 + ${TO_BE_CCPPIZED_SRC_PATH}/ccpp_const_utils.F90 + ${CCPP_SRC_PATH}/ccpp_constituent_prop_mod.F90 + ${CCPP_SRC_PATH}/ccpp_hash_table.F90 + ${CCPP_SRC_PATH}/ccpp_hashable.F90 + ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 +) + +target_link_libraries(test_tuvx_load_species + PRIVATE + musica::musica-fortran +) + +set_target_properties(test_tuvx_load_species + PROPERTIES + LINKER_LANGUAGE Fortran +) + +add_test( + NAME test_tuvx_load_species + COMMAND $ + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} +) + +add_memory_check_test(test_tuvx_load_species $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) \ No newline at end of file diff --git a/test/musica/tuvx/test_tuvx_gas_species.F90 b/test/musica/tuvx/test_tuvx_gas_species.F90 new file mode 100644 index 00000000..c226f6c8 --- /dev/null +++ b/test/musica/tuvx/test_tuvx_gas_species.F90 @@ -0,0 +1,476 @@ +program test_tuvx_gas_species_profile + + use musica_ccpp_tuvx_gas_species + + 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_create_gas_species_profile() + call test_initialize_tuvx_species() + +contains + + !> 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 + + real(kind_phys), dimension(:), intent(out) :: edges + + edges = (/ & + 120.0e-9_kind_phys, & + 121.4e-9_kind_phys, & + 121.9e-9_kind_phys, & + 123.5e-9_kind_phys, & + 124.3e-9_kind_phys, & + 125.5e-9_kind_phys, & + 126.3e-9_kind_phys, & + 127.1e-9_kind_phys, & + 130.1e-9_kind_phys, & + 131.1e-9_kind_phys, & + 135.0e-9_kind_phys, & + 140.0e-9_kind_phys, & + 145.0e-9_kind_phys, & + 150.0e-9_kind_phys, & + 155.0e-9_kind_phys, & + 160.0e-9_kind_phys, & + 165.0e-9_kind_phys, & + 168.0e-9_kind_phys, & + 171.0e-9_kind_phys, & + 173.0e-9_kind_phys, & + 174.4e-9_kind_phys, & + 175.4e-9_kind_phys, & + 177.0e-9_kind_phys, & + 178.6e-9_kind_phys, & + 180.2e-9_kind_phys, & + 181.8e-9_kind_phys, & + 183.5e-9_kind_phys, & + 185.2e-9_kind_phys, & + 186.9e-9_kind_phys, & + 188.7e-9_kind_phys, & + 190.5e-9_kind_phys, & + 192.3e-9_kind_phys, & + 194.2e-9_kind_phys, & + 196.1e-9_kind_phys, & + 198.0e-9_kind_phys, & + 200.0e-9_kind_phys, & + 202.0e-9_kind_phys, & + 204.1e-9_kind_phys, & + 206.2e-9_kind_phys, & + 208.0e-9_kind_phys, & + 211.0e-9_kind_phys, & + 214.0e-9_kind_phys, & + 217.0e-9_kind_phys, & + 220.0e-9_kind_phys, & + 223.0e-9_kind_phys, & + 226.0e-9_kind_phys, & + 229.0e-9_kind_phys, & + 232.0e-9_kind_phys, & + 235.0e-9_kind_phys, & + 238.0e-9_kind_phys, & + 241.0e-9_kind_phys, & + 244.0e-9_kind_phys, & + 247.0e-9_kind_phys, & + 250.0e-9_kind_phys, & + 253.0e-9_kind_phys, & + 256.0e-9_kind_phys, & + 259.0e-9_kind_phys, & + 263.0e-9_kind_phys, & + 267.0e-9_kind_phys, & + 271.0e-9_kind_phys, & + 275.0e-9_kind_phys, & + 279.0e-9_kind_phys, & + 283.0e-9_kind_phys, & + 287.0e-9_kind_phys, & + 291.0e-9_kind_phys, & + 295.0e-9_kind_phys, & + 298.5e-9_kind_phys, & + 302.5e-9_kind_phys, & + 305.5e-9_kind_phys, & + 308.5e-9_kind_phys, & + 311.5e-9_kind_phys, & + 314.5e-9_kind_phys, & + 317.5e-9_kind_phys, & + 322.5e-9_kind_phys, & + 327.5e-9_kind_phys, & + 332.5e-9_kind_phys, & + 337.5e-9_kind_phys, & + 342.5e-9_kind_phys, & + 347.5e-9_kind_phys, & + 350.0e-9_kind_phys, & + 355.0e-9_kind_phys, & + 360.0e-9_kind_phys, & + 365.0e-9_kind_phys, & + 370.0e-9_kind_phys, & + 375.0e-9_kind_phys, & + 380.0e-9_kind_phys, & + 385.0e-9_kind_phys, & + 390.0e-9_kind_phys, & + 395.0e-9_kind_phys, & + 400.0e-9_kind_phys, & + 405.0e-9_kind_phys, & + 410.0e-9_kind_phys, & + 415.0e-9_kind_phys, & + 420.0e-9_kind_phys, & + 430.0e-9_kind_phys, & + 440.0e-9_kind_phys, & + 450.0e-9_kind_phys, & + 500.0e-9_kind_phys, & + 550.0e-9_kind_phys, & + 600.0e-9_kind_phys, & + 650.0e-9_kind_phys, & + 700.0e-9_kind_phys, & + 750.0e-9_kind_phys & + /) + + end subroutine get_wavelength_edges + + ! Calculates the expected values for comparison with the answer + subroutine calculate_gas_species_interfaces_and_densities( & + molar_mass, dry_air_density, constituents, height_deltas, & + is_O3, interfaces, densities) + use ccpp_kinds, only: kind_phys + + real(kind_phys), intent(in) :: molar_mass + real(kind_phys), intent(in) :: dry_air_density(:) + real(kind_phys), intent(in) :: constituents(:) + real(kind_phys), intent(in) :: height_deltas(:) + logical, intent(in) :: is_O3 + real(kind_phys), intent(inout) :: interfaces(size(constituents) + 2) + real(kind_phys), intent(inout) :: densities(size(constituents) + 1) + + ! local variables + real(kind_phys) :: constituent_mol_per_cm_3(size(constituents)) ! mol cm-3 + integer :: num_vertical_levels + + constituent_mol_per_cm_3(:) = constituents(:) * dry_air_density(:) / molar_mass / m_3_to_cm_3 + + num_vertical_levels = size(constituents) + interfaces(1) = constituent_mol_per_cm_3(num_vertical_levels) + interfaces(2:num_vertical_levels+1) = constituent_mol_per_cm_3(num_vertical_levels:1:-1) + interfaces(num_vertical_levels+2) = constituent_mol_per_cm_3(1) + + if ( is_O3 ) then + densities(:) = height_deltas(:) * km_to_cm & + * ( interfaces(1:num_vertical_levels+1) & + + interfaces(2:num_vertical_levels+2) ) * 0.5_kind_phys + else + densities(:) = height_deltas(:) * km_to_cm & + * sqrt(interfaces(1:num_vertical_levels+1)) & + * sqrt(interfaces(2:num_vertical_levels+2)) + end if + + end subroutine calculate_gas_species_interfaces_and_densities + + subroutine test_create_gas_species_profile() + 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_tuvx, only: grid_t, profile_t + use musica_util, only: error_t + use musica_ccpp, only: musica_ccpp_register, musica_ccpp_final + use musica_ccpp_tuvx_height_grid, only: create_height_grid + use musica_ccpp_namelist, only: filename_of_micm_configuration, & + filename_of_tuvx_configuration, & + filename_of_tuvx_micm_mapping_configuration + use musica_ccpp_tuvx_load_species, only: index_dry_air, index_O2, index_O3, & + SCALE_HEIGHT_DRY_AIR, SCALE_HEIGHT_O2, SCALE_HEIGHT_O3, & + MOLAR_MASS_DRY_AIR, MOLAR_MASS_O2, MOLAR_MASS_O3 + use musica_ccpp_species, only: tuvx_species_set, MUSICA_INT_UNASSIGNED + + integer, parameter :: NUM_COLUMNS = 2 + integer, parameter :: NUM_LAYERS = 2 + integer, parameter :: NUM_TUVX_SPECIES = 4 + real, parameter :: ABS_ERROR = 10.0 + real, parameter :: m_to_cm = 100.0_kind_phys + type(ccpp_constituent_properties_t), & + allocatable, target :: constituent_props(:) + type(grid_t), pointer :: height_grid => null() + type(profile_t), pointer :: dry_air_profile => null() + type(profile_t), pointer :: O2_profile => null() + type(profile_t), pointer :: O3_profile => null() + real(kind_phys) :: dry_air_density(NUM_COLUMNS,NUM_LAYERS) ! kg m-3 + real(kind_phys) :: constituents(NUM_COLUMNS,NUM_LAYERS,NUM_TUVX_SPECIES) + real(kind_phys) :: height_deltas(NUM_LAYERS+1) ! km + integer :: errcode + character(len=512) :: errmsg + type(error_t) :: error + character(len=50) :: name, unit + real(kind_phys) :: molar_mass ! kg mol-1 + real(kind_phys) :: scale_height ! km + integer :: index_musica + logical :: tmp_bool + integer :: i_elem, i_col, i, j, k + real(kind_phys) :: dry_air_interfaces(NUM_LAYERS+2), expected_dry_air_interfaces(NUM_COLUMNS,NUM_LAYERS+2) + real(kind_phys) :: O2_interfaces(NUM_LAYERS+2), expected_O2_interfaces(NUM_COLUMNS,NUM_LAYERS+2) + real(kind_phys) :: O3_interfaces(NUM_LAYERS+2), expected_O3_interfaces(NUM_COLUMNS,NUM_LAYERS+2) + real(kind_phys) :: dry_air_densities(NUM_LAYERS+1), expected_dry_air_densities(NUM_COLUMNS,NUM_LAYERS+1) + real(kind_phys) :: O2_densities(NUM_LAYERS+1), expected_O2_densities(NUM_COLUMNS,NUM_LAYERS+1) + real(kind_phys) :: O3_densities(NUM_LAYERS+1), expected_O3_densities(NUM_COLUMNS,NUM_LAYERS+1) + real(kind_phys) :: dry_air_exo_layer_density + real(kind_phys) :: expected_dry_air_exo_layer_density = SCALE_HEIGHT_DRY_AIR + real(kind_phys) :: O2_exo_layer_density + real(kind_phys) :: expected_O2_exo_layer_density = SCALE_HEIGHT_O2 + real(kind_phys) :: O3_exo_layer_density + real(kind_phys) :: expected_O3_exo_layer_density = SCALE_HEIGHT_O3 + + 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' + + dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) + dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) + height_deltas(:) = (/ 0.5_kind_phys, 1.5_kind_phys, 1.0_kind_phys /) + + ! Initialize the constituents array with real values between 0 and 1 + do k = 1, NUM_TUVX_SPECIES + do j = 1, NUM_LAYERS + do i = 1, NUM_COLUMNS + constituents(i, j, k) = (i + j + k) * 0.1_kind_phys + end do + end do + end do + + call musica_ccpp_register(constituent_props, errmsg, errcode) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + + do i_col = 1, NUM_COLUMNS + call calculate_gas_species_interfaces_and_densities( & + MOLAR_MASS_DRY_AIR, dry_air_density(i_col,:), constituents(i_col,:,index_dry_air), & + height_deltas, .false., expected_dry_air_interfaces(i_col,:), expected_dry_air_densities(i_col,:) ) + + call calculate_gas_species_interfaces_and_densities( & + MOLAR_MASS_O2, dry_air_density(i_col,:), constituents(i_col,:,index_O2), & + height_deltas, .false., expected_O2_interfaces(i_col,:), expected_O2_densities(i_col,:) ) + + call calculate_gas_species_interfaces_and_densities( & + MOLAR_MASS_O3, dry_air_density(i_col,:), constituents(i_col,:,index_O3), & + height_deltas, .true., expected_O3_interfaces(i_col,:), expected_O3_densities(i_col,:) ) + write(*,*) " $$$$$$$$$ " + write(*,*) "[F] interfaces: ", expected_dry_air_interfaces(i_col,:) + write(*,*) "[F] densities: ", expected_dry_air_densities(i_col,:) + write(*,*) " $$$$$$$$$ " + end do + + height_grid => create_height_grid( NUM_LAYERS, NUM_LAYERS + 1 , errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(associated(height_grid)) + + dry_air_profile => create_dry_air_profile( height_grid, errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(associated(dry_air_profile)) + + O2_profile => create_O2_profile( height_grid, errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(associated(O2_profile)) + + O3_profile => create_O3_profile( height_grid, errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(associated(O3_profile)) + + do i_col = 1, NUM_COLUMNS + call set_gas_species_values( dry_air_profile, dry_air_density(i_col,:), & + constituents(i_col,:,index_dry_air), height_deltas, index_dry_air, & + errmsg, errcode ) + ASSERT(errcode == 0) + + call set_gas_species_values( O2_profile, dry_air_density(i_col,:), & + constituents(i_col,:,index_O2), height_deltas, index_O2, & + errmsg, errcode ) + ASSERT(errcode == 0) + + call set_gas_species_values( O3_profile, dry_air_density(i_col,:), & + constituents(i_col,:,index_O3), height_deltas, index_O3, & + errmsg, errcode ) + ASSERT(errcode == 0) + + ! Retrieve the values set by the profile + ! Dry air + call dry_air_profile%get_edge_values( dry_air_interfaces, error ) + ASSERT(error%is_success()) + do i = 1, size(dry_air_interfaces) + ASSERT(dry_air_interfaces(i) == expected_dry_air_interfaces(i_col, i)) + end do + call dry_air_profile%get_layer_densities( dry_air_densities, error ) + ASSERT(error%is_success()) + do i = 1, size(dry_air_densities) - 1 + ASSERT(dry_air_densities(i) == expected_dry_air_densities(i_col, i)) + end do + ! the calculate_exo_layer_density call uses the scale height and the density of the uppermost + ! layer to estimate the density above the column top, which affects the uppermost top value + ASSERT_NEAR(dry_air_densities(size(dry_air_densities)), expected_dry_air_densities(i_col, size(dry_air_densities)) * SCALE_HEIGHT_DRY_AIR * m_to_cm, ABS_ERROR) + dry_air_exo_layer_density = dry_air_profile%get_exo_layer_density( error ) + ASSERT(error%is_success()) + expected_dry_air_exo_layer_density = expected_dry_air_densities(i_col, size(dry_air_densities)) * SCALE_HEIGHT_DRY_AIR * m_to_cm + ASSERT_NEAR(dry_air_exo_layer_density, expected_dry_air_exo_layer_density, ABS_ERROR) + + ! O2 + call O2_profile%get_edge_values( O2_interfaces, error ) + ASSERT(error%is_success()) + do i = 1, size(O2_interfaces) + ASSERT(O2_interfaces(i) == expected_O2_interfaces(i_col, i)) + end do + call O2_profile%get_layer_densities( O2_densities, error ) + ASSERT(error%is_success()) + do i = 1, size(O2_densities) - 1 + ASSERT(O2_densities(i) == expected_O2_densities(i_col, i)) + end do + ASSERT_NEAR(O2_densities(size(O2_densities)), expected_O2_densities(i_col, size(O2_densities)) * SCALE_HEIGHT_O2 * m_to_cm, ABS_ERROR) + O2_exo_layer_density = O2_profile%get_exo_layer_density( error ) + ASSERT(error%is_success()) + expected_O2_exo_layer_density = expected_O2_densities(i_col, size(O2_densities)) * SCALE_HEIGHT_O2 * m_to_cm + ASSERT_NEAR(O2_exo_layer_density, expected_O2_exo_layer_density, ABS_ERROR) + + ! O3 + call O3_profile%get_edge_values( O3_interfaces, error ) + ASSERT(error%is_success()) + do i = 1, size(O3_interfaces) + ASSERT(O3_interfaces(i) == expected_O3_interfaces(i_col, i)) + end do + call O3_profile%get_layer_densities( O3_densities, error ) + ASSERT(error%is_success()) + do i = 1, size(O3_densities) - 1 + ASSERT(O3_densities(i) == expected_O3_densities(i_col, i)) + end do + ASSERT_NEAR(O3_densities(size(O3_densities)), expected_O3_densities(i_col, size(O3_densities)) * SCALE_HEIGHT_O3 * m_to_cm, ABS_ERROR) + O3_exo_layer_density = O3_profile%get_exo_layer_density( error ) + ASSERT(error%is_success()) + expected_O3_exo_layer_density = expected_O3_densities(i_col, size(O3_densities)) * SCALE_HEIGHT_O3 * m_to_cm + ASSERT_NEAR(O3_exo_layer_density, expected_O3_exo_layer_density, ABS_ERROR) + end do + + ! The gas species index starts at 2 because index 1 is reserved for cloud liquid + do i_elem = 2, size(tuvx_species_set) + name = tuvx_species_set(i_elem)%name + unit = tuvx_species_set(i_elem)%unit + molar_mass = tuvx_species_set(i_elem)%molar_mass + scale_height = tuvx_species_set(i_elem)%scale_height + index_musica = tuvx_species_set(i_elem)%index_musica_species + tmp_bool = (trim(name) == 'air' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0289644_kind_phys .and. & + scale_height == 8.01_kind_phys .and. index_musica == 2) .or. & + (trim(name) == 'O2' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0319988_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 3) .or. & + (trim(name) == 'O3' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0479982_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 4) + ASSERT(tmp_bool) + end do + + deallocate( dry_air_profile ) + deallocate( O2_profile ) + deallocate( O3_profile ) + deallocate( height_grid ) + + call musica_ccpp_final(errmsg, errcode) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + + end subroutine test_create_gas_species_profile + + subroutine test_initialize_tuvx_species() + 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_tuvx, only: grid_t, profile_t + use musica_ccpp, only: musica_ccpp_register, musica_ccpp_init, musica_ccpp_final + use musica_ccpp_tuvx_height_grid, only: create_height_grid + use musica_ccpp_namelist, only: filename_of_micm_configuration, & + filename_of_tuvx_configuration, & + filename_of_tuvx_micm_mapping_configuration + use musica_ccpp_tuvx_load_species, only: index_dry_air, index_O2, index_O3 + use musica_ccpp_species, only: tuvx_species_set, MUSICA_INT_UNASSIGNED + + integer, parameter :: NUM_COLUMNS = 2 + integer, parameter :: NUM_LAYERS = 2 + integer, parameter :: NUM_WAVELENGTH_BINS = 102 + integer, parameter :: NUM_TUVX_SPECIES = 4 + real(kind_phys) :: photolysis_wavelength_grid_interfaces(NUM_WAVELENGTH_BINS+1) ! m + type(ccpp_constituent_properties_t), allocatable, & + target :: constituent_props(:) + type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:) + type(ccpp_constituent_properties_t), pointer :: const_prop + real(kind_phys) :: dry_air_density(NUM_COLUMNS,NUM_LAYERS) ! kg m-3 + real(kind_phys) :: constituents(NUM_COLUMNS,NUM_LAYERS, NUM_TUVX_SPECIES) + integer :: errcode + character(len=512) :: errmsg + character(len=50) :: name, unit + real(kind_phys) :: molar_mass ! kg mol-1 + real(kind_phys) :: scale_height ! km + integer :: index_musica, index_constituent_props + logical :: tmp_bool, has_profile + integer :: i_elem, i_col, i, j, k + + 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' + + dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) + dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) + call get_wavelength_edges( photolysis_wavelength_grid_interfaces ) + do k = 1, 4 + do j = 1, 2 + do i = 1, 2 + constituents(i, j, k) = (i + j + k) * 0.1_kind_phys + end do + end do + end do + + call musica_ccpp_register( constituent_props, errmsg, errcode ) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + + 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_COLUMNS, NUM_LAYERS, NUM_LAYERS+1, photolysis_wavelength_grid_interfaces, & + constituent_props_ptr, errmsg, errcode ) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + + ! The gas species index starts at 2 because index 1 is reserved for cloud liquid + do i_elem = 2, size(tuvx_species_set) + name = tuvx_species_set(i_elem)%name + unit = tuvx_species_set(i_elem)%unit + molar_mass = tuvx_species_set(i_elem)%molar_mass + scale_height = tuvx_species_set(i_elem)%scale_height + index_musica = tuvx_species_set(i_elem)%index_musica_species + index_constituent_props = tuvx_species_set(i_elem)%index_constituent_props + has_profile = tuvx_species_set(i_elem)%profiled + tmp_bool = (trim(name) == 'air' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0289644_kind_phys .and. & + scale_height == 8.01_kind_phys .and. index_musica == 2 .and. index_constituent_props /= MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'O2' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0319988_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 3 .and. index_constituent_props /= MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'O3' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0479982_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 4 .and. index_constituent_props /= MUSICA_INT_UNASSIGNED & + .and. has_profile) + ASSERT(tmp_bool) + end do + + call musica_ccpp_final(errmsg, errcode) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + + end subroutine test_initialize_tuvx_species + +end program test_tuvx_gas_species_profile \ No newline at end of file diff --git a/test/musica/tuvx/test_tuvx_height_grid.F90 b/test/musica/tuvx/test_tuvx_height_grid.F90 index 6da26960..c18ce338 100644 --- a/test/musica/tuvx/test_tuvx_height_grid.F90 +++ b/test/musica/tuvx/test_tuvx_height_grid.F90 @@ -25,14 +25,19 @@ subroutine test_create_height_grid() real(kind_phys) :: host_interfaces(NUM_HOST_INTERFACES) = [250.3_kind_phys, 150.2_kind_phys, 50.1_kind_phys] real(kind_phys) :: expected_midpoints(NUM_HOST_MIDPOINTS+1) = [(100.6 + 50.1) * 0.5, 150.2, (250.3 + 200.8) * 0.5] real(kind_phys) :: expected_interfaces(NUM_HOST_INTERFACES+1) = [50.1_kind_phys, 100.6_kind_phys, 200.8_kind_phys, 250.3_kind_phys] + real(kind_phys) :: expected_height_deltas(NUM_HOST_INTERFACES) real(kind_phys) :: midpoints(NUM_HOST_MIDPOINTS+1) real(kind_phys) :: interfaces(NUM_HOST_INTERFACES+1) + real(kind_phys) :: height_deltas(NUM_HOST_INTERFACES) type(grid_t), pointer :: height_grid => null() type(error_t) :: error character(len=512) :: errmsg integer :: errcode integer :: i + expected_height_deltas(:) = expected_interfaces(2:size(expected_interfaces)) & + - expected_interfaces(1:size(expected_interfaces)-1) + height_grid => create_height_grid(-1, 0, errmsg, errcode) ASSERT(errcode == 1) ASSERT(.not. associated(height_grid)) @@ -47,12 +52,16 @@ subroutine test_create_height_grid() ASSERT(associated(height_grid)) call set_height_grid_values(height_grid, host_midpoints, host_interfaces, & - errmsg, errcode) + height_deltas, errmsg, errcode) ASSERT(errcode == 0) ASSERT(height_grid%number_of_sections(error) == NUM_HOST_MIDPOINTS + 1) ASSERT(error%is_success()) + do i = 1, size(height_deltas) + ASSERT_NEAR(height_deltas(i), expected_height_deltas(i), ABS_ERROR) + end do + call height_grid%get_midpoints(midpoints, error) ASSERT(error%is_success()) do i = 1, size(midpoints) diff --git a/test/musica/tuvx/test_tuvx_load_species.F90 b/test/musica/tuvx/test_tuvx_load_species.F90 new file mode 100644 index 00000000..5d34f1e5 --- /dev/null +++ b/test/musica/tuvx/test_tuvx_load_species.F90 @@ -0,0 +1,362 @@ +program test_tuvx_load_species + + use musica_ccpp_tuvx_load_species + + 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_configure_shared_gas_species_tuvx_micm() + call test_configure_partial_shared_gas_species() + call test_configure_no_shared_gas_species() + +contains + + subroutine test_configure_shared_gas_species_tuvx_micm() + ! There are three gas species required for TUV-x: dry air, O2, and O3. + ! This test focuses on configuring MUSICA species and constituent properties + ! when the MICM species include all of these. Cloud liquid water content + ! is the only component specific to TUVX. + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_species, only: musica_species_t, MUSICA_INT_UNASSIGNED + use musica_util, only: error_t + + integer, parameter :: NUM_MICM_SPECIES = 6 + integer, parameter :: NUM_SHARED_SPECIES_BETWEEN_MICM_TUVX = 3 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 4 + type(musica_species_t) :: micm_species(NUM_MICM_SPECIES) + type(musica_species_t), allocatable :: tuvx_species(:) + type(ccpp_constituent_properties_t) :: micm_constituent_props(NUM_MICM_SPECIES) + type(ccpp_constituent_properties_t), allocatable :: tuvx_constituent_props(:) + character(len=512) :: errmsg + integer :: errcode + real(kind_phys) :: molar_mass_group(NUM_MICM_SPECIES) = & + [0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys, 0.4_kind_phys, 0.5_kind_phys, 0.6_kind_phys] + integer :: i_species + character(len=512) :: species_names(NUM_MICM_SPECIES) + character(len=512) :: name, unit, species_name + real(kind_phys) :: molar_mass ! kg mol-1 + real(kind_phys) :: scale_height ! km + integer :: index_musica, index_constituent_props + logical :: is_advected, tmp_bool, has_profile + + species_names(1) = 'N2' + species_names(2) = 'O2' ! shared species + species_names(3) = 'FOO' + species_names(4) = 'O1D' + species_names(5) = 'air' ! shared species + species_names(6) = 'O3' ! shared species + + do i_species = 1, NUM_MICM_SPECIES + call micm_constituent_props(i_species)%instantiate( & + std_name = trim(species_names(i_species)), & + long_name = trim(species_names(i_species)), & + units = 'kg kg-1', & + vertical_dim = 'vertical_layer_dimension', & + default_value = 0.0_kind_phys, & + min_value = 0.0_kind_phys, & + molar_mass = molar_mass_group(i_species), & + advected = .true., & + errcode = errcode, & + errmsg = errmsg) + + micm_species(i_species) = musica_species_t( & + name = species_names(i_species), & + unit = 'kg kg-1', & + molar_mass = molar_mass_group(i_species), & + index_musica_species = i_species ) + end do + + call configure_tuvx_species( micm_species, tuvx_species, tuvx_constituent_props, & + errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(allocated(tuvx_constituent_props)) + ASSERT(size(tuvx_constituent_props) == NUM_TUVX_CONSTITUENTS - NUM_SHARED_SPECIES_BETWEEN_MICM_TUVX) + do i_species = 1, size(tuvx_constituent_props) + ASSERT(tuvx_constituent_props(i_species)%is_instantiated(errcode, errmsg)) + call tuvx_constituent_props(i_species)%standard_name(species_name, errcode, errmsg) + ASSERT(errcode == 0) + call tuvx_constituent_props(i_species)%molar_mass(molar_mass, errcode, errmsg) + ASSERT(errcode == 0) + call tuvx_constituent_props(i_species)%is_advected(is_advected, errcode, errmsg) + ASSERT(errcode == 0) + tmp_bool = (trim(species_name) == 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' .and. & + molar_mass == 0.018_kind_phys .and. is_advected) + ASSERT(tmp_bool) + end do + + ASSERT(allocated(tuvx_species)) + ASSERT(size(tuvx_species) == NUM_TUVX_CONSTITUENTS) + do i_species = 1, size(tuvx_species) + name = tuvx_species(i_species)%name + unit = tuvx_species(i_species)%unit + molar_mass = tuvx_species(i_species)%molar_mass + scale_height = tuvx_species(i_species)%scale_height + index_musica = tuvx_species(i_species)%index_musica_species + index_constituent_props = tuvx_species(i_species)%index_constituent_props + has_profile = tuvx_species(i_species)%profiled + tmp_bool = (trim(name) == 'air' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0289644_kind_phys .and. & + scale_height == 8.01_kind_phys .and. index_musica == 2 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'O2' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0319988_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 3 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'O3' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0479982_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 4 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' .and. & + trim(unit) == 'kg kg-1' .and. molar_mass == 0.018_kind_phys .and. & + scale_height == 0.0_kind_phys .and. index_musica == 1 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. .not. has_profile) + ASSERT(tmp_bool) + end do + + call check_tuvx_species_initialization( errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(index_cloud_liquid_water_content == 1) + ASSERT(index_dry_air == 2) + ASSERT(index_O2 == 3) + ASSERT(index_O3 == 4) + + do i_species = 1, size(tuvx_species) + call tuvx_species(i_species)%deallocate() + end do + deallocate(tuvx_species) + deallocate(tuvx_constituent_props) + + end subroutine test_configure_shared_gas_species_tuvx_micm + + subroutine test_configure_partial_shared_gas_species() + ! This test case applies when some gas species are registered in MICM. + ! It checks which species are already registered and which are not, + ! adding only the new species to the constituent properties. + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_species, only: musica_species_t, MUSICA_INT_UNASSIGNED + use musica_util, only: error_t + + integer, parameter :: NUM_MICM_SPECIES = 6 + integer, parameter :: NUM_SHARED_SPECIES_BETWEEN_MICM_TUVX = 2 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 4 + type(musica_species_t) :: micm_species(NUM_MICM_SPECIES) + type(musica_species_t), allocatable :: tuvx_species(:) + type(ccpp_constituent_properties_t) :: micm_constituent_props(NUM_MICM_SPECIES) + type(ccpp_constituent_properties_t), allocatable :: tuvx_constituent_props(:) + character(len=512) :: errmsg + integer :: errcode + real(kind_phys) :: molar_mass_group(NUM_MICM_SPECIES) = & + [0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys, 0.4_kind_phys, 0.5_kind_phys, 0.6_kind_phys] + integer :: i_species + character(len=512) :: species_names(NUM_MICM_SPECIES) + character(len=512) :: name, unit, species_name + real(kind_phys) :: molar_mass ! kg mol-1 + real(kind_phys) :: scale_height ! km + integer :: index_musica, index_constituent_props + logical :: is_advected, tmp_bool, has_profile + + species_names(1) = 'N2' + species_names(2) = 'O2' ! shared species + species_names(3) = 'FOO' + species_names(4) = 'O1D' + species_names(5) = 'BAZ' + species_names(6) = 'O3' ! shared species + + do i_species = 1, NUM_MICM_SPECIES + call micm_constituent_props(i_species)%instantiate( & + std_name = trim(species_names(i_species)), & + long_name = trim(species_names(i_species)), & + units = 'kg kg-1', & + vertical_dim = 'vertical_layer_dimension', & + default_value = 0.0_kind_phys, & + min_value = 0.0_kind_phys, & + molar_mass = molar_mass_group(i_species), & + advected = .true., & + errcode = errcode, & + errmsg = errmsg) + + micm_species(i_species) = musica_species_t( & + name = species_names(i_species), & + unit = 'kg kg-1', & + molar_mass = molar_mass_group(i_species), & + index_musica_species = i_species ) + end do + + call configure_tuvx_species( micm_species, tuvx_species, tuvx_constituent_props, & + errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(allocated(tuvx_constituent_props)) + ASSERT(size(tuvx_constituent_props) == NUM_TUVX_CONSTITUENTS - NUM_SHARED_SPECIES_BETWEEN_MICM_TUVX) + do i_species = 1, size(tuvx_constituent_props) + ASSERT(tuvx_constituent_props(i_species)%is_instantiated(errcode, errmsg)) + call tuvx_constituent_props(i_species)%standard_name(species_name, errcode, errmsg) + ASSERT(errcode == 0) + call tuvx_constituent_props(i_species)%molar_mass(molar_mass, errcode, errmsg) + ASSERT(errcode == 0) + call tuvx_constituent_props(i_species)%is_advected(is_advected, errcode, errmsg) + ASSERT(errcode == 0) + tmp_bool = (trim(species_name) == 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' .and. & + molar_mass == 0.018_kind_phys .and. is_advected) .or. & + (trim(species_name) == 'air' .and. molar_mass == 0.0289644_kind_phys .and. .not. is_advected) + ASSERT(tmp_bool) + end do + + ASSERT(allocated(tuvx_species)) + ASSERT(size(tuvx_species) == NUM_TUVX_CONSTITUENTS) + do i_species = 1, size(tuvx_species) + name = tuvx_species(i_species)%name + unit = tuvx_species(i_species)%unit + molar_mass = tuvx_species(i_species)%molar_mass + scale_height = tuvx_species(i_species)%scale_height + index_musica = tuvx_species(i_species)%index_musica_species + index_constituent_props = tuvx_species(i_species)%index_constituent_props + has_profile = tuvx_species(i_species)%profiled + tmp_bool = (trim(name) == 'air' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0289644_kind_phys .and. & + scale_height == 8.01_kind_phys .and. index_musica == 2 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'O2' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0319988_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 3 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'O3' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0479982_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 4 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' .and. & + trim(unit) == 'kg kg-1' .and. molar_mass == 0.018_kind_phys .and. & + scale_height == 0.0_kind_phys .and. index_musica == 1 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. .not. has_profile) + ASSERT(tmp_bool) + end do + + call check_tuvx_species_initialization( errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(index_cloud_liquid_water_content == 1) + ASSERT(index_dry_air == 2) + ASSERT(index_O2 == 3) + ASSERT(index_O3 == 4) + + do i_species = 1, size(tuvx_species) + call tuvx_species(i_species)%deallocate() + end do + deallocate(tuvx_species) + deallocate(tuvx_constituent_props) + + end subroutine test_configure_partial_shared_gas_species + + subroutine test_configure_no_shared_gas_species() + ! This test case applies when there are no shared species between MICM and TUV-x. + ! All configured components are added to the constituent properties. + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_species, only: musica_species_t, MUSICA_INT_UNASSIGNED + use musica_util, only: error_t + + integer, parameter :: NUM_MICM_SPECIES = 6 + integer, parameter :: NUM_SHARED_SPECIES_BETWEEN_MICM_TUVX = 0 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 4 + type(musica_species_t) :: micm_species(NUM_MICM_SPECIES) + type(musica_species_t), allocatable :: tuvx_species(:) + type(ccpp_constituent_properties_t) :: micm_constituent_props(NUM_MICM_SPECIES) + type(ccpp_constituent_properties_t), allocatable :: tuvx_constituent_props(:) + character(len=512) :: errmsg + integer :: errcode + real(kind_phys) :: molar_mass_group(NUM_MICM_SPECIES) = & + [0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys, 0.4_kind_phys, 0.5_kind_phys, 0.6_kind_phys] + integer :: i_species + character(len=512) :: species_names(NUM_MICM_SPECIES) + character(len=512) :: name, unit, species_name + real(kind_phys) :: molar_mass ! kg mol-1 + real(kind_phys) :: scale_height ! km + integer :: index_musica, index_constituent_props + logical :: is_advected, tmp_bool, has_profile + + species_names(1) = 'N2' + species_names(2) = 'BAR' + species_names(3) = 'FOO' + species_names(4) = 'O1D' + species_names(5) = 'BAZ' + species_names(6) = 'BOB' + + do i_species = 1, NUM_MICM_SPECIES + call micm_constituent_props(i_species)%instantiate( & + std_name = trim(species_names(i_species)), & + long_name = trim(species_names(i_species)), & + units = 'kg kg-1', & + vertical_dim = 'vertical_layer_dimension', & + default_value = 0.0_kind_phys, & + min_value = 0.0_kind_phys, & + molar_mass = molar_mass_group(i_species), & + advected = .true., & + errcode = errcode, & + errmsg = errmsg) + + micm_species(i_species) = musica_species_t( & + name = species_names(i_species), & + unit = 'kg kg-1', & + molar_mass = molar_mass_group(i_species), & + index_musica_species = i_species ) + end do + + call configure_tuvx_species( micm_species, tuvx_species, tuvx_constituent_props, & + errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(allocated(tuvx_constituent_props)) + ASSERT(size(tuvx_constituent_props) == NUM_TUVX_CONSTITUENTS - NUM_SHARED_SPECIES_BETWEEN_MICM_TUVX) + do i_species = 1, size(tuvx_constituent_props) + ASSERT(tuvx_constituent_props(i_species)%is_instantiated(errcode, errmsg)) + call tuvx_constituent_props(i_species)%standard_name(species_name, errcode, errmsg) + ASSERT(errcode == 0) + call tuvx_constituent_props(i_species)%molar_mass(molar_mass, errcode, errmsg) + ASSERT(errcode == 0) + call tuvx_constituent_props(i_species)%is_advected(is_advected, errcode, errmsg) + ASSERT(errcode == 0) + tmp_bool = (trim(species_name) == 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' .and. & + molar_mass == 0.018_kind_phys .and. is_advected) .or. & + (trim(species_name) == 'air' .and. molar_mass == 0.0289644_kind_phys .and. .not. is_advected) .or. & + (trim(species_name) == 'O2' .and. molar_mass == 0.0319988_kind_phys .and. .not. is_advected) .or. & + (trim(species_name) == 'O3' .and. molar_mass == 0.0479982_kind_phys .and. .not. is_advected) + ASSERT(tmp_bool) + end do + + ASSERT(allocated(tuvx_species)) + ASSERT(size(tuvx_species) == NUM_TUVX_CONSTITUENTS) + do i_species = 1, size(tuvx_species) + name = tuvx_species(i_species)%name + unit = tuvx_species(i_species)%unit + molar_mass = tuvx_species(i_species)%molar_mass + scale_height = tuvx_species(i_species)%scale_height + index_musica = tuvx_species(i_species)%index_musica_species + index_constituent_props = tuvx_species(i_species)%index_constituent_props + has_profile = tuvx_species(i_species)%profiled + tmp_bool = (trim(name) == 'air' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0289644_kind_phys .and. & + scale_height == 8.01_kind_phys .and. index_musica == 2 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'O2' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0319988_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 3 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'O3' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0479982_kind_phys .and. & + scale_height == 7.0_kind_phys .and. index_musica == 4 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. has_profile) .or. & + (trim(name) == 'cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water' .and. & + trim(unit) == 'kg kg-1' .and. molar_mass == 0.018_kind_phys .and. & + scale_height == 0.0_kind_phys .and. index_musica == 1 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. .not. has_profile) + ASSERT(tmp_bool) + end do + + call check_tuvx_species_initialization( errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(index_cloud_liquid_water_content == 1) + ASSERT(index_dry_air == 2) + ASSERT(index_O2 == 3) + ASSERT(index_O3 == 4) + + do i_species = 1, size(tuvx_species) + call tuvx_species(i_species)%deallocate() + end do + deallocate(tuvx_species) + deallocate(tuvx_constituent_props) + + end subroutine test_configure_no_shared_gas_species + +end program test_tuvx_load_species \ No newline at end of file diff --git a/to_be_ccppized/ccpp_tuvx_utils.F90 b/to_be_ccppized/ccpp_tuvx_utils.F90 index 837801bc..e7329d56 100644 --- a/to_be_ccppized/ccpp_tuvx_utils.F90 +++ b/to_be_ccppized/ccpp_tuvx_utils.F90 @@ -10,7 +10,7 @@ module ccpp_tuvx_utils contains !> Regrids normalized flux data to match a specified wavelength grid - !! This function is copied from CAM/src/chemistry/utils/mo_util.F90 + ! This function is copied from CAM/src/chemistry/utils/mo_util.F90 subroutine rebin( source_dimension, target_dimension, source_coordinates, & target_coordinates, source, target ) use ccpp_kinds, only: kind_phys