From 61ca4ebadf98a691d1ee590ea1f89bab9d9b1ca9 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Thu, 20 Feb 2025 12:38:29 -0800 Subject: [PATCH] move NO photolysis logic to TUVx modules --- schemes/musica/micm/musica_ccpp_micm.F90 | 16 +- schemes/musica/musica_ccpp.F90 | 35 +- schemes/musica/musica_ccpp_species.F90 | 3 +- schemes/musica/musica_ccpp_util.F90 | 4 +- schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 57 +- .../tuvx/musica_ccpp_tuvx_load_species.F90 | 46 +- ...F90 => musica_ccpp_tuvx_no_photolysis.F90} | 343 ++++++------ test/docker/Dockerfile.musica | 4 +- test/docker/Dockerfile.musica.no_install | 4 +- .../data/NO_photolysis/mapping.nojNO.json | 10 + .../data/NO_photolysis/mapping.valid.json | 10 + test/musica/test_musica_api.F90 | 21 +- test/musica/tuvx/test_tuvx_gas_species.F90 | 6 +- test/musica/tuvx/test_tuvx_load_species.F90 | 490 +++++++++++++++++- test/musica/tuvx/test_tuvx_no_photolysis.F90 | 313 +++-------- 15 files changed, 851 insertions(+), 511 deletions(-) rename schemes/musica/tuvx/{musica_ccpp_tuvx_no_photolysis_rate.F90 => musica_ccpp_tuvx_no_photolysis.F90} (75%) create mode 100644 test/musica/data/NO_photolysis/mapping.nojNO.json create mode 100644 test/musica/data/NO_photolysis/mapping.valid.json diff --git a/schemes/musica/micm/musica_ccpp_micm.F90 b/schemes/musica/micm/musica_ccpp_micm.F90 index d13b10e5..d51b51cc 100644 --- a/schemes/musica/micm/musica_ccpp_micm.F90 +++ b/schemes/musica/micm/musica_ccpp_micm.F90 @@ -19,18 +19,16 @@ module musica_ccpp_micm !> Registers MICM constituent properties with the CCPP subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, & - micm_species, includes_no_photolysis, errmsg, errcode) + micm_species, 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 use iso_c_binding, only: c_int - use musica_ccpp_tuvx_no_photolysis_rate, only : set_NO_photolysis_index 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(:) - logical, intent(out) :: includes_no_photolysis character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -100,18 +98,6 @@ subroutine micm_register(solver_type, number_of_grid_cells, constituent_props, & end do number_of_rate_parameters = micm%user_defined_reaction_rates%size() - includes_no_photolysis = .false. - ! iterate through the user defined reaction rates and print their names and indices - do i = 1, number_of_rate_parameters - associate( map => micm%user_defined_reaction_rates ) - ! check if any rate is named PHOTO.jNO - if (index(map%name(i), "PHOTO.jNO") > 0) then - includes_no_photolysis = .true. - call set_NO_photolysis_index(map%index(i)) - end if - end associate ! map - end do - end subroutine micm_register !> Initializes MICM diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index d2ff53eb..41b4796c 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -18,7 +18,6 @@ subroutine musica_ccpp_register(constituent_props, errmsg, errcode) 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 - use musica_ccpp_tuvx_no_photolysis_rate, only: check_NO_exist type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituent_props(:) character(len=512), intent(out) :: errmsg @@ -28,18 +27,15 @@ subroutine musica_ccpp_register(constituent_props, errmsg, errcode) type(ccpp_constituent_properties_t), allocatable :: constituent_props_subset(:) type(musica_species_t), allocatable :: micm_species(:) type(musica_species_t), allocatable :: tuvx_species(:) - logical :: includes_no_photolysis integer :: number_of_grid_cells - includes_no_photolysis = .false. - ! Temporary fix until the number of grid cells is only needed to create a MICM state ! instead of when the solver is created. ! The number of grid cells is not known at this point, so we set it to 1 and recreate ! 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, & - micm_species, includes_no_photolysis,errmsg, errcode) + micm_species, errmsg, errcode) if (errcode /= 0) return constituent_props = constituent_props_subset deallocate(constituent_props_subset) @@ -53,8 +49,6 @@ subroutine musica_ccpp_register(constituent_props, errmsg, errcode) call check_tuvx_species_initialization(errmsg, errcode) if (errcode /= 0) return - call check_NO_exist(micm_species, includes_no_photolysis) - end subroutine musica_ccpp_register !> \section arg_table_musica_ccpp_init Argument Table @@ -70,8 +64,6 @@ subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, & 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 - use musica_ccpp_tuvx_no_photolysis_rate, & - only: is_NO_photolysis_active, set_NO_index_constituent_props, check_NO_initialization integer, intent(in) :: horizontal_dimension ! (count) integer, intent(in) :: vertical_layer_dimension ! (count) @@ -84,7 +76,6 @@ subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, & ! local variables type(ccpp_constituent_properties_t), allocatable :: constituent_props(:) type(musica_species_t), allocatable :: micm_species(:) - logical :: includes_no_photolysis integer :: number_of_grid_cells ! Temporary fix until the number of grid cells is only needed to create a MICM state @@ -92,7 +83,8 @@ subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, & ! 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, constituent_props, & - micm_species, includes_no_photolysis, errmsg, errcode) + micm_species, errmsg, errcode) + if (errcode /= 0) return call micm_init(errmsg, errcode) if (errcode /= 0) return call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & @@ -107,13 +99,6 @@ subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, & call check_initialization(errmsg, errcode) if (errcode /= 0) return - if (is_NO_photolysis_active) then - call set_NO_index_constituent_props(constituent_props_ptr, errmsg, errcode) - if (errcode /= 0) return - call check_NO_initialization(errmsg, errcode) - if (errcode /= 0) return - end if - end subroutine musica_ccpp_init !> \section arg_table_musica_ccpp_run Argument Table @@ -137,8 +122,6 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co 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 - use musica_ccpp_tuvx_no_photolysis_rate, & - only: is_NO_photolysis_active, num_NO_photolysis_constituents, NO_photolysis_indices_constituent_props real(kind_phys), intent(in) :: time_step ! s real(kind_phys), target, intent(in) :: temperature(:,:) ! K (column, layer) @@ -172,24 +155,14 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co real(kind_phys), dimension(size(constituents, dim=1), & size(constituents, dim=2), & number_of_tuvx_species) :: constituents_tuvx_species ! kg kg-1 - real(kind_phys), dimension(size(constituents, dim=1), & - size(constituents, dim=2), & - num_NO_photolysis_constituents) :: constituents_NO_photolysis ! kg kg-1 call extract_subset_constituents(tuvx_indices_constituent_props, constituents, & constituents_tuvx_species, errmsg, errcode) if (errcode /= 0) return - if (is_NO_photolysis_active) then - call extract_subset_constituents(NO_photolysis_indices_constituent_props, constituents, & - constituents_NO_photolysis, errmsg, errcode) - if (errcode /= 0) return - end if - ! Calculate photolysis rate constants using TUV-x call tuvx_run(temperature, dry_air_density, & constituents_tuvx_species, & - constituents_NO_photolysis, & geopotential_height_wrt_surface_at_midpoint, & geopotential_height_wrt_surface_at_interface, & surface_geopotential, surface_temperature, & @@ -203,6 +176,7 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co earth_sun_distance, & rate_parameters, & errmsg, errcode) + if (errcode /= 0) return call extract_subset_constituents(micm_indices_constituent_props, constituents, & constituents_micm_species, errmsg, errcode) @@ -214,6 +188,7 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co ! Solve chemistry at the current time step call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, & constituents_micm_species, errmsg, errcode) + if (errcode /= 0) return ! Convert MICM unit back to CAM-SIMA unit (mol m-3 -> kg kg-1) call convert_to_mass_mixing_ratio(dry_air_density, micm_molar_mass_array, constituents_micm_species) diff --git a/schemes/musica/musica_ccpp_species.F90 b/schemes/musica/musica_ccpp_species.F90 index 3f746a2f..54b19c24 100644 --- a/schemes/musica/musica_ccpp_species.F90 +++ b/schemes/musica/musica_ccpp_species.F90 @@ -1,5 +1,6 @@ module musica_ccpp_species use ccpp_kinds, only: kind_phys + use musica_ccpp_util, only: MUSICA_INT_UNASSIGNED implicit none private @@ -8,8 +9,6 @@ module musica_ccpp_species 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 diff --git a/schemes/musica/musica_ccpp_util.F90 b/schemes/musica/musica_ccpp_util.F90 index f51fcb2c..1ae6cca9 100644 --- a/schemes/musica/musica_ccpp_util.F90 +++ b/schemes/musica/musica_ccpp_util.F90 @@ -11,7 +11,9 @@ module musica_ccpp_util real(kind_phys), parameter, public :: PI = 3.14159265358979323846_kind_phys real(kind_phys), parameter, public :: DEGREE_TO_RADIAN = PI / 180.0_kind_phys - real(kind_phys), parameter, public :: EARTH_RADIUS_KM = 6378_kind_phys + real(kind_phys), parameter, public :: EARTH_RADIUS_M = 6378000.0_kind_phys + real(kind_phys), parameter, public :: AVOGADRO = 6.02214076e23_kind_phys + integer, parameter, public :: MUSICA_INT_UNASSIGNED = -99999 contains diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index 2c4c508d..e7de83d3 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -133,7 +133,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & 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 + use musica_util, only: error_t, configuration_t, MUSICA_INDEX_MAPPINGS_MAP_ANY use musica_ccpp_namelist, only: filename_of_tuvx_micm_mapping_configuration use musica_ccpp_util, only: PI use musica_ccpp_tuvx_height_grid, & @@ -155,6 +155,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & 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 + use musica_ccpp_tuvx_no_photolysis, only: NO_photolysis_init integer, intent(in) :: vertical_layer_dimension ! (count) integer, intent(in) :: vertical_interface_dimension ! (count) @@ -486,24 +487,30 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & end if photolysis_rate_constants_mapping => & - index_mappings_t( config, photolysis_rate_constants_ordering, & + index_mappings_t( config, MUSICA_INDEX_MAPPINGS_MAP_ANY, & + photolysis_rate_constants_ordering, & micm_rate_parameter_ordering, error ) + deallocate( photolysis_rate_constants_ordering ) if (has_error_occurred( error, errmsg, errcode )) then deallocate( tuvx ) tuvx => null() call cleanup_tuvx_resources() - deallocate( photolysis_rate_constants_ordering ) return end if - deallocate( photolysis_rate_constants_ordering ) + call NO_photolysis_init(config, micm_rate_parameter_ordering, errmsg, errcode) + if (errcode /= 0) then + deallocate( tuvx ) + tuvx => null() + call cleanup_tuvx_resources() + return + endif end subroutine tuvx_init !> Calculates photolysis rate constants for the current model conditions subroutine tuvx_run(temperature, dry_air_density, & constituents, & - constituents_NO_photolysis, & geopotential_height_wrt_surface_at_midpoint, & geopotential_height_wrt_surface_at_interface, & surface_geopotential, surface_temperature, & @@ -520,7 +527,7 @@ subroutine tuvx_run(temperature, dry_air_density, & 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_util, only: has_error_occurred, PI, MUSICA_INT_UNASSIGNED 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 @@ -528,14 +535,11 @@ subroutine tuvx_run(temperature, dry_air_density, & 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 - use musica_ccpp_tuvx_no_photolysis_rate, only: is_NO_photolysis_active, calculate_NO_photolysis_rate, & - index_NO_photolysis_rate - use musica_ccpp_species, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_tuvx_no_photolysis, only: calculate_NO_photolysis_rate_constants 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(:,:,:) ! kg kg-1 (column, layer, constituent) - real(kind_phys), intent(in) :: constituents_NO_photolysis(:,:,:) ! 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 (column) @@ -548,14 +552,15 @@ subroutine tuvx_run(temperature, dry_air_density, & real(kind_phys), intent(in) :: air_pressure_thickness(:,:) ! Pa (column, layer) real(kind_phys), intent(in) :: solar_zenith_angle(:) ! radians (column) real(kind_phys), intent(in) :: earth_sun_distance ! m - real(kind_phys), intent(inout) :: rate_parameters(:,:,:) ! various units (column, layer, reaction) + real(kind_phys), intent(inout) :: rate_parameters(:,:,:) ! various units (column, layer, parameter) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode ! local variables - real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_midpoint, dim = 2)) :: height_midpoints - real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_interface, dim = 2)) :: height_interfaces - real(kind_phys), dimension(size(height_interfaces)) :: height_deltas ! km + real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_midpoint, dim = 2)) :: height_midpoints ! km (layer) + real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_interface, dim = 1), & + size(geopotential_height_wrt_surface_at_interface, dim = 2)) :: height_interfaces ! km (column, interface) + real(kind_phys), dimension(size(height_interfaces, dim=2)) :: height_deltas ! km (layer) 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) @@ -584,9 +589,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(i_col,:) ) - call set_height_grid_values( height_grid, height_midpoints, height_interfaces, & + call set_height_grid_values( height_grid, height_midpoints, height_interfaces(i_col,:), & height_deltas, errmsg, errcode ) if (errcode /= 0) return @@ -634,30 +639,24 @@ subroutine tuvx_run(temperature, dry_air_density, & max( photolysis_rate_constants(:,:), 0.0_kind_phys ) end if ! solar zenith angle check - if (is_NO_photolysis_active) then - ! TODO: How do I ensure that the extraterrestrial flux matches the wavelength grid needed for NO photolysis? - NO_photolysis_rate_constant = calculate_NO_photolysis_rate( & - solar_zenith_angle(i_col), extraterrestrial_flux, height_interfaces, & - dry_air_density(i_col,:), constituents(i_col,:,index_O2), & - constituents(i_col,:,index_O3), constituents_NO_photolysis(i_col,:,:) ) - end if - - ! TODO: Where does jno go into the photolysis_rate_constants array? ! map photolysis rate constants to the host model's rate parameters and vertical grid do i_level = 1, size(rate_parameters, dim=2) call photolysis_rate_constants_mapping%copy_data( & photolysis_rate_constants(size(rate_parameters, dim=2)-i_level+2,:), & rate_parameters(i_col,i_level,:) ) - if (is_NO_photolysis_active) then - rate_parameters(i_col,i_level,index_NO_photolysis_rate) = NO_photolysis_rate_constant(i_level) - end if end do end do + ! TODO: How do I ensure that the extraterrestrial flux matches the wavelength grid needed for NO photolysis? + call calculate_NO_photolysis_rate_constants( solar_zenith_angle, extraterrestrial_flux, & + height_interfaces, dry_air_density, constituents, rate_parameters ) + end subroutine tuvx_run !> Finalizes TUV-x subroutine tuvx_final(errmsg, errcode) + use musica_ccpp_tuvx_no_photolysis, only: NO_photolysis_final + character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -671,6 +670,8 @@ subroutine tuvx_final(errmsg, errcode) tuvx => null() end if + call NO_photolysis_final() + end subroutine tuvx_final end module musica_ccpp_tuvx diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 index 243849ca..9780fbb2 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 @@ -1,6 +1,6 @@ module musica_ccpp_tuvx_load_species use ccpp_kinds, only: kind_phys - use musica_ccpp_species, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_util, only: MUSICA_INT_UNASSIGNED implicit none private @@ -11,6 +11,8 @@ module musica_ccpp_tuvx_load_species 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 + integer, protected, public :: index_NO = MUSICA_INT_UNASSIGNED + integer, protected, public :: index_N2 = MUSICA_INT_UNASSIGNED ! Constants ! Cloud liquid water @@ -24,6 +26,8 @@ module musica_ccpp_tuvx_load_species 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 :: NO_LABEL = 'NO' + character(len=*), parameter, public :: N2_LABEL = 'N2' 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 @@ -32,6 +36,8 @@ module musica_ccpp_tuvx_load_species 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 + real(kind_phys), parameter, public :: MOLAR_MASS_NO = 0.0300061_kind_phys ! kg mol-1 + real(kind_phys), parameter, public :: MOLAR_MASS_N2 = 0.0280134_kind_phys ! kg mol-1 contains @@ -42,7 +48,6 @@ 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(:) @@ -57,15 +62,19 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props, ! 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. + logical :: is_dry_air_registered + logical :: is_O2_registered + logical :: is_O3_registered + logical :: is_NO_registered + logical :: is_N2_registered 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. + is_NO_registered = .false. + is_N2_registered = .false. ! Register cloud liquid water content needed for cloud optics calculations i_new = 1 @@ -99,8 +108,15 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props, is_O3_registered = .true. micm_species(i_species)%profiled = .true. micm_species(i_species)%scale_height = SCALE_HEIGHT_O3 + else if ( micm_species(i_species)%name == NO_LABEL ) then + is_NO_registered = .true. + else if ( micm_species(i_species)%name == N2_LABEL ) then + is_N2_registered = .true. end if end do + if ( is_NO_registered .and. is_N2_registered ) then + num_new_species = num_new_species + 2 + end if if (.not. is_dry_air_registered) then i_new = i_new + 1 @@ -191,6 +207,26 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props, index_musica_species = i_tuvx_species, & profiled = .true., & scale_height = SCALE_HEIGHT_O3 ) + + if ( is_NO_registered .and. is_N2_registered ) then + i_tuvx_species = i_tuvx_species + 1 + index_NO = i_tuvx_species + tuvx_species(i_tuvx_species) = musica_species_t( & + name = NO_LABEL, & + unit = TUVX_GAS_SPECIES_UNITS, & ! TUV-x profile unit, different from molar mass unit + molar_mass = MOLAR_MASS_NO, & ! kg mol-1 + index_musica_species = i_tuvx_species, & + scale_height = 0.0_kind_phys ) + + i_tuvx_species = i_tuvx_species + 1 + index_N2 = i_tuvx_species + tuvx_species(i_tuvx_species) = musica_species_t( & + name = N2_LABEL, & + unit = TUVX_GAS_SPECIES_UNITS, & ! TUV-x profile unit, different from molar mass unit + molar_mass = MOLAR_MASS_N2, & ! kg mol-1 + index_musica_species = i_tuvx_species, & + scale_height = 0.0_kind_phys ) + end if end subroutine configure_tuvx_species diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_no_photolysis_rate.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_no_photolysis.F90 similarity index 75% rename from schemes/musica/tuvx/musica_ccpp_tuvx_no_photolysis_rate.F90 rename to schemes/musica/tuvx/musica_ccpp_tuvx_no_photolysis.F90 index 3373e66e..033b7fca 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_no_photolysis_rate.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_no_photolysis.F90 @@ -1,12 +1,13 @@ -module musica_ccpp_tuvx_no_photolysis_rate +module musica_ccpp_tuvx_no_photolysis use ccpp_kinds, only: dk => kind_phys - use musica_ccpp_species, only: MUSICA_INT_UNASSIGNED + use musica_util, only: index_mappings_t + use musica_ccpp_util, only: MUSICA_INT_UNASSIGNED implicit none !> @file musica_ccpp_tuvx_no_photolysis_rate.F90 - !> @brief Module for calculating the NO photolysis rate + !> @brief Module for calculating the NO photolysis rate constant ! This module calculates the NO photolysis rate using the same method in CAM ! https://github.com/ESCOMP/CAM/blob/ab476f9b7345cbefdc4cf67ff17f0fe85d8c7387/src/chemistry/mozart/mo_jshort.F90#L1796-L1870 ! Much of the code is also based off of a PR which incorporated TUVx into CAM @@ -22,214 +23,200 @@ module musica_ccpp_tuvx_no_photolysis_rate ! MS93: A reference to the Minschwaner and Siskind paper private - public :: check_NO_exist, set_NO_index_constituent_props, check_NO_initialization, & - calculate_NO_photolysis_rate, set_NO_photolysis_index - - character(len=*), parameter, public :: N2_label = 'N2' - character(len=*), parameter, public :: NO_label = 'NO' - integer, parameter, public :: num_NO_photolysis_constituents = 2 - logical, protected, public :: is_NO_photolysis_active = .false. - integer, protected, public :: index_NO_photolysis_rate = MUSICA_INT_UNASSIGNED - integer, protected, public :: index_N2 = MUSICA_INT_UNASSIGNED - integer, protected, public :: index_NO = MUSICA_INT_UNASSIGNED - integer, protected, public :: constituent_props_index_N2 = MUSICA_INT_UNASSIGNED - integer, protected, public :: constituent_props_index_NO = MUSICA_INT_UNASSIGNED - integer, protected, public :: NO_photolysis_indices_constituent_props(2) = & - [MUSICA_INT_UNASSIGNED, MUSICA_INT_UNASSIGNED] - real(dk), protected, public :: MOLAR_MASS_N2 = MUSICA_INT_UNASSIGNED ! kg mol-1 - real(dk), protected, public :: MOLAR_MASS_NO = MUSICA_INT_UNASSIGNED ! kg mol-1 + public :: NO_photolysis_init, calculate_NO_photolysis_rate_constants, NO_photolysis_final -contains + real(dk), parameter :: MAX_SOLAR_ZENITH_ANGLE = 110.0_dk ! degrees + real(dk), parameter :: MIN_SOLAR_ZENITH_ANGLE = 0.0_dk ! degrees - !> Checks for the presence of NO and N2, and sets the flag to indicate - !! whether the NO photolysis calculation is enabled or disabled - subroutine check_NO_exist(micm_species, includes_no_photolysis) - use musica_ccpp_species, only: musica_species_t + character(len=*), parameter, public :: NO_PHOTOLYSIS_LABEL = "jNO" + logical, protected, public :: do_NO_photolysis = .false. + type(index_mappings_t), pointer, protected, public :: & + NO_photolysis_rate_constants_mapping => null() - type(musica_species_t), intent(in) :: micm_species(:) - logical, intent(in) :: includes_no_photolysis +contains - ! local variables - logical :: is_N2_registered = .false. - logical :: is_NO_registered = .false. - integer :: num_micm_species - integer :: i_species - - is_N2_registered = .false. - is_NO_registered = .false. - num_micm_species = size(micm_species) - - do i_species = 1, num_micm_species - if ( is_N2_registered .and. is_NO_registered ) then - is_NO_photolysis_active = .true. - exit - end if + !> Creates a set of mapped indices for the NO photolysis rate constant to + !! the expected rate parameter ordering from the chemistry solver. + !! Also determines whether NO photolysis should be calculated. + subroutine NO_photolysis_init(config, expected_rate_parameter_ordering, & + errmsg, errcode) + use musica_util, only : mappings_t, mapping_t_c, mapping_t, create_string_c, & + to_c_string, copy_mappings, configuration_t, error_t, & + MUSICA_INDEX_MAPPINGS_MAP_ANY, delete_string_c + use musica_ccpp_util, only : has_error_occurred, MUSICA_INT_UNASSIGNED + use musica_ccpp_tuvx_load_species, only : index_NO, index_N2, index_O2, index_O3 + use iso_c_binding, only : c_loc, c_size_t + + type(configuration_t), intent(in) :: config + type(mappings_t), intent(in) :: expected_rate_parameter_ordering + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + type(mapping_t_c), target :: c_mappings(1) + type(mapping_t), allocatable :: f_mappings(:) + type(mappings_t), pointer :: NO_parameter_map + type(error_t) :: error - if ( micm_species(i_species)%name == N2_label ) then - is_N2_registered = .true. - MOLAR_MASS_N2 = micm_species(i_species)%molar_mass - else if ( micm_species(i_species)%name == NO_label ) then - is_NO_registered = .true. - MOLAR_MASS_NO = micm_species(i_species)%molar_mass + errcode = 0 + errmsg = '' + c_mappings(1)%index_ = 0 + c_mappings(1)%name_ = create_string_c( to_c_string( NO_PHOTOLYSIS_LABEL ) ) + call copy_mappings( c_loc( c_mappings ), 1_c_size_t, f_mappings ) + call delete_string_c( c_mappings(1)%name_ ) + NO_parameter_map => mappings_t( f_mappings ) + NO_photolysis_rate_constants_mapping => index_mappings_t( config, & + MUSICA_INDEX_MAPPINGS_MAP_ANY, NO_parameter_map, & + expected_rate_parameter_ordering, error ) + deallocate( no_parameter_map ) + if (has_error_occurred( error, errmsg, errcode )) return + if (NO_photolysis_rate_constants_mapping%size() > 0) then + do_NO_photolysis = .true. + if (index_NO == MUSICA_INT_UNASSIGNED .or. & + index_N2 == MUSICA_INT_UNASSIGNED .or. & + index_O3 == MUSICA_INT_UNASSIGNED .or. & + index_O2 == MUSICA_INT_UNASSIGNED) then + do_NO_photolysis = .false. + errmsg = "[MUSICA Error] The indices for NO and N2 have not "// & + "been initialized, and are needed for jNO calculations." + errcode = 1 + return end if - end do - - if (is_N2_registered .and. is_NO_registered .and. includes_no_photolysis) then - is_NO_photolysis_active = .true. end if - end subroutine check_NO_exist + end subroutine NO_photolysis_init - subroutine set_NO_photolysis_index(no_index) - integer, intent(in) :: no_index + !> Prepares to calculate the photolysis rate of NO (nitiric acid). + !! This function transforms the inputs into the units + !! needed by the calculate_jno routine which actually produces the photolysis rate + subroutine calculate_NO_photolysis_rate_constants(solar_zenith_angle, & + extraterrestrial_flux, height_at_interfaces, dry_air_density, & + constituents, rate_parameters) + use musica_ccpp_tuvx_load_species, only: MOLAR_MASS_O2, MOLAR_MASS_O3, MOLAR_MASS_NO, & + MOLAR_MASS_N2, index_O2, index_O3, index_NO, & + index_N2 + use musica_ccpp_tuvx_gas_species, only: km_to_cm + use musica_ccpp_species, only: tuvx_indices_constituent_props + use musica_ccpp_util, only: PI - index_NO_photolysis_rate = no_index + real(dk), intent(in) :: solar_zenith_angle(:) ! radians (column) + real(dk), intent(in) :: extraterrestrial_flux(:) ! photons cm-2 s-1 nm-1 (wavelength bin) + real(dk), intent(in) :: height_at_interfaces(:,:) ! km (column, layer) + real(dk), intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer) + real(dk), intent(in) :: constituents(:,:,:) ! kg kg-1 (column, layer, species) + real(dk), intent(inout) :: rate_parameters(:,:,:) ! s-1 (column, layer, parameter) - end subroutine set_NO_photolysis_index + ! local variables + integer :: number_of_columns + integer :: number_of_vertical_layers + real(dk) :: N2_densities(size(dry_air_density, dim=2)+1) ! molecule cm-3 + real(dk) :: O2_densities(size(dry_air_density, dim=2)+1) ! molecule cm-3 + real(dk) :: O3_densities(size(dry_air_density, dim=2)+1) ! molecule cm-3 + real(dk) :: NO_densities(size(dry_air_density, dim=2)+1) ! molecule cm-3 + real(dk) :: O2_slant_column_densities(size(dry_air_density, dim=2)+1) ! molecule cm-2 + real(dk) :: O3_slant_column_densities(size(dry_air_density, dim=2)+1) ! molecule cm-2 + real(dk) :: NO_slant_column_densities(size(dry_air_density, dim=2)+1) ! molecule cm-2 - !> Gets the indices of NO constituents from the CCPP constituent properties, - !! and stores them in the indices array along with their relative positions - !! in the NO constituents array - subroutine set_NO_index_constituent_props(constituent_props, errmsg, errcode) - use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t - use ccpp_const_utils, only: ccpp_const_get_idx + ! parameters needed to calculate slant column densities + ! (see sphers routine description for details) + integer :: number_of_crossed_layers(size(dry_air_density, dim=2)+1) + real(dk) :: slant_path(0:size(dry_air_density, dim=2)+1, size(dry_air_density, dim=2)+1) + real(dk) :: delta_z(size(dry_air_density, dim=2)+1) ! layer thickness (cm) + real(dk) :: sza_degrees + integer :: i_col, i_level - type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode + ! final photolysis rate (column, layer) + real(dk) :: jNO(size(dry_air_density, dim=2), 1) - call ccpp_const_get_idx( constituent_props, N2_label, & - constituent_props_index_N2, errmsg, errcode ) - if (errcode /= 0) return - index_N2 = 1 + if (.not. do_NO_photolysis) return - call ccpp_const_get_idx( constituent_props, NO_label, & - constituent_props_index_NO, errmsg, errcode ) - if (errcode /= 0) return - index_NO = index_N2 + 1 + number_of_columns = size(dry_air_density, dim=1) + number_of_vertical_layers = size(dry_air_density, dim=2) - NO_photolysis_indices_constituent_props = & - [constituent_props_index_N2, constituent_props_index_NO] + do i_col = 1, number_of_columns - end subroutine set_NO_index_constituent_props + sza_degrees = solar_zenith_angle(i_col) * 180.0_dk / PI + if (sza_degrees > MAX_SOLAR_ZENITH_ANGLE .or. & + sza_degrees < MIN_SOLAR_ZENITH_ANGLE) then + do i_level = 1, number_of_vertical_layers + call NO_photolysis_rate_constants_mapping%copy_data( & + (/ 0.0_dk /), rate_parameters(i_col, i_level, :) ) + end do + cycle + end if - !> Checks if the variables associated with NO photolysis constituents - !! have been initialized - subroutine check_NO_initialization(errmsg, errcode) - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode + ! TODO: what are these constants? scale heights? + ! TODO: the values at index 1 appear to be for values above the model top in CAM, but how does that affect cam sima? + call convert_mixing_ratio_to_molecule_cm3(constituents(i_col,:,tuvx_indices_constituent_props(index_N2)), & + dry_air_density(i_col,:), MOLAR_MASS_N2, N2_densities(2:)) + N2_densities(1) = N2_densities(2) * 0.9_dk + + call convert_mixing_ratio_to_molecule_cm3(constituents(i_col,:,tuvx_indices_constituent_props(index_O2)), & + dry_air_density(i_col,:), MOLAR_MASS_O2, O2_densities(2:)) + O2_densities(1) = O2_densities(2) * 7.0_dk / ( height_at_interfaces(i_col,1) - height_at_interfaces(i_col,2) ) + + call convert_mixing_ratio_to_molecule_cm3(constituents(i_col,:,tuvx_indices_constituent_props(index_O3)), & + dry_air_density(i_col,:), MOLAR_MASS_O3, O3_densities(2:)) + O3_densities(1) = O3_densities(2) * 7.0_dk / ( height_at_interfaces(i_col,1) - height_at_interfaces(i_col,2) ) + + call convert_mixing_ratio_to_molecule_cm3(constituents(i_col,:,tuvx_indices_constituent_props(index_NO)), & + dry_air_density(i_col,:), MOLAR_MASS_NO, NO_densities(2:)) + NO_densities(1) = NO_densities(2) * 0.9_dk + + ! calculate slant column densities + call calculate_slant_path( number_of_vertical_layers, height_at_interfaces(i_col,:), & + solar_zenith_angle(i_col), slant_path, number_of_crossed_layers ) + delta_z(1:number_of_vertical_layers) = km_to_cm * & + ( height_at_interfaces(i_col, 1:number_of_vertical_layers) & + - height_at_interfaces(i_col, 2:number_of_vertical_layers+1) ) + call calculate_slant_column_density( number_of_vertical_layers+1, delta_z,& + slant_path, number_of_crossed_layers, O2_densities, O2_slant_column_densities ) + call calculate_slant_column_density( number_of_vertical_layers+1, delta_z, & + slant_path, number_of_crossed_layers, O3_densities, O3_slant_column_densities ) + call calculate_slant_column_density( number_of_vertical_layers+1, delta_z, & + slant_path, number_of_crossed_layers, NO_densities, NO_slant_column_densities ) + + jNO(:,1) = calculate_jno(number_of_vertical_layers, extraterrestrial_flux, & + N2_densities, O2_slant_column_densities, O3_slant_column_densities, & + NO_slant_column_densities) + + do i_level = 1, number_of_vertical_layers + call NO_photolysis_rate_constants_mapping%copy_data( & + jNO(i_level,:), rate_parameters(i_col, i_level, :) ) + end do + end do - errmsg = '' - errcode = 0 + end subroutine calculate_NO_photolysis_rate_constants - if ((index_N2 == MUSICA_INT_UNASSIGNED) .or. & - (index_NO == MUSICA_INT_UNASSIGNED) .or. & - (constituent_props_index_N2 == MUSICA_INT_UNASSIGNED) .or. & - (constituent_props_index_NO == MUSICA_INT_UNASSIGNED) .or. & - (NO_photolysis_indices_constituent_props(1) == MUSICA_INT_UNASSIGNED) .or. & - (NO_photolysis_indices_constituent_props(2) == MUSICA_INT_UNASSIGNED) .or. & - (MOLAR_MASS_N2 == MUSICA_INT_UNASSIGNED) .or. & - (MOLAR_MASS_NO == MUSICA_INT_UNASSIGNED)) then - errmsg = "[MUSICA Error] The index (or indices) or molar mass values & - for NO, N2 have not been initialized." - errcode = 1 - return + !> Finalizes the NO photolysis rate constant calculation + subroutine NO_photolysis_final() + if (associated(no_photolysis_rate_constants_mapping)) then + deallocate( no_photolysis_rate_constants_mapping ) end if + end subroutine NO_photolysis_final - end subroutine check_NO_initialization +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! Internal routines +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Converts mixing ratio to molecule cm-3 subroutine convert_mixing_ratio_to_molecule_cm3(mixing_ratio, & dry_air_density, molar_mass, molecule_cm3) use musica_ccpp_tuvx_gas_species, only: m_3_to_cm_3 + use musica_ccpp_util, only: AVOGADRO real(dk), intent(in) :: mixing_ratio(:) ! kg kg-1 real(dk), intent(in) :: dry_air_density(:) ! kg m-3 real(dk), intent(in) :: molar_mass ! kg mol-1 real(dk), intent(out) :: molecule_cm3(:) ! molecule cm-3 - ! TODO: get from CAM-SIMA - real(dk) :: avogadro = 6.022e23_dk - - molecule_cm3 = mixing_ratio * dry_air_density / molar_mass / m_3_to_cm_3 * avogadro + molecule_cm3 = mixing_ratio * dry_air_density / molar_mass / m_3_to_cm_3 * AVOGADRO end subroutine convert_mixing_ratio_to_molecule_cm3 - !> Prepares to calculate the photolysis rate of NO (nitiric acid). - !! This function transforms the inputs into the units - !! needed by the calculate_jno routine which actually produces the photolysis rate - function calculate_NO_photolysis_rate(solar_zenith_angle, extraterrestrial_flux, & - height_at_interfaces, dry_air_density, constituents_O2, constituents_O3, & - constituents_NO_photolysis) result(jNO) - use musica_ccpp_tuvx_load_species, only: MOLAR_MASS_O2, MOLAR_MASS_O3 - use musica_ccpp_tuvx_gas_species, only: km_to_cm - - real(dk), intent(in) :: solar_zenith_angle ! degrees - real(dk), intent(in) :: extraterrestrial_flux(:) ! photons cm-2 s-1 nm-1 (layer) - real(dk), intent(in) :: height_at_interfaces(:) ! km (layer) - real(dk), intent(in) :: dry_air_density(:) ! kg m-3 (layer) - real(dk), intent(in) :: constituents_O2(:) ! kg kg-1 (layer) - real(dk), intent(in) :: constituents_O3(:) ! kg kg-1 (layer) - real(dk), intent(in) :: constituents_NO_photolysis(:,:) ! kg kg-1 (layer, constituent) - - ! local variables - integer :: number_of_vertical_layers - real(dk) :: N2_densities(size(dry_air_density)+1) ! mol cm-3 - real(dk) :: O2_densities(size(dry_air_density)+1) ! mol cm-3 - real(dk) :: O3_densities(size(dry_air_density)+1) ! mol cm-3 - real(dk) :: NO_densities(size(dry_air_density)+1) ! mol cm-3 - real(dk) :: O2_slant_colunm_densities(size(dry_air_density)+1) ! mol cm-2 - real(dk) :: O3_slant_colunm_densities(size(dry_air_density)+1) ! mol cm-2 - real(dk) :: NO_slant_colunm_densities(size(dry_air_density)+1) ! mol cm-2 - real(dk) :: current_jNO(size(dry_air_density)+1) - - ! parameters needed to calculate slant column densities - ! (see sphers routine description for details) - integer :: number_of_crossed_layers(size(dry_air_density)+1) - real(dk) :: slant_path(0:size(dry_air_density)+1, size(dry_air_density)+1) - real(dk) :: delta_z(size(dry_air_density)+1) ! layer thickness (cm) - real(dk) :: jNO(size(dry_air_density)) ! final photolysis rate - - number_of_vertical_layers = size(dry_air_density) - ! TODO: what are these constants? scale heights? - ! TODO: the values at index 1 appear to be for values above the model top in CAM, but how does that affect cam sima? - call convert_mixing_ratio_to_molecule_cm3(constituents_NO_photolysis(:,index_N2), & - dry_air_density, MOLAR_MASS_N2, N2_densities(2:)) - N2_densities(1) = N2_densities(2) * 0.9_dk - - call convert_mixing_ratio_to_molecule_cm3(constituents_O2, & - dry_air_density, MOLAR_MASS_O2, O2_densities(2:)) - O2_densities(1) = O2_densities(2) * 7.0_dk / ( height_at_interfaces(1) - height_at_interfaces(2) ) - - call convert_mixing_ratio_to_molecule_cm3(constituents_O3, & - dry_air_density, MOLAR_MASS_O3, O3_densities(2:)) - O3_densities(1) = O3_densities(2) * 7.0_dk / ( height_at_interfaces(1) - height_at_interfaces(2) ) - - call convert_mixing_ratio_to_molecule_cm3(constituents_NO_photolysis(:,index_NO), & - dry_air_density, MOLAR_MASS_NO, NO_densities(2:)) - NO_densities(1) = NO_densities(2) * 0.9_dk - - ! calculate slant column densities - call calculate_slant_path( number_of_vertical_layers, height_at_interfaces, & - solar_zenith_angle, slant_path, number_of_crossed_layers ) - delta_z(1:number_of_vertical_layers) = km_to_cm * & - (height_at_interfaces(1:number_of_vertical_layers) & - - height_at_interfaces(2:number_of_vertical_layers+1) ) - call calculate_slant_column_density( number_of_vertical_layers+1, delta_z,& - slant_path, number_of_crossed_layers, O2_densities, O2_slant_colunm_densities ) - call calculate_slant_column_density( number_of_vertical_layers+1, delta_z, & - slant_path, number_of_crossed_layers, O3_densities, O3_slant_colunm_densities ) - call calculate_slant_column_density( number_of_vertical_layers+1, delta_z, & - slant_path, number_of_crossed_layers, NO_densities, NO_slant_colunm_densities ) - - jNO = calculate_jno(number_of_vertical_layers, extraterrestrial_flux, N2_densities, & - O2_slant_colunm_densities, O3_slant_colunm_densities, & - NO_slant_colunm_densities, current_jNO) - - end function calculate_NO_photolysis_rate - !> Calculate the photolysis rate of NO (nitric acid) - function calculate_jno(num_vertical_layers, extraterrestrial_flux, n2_dens, o2_slant, o3_slant, no_slant, work_jno) & + function calculate_jno(num_vertical_layers, extraterrestrial_flux, n2_dens, o2_slant, o3_slant, no_slant) & result(jno) ! inputs @@ -239,12 +226,12 @@ function calculate_jno(num_vertical_layers, extraterrestrial_flux, n2_dens, o2_s real(dk), intent(in) :: o2_slant(:) ! molecule cm-2 real(dk), intent(in) :: o3_slant(:) ! molecule cm-2 real(dk), intent(in) :: no_slant(:) ! molecule cm-2 - real(dk), intent(inout) :: work_jno(:) ! various ! local variables ! NO in an excited electronic state may result in an emission of NO or be quenched by an interaction ! with N2. The predissociation faction is the probability that a precissociation of NO will occur from ! an excited electronic state + real(dk) :: work_jno(size(n2_dens)) ! various real(dk) :: o3_transmission_factor(size(o3_slant), 4) ! Define the correct size real(dk) :: jno(num_vertical_layers) real(dk) :: jno100, jno90, jno50 @@ -507,7 +494,7 @@ subroutine calculate_slant_path( nlev, altitude, zenith_angle_degrees, slant_pat ! ... Dummy arguments !------------------------------------------------------------------------------ - use musica_ccpp_util, only: DEGREE_TO_RADIAN, EARTH_RADIUS_KM + use musica_ccpp_util, only: DEGREE_TO_RADIAN, EARTH_RADIUS_M integer, intent(in) :: nlev ! number model vertical levels integer, intent(out) :: number_of_crossed_layers(0:nlev) ! see above @@ -552,7 +539,7 @@ subroutine calculate_slant_path( nlev, altitude, zenith_angle_degrees, slant_pat !------------------------------------------------------------------------------ ! ... include the elevation above sea level to the radius of the earth: !------------------------------------------------------------------------------ - radius_at_current_height = EARTH_RADIUS_KM + altitude(nlev) + radius_at_current_height = EARTH_RADIUS_M * 0.001_dk + altitude(nlev) !------------------------------------------------------------------------------ ! ... inverse coordinate of z @@ -715,4 +702,4 @@ subroutine calculate_slant_column_density( nlev, delta_z, slant_path, number_of_ slant_column(nlev) = .95_dk*slant_column(nlev-1) end subroutine calculate_slant_column_density -end module musica_ccpp_tuvx_no_photolysis_rate \ No newline at end of file +end module musica_ccpp_tuvx_no_photolysis \ No newline at end of file diff --git a/test/docker/Dockerfile.musica b/test/docker/Dockerfile.musica index d3be2949..c16dbf9a 100644 --- a/test/docker/Dockerfile.musica +++ b/test/docker/Dockerfile.musica @@ -5,8 +5,8 @@ FROM ubuntu:22.04 -ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8 -ARG CAM_SIMA_CHEMISTRY_DATA_TAG=2b58f2410ec7a565bcf80dee16ec20f6bc35d78b +ARG MUSICA_GIT_TAG=869f24102360c9fef54d710101cc62f5c2513745 +ARG CAM_SIMA_CHEMISTRY_DATA_TAG=e9dfa46cdd718c3d0110edfb33f29c4884794f49 ARG BUILD_TYPE=Debug RUN apt update \ diff --git a/test/docker/Dockerfile.musica.no_install b/test/docker/Dockerfile.musica.no_install index 2c9f9714..71c2f8a9 100644 --- a/test/docker/Dockerfile.musica.no_install +++ b/test/docker/Dockerfile.musica.no_install @@ -8,8 +8,8 @@ FROM ubuntu:22.04 -ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8 -ARG CAM_SIMA_CHEMISTRY_DATA_TAG=2b58f2410ec7a565bcf80dee16ec20f6bc35d78b +ARG MUSICA_GIT_TAG=869f24102360c9fef54d710101cc62f5c2513745 +ARG CAM_SIMA_CHEMISTRY_DATA_TAG=e9dfa46cdd718c3d0110edfb33f29c4884794f49 ARG BUILD_TYPE=Debug RUN apt update \ diff --git a/test/musica/data/NO_photolysis/mapping.nojNO.json b/test/musica/data/NO_photolysis/mapping.nojNO.json new file mode 100644 index 00000000..49f3c3da --- /dev/null +++ b/test/musica/data/NO_photolysis/mapping.nojNO.json @@ -0,0 +1,10 @@ +[ + { + "source": "jNO2", + "target": "PHOTO.NO2" + }, + { + "source": "jOther", + "target": "PHOTO.jOther" + } +] \ No newline at end of file diff --git a/test/musica/data/NO_photolysis/mapping.valid.json b/test/musica/data/NO_photolysis/mapping.valid.json new file mode 100644 index 00000000..fc721fa6 --- /dev/null +++ b/test/musica/data/NO_photolysis/mapping.valid.json @@ -0,0 +1,10 @@ +[ + { + "source": "jNO", + "target": "PHOTO.NO" + }, + { + "source": "jOther", + "target": "PHOTO.jOther" + } +] \ No newline at end of file diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index 1cc22d3e..4994b572 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -506,11 +506,13 @@ subroutine test_chapman_with_no_photolysis() use musica_ccpp_namelist, only: filename_of_micm_configuration, & filename_of_tuvx_configuration, & filename_of_tuvx_micm_mapping_configuration + use musica_ccpp_util, only: PI implicit none - integer, parameter :: NUM_SPECIES = 11 - integer, parameter :: NUM_TUVX_CONSTITUENTS = 1 + integer, parameter :: NUM_SPECIES = 10 + 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. @@ -537,9 +539,9 @@ subroutine test_chapman_with_no_photolysis() 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(:) @@ -566,8 +568,8 @@ subroutine test_chapman_with_no_photolysis() temperature(:,2) = (/ 300._kind_phys, 400._kind_phys /) pressure(:,1) = (/ 6000.04_kind_phys, 7000.04_kind_phys /) pressure(:,2) = (/ 8000.04_kind_phys, 9000.04_kind_phys /) - dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) - dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) + dry_air_density(:,1) = (/ 0.95_kind_phys, 0.93_kind_phys /) + dry_air_density(:,2) = (/ 1.225_kind_phys, 1.15_kind_phys /) flux_data_photolysis_wavelength_interfaces(:) = & (/ 200.0_kind_phys, 210.0_kind_phys, 220.0_kind_phys, 230.0_kind_phys, & 240.0_kind_phys, 250.0_kind_phys, 260.0_kind_phys, 270.0_kind_phys, 280.0_kind_phys /) @@ -591,7 +593,7 @@ subroutine test_chapman_with_no_photolysis() 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) allocate(constituent_props_ptr(size(constituent_props))) do i = 1, size(constituent_props) @@ -606,6 +608,8 @@ subroutine test_chapman_with_no_photolysis() stop 3 endif + ASSERT(NUM_SPECIES == micm%species_ordering%size()) + constituents = 0.0_kind_phys do i = 1, micm%species_ordering%size() micm_species_name = micm%species_ordering%name(i) @@ -649,7 +653,7 @@ subroutine test_chapman_with_no_photolysis() endif write(*,*) "[MUSICA INFO] Solved Concentrations" - ! write(*,fmt="(4(3x,e13.6))") constituents + write(*,fmt="(4(3x,e13.6))") constituents call musica_ccpp_final(errmsg, errcode) @@ -662,6 +666,7 @@ subroutine test_chapman_with_no_photolysis() do i = 1, NUM_COLUMNS do j = 1, NUM_LAYERS ! NO should be different + ! TODO: find a more accurate way to check this ASSERT(constituents(i,j,NO_index) .ne. initial_constituents(i,j,NO_index)) end do end do diff --git a/test/musica/tuvx/test_tuvx_gas_species.F90 b/test/musica/tuvx/test_tuvx_gas_species.F90 index ec2374ee..bf1f540e 100644 --- a/test/musica/tuvx/test_tuvx_gas_species.F90 +++ b/test/musica/tuvx/test_tuvx_gas_species.F90 @@ -63,7 +63,8 @@ subroutine test_create_gas_species_profile() 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 + use musica_ccpp_util, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_species, only: tuvx_species_set integer, parameter :: NUM_COLUMNS = 2 integer, parameter :: NUM_LAYERS = 2 @@ -261,11 +262,12 @@ subroutine test_initialize_tuvx_species() 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_util, only: MUSICA_INT_UNASSIGNED 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 + use musica_ccpp_species, only: tuvx_species_set use test_musica_data, only: get_wavelength_edges integer, parameter :: NUM_COLUMNS = 2 diff --git a/test/musica/tuvx/test_tuvx_load_species.F90 b/test/musica/tuvx/test_tuvx_load_species.F90 index 5d34f1e5..df30f000 100644 --- a/test/musica/tuvx/test_tuvx_load_species.F90 +++ b/test/musica/tuvx/test_tuvx_load_species.F90 @@ -10,6 +10,10 @@ program test_tuvx_load_species call test_configure_shared_gas_species_tuvx_micm() call test_configure_partial_shared_gas_species() call test_configure_no_shared_gas_species() + call test_NO_and_N2_absence() + call test_only_NO_exist() + call test_only_N2_exist() + call test_NO_and_N2_exist() contains @@ -20,7 +24,8 @@ subroutine test_configure_shared_gas_species_tuvx_micm() ! 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_ccpp_util, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_species, only: musica_species_t use musica_util, only: error_t integer, parameter :: NUM_MICM_SPECIES = 6 @@ -134,7 +139,8 @@ subroutine test_configure_partial_shared_gas_species() ! 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_ccpp_util, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_species, only: musica_species_t use musica_util, only: error_t integer, parameter :: NUM_MICM_SPECIES = 6 @@ -248,7 +254,8 @@ subroutine test_configure_no_shared_gas_species() ! 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_ccpp_util, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_species, only: musica_species_t use musica_util, only: error_t integer, parameter :: NUM_MICM_SPECIES = 6 @@ -359,4 +366,481 @@ subroutine test_configure_no_shared_gas_species() end subroutine test_configure_no_shared_gas_species + subroutine test_NO_and_N2_absence() + ! This test case applies when NO and N2 are not present in the MICM species. + ! It ensures that the TUV-x species are correctly configured. + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_util, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_species, only: musica_species_t + 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) = 'Ar' + 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) + + ASSERT(index_NO == MUSICA_INT_UNASSIGNED) + ASSERT(index_N2 == MUSICA_INT_UNASSIGNED) + + end subroutine test_NO_and_N2_absence + + subroutine test_only_NO_exist() + ! This test case applies when NO (but not N2) is present in the MICM species. + ! It ensures that the TUV-x species are correctly configured. + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_util, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_species, only: musica_species_t + 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) = 'Ar' + species_names(2) = 'O2' ! shared species + species_names(3) = 'NO' + 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) + + ASSERT(index_NO == MUSICA_INT_UNASSIGNED) + ASSERT(index_N2 == MUSICA_INT_UNASSIGNED) + + end subroutine test_only_NO_exist + + subroutine test_only_N2_exist() + ! This test case applies when N2 (but not NO) is present in the MICM species. + ! It ensures that the TUV-x species are correctly configured. + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_util, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_species, only: musica_species_t + 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) + + ASSERT(index_NO == MUSICA_INT_UNASSIGNED) + ASSERT(index_N2 == MUSICA_INT_UNASSIGNED) + + end subroutine test_only_N2_exist + + subroutine test_NO_and_N2_exist() + ! This test case applies when NO and N2 are present in the MICM species. + ! It ensures that the TUV-x species are correctly configured. + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_util, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_species, only: musica_species_t + use musica_util, only: error_t + + integer, parameter :: NUM_MICM_SPECIES = 6 + integer, parameter :: NUM_SHARED_SPECIES_BETWEEN_MICM_TUVX = 4 + integer, parameter :: NUM_TUVX_CONSTITUENTS = 6 + 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) = 'NO' + 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) .or. & + (trim(name) == 'NO' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0300061_kind_phys .and. & + scale_height == 0.0_kind_phys .and. index_musica == 5 .and. index_constituent_props == MUSICA_INT_UNASSIGNED & + .and. .not. has_profile) .or. & + (trim(name) == 'N2' .and. trim(unit) == 'molecule cm-3' .and. molar_mass == 0.0280134_kind_phys .and. & + scale_height == 0.0_kind_phys .and. index_musica == 6 .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) + ASSERT(index_NO == 5) + ASSERT(index_N2 == 6) + + 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_NO_and_N2_exist + end program test_tuvx_load_species \ No newline at end of file diff --git a/test/musica/tuvx/test_tuvx_no_photolysis.F90 b/test/musica/tuvx/test_tuvx_no_photolysis.F90 index 4513e115..b17cc02e 100644 --- a/test/musica/tuvx/test_tuvx_no_photolysis.F90 +++ b/test/musica/tuvx/test_tuvx_no_photolysis.F90 @@ -1,202 +1,34 @@ program test_tuvx_no_photolysis - use musica_ccpp_tuvx_no_photolysis_rate + use musica_ccpp_tuvx_no_photolysis 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_NO_and_N2_absence() - call test_only_NO_exist() - call test_only_N2_exist() - call test_NO_and_N2_exist() - call test_initialize_NO_constituents_indices_and_molar_mass() call test_calculate_NO_photolysis_rate() contains - subroutine test_NO_and_N2_absence() - use ccpp_kinds, only: kind_phys - use musica_ccpp_species, only: musica_species_t - - integer, parameter :: NUM_MICM_SPECIES = 3 - type(musica_species_t) :: micm_species(NUM_MICM_SPECIES) - real(kind_phys) :: molar_mass_group(NUM_MICM_SPECIES) = & - [0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys] - integer :: i_species - character(len=512) :: species_names(NUM_MICM_SPECIES) - - species_names(1) = 'O2' - species_names(2) = 'O1D' - species_names(3) = 'O3' - - do i_species = 1, NUM_MICM_SPECIES - 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 check_NO_exist( micm_species, .true. ) - ASSERT(.not. is_NO_photolysis_active) - - end subroutine test_NO_and_N2_absence - - subroutine test_only_NO_exist() - use ccpp_kinds, only: kind_phys - use musica_ccpp_species, only: musica_species_t - - integer, parameter :: NUM_MICM_SPECIES = 3 - type(musica_species_t) :: micm_species(NUM_MICM_SPECIES) - real(kind_phys) :: molar_mass_group(NUM_MICM_SPECIES) = & - [0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys] - integer :: i_species - character(len=512) :: species_names(NUM_MICM_SPECIES) - - species_names(1) = 'NO' - species_names(2) = 'O1D' - species_names(3) = 'O3' - - do i_species = 1, NUM_MICM_SPECIES - 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 check_NO_exist( micm_species, .true. ) - ASSERT(.not. is_NO_photolysis_active) - - end subroutine test_only_NO_exist - - subroutine test_only_N2_exist() - use ccpp_kinds, only: kind_phys - use musica_ccpp_species, only: musica_species_t - - integer, parameter :: NUM_MICM_SPECIES = 3 - type(musica_species_t) :: micm_species(NUM_MICM_SPECIES) - real(kind_phys) :: molar_mass_group(NUM_MICM_SPECIES) = & - [0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys] - integer :: i_species - character(len=512) :: species_names(NUM_MICM_SPECIES) - - species_names(1) = 'N2' - species_names(2) = 'O1D' - species_names(3) = 'O3' - - do i_species = 1, NUM_MICM_SPECIES - 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 check_NO_exist( micm_species, .true. ) - ASSERT(.not. is_NO_photolysis_active) - - end subroutine test_only_N2_exist - - subroutine test_NO_and_N2_exist() - use ccpp_kinds, only: kind_phys - use musica_ccpp_species, only: musica_species_t - - integer, parameter :: NUM_MICM_SPECIES = 3 - type(musica_species_t) :: micm_species(NUM_MICM_SPECIES) - real(kind_phys) :: molar_mass_group(NUM_MICM_SPECIES) = & - [0.0280134_kind_phys, 0.2_kind_phys, 0.0300100_kind_phys] - integer :: i_species - character(len=512) :: species_names(NUM_MICM_SPECIES) - - species_names(1) = 'N2' - species_names(2) = 'O2' - species_names(3) = 'NO' - - do i_species = 1, NUM_MICM_SPECIES - 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 check_NO_exist( micm_species, .true. ) - ASSERT(is_NO_photolysis_active) - ASSERT(MOLAR_MASS_N2 == 0.0280134_kind_phys) - ASSERT(MOLAR_MASS_NO == 0.0300100_kind_phys) - - end subroutine test_NO_and_N2_exist - - subroutine test_initialize_NO_constituents_indices_and_molar_mass() - use ccpp_kinds, only: kind_phys - use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t, ccpp_constituent_prop_ptr_t - use musica_ccpp, only: musica_ccpp_register, musica_ccpp_final - use musica_ccpp_namelist, only: filename_of_micm_configuration, & - filename_of_tuvx_configuration, & - filename_of_tuvx_micm_mapping_configuration - - integer, parameter :: expected_index_N2 = 1 - integer, parameter :: expected_index_NO = 2 - real(kind_phys), parameter :: expected_MOLAR_MASS_N2 = 0.0280134_kind_phys - real(kind_phys), parameter :: expected_MOLAR_MASS_NO = 0.0300100_kind_phys - 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/NO_photolysis/micm/config.json' - filename_of_tuvx_configuration = 'musica_configurations/NO_photolysis/tuvx/config.json' - filename_of_tuvx_micm_mapping_configuration = 'musica_configurations/NO_photolysis/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)) - ASSERT(is_NO_photolysis_active) - ASSERT(MOLAR_MASS_N2 == expected_MOLAR_MASS_N2) - ASSERT(MOLAR_MASS_NO == expected_MOLAR_MASS_NO) - - 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 set_NO_index_constituent_props(constituent_props_ptr, errmsg, errcode) - ASSERT(errcode == 0) - ASSERT(index_N2 == expected_index_N2) - ASSERT(index_NO == expected_index_NO) - ASSERT(NO_photolysis_indices_constituent_props(1) == constituent_props_index_N2) - ASSERT(NO_photolysis_indices_constituent_props(2) == constituent_props_index_NO) - - call check_NO_initialization(errmsg, errcode) - ASSERT(errcode == 0) - - call musica_ccpp_final(errmsg, errcode) - ASSERT(errcode == 0) - - end subroutine test_initialize_NO_constituents_indices_and_molar_mass - subroutine test_calculate_NO_photolysis_rate() use ccpp_kinds, only: kind_phys use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t - use musica_ccpp_micm, only: micm - use musica_ccpp, only: musica_ccpp_register, musica_ccpp_init, musica_ccpp_final + use musica_ccpp, only: musica_ccpp_register, musica_ccpp_init, musica_ccpp_final 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_O2, index_O3 - use musica_ccpp_species, only: tuvx_indices_constituent_props, extract_subset_constituents - use musica_ccpp_tuvx_no_photolysis_rate, only: NO_photolysis_indices_constituent_props + use musica_ccpp_tuvx_load_species, only: index_cloud_liquid_water_content, index_dry_air, & + tuvx_index_O2 => index_O2, & + tuvx_index_O3 => index_O3, & + tuvx_index_NO => index_NO, & + tuvx_index_N2 => index_N2 + use musica_ccpp_species, only: tuvx_indices_constituent_props, extract_subset_constituents, & + initialize_musica_species_indices, initialize_molar_mass_array, & + check_initialization use test_musica_data, only: get_wavelength_edges + use musica_ccpp_util, only: PI implicit none @@ -211,6 +43,8 @@ subroutine test_calculate_NO_photolysis_rate() integer, parameter :: NUM_COLUMNS = 2 integer, parameter :: NUM_LAYERS = 2 integer, parameter :: NUM_WAVELENGTH_BINS = 102 + integer, parameter :: NUM_MICM_RATE_PARAMETERS = 4 + integer, parameter :: INDEX_MICM_JNO = 3 integer :: errcode character(len=512) :: errmsg real(kind_phys), dimension(NUM_WAVELENGTH_BINS+1) :: photolysis_wavelength_grid_interfaces ! m @@ -219,34 +53,33 @@ subroutine test_calculate_NO_photolysis_rate() real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: dry_air_density ! kg m-3 real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & NUM_SPECIES+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) :: constituents ! kg kg-1 - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & - NUM_TUVX_CONSTITUENTS_INCLUDE_SHARED_SPECIES) :: constituents_tuvx_species ! kg kg-1 - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & - NUM_NO_CONSTITUENTS) :: constituents_NO_photolysis ! kg kg-1 real(kind_phys), dimension(NUM_COLUMNS) :: solar_zenith_angle ! radians - real(kind_phys), dimension(NUM_LAYERS+1) :: height_at_interfaces ! km - real(kind_phys), dimension(NUM_LAYERS) :: NO_photolysis_rate_constant(NUM_LAYERS) + real(kind_phys), dimension(NUM_COLUMNS, NUM_LAYERS+1) :: height_at_interfaces ! km + real(kind_phys), dimension(NUM_COLUMNS, NUM_LAYERS, & + NUM_MICM_RATE_PARAMETERS) :: rate_parameters type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:) type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:) type(ccpp_constituent_properties_t), pointer :: const_prop real(kind_phys) :: base_conc - character(len=:), allocatable :: micm_species_name + character(len=256) :: constituent_name integer :: i, j, k, i_col - integer :: O2_index, Ar_index, CO2_index, H2O_index, N_index - integer :: NO_index, O_index, O1D_index, O3_index, N2_index + integer :: Ar_index, CO2_index, H2O_index, N_index + integer :: O_index, O1D_index, O2_index, O3_index + integer :: N2_index, NO_index, cloud_index, air_index filename_of_micm_configuration = 'musica_configurations/NO_photolysis/micm/config.json' filename_of_tuvx_configuration = 'musica_configurations/NO_photolysis/tuvx/config.json' filename_of_tuvx_micm_mapping_configuration = 'musica_configurations/NO_photolysis/tuvx_micm_mapping.json' call get_wavelength_edges(photolysis_wavelength_grid_interfaces) - dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) - dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) + dry_air_density(:,1) = (/ 0.95_kind_phys, 0.93_kind_phys /) + dry_air_density(:,2) = (/ 1.2255_kind_phys, 1.15_kind_phys /) extraterrestrial_flux(:) = & (/ 1.5e13_kind_phys, 1.5e13_kind_phys, 1.4e13_kind_phys, 1.4e13_kind_phys, & 1.3e13_kind_phys, 1.2e13_kind_phys, 1.1e13_kind_phys, 1.0e13_kind_phys /) solar_zenith_angle(:) = (/ 0.0_kind_phys, 2.1_kind_phys /) - height_at_interfaces(:) = (/ 3.01_kind_phys, 1.01_kind_phys, 0.01_kind_phys /) + height_at_interfaces(1,:) = (/ 3.01_kind_phys, 1.01_kind_phys, 0.01_kind_phys /) + height_at_interfaces(2,:) = (/ 6.01_kind_phys, 1.01_kind_phys, 0.51_kind_phys /) call musica_ccpp_register( constituent_props, errmsg, errcode ) ASSERT(errcode == 0) @@ -258,83 +91,93 @@ subroutine test_calculate_NO_photolysis_rate() ASSERT(errcode == 0) end do - call musica_ccpp_init( NUM_COLUMNS, NUM_LAYERS, NUM_LAYERS+1, photolysis_wavelength_grid_interfaces, & - constituent_props_ptr, errmsg, errcode ) - ASSERT(errcode == 0) - - do i = 1, micm%species_ordering%size() - micm_species_name = micm%species_ordering%name(i) - if (micm_species_name == "O2") then + do i = 1, size(constituent_props) + call constituent_props(i)%standard_name(constituent_name) + if (constituent_name == "O2") then O2_index = i base_conc = 0.21_kind_phys - else if (micm_species_name == "Ar") then + else if (constituent_name == "Ar") then Ar_index = i base_conc = 1.0e-9_kind_phys - else if (micm_species_name == "CO2") then + else if (constituent_name == "CO2") then CO2_index = i base_conc = 1.0e-9_kind_phys - else if (micm_species_name == "H2O") then + else if (constituent_name == "H2O") then H2O_index = i base_conc = 1.0e-9_kind_phys - else if (micm_species_name == "N") then + else if (constituent_name == "N") then N_index = i base_conc = 1.0e-9_kind_phys - else if (micm_species_name == "NO") then + else if (constituent_name == "NO") then NO_index = i base_conc = 1.0e-9_kind_phys - else if (micm_species_name == "O") then + else if (constituent_name == "O") then O_index = i base_conc = 1.0e-9_kind_phys - else if (micm_species_name == "O1D") then + else if (constituent_name == "O1D") then O1D_index = i base_conc = 1.0e-9_kind_phys - else if (micm_species_name == "O3") then + else if (constituent_name == "O3") then O3_index = i base_conc = 1.0e-4_kind_phys - else if (micm_species_name == "N2") then + else if (constituent_name == "N2") then N2_index = i base_conc = 0.79_kind_phys + else if (constituent_name == "cloud_liquid_water_mixing_ratio_wrt_moist_air_and_condensed_water") then + cloud_index = i + base_conc = 1.0e-3_kind_phys + else if (constituent_name == "air") then + air_index = i + base_conc = 1.2_kind_phys else - write(*,*) "Unknown species: ", micm_species_name + write(*,*) "Unknown species: ", constituent_name stop 3 endif + end do + + call musica_ccpp_init( NUM_COLUMNS, NUM_LAYERS, NUM_LAYERS+1, photolysis_wavelength_grid_interfaces, & + constituent_props_ptr, errmsg, errcode ) + ASSERT(errcode == 0) + + ! check species indices + ASSERT( O2_index == tuvx_indices_constituent_props( tuvx_index_O2 ) ) + ASSERT( O3_index == tuvx_indices_constituent_props( tuvx_index_O3 ) ) + ASSERT( NO_index == tuvx_indices_constituent_props( tuvx_index_NO ) ) + ASSERT( N2_index == tuvx_indices_constituent_props( tuvx_index_N2 ) ) + ASSERT( cloud_index == tuvx_indices_constituent_props( index_cloud_liquid_water_content ) ) + ASSERT( air_index == tuvx_indices_constituent_props( index_dry_air ) ) + + ! Set initial species concentrations + do i = 1, size(constituent_props) do j = 1, NUM_COLUMNS do k = 1, NUM_LAYERS constituents(j,k,i) = base_conc * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) end do end do end do - ! set initial cloud liquid water mixing ratio to ~1e-3 kg kg-1 - do j = 1, NUM_COLUMNS - do k = 1, NUM_LAYERS - constituents(j,k,NUM_SPECIES+NUM_TUVX_CONSTITUENTS) = 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+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) = & - 1.0e-3_kind_phys * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) - end do - end do - - call extract_subset_constituents(tuvx_indices_constituent_props, constituents, & - constituents_tuvx_species, errmsg, errcode) - ASSERT(errcode == 0) - - call extract_subset_constituents(NO_photolysis_indices_constituent_props, constituents, & - constituents_NO_photolysis, errmsg, errcode) - ASSERT(errcode == 0) + rate_parameters = 0.0_kind_phys + call calculate_NO_photolysis_rate_constants( & + solar_zenith_angle, extraterrestrial_flux, & + height_at_interfaces, dry_air_density, & + constituents, rate_parameters ) do i_col = 1, size(dry_air_density, dim=1) - NO_photolysis_rate_constant = calculate_NO_photolysis_rate( & - solar_zenith_angle(i_col), extraterrestrial_flux, & - height_at_interfaces, dry_air_density(i_col,:), & - constituents_tuvx_species(i_col,:,index_O2), & - constituents_tuvx_species(i_col,:,index_O3), & - constituents_NO_photolysis(i_col,:,:) ) + ! TODO: figure out why the uppermost layer has a jNO of 0.0 + do i = 1, NUM_LAYERS-1 + do j = 1, NUM_MICM_RATE_PARAMETERS + if (j == INDEX_MICM_JNO) then + if (solar_zenith_angle(i_col) <= 110.0_kind_phys / 180.0_kind_phys * PI) then + ASSERT( rate_parameters(i_col,i,j) > 0.0_kind_phys ) + else + ASSERT( rate_parameters(i_col,i,j) == 0.0_kind_phys ) + end if + else + ASSERT( rate_parameters(i_col,i,j) == 0.0_kind_phys ) + end if + end do + end do end do - call musica_ccpp_final(errmsg, errcode) ASSERT(errcode == 0)