diff --git a/examples/ts1_tsmlt.json b/examples/ts1_tsmlt.json index 5cd0e33e..7ea1806e 100644 --- a/examples/ts1_tsmlt.json +++ b/examples/ts1_tsmlt.json @@ -199,6 +199,9 @@ "value": 1.0 } ] + }, + "heating" : { + "energy term": 175.05 } }, { @@ -228,6 +231,9 @@ "value": 0.0 } ] + }, + "heating" : { + "energy term": 242.37 } }, { @@ -244,6 +250,9 @@ }, "quantum yield": { "type": "O3+hv->O2+O(1D)" + }, + "heating" : { + "energy term": 310.32 } }, { @@ -260,6 +269,9 @@ }, "quantum yield": { "type": "O3+hv->O2+O(3P)" + }, + "heating" : { + "energy term": 1179.87 } }, { diff --git a/examples/ts1_tsmlt.yml b/examples/ts1_tsmlt.yml index 12de351d..c6d1e4e8 100644 --- a/examples/ts1_tsmlt.yml +++ b/examples/ts1_tsmlt.yml @@ -210,6 +210,8 @@ photolysis: lower extrapolation: type: boundary type: base + heating: + energy term: 175.05 name: jo2_a quantum yield: constant value: 0 @@ -229,6 +231,8 @@ photolysis: lower extrapolation: type: boundary type: base + heating: + energy term: 242.37 name: jo2_b quantum yield: constant value: 1.0 @@ -246,6 +250,8 @@ photolysis: - file path: data/cross_sections/O3_3.nc - file path: data/cross_sections/O3_4.nc type: O3 + heating: + energy term: 310.32 name: jo3_a quantum yield: type: O3+hv->O2+O(1D) @@ -257,6 +263,8 @@ photolysis: - file path: data/cross_sections/O3_3.nc - file path: data/cross_sections/O3_4.nc type: O3 + heating: + energy term: 1179.87 name: jo3_b quantum yield: type: O3+hv->O2+O(3P) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 37b64f99..bb185d68 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -46,6 +46,7 @@ target_sources(tuvx_object grid.F90 grid_factory.F90 grid_warehouse.F90 + heating_rates.F90 interpolate.F90 la_sr_bands.F90 linear_algebra.F90 diff --git a/src/core.F90 b/src/core.F90 index 43c7e482..72828dda 100644 --- a/src/core.F90 +++ b/src/core.F90 @@ -9,6 +9,7 @@ module tuvx_core use musica_constants, only : dk => musica_dk use tuvx_dose_rates, only : dose_rates_t use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_heating_rates, only : heating_rates_t use tuvx_la_sr_bands, only : la_sr_bands_t use tuvx_photolysis_rates, only : photolysis_rates_t use tuvx_profile_warehouse, only : profile_warehouse_t @@ -32,10 +33,11 @@ module tuvx_core type(radiative_transfer_t), pointer :: radiative_transfer_ => null() type(photolysis_rates_t), pointer :: photolysis_rates_ => null() type(dose_rates_t), pointer :: dose_rates_ => null() + type(heating_rates_t), pointer :: heating_rates_ => null() type(radiation_field_t), pointer :: radiation_field_ => null() logical :: enable_diagnostics_ ! determines if diagnostic output is written or not contains - ! Calculate photolysis rate constants and dose rates + ! Calculate photolysis rate constants, dose rates, and heating rates procedure :: run ! Returns a grid from the warehouse procedure :: get_grid @@ -50,10 +52,14 @@ module tuvx_core procedure :: number_of_photolysis_reactions ! Returns the number of dose rates procedure :: number_of_dose_rates + ! Returns the number of heating rates + procedure :: number_of_heating_rates ! Returns the set of photolysis reaction labels procedure :: photolysis_reaction_labels ! Returns the set of dose rate labels procedure :: dose_rate_labels + ! Returns the set of heating rate labels + procedure :: heating_rate_labels ! Returns the photolysis reaction cross section for the current conditions procedure :: get_photolysis_cross_section ! Returns the photolysis reaction quantum yield for the current conditions @@ -165,6 +171,9 @@ function constructor( config, grids, profiles, radiators ) result( new_core ) photolysis_rates_t( child_config, & new_core%grid_warehouse_, & new_core%profile_warehouse_ ) + new_core%heating_rates_ => heating_rates_t( child_config, & + new_core%grid_warehouse_, & + new_core%profile_warehouse_ ) end if ! dose rates @@ -191,7 +200,7 @@ end function constructor !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine run( this, solar_zenith_angle, earth_sun_distance, & - photolysis_rate_constants, dose_rates, diagnostic_label ) + photolysis_rate_constants, dose_rates, heating_rates, diagnostic_label ) ! Performs calculations for specified photolysis and dose rates for a ! given set of conditions @@ -205,6 +214,7 @@ subroutine run( this, solar_zenith_angle, earth_sun_distance, & real(dk), intent(in) :: earth_sun_distance ! [AU] real(dk), optional, intent(out) :: photolysis_rate_constants(:,:) ! (vertical level, reaction) [s-1] real(dk), optional, intent(out) :: dose_rates(:,:) ! (vertical level, reaction) [s-1] + real(dk), optional, intent(out) :: heating_rates(:,:) ! (vertical level, reaction) [J s-1] character(len=*), optional, intent(in) :: diagnostic_label ! label used in diagnostic file names ! Local variables @@ -247,6 +257,14 @@ subroutine run( this, solar_zenith_angle, earth_sun_distance, & photolysis_rate_constants, & diag_label ) end if + if( associated( this%heating_rates_ ) .and. present( heating_rates ) ) then + call this%heating_rates_%get( this%la_sr_bands_, & + this%spherical_geometry_, & + this%grid_warehouse_, & + this%profile_warehouse_, & + this%radiation_field_, & + heating_rates ) + end if if( associated( this%dose_rates_ ) .and. present( dose_rates ) ) then call this%dose_rates_%get( this%grid_warehouse_, & this%profile_warehouse_, & @@ -410,6 +428,20 @@ integer function number_of_dose_rates( this ) end function number_of_dose_rates +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + integer function number_of_heating_rates( this ) + ! Returns the number of heating rates + + class(core_t), intent(in) :: this + + number_of_heating_rates = 0 + if( associated( this%heating_rates_ ) ) then + number_of_heating_rates = this%heating_rates_%size( ) + end if + + end function number_of_heating_rates + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function photolysis_reaction_labels( this ) result( labels ) @@ -442,6 +474,22 @@ function dose_rate_labels( this ) result( labels ) end function dose_rate_labels +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function heating_rate_labels( this ) result( labels ) + ! Returns the set of heating rate labels + + class(core_t), intent(in) :: this + type(string_t), allocatable :: labels(:) + + if( associated( this%heating_rates_ ) ) then + labels = this%heating_rates_%labels( ) + else + allocate( labels( 0 ) ) + end if + + end function heating_rate_labels + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function get_photolysis_cross_section( this, reaction_label, found ) & @@ -558,6 +606,11 @@ integer function pack_size( this, comm ) if( associated( this%dose_rates_ ) ) then pack_size = pack_size + this%dose_rates_%pack_size( comm ) end if + pack_size = pack_size + & + musica_mpi_pack_size( associated( this%heating_rates_ ), comm ) + if( associated( this%heating_rates_ ) ) then + pack_size = pack_size + this%heating_rates_%pack_size( comm ) + end if #else pack_size = 0 #endif @@ -617,6 +670,11 @@ subroutine mpi_pack( this, buffer, position, comm ) if( associated( this%dose_rates_ ) ) then call this%dose_rates_%mpi_pack( buffer, position, comm ) end if + call musica_mpi_pack( buffer, position, & + associated( this%heating_rates_ ), comm ) + if( associated( this%heating_rates_ ) ) then + call this%heating_rates_%mpi_pack( buffer, position, comm ) + end if call assert( 332208077, position - prev_pos <= this%pack_size( comm ) ) #endif @@ -676,6 +734,11 @@ subroutine mpi_unpack( this, buffer, position, comm ) allocate( this%dose_rates_ ) call this%dose_rates_%mpi_unpack( buffer, position, comm ) end if + call musica_mpi_unpack( buffer, position, alloced, comm ) + if( alloced ) then + allocate( this%heating_rates_ ) + call this%heating_rates_%mpi_unpack( buffer, position, comm ) + end if call assert( 332208077, position - prev_pos <= this%pack_size( comm ) ) #endif @@ -713,6 +776,9 @@ subroutine finalize( this ) if( associated( this%radiation_field_ ) ) then deallocate( this%radiation_field_ ) end if + if( associated( this%heating_rates_ ) ) then + deallocate( this%heating_rates_ ) + end if end subroutine finalize diff --git a/src/grid_warehouse.F90 b/src/grid_warehouse.F90 index 04f68d8f..fa14251b 100644 --- a/src/grid_warehouse.F90 +++ b/src/grid_warehouse.F90 @@ -152,10 +152,10 @@ function get_grid_string( this, name, units ) result( a_grid_ptr ) use musica_string, only : string_t use tuvx_grid, only : grid_t - class(grid_warehouse_t), intent(inout) :: this ! This :f:type:`~tuvx_grid_warehouse/grid_warehouse_t` - type(string_t), intent(in) :: name ! The name of a grid, see :ref:`configuration-grids` for grid names - type(string_t), intent(in) :: units ! The units of the grid - class(grid_t), pointer :: a_grid_ptr ! The :f:type:`~tuvx_grid/grid_t` which matches the name passed in + class(grid_warehouse_t), intent(in) :: this ! This :f:type:`~tuvx_grid_warehouse/grid_warehouse_t` + type(string_t), intent(in) :: name ! The name of a grid, see :ref:`configuration-grids` for grid names + type(string_t), intent(in) :: units ! The units of the grid + class(grid_t), pointer :: a_grid_ptr ! The :f:type:`~tuvx_grid/grid_t` which matches the name passed in a_grid_ptr => this%get_grid_char( name%to_char( ), units%to_char( ) ) @@ -169,9 +169,9 @@ function get_grid_ptr( this, ptr ) result( grid ) use musica_assert, only : assert_msg use tuvx_grid, only : grid_t - class(grid_warehouse_t), intent(inout) :: this ! This grid warehouse - type(grid_warehouse_ptr), intent(in) :: ptr ! Pointer to a grid in the warehouse - class(grid_t), pointer :: grid + class(grid_warehouse_t), intent(in) :: this ! This grid warehouse + type(grid_warehouse_ptr), intent(in) :: ptr ! Pointer to a grid in the warehouse + class(grid_t), pointer :: grid call assert_msg( 870082797, ptr%index_ > 0, "Invalid grid pointer" ) allocate( grid, source = this%grids_( ptr%index_ )%val_ ) diff --git a/src/heating_rates.F90 b/src/heating_rates.F90 new file mode 100644 index 00000000..c3d4ab62 --- /dev/null +++ b/src/heating_rates.F90 @@ -0,0 +1,614 @@ +! Copyright (C) 2024 National Center for Atmospheric Research +! SPDX-License-Identifier: Apache-2.0 +! +module tuvx_heating_rates + ! The chemical potential heating rates type heating_rates_t and related functions + + use musica_constants, only : dk => musica_dk + use musica_string, only : string_t + use tuvx_cross_section, only : cross_section_ptr + use tuvx_grid_warehouse, only : grid_warehouse_ptr + use tuvx_profile_warehouse, only : profile_warehouse_ptr + use tuvx_quantum_yield, only : quantum_yield_ptr + + implicit none + + private + public :: heating_rates_t + + type :: heating_parameters_t + ! Heating parameters for a single photolyzing species + type(string_t) :: label_ ! label for the heating rate + type(cross_section_ptr) :: cross_section_ ! cross section + type(quantum_yield_ptr) :: quantum_yield_ ! quantum yield + real(kind=dk) :: scaling_factor_ ! scaling factor for the heating rate + real(kind=dk), allocatable :: energy_(:) ! wavelength resolved bond-dissociation energy [J] + contains + !> Returns the size of a character buffer needed to pack the heating parameters + procedure :: pack_size => heating_parameters_pack_size + !> Packs the heating parameters into a character buffer + procedure :: mpi_pack => heating_parameters_mpi_pack + !> Unpacks the heating parameters from a character buffer + procedure :: mpi_unpack => heating_parameters_mpi_unpack + end type heating_parameters_t + + !> heating_parameters_t constructor + interface heating_parameters_t + module procedure :: heating_parameters_constructor + end interface heating_parameters_t + + type, public :: heating_rates_t + type(heating_parameters_t), allocatable :: heating_parameters_(:) ! heating parameters for each photolyzing species + type(grid_warehouse_ptr) :: height_grid_ ! height grid + type(grid_warehouse_ptr) :: wavelength_grid_ ! wavelength grid + type(profile_warehouse_ptr) :: etfl_profile_ ! Extraterrestrial flux profile + type(profile_warehouse_ptr) :: air_profile_ ! Air profile + integer, allocatable :: o2_rate_indices_(:) ! indices in the heating rates array where O2 + ! corrections to the cross-section in the + ! Lyman-Alpha and Schumann-Runge bands should + ! be applied + contains + !> Calulates the heating rates + procedure :: get + !> Returns the names of each photolysis reaction with a heating rate + procedure :: labels + !> Returns the number of heating rates + procedure :: size => get_number + !> Returns the size of a character buffer needed to pack the heating rates + procedure :: pack_size + !> Packs the heating rates into a character buffer + procedure :: mpi_pack + !> Unpacks the heating rates from a character buffer + procedure :: mpi_unpack + !> Cleans up memory + final :: destructor + end type heating_rates_t + + !> heating_rates_t constructor + interface heating_rates_t + module procedure :: constructor + end interface heating_rates_t + +contains +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> heating_rates_t constructor + function constructor( config, grids, profiles ) result( this ) + + use musica_assert, only : assert, assert_msg + use musica_config, only : config_t + use musica_iterator, only : iterator_t + use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_profile_warehouse, only : profile_warehouse_t + + !> Heating rate collection + type(heating_rates_t), pointer :: this + !> Configuration + type(config_t), intent(inout) :: config + !> Grids + type(grid_warehouse_t), intent(inout) :: grids + !> Profiles + type(profile_warehouse_t), intent(inout) :: profiles + + character(len=*), parameter :: Iam = 'heating rates constructor' + type(config_t) :: reaction_set, reaction_config, heating_config + class(iterator_t), pointer :: iter + type(string_t) :: label + type(string_t) :: required_keys(1), optional_keys(1) + logical :: found, do_apply_bands + integer :: n_hr, i_hr, n_O2, i_O2 + + required_keys(1) = "reactions" + optional_keys(1) = "enable diagnostics" + + call assert_msg( 310567326, & + config%validate( required_keys, optional_keys ), & + "Invalid configuration for heating rates" ) + + allocate( this ) + this%height_grid_ = grids%get_ptr( "height", "km" ) + this%wavelength_grid_ = grids%get_ptr( "wavelength", "nm" ) + this%etfl_profile_ = profiles%get_ptr( "extraterrestrial flux", & + "photon cm-2 s-1" ) + this%air_profile_ = profiles%get_ptr( "air", "molecule cm-3" ) + + ! iterate over photolysis reactions looking for those with + ! heating rate parameters + allocate( this%o2_rate_indices_( 0 ) ) + call config%get( "reactions", reaction_set, Iam ) + iter => reaction_set%get_iterator( ) + n_hr = 0 + n_O2 = 0 + do while( iter%next( ) ) + call reaction_set%get( iter, reaction_config, Iam ) + call reaction_config%get( "heating", heating_config, Iam, found = found ) + if( found ) then + n_hr = n_hr + 1 + call reaction_config%get( "apply O2 bands", do_apply_bands, Iam, & + default = .false. ) + if( do_apply_bands ) n_O2 = n_O2 + 1 + end if + end do + allocate( this%heating_parameters_( n_hr ) ) + call iter%reset( ) + i_hr = 0 + i_O2 = 0 + do while( iter%next( ) ) + call reaction_set%get( iter, reaction_config, Iam ) + call reaction_config%get( "heating", heating_config, Iam, found = found ) + if( found ) then + i_hr = i_hr + 1 + call reaction_config%get( "name", label, Iam ) + this%heating_parameters_( i_hr ) = & + heating_parameters_constructor( reaction_config, grids, profiles ) + call reaction_config%get( "apply O2 bands", do_apply_bands, Iam, & + default = .false. ) + if( do_apply_bands ) then + i_O2 = i_O2 + 1 + this%o2_rate_indices_( i_O2 ) = i_hr + end if + end if + end do + call assert( 357615745, i_hr .eq. n_hr ) + call assert( 336635308, i_O2 .eq. n_O2 ) + deallocate( iter ) + + end function constructor + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> heating_parameters_t constructor + function heating_parameters_constructor( config, grids, profiles ) & + result( this ) + + use musica_assert, only : assert_msg + use musica_config, only : config_t + use tuvx_constants, only : hc + use tuvx_cross_section_factory, only : cross_section_builder + use tuvx_grid, only : grid_t + use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_profile_warehouse, only : profile_warehouse_t + use tuvx_quantum_yield_factory, only : quantum_yield_builder + + !> Heating parameters for a single photolyzing species + type(heating_parameters_t) :: this + !> Configuration for the photolysis reaction + type(config_t), intent(inout) :: config + !> Grids + type(grid_warehouse_t), intent(inout) :: grids + !> Profiles + type(profile_warehouse_t), intent(inout) :: profiles + + character(len=*), parameter :: Iam = 'heating parameters constructor' + type(config_t) :: heating_config, cs_config, qy_config + class(grid_t), pointer :: wavelengths + real(kind=dk) :: energy_term + type(string_t) :: required_keys(4), optional_keys(1) + type(string_t) :: heating_required_keys(1), heating_optional_keys(0) + + required_keys(1) = "name" + required_keys(2) = "cross section" + required_keys(3) = "quantum yield" + required_keys(4) = "heating" + optional_keys(1) = "scaling factor" + + call assert_msg( 316144353, & + config%validate( required_keys, optional_keys ), & + "Invalid configuration for photolysis reactions with "// & + "heating parameters" ) + + call config%get( "heating", heating_config, Iam ) + + heating_required_keys(1) = "energy term" + + call assert_msg( 316144354, & + heating_config%validate( heating_required_keys, & + heating_optional_keys ), & + "Invalid configuration for heating parameters" ) + + call config%get( "name", this%label_, Iam ) + call config%get( "cross section", cs_config, Iam ) + this%cross_section_%val_ => cross_section_builder( cs_config, grids, & + profiles ) + call config%get( "quantum yield", qy_config, Iam ) + this%quantum_yield_%val_ => quantum_yield_builder( qy_config, grids, & + profiles ) + call config%get( "scaling factor", this%scaling_factor_, Iam, & + default = 1.0_dk ) + call heating_config%get( "energy term", energy_term, Iam ) + wavelengths => grids%get_grid( "wavelength", "nm" ) + allocate( this%energy_( wavelengths%ncells_ ) ) + this%energy_(:) = & + max( 0.0_dk, hc * 1.0e9_dk * ( energy_term - wavelengths%mid_(:) ) / & + ( energy_term * wavelengths%mid_(:) ) ) + deallocate( wavelengths ) + + end function heating_parameters_constructor + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> calculate heating rates + subroutine get( this, la_srb, spherical_geometry, grids, profiles, & + radiation_field, heating_rates ) + + use musica_assert, only : assert + use tuvx_grid, only : grid_t + use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_la_sr_bands, only : la_sr_bands_t + use tuvx_profile, only : profile_t + use tuvx_profile_warehouse, only : profile_warehouse_t + use tuvx_solver, only : radiation_field_t + use tuvx_spherical_geometry, only : spherical_geometry_t + + !> Heating rate collection + class(heating_rates_t), intent(in) :: this + !> Lyman Alpha and Schumann-Runge bands + class(la_sr_bands_t), intent(inout) :: la_srb + !> Spherical geometry + class(spherical_geometry_t), intent(inout) :: spherical_geometry + !> Grids + class(grid_warehouse_t), intent(inout) :: grids + !> Profiles + class(profile_warehouse_t), intent(inout) :: profiles + !> Radiation field + class(radiation_field_t), intent(in) :: radiation_field + !> Heating rates (vertical interface, reaction) [J s-1] + real(kind=dk), intent(inout) :: heating_rates(:,:) + + character(len=*), parameter :: Iam = 'heating rates get' + class(grid_t), pointer :: heights, wavelengths + class(profile_t), pointer :: etfl, air + real(kind=dk), allocatable :: actinic_flux(:,:), xsqy(:,:) + real(kind=dk), allocatable :: cross_section(:,:), quantum_yield(:,:) + real(kind=dk), allocatable :: air_vertical_column(:), air_slant_column(:) + integer :: i_rate, n_rates, i_height + + heights => grids%get_grid( this%height_grid_ ) + wavelengths => grids%get_grid( this%wavelength_grid_ ) + etfl => profiles%get_profile( this%etfl_profile_ ) + air => profiles%get_profile( this%air_profile_ ) + + n_rates = size( this%heating_parameters_ ) + call assert( 966855732, & + size( heating_rates, 1 ) .eq. heights%ncells_ + 1 .and. & + size( heating_rates, 2 ) .eq. n_rates ) + + actinic_flux = transpose( radiation_field%fdr_ + radiation_field%fup_ + & + radiation_field%fdn_ ) + do i_height = 1, heights%ncells_ + 1 + actinic_flux( :, i_height ) = actinic_flux( :, i_height ) * etfl%mid_val_ + end do + where( actinic_flux < 0.0_dk ) + actinic_flux = 0.0_dk + end where + + do i_rate = 1, n_rates + associate( params => this%heating_parameters_( i_rate ) ) + cross_section = params%cross_section_%val_%calculate( grids, profiles ) + quantum_yield = params%quantum_yield_%val_%calculate( grids, profiles ) + + ! O2 photolysis can have special la & srb band handling + if( any( this%o2_rate_indices_ == i_rate ) ) then + allocate( air_vertical_column( air%ncells_ ), & + air_slant_column( air%ncells_ + 1 ) ) + call spherical_geometry%air_mass( air%exo_layer_dens_, & + air_vertical_column, & + air_slant_column ) + call la_srb%cross_section( grids, profiles, air_vertical_column, & + air_slant_column, cross_section, & + spherical_geometry ) + deallocate( air_vertical_column, air_slant_column ) + end if + + xsqy = transpose( cross_section * quantum_yield ) + do i_height = 1, heights%ncells_ + 1 + heating_rates( i_height, i_rate ) = & + dot_product( actinic_flux( :, i_height ), & + params%energy_(:) * xsqy( :, i_height ) ) * & + params%scaling_factor_ + end do + end associate + end do + + deallocate( heights ) + deallocate( wavelengths ) + deallocate( etfl ) + deallocate( air ) + + end subroutine get + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the names of each photolysis reaction with a heating rate + function labels( this ) + + !> Photolysis reaction labels + type(string_t), allocatable :: labels(:) + !> Heating rate collection + class(heating_rates_t), intent(in) :: this + + allocate( labels( size( this%heating_parameters_ ) ) ) + labels(:) = this%heating_parameters_(:)%label_ + + end function labels + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the number of heating rates + function get_number( this ) result( n_rates ) + + !> Number of heating rates + integer :: n_rates + !> Heating rate collection + class(heating_rates_t), intent(in) :: this + + n_rates = 0 + if( allocated( this%heating_parameters_ ) ) then + n_rates = size( this%heating_parameters_ ) + end if + + end function get_number + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the size of a character buffer needed to pack the heating rates + function pack_size( this, comm ) + + use musica_mpi, only : musica_mpi_pack_size + + !> Heating rate collection + class(heating_rates_t), intent(in) :: this + !> MPI communicator + integer, intent(in) :: comm + !> Size of the character buffer + integer :: pack_size + +#ifdef MUSICA_USE_MPI + integer :: i_elem + + pack_size = musica_mpi_pack_size( allocated( this%heating_parameters_ ), & + comm ) + if( allocated( this%heating_parameters_ ) ) then + pack_size = pack_size + & + musica_mpi_pack_size( size( this%heating_parameters_ ), comm ) + do i_elem = 1, size( this%heating_parameters_ ) + pack_size = pack_size + & + this%heating_parameters_( i_elem )%pack_size( comm ) + end do + end if + pack_size = pack_size + & + this%height_grid_%pack_size( comm ) + & + this%wavelength_grid_%pack_size( comm ) + & + this%etfl_profile_%pack_size( comm ) + & + this%air_profile_%pack_size( comm ) + & + musica_mpi_pack_size( this%o2_rate_indices_, comm ) +#else + pack_size = 0 +#endif + + end function pack_size + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Packs the heating rates into a character buffer + subroutine mpi_pack( this, buffer, position, comm ) + + use musica_assert, only : assert + use musica_mpi, only : musica_mpi_pack + + !> Heating rate collection + class(heating_rates_t), intent(in) :: this + !> Character buffer + character, intent(inout) :: buffer(:) + !> Position in the buffer + integer, intent(inout) :: position + !> MPI communicator + integer, intent(in) :: comm + +#ifdef MUSICA_USE_MPI + integer :: prev_pos, i_elem + + prev_pos = position + call musica_mpi_pack( buffer, position, & + allocated( this%heating_parameters_ ), comm ) + if( allocated( this%heating_parameters_ ) ) then + call musica_mpi_pack( buffer, position, & + size( this%heating_parameters_ ), comm ) + do i_elem = 1, size( this%heating_parameters_ ) + call this%heating_parameters_( i_elem )%mpi_pack( buffer, position, & + comm ) + end do + end if + call this%height_grid_%mpi_pack( buffer, position, comm ) + call this%wavelength_grid_%mpi_pack( buffer, position, comm ) + call this%etfl_profile_%mpi_pack( buffer, position, comm ) + call this%air_profile_%mpi_pack( buffer, position, comm ) + call musica_mpi_pack( buffer, position, this%o2_rate_indices_, comm ) + call assert( 247051769, position - prev_pos <= this%pack_size( comm ) ) +#endif + + end subroutine mpi_pack + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Unpacks the heating rates from a character buffer + subroutine mpi_unpack( this, buffer, position, comm ) + + use musica_assert, only : assert + use musica_mpi, only : musica_mpi_unpack + + !> Heating rate collection + class(heating_rates_t), intent(out) :: this + !> Character buffer + character, intent(inout) :: buffer(:) + !> Position in the buffer + integer, intent(inout) :: position + !> MPI communicator + integer, intent(in) :: comm + +#ifdef MUSICA_USE_MPI + integer :: prev_pos, i_elem, n_elems + logical :: is_allocated + + prev_pos = position + call musica_mpi_unpack( buffer, position, is_allocated, comm ) + if( is_allocated ) then + call musica_mpi_unpack( buffer, position, n_elems, comm ) + allocate( this%heating_parameters_( n_elems ) ) + do i_elem = 1, n_elems + call this%heating_parameters_( i_elem )%mpi_unpack( buffer, position, & + comm ) + end do + end if + call this%height_grid_%mpi_unpack( buffer, position, comm ) + call this%wavelength_grid_%mpi_unpack( buffer, position, comm ) + call this%etfl_profile_%mpi_unpack( buffer, position, comm ) + call this%air_profile_%mpi_unpack( buffer, position, comm ) + call musica_mpi_unpack( buffer, position, this%o2_rate_indices_, comm ) + call assert( 631316749, position - prev_pos <= this%pack_size( comm ) ) +#endif + + end subroutine mpi_unpack + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Returns the size of a character buffer needed to pack the heating + !! parameters + function heating_parameters_pack_size( this, comm ) result( pack_size ) + + use musica_mpi, only : musica_mpi_pack_size + use tuvx_cross_section_factory, only : cross_section_type_name + use tuvx_quantum_yield_factory, only : quantum_yield_type_name + + !> Heating parameters for a single photolyzing species + class(heating_parameters_t), intent(in) :: this + !> MPI communicator + integer, intent(in) :: comm + !> Size of the character buffer + integer :: pack_size + +#ifdef MUSICA_USE_MPI + type(string_t) :: cs_type_name, qy_type_name + + cs_type_name = cross_section_type_name( this%cross_section_%val_ ) + qy_type_name = quantum_yield_type_name( this%quantum_yield_%val_ ) + pack_size = this%label_%pack_size( comm ) + & + cs_type_name%pack_size( comm ) + & + this%cross_section_%val_%pack_size( comm ) + & + qy_type_name%pack_size( comm ) + & + this%quantum_yield_%val_%pack_size( comm ) + & + musica_mpi_pack_size( this%scaling_factor_, comm ) + & + musica_mpi_pack_size( this%energy_, comm ) +#else + pack_size = 0 +#endif + + end function heating_parameters_pack_size + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Packs the heating parameters into a character buffer + subroutine heating_parameters_mpi_pack( this, buffer, position, comm ) + + use musica_assert, only : assert + use musica_mpi, only : musica_mpi_pack + use tuvx_cross_section_factory, only : cross_section_type_name + use tuvx_quantum_yield_factory, only : quantum_yield_type_name + + !> Heating parameters for a single photolyzing species + class(heating_parameters_t), intent(in) :: this + !> Character buffer + character, intent(inout) :: buffer(:) + !> Position in the buffer + integer, intent(inout) :: position + !> MPI communicator + integer, intent(in) :: comm + +#ifdef MUSICA_USE_MPI + integer :: prev_pos + type(string_t) :: cs_type_name, qy_type_name + + prev_pos = position + cs_type_name = cross_section_type_name( this%cross_section_%val_ ) + qy_type_name = quantum_yield_type_name( this%quantum_yield_%val_ ) + call this%label_%mpi_pack( buffer, position, comm ) + call cs_type_name%mpi_pack( buffer, position, comm ) + call this%cross_section_%val_%mpi_pack( buffer, position, comm ) + call qy_type_name%mpi_pack( buffer, position, comm ) + call this%quantum_yield_%val_%mpi_pack( buffer, position, comm ) + call musica_mpi_pack( buffer, position, this%scaling_factor_, comm ) + call musica_mpi_pack( buffer, position, this%energy_, comm ) + call assert( 243240701, position - prev_pos <= this%pack_size( comm ) ) +#endif + + end subroutine heating_parameters_mpi_pack + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Unpacks the heating parameters from a character buffer + subroutine heating_parameters_mpi_unpack( this, buffer, position, comm ) + + use musica_assert, only : assert + use musica_mpi, only : musica_mpi_unpack + use tuvx_cross_section_factory, only : cross_section_allocate + use tuvx_quantum_yield_factory, only : quantum_yield_allocate + + !> Heating parameters for a single photolyzing species + class(heating_parameters_t), intent(out) :: this + !> Character buffer + character, intent(inout) :: buffer(:) + !> Position in the buffer + integer, intent(inout) :: position + !> MPI communicator + integer, intent(in) :: comm + +#ifdef MUSICA_USE_MPI + integer :: prev_pos + type(string_t) :: cs_type_name, qy_type_name + + prev_pos = position + call this%label_%mpi_unpack( buffer, position, comm ) + call cs_type_name%mpi_unpack( buffer, position, comm ) + this%cross_section_%val_ => cross_section_allocate( cs_type_name ) + call this%cross_section_%val_%mpi_unpack( buffer, position, comm ) + call qy_type_name%mpi_unpack( buffer, position, comm ) + this%quantum_yield_%val_ => quantum_yield_allocate( qy_type_name ) + call this%quantum_yield_%val_%mpi_unpack( buffer, position, comm ) + call musica_mpi_unpack( buffer, position, this%scaling_factor_, comm ) + call musica_mpi_unpack( buffer, position, this%energy_, comm ) + call assert( 243240702, position - prev_pos <= this%pack_size( comm ) ) +#endif + + end subroutine heating_parameters_mpi_unpack + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> Cleans up memory + elemental subroutine destructor( this ) + + !> Heating rates + type(heating_rates_t), intent(inout) :: this + + integer :: i_rate + + if( allocated( this%heating_parameters_ ) ) then + do i_rate = 1, size( this%heating_parameters_ ) + associate( params => this%heating_parameters_( i_rate ) ) + if( associated( params%cross_section_%val_ ) ) then + deallocate( params%cross_section_%val_ ) + nullify( params%cross_section_%val_ ) + end if + if( associated( params%quantum_yield_%val_ ) ) then + deallocate( params%quantum_yield_%val_ ) + nullify( params%quantum_yield_%val_ ) + end if + end associate + end do + deallocate( this%heating_parameters_ ) + end if + + end subroutine destructor + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +end module tuvx_heating_rates \ No newline at end of file diff --git a/src/output.F90 b/src/output.F90 index faf041a7..e2561a93 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -23,9 +23,11 @@ module tuvx_output class(io_t), pointer :: file_ => null( ) logical :: do_photo_ = .false. logical :: do_dose_ = .false. + logical :: do_heating_ = .false. logical :: do_radiation_ = .false. type(string_t), allocatable :: photo_labels_(:) type(string_t), allocatable :: dose_labels_(:) + type(string_t), allocatable :: heating_labels_(:) type(string_t), allocatable :: photo_cross_sections_(:) type(string_t), allocatable :: photo_quantum_yields_(:) contains @@ -58,13 +60,14 @@ function constructor( config, core ) result( this ) character(len=*), parameter :: Iam = "output writer" integer :: stat type(string_t) :: file_path - type(string_t) :: required_keys(2), optional_keys(3) + type(string_t) :: required_keys(2), optional_keys(4) type(config_t) :: tuvx_config, rad_config required_keys(1) = "file path" required_keys(2) = "tuv-x configuration" optional_keys(1) = "include photolysis" optional_keys(2) = "include dose rates" + optional_keys(3) = "include heating rates" call assert_msg( 215370625, & config%validate( required_keys, optional_keys ), & @@ -93,6 +96,10 @@ function constructor( config, core ) result( this ) default = .false. ) if( this%do_dose_ ) this%dose_labels_ = core%dose_rate_labels( ) + call config%get( "include heating rates", this%do_heating_, Iam, & + default = .false. ) + if( this%do_heating_ ) this%heating_labels_ = core%heating_rate_labels( ) + ! Add custom diagnostics call config%get( "tuv-x configuration", tuvx_config, Iam ) call this%add_photolysis_diagnostics( tuvx_config ) @@ -106,7 +113,7 @@ end function constructor !> Outputs results subroutine output( this, step, core, photolysis_rate_constants, dose_rates, & - time, solar_zenith_angle, earth_sun_distance ) + heating_rates, time, solar_zenith_angle, earth_sun_distance ) use musica_assert, only : assert_msg use musica_constants, only : dk => musica_dk @@ -124,6 +131,8 @@ subroutine output( this, step, core, photolysis_rate_constants, dose_rates, & real(dk), optional, intent(in) :: photolysis_rate_constants(:,:) !> Dose rates (vertical level, dose rate type) real(dk), optional, intent(in) :: dose_rates(:,:) + !> Heating rates (vertical level, reaction) + real(dk), optional, intent(in) :: heating_rates(:,:) !> Time [hours] real(dk), optional, intent(in) :: time !> Solar zenith angle [degrees] @@ -209,6 +218,18 @@ subroutine output( this, step, core, photolysis_rate_constants, dose_rates, & end do end if + if( present( heating_rates ) ) then + call assert_msg( 935671025, this%do_heating_, "Heating rates are not " & + //"configured to be output" ) + dim_names(1) = "vertical_level" + units = "J s-1" + do i_rate = 1, size( this%heating_labels_ ) + var_name = clean_string( this%heating_labels_( i_rate ) ) + call this%file_%append( var_name, units, append_dim, step, & + dim_names(1), heating_rates( :, i_rate ), Iam ) + end do + end if + dim_names(1) = "vertical_level" dim_names(2) = "wavelength" units = "cm2 molecule-1" diff --git a/src/photolysis_rates.F90 b/src/photolysis_rates.F90 index c9d3493a..cdfbfa8f 100644 --- a/src/photolysis_rates.F90 +++ b/src/photolysis_rates.F90 @@ -178,7 +178,7 @@ subroutine add( this, config, grid_warehouse, profile_warehouse ) real(dk) :: scale_factor type(string_t) :: reaction_key logical :: do_apply_bands, found - type(string_t) :: required_keys(3), optional_keys(1) + type(string_t) :: required_keys(3), optional_keys(2) type(cross_section_ptr), allocatable :: temp_cs(:) type(quantum_yield_ptr), allocatable :: temp_qy(:) type(string_t), allocatable :: temp_handle(:) @@ -190,6 +190,7 @@ subroutine add( this, config, grid_warehouse, profile_warehouse ) required_keys(2) = "cross section" required_keys(3) = "quantum yield" optional_keys(1) = "scaling factor" + optional_keys(2) = "heating" call assert_msg( 780273355, & config%validate( required_keys, optional_keys ), & @@ -372,7 +373,8 @@ subroutine get( this, la_srb, spherical_geometry, grid_warehouse, & xsqy = transpose( cross_section * quantum_yield ) do vertNdx = 1, zGrid%ncells_ + 1 photolysis_rates( vertNdx, rateNdx ) = & - dot_product( actinicFlux( :, vertNdx ), xsqy( :, vertNdx ) ) + dot_product( actinicFlux( :, vertNdx ), xsqy( :, vertNdx ) ) * & + this%scaling_factors_( rateNdx ) enddo if( allocated( cross_section ) ) deallocate( cross_section ) if( allocated( quantum_yield ) ) deallocate( quantum_yield ) diff --git a/src/tuvx.F90 b/src/tuvx.F90 index e0e63126..b56a6602 100644 --- a/src/tuvx.F90 +++ b/src/tuvx.F90 @@ -139,11 +139,13 @@ subroutine run_tuvx( ) class(profile_t), pointer :: earth_sun_distance ! [AU] real(dk), allocatable :: photo_rates(:,:,:) ! (time, vertical level, reaction) [s-1] real(dk), allocatable :: dose_rates(:,:,:) ! (time, vertical level, dose rate) [?] + real(dk), allocatable :: heating_rates(:,:,:) ! (vertical level, reaction, thread) [K s-1] real(dk), allocatable :: thread_photo_rates(:,:,:) ! (vertical level, reaction, thread) [s-1] real(dk), allocatable :: thread_dose_rates(:,:,:) ! (vertical level, dose rate, thread) [?] + real(dk), allocatable :: thread_heating_rates(:,:,:)! (vertical level, reaction, thread) [K s-1] type(string_t) :: file_path character(len=2) :: diagnostic_label - class(output_t), pointer :: photo_output, dose_output + class(output_t), pointer :: photo_output, dose_output, heating_output type(config_t) :: config height => core%get_grid( "height", "km" ) @@ -161,10 +163,15 @@ subroutine run_tuvx( ) allocate( dose_rates( sza%ncells_ + 1, & height%ncells_ + 1, & core%number_of_dose_rates( ) ) ) + allocate( heating_rates( sza%ncells_ + 1, & + height%ncells_ + 1, & + core%number_of_heating_rates( ) ) ) + ! set up output files nullify( photo_output ) nullify( dose_output ) + nullify( heating_output ) if( core%number_of_photolysis_reactions( ) > 0 ) then call config%empty( ) call config%add( "file path", "photolysis_rate_constants.nc", Iam ) @@ -179,6 +186,13 @@ subroutine run_tuvx( ) call config%add( "tuv-x configuration", tuvx_config, Iam ) dose_output => output_t( config, core ) end if + if( core%number_of_heating_rates( ) > 0 ) then + call config%empty( ) + call config%add( "file path", "heating_rates.nc", Iam ) + call config%add( "include heating rates", .true., Iam ) + call config%add( "tuv-x configuration", tuvx_config, Iam ) + heating_output => output_t( config, core ) + end if ! calculate photolysis and dose rates do i_sza = 1, sza%ncells_ + 1 @@ -187,6 +201,7 @@ subroutine run_tuvx( ) earth_sun_distance%edge_val_( i_sza ), & photolysis_rate_constants = photo_rates( i_sza, :, : ), & dose_rates = dose_rates( i_sza, :, : ), & + heating_rates = heating_rates( i_sza, :, : ), & diagnostic_label = diagnostic_label ) ! output results @@ -204,6 +219,13 @@ subroutine run_tuvx( ) solar_zenith_angle = sza%edge_val_( i_sza ), & earth_sun_distance = earth_sun_distance%edge_val_( i_sza ) ) end if + if( associated( heating_output ) ) then + call heating_output%output( i_sza, core, & + heating_rates = heating_rates( i_sza, : , : ), & + time = time%edge_( i_sza ), & + solar_zenith_angle = sza%edge_val_( i_sza ), & + earth_sun_distance = earth_sun_distance%edge_val_( i_sza ) ) + end if end do deallocate( height ) @@ -212,6 +234,7 @@ subroutine run_tuvx( ) deallocate( earth_sun_distance ) if( associated( photo_output ) ) deallocate( photo_output ) if( associated( dose_output ) ) deallocate( dose_output ) + if( associated( heating_output ) ) deallocate( heating_output ) #if MUSICA_USE_OPENMP ! Compare results from threads for fixed solar zenith angle @@ -221,15 +244,20 @@ subroutine run_tuvx( ) allocate( thread_dose_rates( size( dose_rates, 2 ), & size( dose_rates, 3 ), & omp_get_max_threads( ) ) ) + allocate( thread_heating_rates( size( heating_rates, 2 ), & + size( heating_rates, 3 ), & + omp_get_max_threads( ) ) ) !$omp parallel & !$omp shared( threads, thread_photo_rates, thread_dose_rates ) associate( thread => threads( omp_get_thread_num( ) + 1 ), & photos => thread_photo_rates(:,:, omp_get_thread_num( ) + 1 ), & - doses => thread_dose_rates( :,:, omp_get_thread_num( ) + 1 ) ) + doses => thread_dose_rates( :,:, omp_get_thread_num( ) + 1 ), & + heat => thread_heating_rates( :,:, omp_get_thread_num( ) + 1 ) ) call thread%core_%run( 40.0_dk, & 1.0_dk, & photolysis_rate_constants = photos, & - dose_rates = doses ) + dose_rates = doses, & + heating_rates = heat ) end associate !$omp end parallel @@ -249,6 +277,13 @@ subroutine run_tuvx( ) "Thread result mismatch for thread "// & to_char( i_thread ) ) end do + do i_photo = 1, size( thread_heating_rates, 2 ) + call assert_msg( 389419926, & + thread_heating_rates( i_level, i_photo, i_thread ) & + .eq. thread_heating_rates( i_level, i_photo, 1 ), & + "Thread result mismatch for thread "// & + to_char( i_thread ) ) + end do end do end do #endif diff --git a/test/data/heating_rates.json b/test/data/heating_rates.json new file mode 100644 index 00000000..7ff2a056 --- /dev/null +++ b/test/data/heating_rates.json @@ -0,0 +1,100 @@ +{ + "grids" : [ + { + "name": "height", + "type": "equal interval", + "units": "km", + "begins at": 1.0, + "ends at": 5.0, + "cell delta": 1.0 + }, + { + "name": "wavelength", + "type": "equal interval", + "units": "nm", + "begins at": 400.0, + "ends at": 700.0, + "cell delta": 50.0 + } + ], + "profiles": [ + { + "name": "temperature", + "type": "from config file", + "units": "K", + "grid": { + "name": "height", + "units": "km" + }, + "values": [ 200.0, 250.0, 300.0, 350.0, 400.0 ] + }, + { + "name": "extraterrestrial flux", + "type": "from config file", + "units": "photon cm-2 s-1", + "grid": { + "name": "wavelength", + "units": "nm" + }, + "values": [ 1.0e+4, 1.0e+5, 1.0e+6, 1.0e+7, 1.0e+8, 1.0e+9, 1.0e+10 ] + }, + { + "name": "air", + "type": "from config file", + "units": "molecule cm-3", + "grid": { + "name": "height", + "units": "km" + }, + "values": [ 2.5e+19, 2.0e+19, 1.5e+19, 1.0e+19, 5.0e+18 ] + } + ], + "reactions": [ + { + "name": "jfoo", + "cross section": { + "type": "base", + "data": { + "default value": 12.3 + } + }, + "quantum yield": { + "type": "base", + "constant value": 0.75 + }, + "heating": { + "energy term": 2.0 + } + }, + { + "name": "jbar", + "cross section": { + "type": "base", + "data": { + "default value": 45.6 + } + }, + "quantum yield": { + "type": "base", + "constant value": 0.25 + } + }, + { + "name": "jbaz", + "cross section": { + "type": "base", + "data": { + "default value": 78.9 + } + }, + "quantum yield": { + "type": "base", + "constant value": 0.5 + }, + "scaling factor": 1.1, + "heating": { + "energy term": 3000.0 + } + } + ] +} \ No newline at end of file diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index e2e46248..296c17aa 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -19,6 +19,7 @@ add_subdirectory(tuv_doug) # TUV-x tests create_standard_test(NAME grid_warehouse SOURCES grid_warehouse.F90) +create_standard_test(NAME heating_rates SOURCES heating_rates.F90) create_standard_test(NAME la_sr_bands SOURCES la_sr_bands.F90 ) create_standard_test(NAME spherical_geometry SOURCES spherical_geometry.F90 ) diff --git a/test/unit/heating_rates.F90 b/test/unit/heating_rates.F90 new file mode 100644 index 00000000..002821a3 --- /dev/null +++ b/test/unit/heating_rates.F90 @@ -0,0 +1,138 @@ +! Copyright (C) 2024 National Center for Atmospheric Research +! SPDX-License-Identifier: Apache-2.0 +! +program test_heating_rates + + use musica_mpi, only : musica_mpi_init, & + musica_mpi_finalize + use tuvx_heating_rates + + implicit none + + call musica_mpi_init( ) + call test_heating_rates_t( ) + call musica_mpi_finalize( ) + +contains + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !> @brief Test the heating rates + subroutine test_heating_rates_t( ) + + use musica_assert, only : assert + use musica_config, only : config_t + use musica_constants, only : dk => musica_dk + use musica_mpi, only : musica_mpi_bcast, & + musica_mpi_rank, & + MPI_COMM_WORLD + use musica_string, only : string_t + use tuvx_grid_warehouse, only : grid_warehouse_t + use tuvx_la_sr_bands, only : la_sr_bands_t + use tuvx_profile_warehouse, only : profile_warehouse_t + use tuvx_solver, only : radiation_field_t + use tuvx_spherical_geometry, only : spherical_geometry_t + use tuvx_test_utils, only : check_values + + type(heating_rates_t), pointer :: heating_rates + class(grid_warehouse_t), pointer :: grids + class(profile_warehouse_t), pointer :: profiles + + character(len=*), parameter :: Iam = "heating_rates_t tests" + type(config_t) :: config, sub_config, reactions_config + type(string_t), allocatable :: labels(:) + character, allocatable :: buffer(:) + integer :: pos, pack_size, i_height, i_wavelength + integer, parameter :: comm = MPI_COMM_WORLD + real(dk), parameter :: hc = 6.62608e-34_dk * 2.9979e8_dk / 1.e-9_dk + real(dk) :: bde(6,2), actinic_flux(5,6), etfl(6) + real(dk) :: wc(6) = (/ 425.0_dk, 475.0_dk, 525.0_dk, 575.0_dk, 625.0_dk, & + 675.0_dk /) + type(radiation_field_t) :: radiation_field + real(dk) :: calculated_rates(5,2), expected_rates(5,2) + type(la_sr_bands_t) :: la_srb + type(spherical_geometry_t) :: spherical_geometry + + call config%from_file( "test/data/heating_rates.json" ) + call config%get( "grids", sub_config, Iam ) + grids => grid_warehouse_t( sub_config ) + call config%get( "profiles", sub_config, Iam ) + profiles => profile_warehouse_t( sub_config, grids ) + + if( musica_mpi_rank( comm ) == 0 ) then + call config%get( "reactions", reactions_config, Iam ) + call sub_config%empty( ) + call sub_config%add( "reactions", reactions_config, Iam ) + heating_rates => heating_rates_t( sub_config, grids, profiles ) + pack_size = heating_rates%pack_size( comm ) + allocate( buffer( pack_size ) ) + pos = 0 + call heating_rates%mpi_pack( buffer, pos, comm ) + call assert( 534250649, pos <= pack_size ) + end if + + call musica_mpi_bcast( pack_size, comm ) + if( musica_mpi_rank( comm ) .ne. 0 ) allocate( buffer( pack_size ) ) + call musica_mpi_bcast( buffer, comm ) + + if( musica_mpi_rank( comm ) .ne. 0 ) then + pos = 0 + allocate( heating_rates ) + call heating_rates%mpi_unpack( buffer, pos, comm ) + call assert( 192483602, pos <= pack_size ) + end if + deallocate( buffer ) + + ! check labels + labels = heating_rates%labels( ) + call assert( 152892147, size(labels) == 2 ) + call assert( 437272930, labels(1) == "jfoo" ) + call assert( 884640776, labels(2) == "jbaz" ) + + ! check bond dissociation energies + call assert( 613305591, size( heating_rates%heating_parameters_ ) == 2 ) + bde(:,1) = max( 0.0_dk, hc * ( 2.0_dk - wc(:) ) / ( 2.0_dk * wc(:) ) ) + call check_values( heating_rates%heating_parameters_(1)%energy_, & + bde(:,1), 1.0e-4_dk ) + bde(:,2) = max( 0.0_dk, hc * ( 3000.0_dk - wc(:) ) / ( 3000.0_dk * wc(:) ) ) + call check_values( heating_rates%heating_parameters_(2)%energy_, & + bde(:,2), 1.0e-4_dk ) + + ! check calculated heating rates + calculated_rates(:,:) = 0.0_dk + allocate( radiation_field%fdr_(5,6), radiation_field%fdn_(5,6), & + radiation_field%fup_(5,6) ) + do i_wavelength = 1, 6 + etfl(i_wavelength) = 0.5_dk * ( 1.0e3_dk * 10.0_dk**i_wavelength + & + 1.0e4_dk * 10.0_dk**i_wavelength ) + end do + do i_height = 1, 5 + do i_wavelength = 1, 6 + radiation_field%fdr_( i_height, i_wavelength ) = & + 1.0_dk * i_height * i_wavelength + radiation_field%fdn_( i_height, i_wavelength ) = & + 2.0_dk * i_height * i_wavelength + radiation_field%fup_( i_height, i_wavelength ) = & + 3.0_dk * i_height * i_wavelength + actinic_flux( i_height, i_wavelength ) = & + 6.0_dk * i_height * i_wavelength * etfl( i_wavelength ) + expected_rates( i_height, 1 ) = & + dot_product( actinic_flux(i_height,:), bde(:,1) * 12.3_dk * 0.75_dk ) + expected_rates( i_height, 2 ) = 1.1_dk * & + dot_product( actinic_flux(i_height,:), bde(:,2) * 78.9_dk * 0.5_dk ) + end do + end do + call heating_rates%get( la_srb, spherical_geometry, grids, profiles, & + radiation_field, calculated_rates ) + + call check_values( calculated_rates, expected_rates, 1.0e-4_dk ) + + deallocate( grids ) + deallocate( profiles ) + deallocate( heating_rates ) + + end subroutine test_heating_rates_t + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +end program test_heating_rates