From 4d977b45497e56539df6b16910b2adf5c3c782e0 Mon Sep 17 00:00:00 2001 From: Jiwon Gim Date: Wed, 22 Jan 2025 20:27:43 -0700 Subject: [PATCH] add function for no photolysis calculation --- schemes/musica/musica_ccpp.F90 | 42 +++- schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 20 +- .../tuvx/musica_ccpp_tuvx_load_species.F90 | 68 +++---- .../musica_ccpp_tuvx_no_photolysis_rate.F90 | 185 ++++++++++++++++-- test/musica/test_musica_species.F90 | 2 +- test/musica/tuvx/CMakeLists.txt | 9 +- test/musica/tuvx/test_tuvx_gas_species.F90 | 4 - test/musica/tuvx/test_tuvx_no_photolysis.F90 | 116 ++++++----- 8 files changed, 322 insertions(+), 124 deletions(-) diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index 3a9c6512..8d6ccb2d 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -14,10 +14,11 @@ 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 - use musica_ccpp_species, only: musica_species_t, register_musica_species - use musica_ccpp_tuvx_load_species, only: check_tuvx_species_initialization + 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 + 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 @@ -49,6 +50,8 @@ 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) + end subroutine musica_ccpp_register !> \section arg_table_musica_ccpp_init Argument Table @@ -64,6 +67,8 @@ 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_present, set_NO_index_constituent_props, check_NO_initialization integer, intent(in) :: horizontal_dimension ! (count) integer, intent(in) :: vertical_layer_dimension ! (count) @@ -98,6 +103,13 @@ subroutine musica_ccpp_init(horizontal_dimension, vertical_layer_dimension, & call check_initialization(errmsg, errcode) if (errcode /= 0) return + if (is_NO_present) 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 @@ -121,6 +133,8 @@ 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_present, 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 @@ -154,14 +168,24 @@ 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_present) 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, & @@ -176,9 +200,19 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co rate_parameters, & errmsg, errcode) + ! TODO (jiwon) - update constituents doens't have an effect since tuvx_run takes it as intent in + ! is this correct? call update_constituents(tuvx_indices_constituent_props, constituents_tuvx_species, & constituents, errmsg, errcode) if (errcode /= 0) return + + if (is_NO_present) then + ! TODO(jiwon) - does it need to update NO constituents? + call update_constituents(NO_photolysis_indices_constituent_props, constituents_NO_photolysis, & + constituents, errmsg, errcode) + if (errcode /= 0) return + end if + call extract_subset_constituents(micm_indices_constituent_props, constituents, & constituents_micm_species, errmsg, errcode) if (errcode /= 0) return diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index c6c2a419..a9f58329 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -503,6 +503,7 @@ 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, & @@ -525,16 +526,16 @@ subroutine tuvx_run(temperature, dry_air_density, & 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, & - index_N2, index_NO, MOLAR_MASS_N2, & - MOLAR_MASS_O2, MOLAR_MASS_O3, MOLAR_MASS_NO + 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: calculate_NO_photolysis_rate use musica_ccpp_species, only: MUSICA_INT_UNASSIGNED + use musica_ccpp_tuvx_no_photolysis_rate, only: is_NO_present 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 @@ -562,7 +563,7 @@ subroutine tuvx_run(temperature, dry_air_density, & real(kind_phys) :: solar_zenith_angle_degrees type(error_t) :: error integer :: i_col, i_level - real(kind_phys) :: jno(size(constituents, dim=2)) ! s-1 + real(kind_phys) :: jno(size(constituents, dim=2)) ! s-1 reciprocal_of_gravitational_acceleration = 1.0_kind_phys / standard_gravitational_acceleration @@ -631,11 +632,14 @@ subroutine tuvx_run(temperature, dry_air_density, & max( photolysis_rate_constants(:,:), 0.0_kind_phys ) end if ! solar zenith angle check - if (index_N2 .ne. MUSICA_INT_UNASSIGNED .and. index_O2 .ne. MUSICA_INT_UNASSIGNED .and. & - index_O3 .ne. MUSICA_INT_UNASSIGNED .and. index_NO .ne. MUSICA_INT_UNASSIGNED) then + if (is_NO_present) then ! TODO: How do I ensure that the extraterrestrial flux matches the wavelength grid needed for NO photolysis? - jno = calculate_NO_photolysis_rate(size(constituents, dim=2), solar_zenith_angle(i_col), extraterrestrial_flux, constituents(i_col,:,:), height_interfaces, & - dry_air_density(i_col,:), index_N2, index_O2, index_O3, index_NO, MOLAR_MASS_N2, MOLAR_MASS_O2, MOLAR_MASS_O3, MOLAR_MASS_NO) + ! jno = calculate_NO_photolysis_rate(size(constituents, dim=2), solar_zenith_angle(i_col), extraterrestrial_flux, constituents(i_col,:,:), height_interfaces, & + ! dry_air_density(i_col,:), index_N2, index_O2, index_O3, index_NO, MOLAR_MASS_N2, MOLAR_MASS_O2, MOLAR_MASS_O3, MOLAR_MASS_NO) + jno = calculate_NO_photolysis_rate( size(constituents, dim=2), & + 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? diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 index 7790f941..4e1fab78 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_load_species.F90 @@ -11,8 +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_N2 = MUSICA_INT_UNASSIGNED - integer, protected, public :: index_NO = MUSICA_INT_UNASSIGNED + ! integer, protected, public :: index_N2 = MUSICA_INT_UNASSIGNED + ! integer, protected, public :: index_NO = MUSICA_INT_UNASSIGNED ! Constants ! Cloud liquid water @@ -34,9 +34,6 @@ 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_N2 = 0.0280134_kind_phys ! kg mol-1 - real(kind_phys), parameter, public :: MOLAR_MASS_NO = 0.0300100_kind_phys ! kg mol-1 - contains @@ -65,8 +62,8 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props, logical :: is_dry_air_registered = .false. logical :: is_O2_registered = .false. logical :: is_O3_registered = .false. - logical :: is_N2_registered = .false. - logical :: is_NO_registered = .false. + ! logical :: is_N2_registered = .false. + ! logical :: is_NO_registered = .false. integer :: i_new, i_species, i_tuvx_species num_micm_species = size(micm_species) @@ -92,8 +89,7 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props, ! 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 .and. & - is_N2_registered .and. is_NO_registered) exit + 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. @@ -107,12 +103,12 @@ 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 == 'N2' ) then - num_new_species = num_new_species + 1 - is_N2_registered = .true. - else if ( micm_species(i_species)%name == 'NO' ) then - num_new_species = num_new_species + 1 - is_NO_registered = .true. + ! else if ( micm_species(i_species)%name == 'N2' ) then + ! num_new_species = num_new_species + 1 + ! is_N2_registered = .true. + ! else if ( micm_species(i_species)%name == 'NO' ) then + ! num_new_species = num_new_species + 1 + ! is_NO_registered = .true. end if end do @@ -206,27 +202,27 @@ subroutine configure_tuvx_species(micm_species, tuvx_species, constituent_props, profiled = .true., & scale_height = SCALE_HEIGHT_O3 ) - if (is_NO_registered) then - i_tuvx_species = i_tuvx_species + 1 - index_NO = i_tuvx_species - tuvx_species(i_tuvx_species) = musica_species_t( & - name = 'NO', & - unit = TUVX_GAS_SPECIES_UNITS, & - molar_mass = MOLAR_MASS_NO, & ! kg mol-1 - index_musica_species = i_tuvx_species, & - profiled = .true.) - end if - - if (is_N2_registered) then - i_tuvx_species = i_tuvx_species + 1 - index_N2 = i_tuvx_species - tuvx_species(i_tuvx_species) = musica_species_t( & - name = 'N2', & - unit = TUVX_GAS_SPECIES_UNITS, & - molar_mass = MOLAR_MASS_N2, & ! kg mol-1 - index_musica_species = i_tuvx_species, & - profiled = .true.) - end if + ! if (is_NO_registered) then + ! i_tuvx_species = i_tuvx_species + 1 + ! index_NO = i_tuvx_species + ! tuvx_species(i_tuvx_species) = musica_species_t( & + ! name = 'NO', & + ! unit = TUVX_GAS_SPECIES_UNITS, & + ! molar_mass = MOLAR_MASS_NO, & ! kg mol-1 + ! index_musica_species = i_tuvx_species, & + ! profiled = .true.) + ! end if + + ! if (is_N2_registered) then + ! i_tuvx_species = i_tuvx_species + 1 + ! index_N2 = i_tuvx_species + ! tuvx_species(i_tuvx_species) = musica_species_t( & + ! name = 'N2', & + ! unit = TUVX_GAS_SPECIES_UNITS, & + ! molar_mass = MOLAR_MASS_N2, & ! kg mol-1 + ! index_musica_species = i_tuvx_species, & + ! profiled = .true.) + ! 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_rate.F90 index 100be1bf..d3d6439c 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_no_photolysis_rate.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_no_photolysis_rate.F90 @@ -1,6 +1,7 @@ module musica_ccpp_tuvx_no_photolysis_rate use ccpp_kinds, only: dk => kind_phys + use musica_ccpp_species, only: MUSICA_INT_UNASSIGNED implicit none @@ -21,10 +22,103 @@ module musica_ccpp_tuvx_no_photolysis_rate ! MS93: A reference to the Minschwaner and Siskind paper private - public :: calculate_NO_photolysis_rate + public :: check_NO_exist, set_NO_index_constituent_props, check_NO_initialization, & + calculate_NO_photolysis_rate + + character(len=*), parameter, public :: N2_label = 'N2' + character(len=*), parameter, public :: NO_label = 'NO' + logical, protected, public :: is_NO_present = .false. + integer, parameter, public :: num_NO_photolysis_constituents = 2 + 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 contains + subroutine check_NO_exist(micm_species) + use musica_ccpp_species, only: musica_species_t + + type(musica_species_t), intent(in) :: micm_species(:) + + ! local variables + integer :: num_micm_species + logical :: is_N2_registered = .false. + logical :: is_NO_registered = .false. + integer :: i_species + + num_micm_species = size(micm_species) + + do i_species = 1, num_micm_species + if ( is_N2_registered .and. is_NO_registered ) then + is_NO_present = .true. + exit + end if + + if ( micm_species(i_species)%name == 'N2' ) then + is_N2_registered = .true. + MOLAR_MASS_N2 = micm_species(i_species)%molar_mass + else if ( micm_species(i_species)%name == 'NO' ) then + is_NO_registered = .true. + MOLAR_MASS_NO = micm_species(i_species)%molar_mass + end if + end do + + if (is_N2_registered .and. is_NO_registered) then + is_NO_present = .true. + end if + + end subroutine check_NO_exist + + 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 + + type(ccpp_constituent_prop_ptr_t), intent(in) :: constituent_props(:) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + call ccpp_const_get_idx( constituent_props, N2_label, & + constituent_props_index_N2, errmsg, errcode ) + if (errcode /= 0) return + index_N2 = 1 + + call ccpp_const_get_idx( constituent_props, NO_label, & + constituent_props_index_NO, errmsg, errcode ) + if (errcode /= 0) return + index_NO = index_N2 + 1 + + NO_photolysis_indices_constituent_props = [constituent_props_index_N2, constituent_props_index_NO] + + end subroutine set_NO_index_constituent_props + + subroutine check_NO_initialization(errmsg, errcode) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + errmsg = '' + errcode = 0 + + 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 + end if + + end subroutine check_NO_initialization + !> Convert mixing ratio to molecule cm-3 subroutine convert_mixing_ratio_to_molecule_cm3(mixing_ratio, dry_air_density, molar_mass, molecule_cm3) real(dk), intent(in) :: mixing_ratio(:) ! kg kg-1 @@ -41,18 +135,77 @@ end subroutine convert_mixing_ratio_to_molecule_cm3 !> Prepare 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(number_of_vertical_layers, solar_zenith_angle, extraterrestrial_flux, constituents, height_at_interfaces, & - dry_air_density, N2_index, O2_index, O3_index, NO_index, molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO) & - result(jNO) + ! function calculate_NO_photolysis_rate(number_of_vertical_layers, solar_zenith_angle, extraterrestrial_flux, constituents, height_at_interfaces, & + ! dry_air_density, N2_index, O2_index, O3_index, NO_index, molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO) & + ! result(jNO) + ! ! inputs + ! integer, intent(in) :: number_of_vertical_layers! unitless + ! real(dk), intent(in) :: solar_zenith_angle ! degrees + ! real(dk), intent(in) :: extraterrestrial_flux(:) ! photons cm-2 s-1 nm-1 (layer) + ! real(dk), target, intent(in) :: constituents(:,:) ! various, but should be mixing ratio (kg kg-1)? (layer, constituent) + ! real(dk), intent(in) :: height_at_interfaces(:) ! km (layer) + ! real(dk), intent(in) :: dry_air_density(:) ! kg m-3 (layer) + ! integer, intent(in) :: N2_index, O2_index, O3_index, NO_index ! position of these species in the constituent arrays + ! real(dk), intent(in) :: molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO ! kg mol-1 + + ! ! local variables + ! real(dk) :: n2_dens(number_of_vertical_layers+1), o2_dens(number_of_vertical_layers+1) + ! real(dk) :: o3_dens(number_of_vertical_layers+1), no_dens(number_of_vertical_layers+1) + ! ! species slant column densities (molecule cm-2) + ! real(dk) :: o2_slant(number_of_vertical_layers+1), o3_slant(number_of_vertical_layers+1) + ! real(dk) :: no_slant(number_of_vertical_layers+1) + ! real(dk) :: work_jno(number_of_vertical_layers+1) ! working photo rate array + + ! ! parameters needed to calculate slant column densities + ! ! (see sphers routine description for details) + ! integer :: number_of_crossed_layers(number_of_vertical_layers+1) + ! real(dk) :: slant_path(0:number_of_vertical_layers+1,number_of_vertical_layers+1) + ! real(dk) :: delta_z(number_of_vertical_layers+1) ! layer thickness (cm) + ! real(dk), parameter :: km2cm = 1.0e5_dk ! conversion from km to cm + ! real(dk) :: jNO(number_of_vertical_layers) ! final photolysis rate + + ! ! 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(:,N2_index), dry_air_density, molar_mass_N2, n2_dens(2:)) + ! n2_dens(1) = n2_dens(2) * 0.9_dk + + ! call convert_mixing_ratio_to_molecule_cm3(constituents(:,O2_index), dry_air_density, molar_mass_O2, o2_dens(2:)) + ! o2_dens(1) = o2_dens(2) * 7.0_dk / ( height_at_interfaces(1) - height_at_interfaces(2) ) + + ! call convert_mixing_ratio_to_molecule_cm3(constituents(:,O3_index), dry_air_density, molar_mass_O3, o3_dens(2:)) + ! o3_dens(1) = o3_dens(2) * 7.0_dk / ( height_at_interfaces(1) - height_at_interfaces(2) ) + + ! call convert_mixing_ratio_to_molecule_cm3(constituents(:,NO_index), dry_air_density, molar_mass_NO, no_dens(2:)) + ! no_dens(1) = no_dens(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) = km2cm * ( 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_dens, o2_slant ) + ! call calculate_slant_column_density( number_of_vertical_layers+1, delta_z, slant_path, number_of_crossed_layers, o3_dens, o3_slant ) + ! call calculate_slant_column_density( number_of_vertical_layers+1, delta_z, slant_path, number_of_crossed_layers, no_dens, no_slant ) + + ! jNO = calculate_jno(number_of_vertical_layers, extraterrestrial_flux, n2_dens, o2_slant, o3_slant, no_slant, work_jno) + + ! end function calculate_NO_photolysis_rate + function calculate_NO_photolysis_rate(number_of_vertical_layers, & + 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 + ! inputs - integer, intent(in) :: number_of_vertical_layers! unitless - real(dk), intent(in) :: solar_zenith_angle ! degrees - real(dk), intent(in) :: extraterrestrial_flux(:) ! photons cm-2 s-1 nm-1 (layer) - real(dk), target, intent(in) :: constituents(:,:) ! various, but should be mixing ratio (kg kg-1)? (layer, constituent) - real(dk), intent(in) :: height_at_interfaces(:) ! km (layer) - real(dk), intent(in) :: dry_air_density(:) ! kg m-3 (layer) - integer, intent(in) :: N2_index, O2_index, O3_index, NO_index ! position of these species in the constituent arrays - real(dk), intent(in) :: molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO ! kg mol-1 + integer, intent(in) :: number_of_vertical_layers ! unitless + 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 + real(dk), intent(in) :: constituents_O3(:) ! kg kg-1 + real(dk), intent(in) :: constituents_NO_photolysis(:,:) ! various, but should be mixing ratio (kg kg-1)? (layer, constituent) ! local variables real(dk) :: n2_dens(number_of_vertical_layers+1), o2_dens(number_of_vertical_layers+1) @@ -72,16 +225,16 @@ function calculate_NO_photolysis_rate(number_of_vertical_layers, solar_zenith_an ! 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(:,N2_index), dry_air_density, molar_mass_N2, n2_dens(2:)) + call convert_mixing_ratio_to_molecule_cm3(constituents_NO_photolysis(:,index_N2), dry_air_density, MOLAR_MASS_N2, n2_dens(2:)) n2_dens(1) = n2_dens(2) * 0.9_dk - call convert_mixing_ratio_to_molecule_cm3(constituents(:,O2_index), dry_air_density, molar_mass_O2, o2_dens(2:)) + call convert_mixing_ratio_to_molecule_cm3(constituents_O2, dry_air_density, MOLAR_MASS_O2, o2_dens(2:)) o2_dens(1) = o2_dens(2) * 7.0_dk / ( height_at_interfaces(1) - height_at_interfaces(2) ) - call convert_mixing_ratio_to_molecule_cm3(constituents(:,O3_index), dry_air_density, molar_mass_O3, o3_dens(2:)) + call convert_mixing_ratio_to_molecule_cm3(constituents_O3, dry_air_density, MOLAR_MASS_O3, o3_dens(2:)) o3_dens(1) = o3_dens(2) * 7.0_dk / ( height_at_interfaces(1) - height_at_interfaces(2) ) - call convert_mixing_ratio_to_molecule_cm3(constituents(:,NO_index), dry_air_density, molar_mass_NO, no_dens(2:)) + call convert_mixing_ratio_to_molecule_cm3(constituents_NO_photolysis(:,index_NO), dry_air_density, MOLAR_MASS_NO, no_dens(2:)) no_dens(1) = no_dens(2) * 0.9_dk ! ================================ diff --git a/test/musica/test_musica_species.F90 b/test/musica/test_musica_species.F90 index 2b5410bd..579031de 100644 --- a/test/musica/test_musica_species.F90 +++ b/test/musica/test_musica_species.F90 @@ -21,7 +21,7 @@ subroutine test_register_musica_species() check_tuvx_species_initialization integer, parameter :: NUM_MICM_SPECIES = 6 - integer, parameter :: NUM_TUVX_CONSTITUENTS = 5 + 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) diff --git a/test/musica/tuvx/CMakeLists.txt b/test/musica/tuvx/CMakeLists.txt index d041570b..a096d2c4 100644 --- a/test/musica/tuvx/CMakeLists.txt +++ b/test/musica/tuvx/CMakeLists.txt @@ -270,13 +270,20 @@ add_test( add_memory_check_test(test_tuvx_load_species $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) - # NO photolysis rate add_executable(test_tuvx_no_photolysis test_tuvx_no_photolysis.F90) target_sources(test_tuvx_no_photolysis PUBLIC ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_no_photolysis_rate.F90 + ${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 ) diff --git a/test/musica/tuvx/test_tuvx_gas_species.F90 b/test/musica/tuvx/test_tuvx_gas_species.F90 index c226f6c8..9b4e0be3 100644 --- a/test/musica/tuvx/test_tuvx_gas_species.F90 +++ b/test/musica/tuvx/test_tuvx_gas_species.F90 @@ -255,10 +255,6 @@ subroutine test_create_gas_species_profile() 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 ) diff --git a/test/musica/tuvx/test_tuvx_no_photolysis.F90 b/test/musica/tuvx/test_tuvx_no_photolysis.F90 index 05729a23..d15b68ea 100644 --- a/test/musica/tuvx/test_tuvx_no_photolysis.F90 +++ b/test/musica/tuvx/test_tuvx_no_photolysis.F90 @@ -1,63 +1,71 @@ -program test_tuvx_surface_albedo +program test_tuvx_no_photolysis - use musica_ccpp_tuvx_no_photolysis_rate + use musica_ccpp_tuvx_no_photolysis_rate - implicit none + 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_calculate_NO_photolysis_rate() + call test_calculate_NO_photolysis_rate() contains - subroutine test_calculate_NO_photolysis_rate() - use ccpp_kinds, only: kind_phys - implicit none - - integer, parameter :: number_of_vertical_layers = 2 - - real(kind_phys) :: solar_zenith_angle - real(kind_phys), dimension(4) :: extraterrestrial_flux - real(kind_phys), dimension(1,number_of_vertical_layers,4) :: constituents ! (column, layers, constituents) kg kg-1 - real(kind_phys), dimension(3) :: height_at_interfaces ! (layers + 1) km - real(kind_phys), dimension(1,number_of_vertical_layers) :: dry_air_density ! (column, layers) kg m-3 - integer :: N2_index, O2_index, O3_index, NO_index - real(kind_phys) :: molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO ! kg mol-1 - real(kind_phys), dimension(number_of_vertical_layers) :: jNO - - ! Some of this initialization data corresponds to values listed in the Minschwaner and Siskind (1993) paper - ! see the musica_ccpp_tuvx_no_photolysis_rate.F90 file for the full citation - - ! Initialize test data - solar_zenith_angle = 60.0_kind_phys - extraterrestrial_flux = (/ 1.5e13_kind_phys, 1.4e13_kind_phys, 1.3e13_kind_phys, 1.2e13_kind_phys /) - height_at_interfaces = reshape([40.0_kind_phys, 30.0_kind_phys, 20.0_kind_phys], shape(height_at_interfaces)) - dry_air_density = reshape([1.2_kind_phys, 1.1_kind_phys], shape(dry_air_density)) - N2_index = 1 - O2_index = 2 - O3_index = 3 - NO_index = 4 - molar_mass_N2 = 0.0280134_kind_phys - molar_mass_O2 = 0.0319988_kind_phys - molar_mass_O3 = 0.0479982_kind_phys - molar_mass_NO = 0.3000610_kind_phys - - constituents(1,:,N2_index) = 0.78_kind_phys - constituents(1,:,O2_index) = 0.21_kind_phys - constituents(1,:,O3_index) = 20.0e-9_kind_phys - constituents(1,:,NO_index) = 0.3e-9_kind_phys - - ! Call the function to test - jNO = calculate_NO_photolysis_rate(size(constituents, dim=2), solar_zenith_angle, extraterrestrial_flux, constituents(1,:,:), height_at_interfaces, & - dry_air_density(1,:), N2_index, O2_index, O3_index, NO_index, molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO) - - ! Validate the results - print *, jNO - ASSERT(jNO(1) .ne. 0.0_kind_phys) - ASSERT(jNO(2) .ne. 0.0_kind_phys) - ASSERT(jNO(1) .lt. 1.0_kind_phys) - ASSERT(jNO(2) .lt. 1.0_kind_phys) - end subroutine test_calculate_NO_photolysis_rate - -end program test_tuvx_surface_albedo + subroutine test_calculate_NO_photolysis_rate() + use ccpp_kinds, only: kind_phys + implicit none + + integer, parameter :: number_of_vertical_layers = 2 + + real(kind_phys) :: solar_zenith_angle + real(kind_phys), dimension(4) :: extraterrestrial_flux + real(kind_phys), dimension(1,number_of_vertical_layers,4) :: constituents ! (column, layers, constituents) kg kg-1 + real(kind_phys), dimension(1,number_of_vertical_layers,2) :: constituents_NO_photolysis + real(kind_phys), dimension(3) :: height_at_interfaces ! (layers + 1) km + real(kind_phys), dimension(1,number_of_vertical_layers) :: dry_air_density ! (column, layers) kg m-3 + integer :: N2_index, O2_index, O3_index, NO_index + real(kind_phys) :: molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO ! kg mol-1 + real(kind_phys), dimension(number_of_vertical_layers) :: jNO + ! Some of this initialization data corresponds to values listed in the Minschwaner and Siskind (1993) paper + ! see the musica_ccpp_tuvx_no_photolysis_rate.F90 file for the full citation + + ! Initialize test data + solar_zenith_angle = 60.0_kind_phys + extraterrestrial_flux = (/ 1.5e13_kind_phys, 1.4e13_kind_phys, 1.3e13_kind_phys, 1.2e13_kind_phys /) + height_at_interfaces = reshape([40.0_kind_phys, 30.0_kind_phys, 20.0_kind_phys], shape(height_at_interfaces)) + dry_air_density = reshape([1.2_kind_phys, 1.1_kind_phys], shape(dry_air_density)) + N2_index = 1 + O2_index = 2 + O3_index = 3 + NO_index = 4 + molar_mass_N2 = 0.0280134_kind_phys + molar_mass_O2 = 0.0319988_kind_phys + molar_mass_O3 = 0.0479982_kind_phys + molar_mass_NO = 0.3000610_kind_phys + + constituents(1,:,N2_index) = 0.78_kind_phys + constituents(1,:,O2_index) = 0.21_kind_phys + constituents(1,:,O3_index) = 20.0e-9_kind_phys + constituents(1,:,NO_index) = 0.3e-9_kind_phys + constituents_NO_photolysis(1,:,1) = constituents(1,:,N2_index) + constituents_NO_photolysis(1,:,2) = constituents(1,:,NO_index) + + ! Call the function to test + ! jNO = calculate_NO_photolysis_rate(size(constituents, dim=2), solar_zenith_angle, extraterrestrial_flux, constituents(1,:,:), height_at_interfaces, & + ! dry_air_density(1,:), N2_index, O2_index, O3_index, NO_index, molar_mass_N2, molar_mass_O2, molar_mass_O3, molar_mass_NO) + + jNO = calculate_NO_photolysis_rate(size(constituents, dim=2), & + solar_zenith_angle, extraterrestrial_flux, height_at_interfaces, & + dry_air_density(1,:), constituents(1,:,O2_index), & + constituents(1,:,O3_index), constituents_NO_photolysis(1,:,:)) + + ! Validate the results + print *, jNO + ASSERT(jNO(1) .ne. 0.0_kind_phys) + ASSERT(jNO(2) .ne. 0.0_kind_phys) + ASSERT(jNO(1) .lt. 1.0_kind_phys) + ASSERT(jNO(2) .lt. 1.0_kind_phys) + + end subroutine test_calculate_NO_photolysis_rate + +end program test_tuvx_no_photolysis \ No newline at end of file