From 903e70c2d223ab52ae69b1b929a86c6c4ee1db47 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Tue, 16 Jan 2024 11:52:39 -0800 Subject: [PATCH] reorganize temperature based cross section --- src/cross_sections/CMakeLists.txt | 2 + src/cross_sections/temperature_based.F90 | 556 +----------------- src/cross_sections/util/CMakeLists.txt | 10 + .../util/temperature_parameterization.F90 | 438 ++++++++++++++ src/cross_sections/util/temperature_range.F90 | 159 +++++ test/unit/tuv_doug/JCALC/CMakeLists.txt | 2 + test/unit/tuv_doug/JCALC/XSQY_BRONO2.f | 176 ++++++ test/unit/tuv_doug/JCALC/XSQY_HO2NO2.f | 210 +++++++ test/unit/tuv_doug/driver.F90 | 6 + 9 files changed, 1006 insertions(+), 553 deletions(-) create mode 100644 src/cross_sections/util/CMakeLists.txt create mode 100644 src/cross_sections/util/temperature_parameterization.F90 create mode 100644 src/cross_sections/util/temperature_range.F90 create mode 100644 test/unit/tuv_doug/JCALC/XSQY_BRONO2.f create mode 100644 test/unit/tuv_doug/JCALC/XSQY_HO2NO2.f diff --git a/src/cross_sections/CMakeLists.txt b/src/cross_sections/CMakeLists.txt index 64db6ead..66ab8245 100644 --- a/src/cross_sections/CMakeLists.txt +++ b/src/cross_sections/CMakeLists.txt @@ -31,4 +31,6 @@ target_sources(tuvx_object rayliegh.F90 ) +add_subdirectory(util) + ################################################################################ diff --git a/src/cross_sections/temperature_based.F90 b/src/cross_sections/temperature_based.F90 index 843fdfb0..f7198727 100644 --- a/src/cross_sections/temperature_based.F90 +++ b/src/cross_sections/temperature_based.F90 @@ -1,4 +1,4 @@ -! Copyright (C) 2020 National Center for Atmospheric Research +! Copyright (C) 2020-4 National Center for Atmospheric Research ! SPDX-License-Identifier: Apache-2.0 module tuvx_cross_section_temperature_based @@ -10,98 +10,14 @@ module tuvx_cross_section_temperature_based use musica_constants, only : dk => musica_dk use tuvx_cross_section, only : cross_section_t use tuvx_interpolate, only : interpolator_conserving_t + use tuvx_temperature_parameterization, & + only : temperature_parameterization_t implicit none private public :: cross_section_temperature_based_t - !> Range for temperature-based calculations - type :: temperature_range_t - !> Minimum temperature [K] for inclusion in range - real(kind=dk) :: min_temperature_ = 0.0_dk - !> Maximum temperature [K] for include in range - real(kind=dk) :: max_temperature_ = huge(1.0_dk) - !> Indicates whether to use a fixed temperature for the - !! parameterization calculation. If FALSE, the actual - !! temperature is used. - logical :: is_fixed_ = .false. - !> Fixed temperature [K] to use in paramterization calculation - !! - !! Is only used if is_fixed == TRUE - real(kind=dk) :: fixed_temperature_ = 0.0_dk - contains - !> Returns the number of bytes required to pack the range onto a - !! character buffer - procedure :: pack_size => temperature_range_pack_size - !> Packs the range onto a character buffer - procedure :: mpi_pack => temperature_range_mpi_pack - !> Unpacks a range from a character buffer - procedure :: mpi_unpack => temperature_range_mpi_unpack - end type temperature_range_t - - !> Constructor for temperature_range_t - interface temperature_range_t - module procedure :: temperature_range_constructor - end interface temperature_range_t - - !> Parameters for calculating cross section values based on - !! temperature - !! - !! Cross section elements are calculated as: - !! - !! \f[ - !! 10^{\sum_i{(AA_i + (T-273)*BB_i)*\lambda^{lp_i}}} - !! \f] - !! - !! where \f$\lambda\f$ is the wavelength [nm] and - !! \f$T\f$ is the temperature [K]. - type :: temperature_parameterization_t - integer :: n_sets_ = 0 - real(kind=dk), allocatable :: AA_(:) - real(kind=dk), allocatable :: BB_(:) - real(kind=dk), allocatable :: lp_(:) - !> Base temperature [K] to use in calculations - real(kind=dk) :: base_temperature_ - !> Base wavelength [nm] to use in calcuations - real(kind=dk) :: base_wavelength_ - !> Flag indicating whether cross section algorithm is base 10 (true) - !! or base e (false) - logical :: is_base_10_ - !> Flad indicating whether to subtract base temperature from - !! actual temperature (false) or to subtract actual temperature - !! from base temperature (true) - logical :: is_temperature_inverted_ - !> Minimum wavelength [nm] to calculate values for - real(kind=dk) :: min_wavelength_ - !> Maximum wavelength [nm] to calculate values for - real(kind=dk) :: max_wavelength_ - !> Index of minimum wavelength [nm] to calculate values for - integer :: min_wavelength_index_ - !> Index of maximum wavelength to calculate values for - integer :: max_wavelength_index_ - !> Temperature ranges used in parameterization - type(temperature_range_t), allocatable :: ranges_(:) - contains - !> Merges NetCDF wavelength grid with parameterization grid - procedure :: merge_wavelength_grids - !> Calculate the cross section value for a specific temperature - !! and wavelength - procedure :: calculate => temperature_parameterization_calculate - !> Returns the number of bytes required to pack the parameterization - !! onto a character buffer - procedure :: pack_size => temperature_parameterization_pack_size - !> Packs the parameterization onto a character buffer - procedure :: mpi_pack => temperature_parameterization_mpi_pack - !> Unpacks the parameterization from a character buffer - procedure :: mpi_unpack => temperature_parameterization_mpi_unpack - end type temperature_parameterization_t - - !> Constructor for temperature_parameterization_t - interface temperature_parameterization_t - module procedure :: temperature_parameterization_constructor - end interface temperature_parameterization_t - !> Calculator for temperature-based cross sections type, extends(cross_section_t) :: cross_section_temperature_based_t real(kind=dk), allocatable :: raw_wavelengths_(:) ! [nm] @@ -374,472 +290,6 @@ subroutine mpi_unpack( this, buffer, position, comm ) end subroutine mpi_unpack -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function temperature_range_constructor( config ) result( this ) - ! Constructs temperature range objects - - use musica_assert, only : assert_msg - use musica_config, only : config_t - use musica_string, only : string_t - - type(temperature_range_t) :: this - type(config_t), intent(inout) :: config - - character(len=*), parameter :: my_name = "temperature range constructor" - type(string_t) :: required_keys(0), optional_keys(3) - logical :: found - - optional_keys(1) = "minimum" - optional_keys(2) = "maximum" - optional_keys(3) = "fixed value" - call assert_msg( 355912601, & - config%validate( required_keys, optional_keys ), & - "Bad configuration for temperature range" ) - - call config%get( "minimum", this%min_temperature_, my_name, & - default = 0.0_dk ) - call config%get( "maximum", this%max_temperature_, my_name, & - default = huge(1.0_dk) ) - call config%get( "fixed value", this%fixed_temperature_, my_name, & - found = found ) - this%is_fixed_ = found - - end function temperature_range_constructor - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - integer function temperature_range_pack_size( this, comm ) & - result( pack_size ) - ! Returns the size of a character buffer required to pack the range - - use musica_mpi, only : musica_mpi_pack_size - - class(temperature_range_t), intent(in) :: this ! temperature range to be packed - integer, intent(in) :: comm ! MPI communicator - -#ifdef MUSICA_USE_MPI - pack_size = musica_mpi_pack_size( this%min_temperature_, comm ) + & - musica_mpi_pack_size( this%max_temperature_, comm ) + & - musica_mpi_pack_size( this%is_fixed_, comm ) + & - musica_mpi_pack_size( this%fixed_temperature_, comm ) -#else - pack_size = 0 -#endif - - end function temperature_range_pack_size - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine temperature_range_mpi_pack( this, buffer, position, comm ) - ! Packs the temperature range onto a character buffer - - use musica_assert, only : assert - use musica_mpi, only : musica_mpi_pack - - class(temperature_range_t), intent(in) :: this ! temperature range to be packed - character, intent(inout) :: buffer(:) ! memory buffer - integer, intent(inout) :: position ! current buffer position - integer, intent(in) :: comm ! MPI communicator - -#ifdef MUSICA_USE_MPI - integer :: prev_pos - - prev_pos = position - call musica_mpi_pack( buffer, position, this%min_temperature_, comm ) - call musica_mpi_pack( buffer, position, this%max_temperature_, comm ) - call musica_mpi_pack( buffer, position, this%is_fixed_, comm ) - call musica_mpi_pack( buffer, position, this%fixed_temperature_, comm ) - call assert( 409699380, position - prev_pos <= this%pack_size( comm ) ) -#endif - - end subroutine temperature_range_mpi_pack - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine temperature_range_mpi_unpack( this, buffer, position, comm ) - ! Unpacks a temperature range from a character buffer - - use musica_assert, only : assert - use musica_mpi, only : musica_mpi_unpack - - class(temperature_range_t), intent(out) :: this ! temperature range to be unpacked - character, intent(inout) :: buffer(:) ! memory buffer - integer, intent(inout) :: position ! current buffer position - integer, intent(in) :: comm ! MPI communicator - -#ifdef MUSICA_USE_MPI - integer :: prev_pos - - prev_pos = position - call musica_mpi_unpack( buffer, position, this%min_temperature_, comm ) - call musica_mpi_unpack( buffer, position, this%max_temperature_, comm ) - call musica_mpi_unpack( buffer, position, this%is_fixed_, comm ) - call musica_mpi_unpack( buffer, position, this%fixed_temperature_, comm ) - call assert( 164457757, position - prev_pos <= this%pack_size( comm ) ) -#endif - - end subroutine temperature_range_mpi_unpack - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function temperature_parameterization_constructor( config, wavelengths ) & - result( this ) - ! Constructs temperature_parameterization_t objects - - use musica_assert, only : assert_msg, die_msg - use musica_config, only : config_t - use musica_iterator, only : iterator_t - use musica_string, only : string_t - use tuvx_grid, only : grid_t - - type(temperature_parameterization_t) :: this - type(config_t), intent(inout) :: config - class(grid_t), intent(in) :: wavelengths - - character(len=*), parameter :: my_name = & - "temperature parameterization constructor" - type(string_t) :: required_keys(6), optional_keys(4), exp_base - type(config_t) :: temp_ranges, temp_range - class(iterator_t), pointer :: iter - integer :: i_range - logical :: found - - required_keys(1) = "AA" - required_keys(2) = "BB" - required_keys(3) = "lp" - required_keys(4) = "base temperature" - required_keys(5) = "base wavelength" - required_keys(6) = "logarithm" - optional_keys(1) = "minimum wavelength" - optional_keys(2) = "maximum wavelength" - optional_keys(3) = "temperature ranges" - optional_keys(4) = "invert temperature offset" - call assert_msg( 256315527, & - config%validate( required_keys, optional_keys ), & - "Bad configuration for temperature parameterization." ) - - call config%get( "AA", this%AA_, my_name ) - call config%get( "BB", this%BB_, my_name ) - call config%get( "lp", this%lp_, my_name ) - call config%get( "base temperature", this%base_temperature_, my_name ) - call config%get( "base wavelength", this%base_wavelength_, my_name ) - call config%get( "logarithm", exp_base, my_name ) - call config%get( "invert temperature offset", & - this%is_temperature_inverted_, my_name, default = .false.) - if( exp_base == "base 10" ) then - this%is_base_10_ = .true. - else if( exp_base == "natural" ) then - this%is_base_10_ = .false. - else - call die_msg( 104603249, "Invalid logarithm type in temperature-based"//& - " cross section: '"//exp_base//"'" ) - end if - call assert_msg( 467090427, size( this%AA_ ) == size( this%BB_ ) .and. & - size( this%AA_ ) == size( this%lp_ ), & - "Arrays AA, BB, and lp must be the same size for "// & - "temperature-based cross sections." ) - call config%get( "minimum wavelength", this%min_wavelength_, my_name, & - default = 0.0_dk ) - call config%get( "maximum wavelength", this%max_wavelength_, my_name, & - default = huge(1.0_dk) ) - this%min_wavelength_index_ = 1 - do while( wavelengths%mid_( this%min_wavelength_index_ ) & - < this%min_wavelength_ & - .and. this%min_wavelength_index_ <= wavelengths%ncells_ ) - this%min_wavelength_index_ = this%min_wavelength_index_ + 1 - end do - call assert_msg( 286143383, & - wavelengths%mid_( this%min_wavelength_index_ ) & - >= this%min_wavelength_, & - "Minimum wavelength for temperature-based cross section is "// & - "outside the bounds of the wavelength grid." ) - this%max_wavelength_index_ = wavelengths%ncells_ - do while( wavelengths%mid_( this%max_wavelength_index_ ) & - > this%max_wavelength_ & - .and. this%max_wavelength_index_ >= 1 ) - this%max_wavelength_index_ = this%max_wavelength_index_ - 1 - end do - call assert_msg( 490175140, & - wavelengths%mid_( this%max_wavelength_index_ ) & - <= this%max_wavelength_, & - "Maximum wavelength for temperature-based cross section is "// & - "outside the bounds of the wavelength grid." ) - call config%get( "temperature ranges", temp_ranges, my_name, & - found = found ) - if( .not. found ) then - allocate( this%ranges_( 1 ) ) - return - end if - allocate( this%ranges_( temp_ranges%number_of_children( ) ) ) - iter => temp_ranges%get_iterator( ) - i_range = 0 - do while( iter%next( ) ) - i_range = i_range + 1 - call temp_ranges%get( iter, temp_range, my_name ) - this%ranges_( i_range ) = temperature_range_t( temp_range ) - end do - deallocate( iter ) - - end function temperature_parameterization_constructor - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - function merge_wavelength_grids( this, input_grid, tuv_grid ) & - result( merged_grid ) - ! Merges wavelength grid from NetCDF input data with parameterization - ! grid (same as the TUV-x grid). - ! Where they overlap, the parameterization is used. - ! Updates the parameterization wavelength indices for new grid. - ! Returns merged wavelength grid. - ! - ! NOTE: Uses mid-points on the TUV-x wavelength grid - - use musica_assert, only : assert - use tuvx_grid, only : grid_t - - class(temperature_parameterization_t), intent(inout) :: this - real(kind=dk), intent(in) :: input_grid(:) - class(grid_t), intent(in) :: tuv_grid - real(kind=dk), allocatable :: merged_grid(:) - - logical :: found_min - integer :: i_wl, n_wl, i_input_wl, i_tuv_wl, n_tuv_wl - - if( size( input_grid ) == 0 ) then - merged_grid = tuv_grid%mid_ - return - end if - - associate( wl_min_index => this%min_wavelength_index_, & - wl_max_index => this%max_wavelength_index_, & - min_wl => this%min_wavelength_, & - max_wl => this%max_wavelength_ ) - n_wl = 0 - do i_input_wl = 1, size( input_grid(:) ) - if( min_wl > input_grid( i_input_wl ) .or. & - max_wl < input_grid( i_input_wl ) ) n_wl = n_wl + 1 - end do - i_tuv_wl = wl_min_index - n_tuv_wl = wl_max_index - n_wl = n_wl + ( n_tuv_wl - i_tuv_wl + 1 ) - allocate( merged_grid( n_wl ) ) - i_input_wl = 1 - i_wl = 1 - found_min = .false. - do - if( i_wl > n_wl ) then - ! end of merged grid - exit - else if( i_tuv_wl > n_tuv_wl .and. & - input_grid( i_input_wl ) <= max_wl ) then - ! skipping input data wavelengths in parameterization range - i_input_wl = i_input_wl + 1 - else if( .not. ( min_wl <= input_grid( i_input_wl ) .and. & - max_wl >= input_grid( i_input_wl ) ) ) then - ! adding input data wavelengths outside of parameterization range - merged_grid( i_wl ) = input_grid( i_input_wl ) - i_input_wl = i_input_wl + 1 - i_wl = i_wl + 1 - else if( i_tuv_wl <= n_tuv_wl ) then - ! adding TUV-x wavelengths in parameterization range - ! - ! TODO This follows logic from original TUV, but perhaps should - ! be modified to assign TUV-x wavelength edges - merged_grid( i_wl ) = tuv_grid%mid_( i_tuv_wl ) - if( .not. found_min ) then - found_min = .true. - wl_min_index = i_wl - end if - wl_max_index = i_wl - i_tuv_wl = i_tuv_wl + 1 - i_wl = i_wl + 1 - end if - end do - call assert( 265861594, i_tuv_wl == n_tuv_wl + 1 ) - call assert( 537808229, i_input_wl <= size( input_grid ) + 1 ) - call assert( 422870529, i_wl == n_wl + 1 ) - end associate - - end function merge_wavelength_grids - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine temperature_parameterization_calculate( this, temperature, & - wavelengths, cross_section ) - - use tuvx_profile, only : profile_t - - class(temperature_parameterization_t), intent(in) :: this - real(kind=dk), intent(in) :: temperature - real(kind=dk), intent(in) :: wavelengths(:) - real(kind=dk), intent(inout) :: cross_section(:) - - ! local variables - real(kind=dk) :: temp, temp_xs( size( cross_section ) ) - integer :: i_lp, i_range, w_min, w_max - - w_min = this%min_wavelength_index_ - w_max = this%max_wavelength_index_ - do i_range = 1, size( this%ranges_ ) - associate( temp_range => this%ranges_( i_range ) ) - if( temperature < temp_range%min_temperature_ .or. & - temperature > temp_range%max_temperature_ ) cycle - if( temp_range%is_fixed_ ) then - temp = temp_range%fixed_temperature_ - else - temp = temperature - end if - if ( this%is_temperature_inverted_ ) then - temp = this%base_temperature_ - temp - else - temp = temp - this%base_temperature_ - end if - temp_xs(:) = 0.0_dk - do i_lp = 1, size( this%lp_ ) - temp_xs( w_min:w_max ) = temp_xs( w_min:w_max ) + & - ( this%AA_( i_lp ) + temp * this%BB_( i_lp ) ) * & - ( wavelengths( w_min:w_max ) & - - this%base_wavelength_ )**this%lp_( i_lp ) - end do - if (this%is_base_10_) then - cross_section( w_min:w_max ) = cross_section( w_min:w_max ) & - + 10**temp_xs( w_min:w_max ) - else - cross_section( w_min:w_max ) = cross_section( w_min:w_max ) & - + exp( temp_xs( w_min:w_max ) ) - end if - end associate - end do - - end subroutine temperature_parameterization_calculate - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - integer function temperature_parameterization_pack_size( this, comm ) & - result( pack_size ) - ! Returns the size of a character buffer required to pack the - ! parameterization - - use musica_mpi, only : musica_mpi_pack_size - - class(temperature_parameterization_t), intent(in) :: this ! parameterization to be packed - integer, intent(in) :: comm ! MPI communicator - -#ifdef MUSICA_USE_MPI - integer :: i_range - - pack_size = musica_mpi_pack_size( this%AA_, comm ) + & - musica_mpi_pack_size( this%BB_, comm ) + & - musica_mpi_pack_size( this%lp_, comm ) + & - musica_mpi_pack_size( this%base_temperature_, comm ) + & - musica_mpi_pack_size( this%base_wavelength_, comm ) + & - musica_mpi_pack_size( this%is_base_10_, comm ) + & - musica_mpi_pack_size( this%is_temperature_inverted_, comm ) + & - musica_mpi_pack_size( this%min_wavelength_, comm ) + & - musica_mpi_pack_size( this%max_wavelength_, comm ) + & - musica_mpi_pack_size( this%min_wavelength_index_, comm ) + & - musica_mpi_pack_size( this%max_wavelength_index_, comm ) + & - musica_mpi_pack_size( allocated( this%ranges_ ), comm ) - if( allocated( this%ranges_ ) ) then - pack_size = pack_size + & - musica_mpi_pack_size( size( this%ranges_ ), comm ) - do i_range = 1, size( this%ranges_ ) - pack_size = pack_size + this%ranges_( i_range )%pack_size( comm ) - end do - end if -#else - pack_size = 0 -#endif - - end function temperature_parameterization_pack_size - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine temperature_parameterization_mpi_pack( this, buffer, position, & - comm ) - ! Packs the parameterization onto a character buffer - - use musica_assert, only : assert - use musica_mpi, only : musica_mpi_pack - - class(temperature_parameterization_t), intent(in) :: this ! parameterization to be packed - character, intent(inout) :: buffer(:) ! memory buffer - integer, intent(inout) :: position ! current buffer position - integer, intent(in) :: comm ! MPI communicator - -#ifdef MUSICA_USE_MPI - integer :: prev_pos, i_range - - prev_pos = position - call musica_mpi_pack( buffer, position, this%AA_, comm ) - call musica_mpi_pack( buffer, position, this%BB_, comm ) - call musica_mpi_pack( buffer, position, this%lp_, comm ) - call musica_mpi_pack( buffer, position, this%base_temperature_, comm ) - call musica_mpi_pack( buffer, position, this%base_wavelength_, comm ) - call musica_mpi_pack( buffer, position, this%is_base_10_, comm ) - call musica_mpi_pack( buffer, position, this%is_temperature_inverted_, & - comm ) - call musica_mpi_pack( buffer, position, this%min_wavelength_, comm ) - call musica_mpi_pack( buffer, position, this%max_wavelength_, comm ) - call musica_mpi_pack( buffer, position, this%min_wavelength_index_, comm ) - call musica_mpi_pack( buffer, position, this%max_wavelength_index_, comm ) - call musica_mpi_pack( buffer, position, allocated( this%ranges_ ), comm ) - if( allocated( this%ranges_ ) ) then - call musica_mpi_pack( buffer, position, size( this%ranges_ ), comm ) - do i_range = 1, size( this%ranges_ ) - call this%ranges_( i_range )%mpi_pack( buffer, position, comm ) - end do - end if - call assert( 267439201, position - prev_pos <= this%pack_size( comm ) ) -#endif - - end subroutine temperature_parameterization_mpi_pack - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - subroutine temperature_parameterization_mpi_unpack( this, buffer, position, & - comm ) - ! Unpacks a parameterization from a character buffer - - use musica_assert, only : assert - use musica_mpi, only : musica_mpi_unpack - - class(temperature_parameterization_t), intent(out) :: this ! parameterization to be unpacked - character, intent(inout) :: buffer(:) ! memory buffer - integer, intent(inout) :: position ! current buffer position - integer, intent(in) :: comm ! MPI communicator - -#ifdef MUSICA_USE_MPI - integer :: prev_pos, i_range, n_ranges - logical :: alloced - - prev_pos = position - call musica_mpi_unpack( buffer, position, this%AA_, comm ) - call musica_mpi_unpack( buffer, position, this%BB_, comm ) - call musica_mpi_unpack( buffer, position, this%lp_, comm ) - call musica_mpi_unpack( buffer, position, this%base_temperature_, comm ) - call musica_mpi_unpack( buffer, position, this%base_wavelength_, comm ) - call musica_mpi_unpack( buffer, position, this%is_base_10_, comm ) - call musica_mpi_unpack( buffer, position, this%is_temperature_inverted_, & - comm ) - call musica_mpi_unpack( buffer, position, this%min_wavelength_, comm ) - call musica_mpi_unpack( buffer, position, this%max_wavelength_, comm ) - call musica_mpi_unpack( buffer, position, this%min_wavelength_index_,comm ) - call musica_mpi_unpack( buffer, position, this%max_wavelength_index_,comm ) - call musica_mpi_unpack( buffer, position, alloced, comm ) - if( alloced ) then - call musica_mpi_unpack( buffer, position, n_ranges, comm ) - allocate( this%ranges_( n_ranges ) ) - do i_range = 1, size( this%ranges_ ) - call this%ranges_( i_range )%mpi_unpack( buffer, position, comm ) - end do - end if - call assert( 483905106, position - prev_pos <= this%pack_size( comm ) ) -#endif - - end subroutine temperature_parameterization_mpi_unpack - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end module tuvx_cross_section_temperature_based diff --git a/src/cross_sections/util/CMakeLists.txt b/src/cross_sections/util/CMakeLists.txt new file mode 100644 index 00000000..7af19114 --- /dev/null +++ b/src/cross_sections/util/CMakeLists.txt @@ -0,0 +1,10 @@ +################################################################################ +# utilities for cross section parameterizations + +target_sources(tuvx_object + PRIVATE + temperature_parameterization.F90 + temperature_range.F90 +) + +################################################################################ \ No newline at end of file diff --git a/src/cross_sections/util/temperature_parameterization.F90 b/src/cross_sections/util/temperature_parameterization.F90 new file mode 100644 index 00000000..7d44bf72 --- /dev/null +++ b/src/cross_sections/util/temperature_parameterization.F90 @@ -0,0 +1,438 @@ +! Copyright (C) 2020-4 National Center for Atmospheric Research +! SPDX-License-Identifier: Apache-2.0 + +module tuvx_temperature_parameterization +! Calculates cross-section elements based on a temperature parameterization + + ! Including musica_config at the module level to avoid an ICE + ! with Intel 2022.1 compiler + use musica_config, only : config_t + use musica_constants, only : dk => musica_dk + use tuvx_temperature_range, only : temperature_range_t + + implicit none + + private + public :: temperature_parameterization_t + + !> Parameters for calculating cross section values based on + !! temperature + !! + !! Cross section elements are calculated as: + !! + !! \f[ + !! 10^{\sum_i{(AA_i + (T-273)*BB_i)*\lambda^{lp_i}}} + !! \f] + !! + !! where \f$\lambda\f$ is the wavelength [nm] and + !! \f$T\f$ is the temperature [K]. + type :: temperature_parameterization_t + integer :: n_sets_ = 0 + real(kind=dk), allocatable :: AA_(:) + real(kind=dk), allocatable :: BB_(:) + real(kind=dk), allocatable :: lp_(:) + !> Base temperature [K] to use in calculations + real(kind=dk) :: base_temperature_ + !> Base wavelength [nm] to use in calcuations + real(kind=dk) :: base_wavelength_ + !> Flag indicating whether cross section algorithm is base 10 (true) + !! or base e (false) + logical :: is_base_10_ + !> Flad indicating whether to subtract base temperature from + !! actual temperature (false) or to subtract actual temperature + !! from base temperature (true) + logical :: is_temperature_inverted_ + !> Minimum wavelength [nm] to calculate values for + real(kind=dk) :: min_wavelength_ + !> Maximum wavelength [nm] to calculate values for + real(kind=dk) :: max_wavelength_ + !> Index of minimum wavelength [nm] to calculate values for + integer :: min_wavelength_index_ + !> Index of maximum wavelength to calculate values for + integer :: max_wavelength_index_ + !> Temperature ranges used in parameterization + type(temperature_range_t), allocatable :: ranges_(:) + contains + !> Merges NetCDF wavelength grid with parameterization grid + procedure :: merge_wavelength_grids + !> Calculate the cross section value for a specific temperature + !! and wavelength + procedure :: calculate => temperature_parameterization_calculate + !> Returns the number of bytes required to pack the parameterization + !! onto a character buffer + procedure :: pack_size => temperature_parameterization_pack_size + !> Packs the parameterization onto a character buffer + procedure :: mpi_pack => temperature_parameterization_mpi_pack + !> Unpacks the parameterization from a character buffer + procedure :: mpi_unpack => temperature_parameterization_mpi_unpack + end type temperature_parameterization_t + + !> Constructor for temperature_parameterization_t + interface temperature_parameterization_t + module procedure :: temperature_parameterization_constructor + end interface temperature_parameterization_t + +contains + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function temperature_parameterization_constructor( config, wavelengths ) & + result( this ) + ! Constructs temperature_parameterization_t objects + + use musica_assert, only : assert_msg, die_msg + use musica_config, only : config_t + use musica_iterator, only : iterator_t + use musica_string, only : string_t + use tuvx_grid, only : grid_t + + type(temperature_parameterization_t) :: this + type(config_t), intent(inout) :: config + class(grid_t), intent(in) :: wavelengths + + character(len=*), parameter :: my_name = & + "temperature parameterization constructor" + type(string_t) :: required_keys(6), optional_keys(4), exp_base + type(config_t) :: temp_ranges, temp_range + class(iterator_t), pointer :: iter + integer :: i_range + logical :: found + + required_keys(1) = "AA" + required_keys(2) = "BB" + required_keys(3) = "lp" + required_keys(4) = "base temperature" + required_keys(5) = "base wavelength" + required_keys(6) = "logarithm" + optional_keys(1) = "minimum wavelength" + optional_keys(2) = "maximum wavelength" + optional_keys(3) = "temperature ranges" + optional_keys(4) = "invert temperature offset" + call assert_msg( 256315527, & + config%validate( required_keys, optional_keys ), & + "Bad configuration for temperature parameterization." ) + + call config%get( "AA", this%AA_, my_name ) + call config%get( "BB", this%BB_, my_name ) + call config%get( "lp", this%lp_, my_name ) + call config%get( "base temperature", this%base_temperature_, my_name ) + call config%get( "base wavelength", this%base_wavelength_, my_name ) + call config%get( "logarithm", exp_base, my_name ) + call config%get( "invert temperature offset", & + this%is_temperature_inverted_, my_name, default = .false.) + if( exp_base == "base 10" ) then + this%is_base_10_ = .true. + else if( exp_base == "natural" ) then + this%is_base_10_ = .false. + else + call die_msg( 104603249, "Invalid logarithm type in temperature-based"//& + " cross section: '"//exp_base//"'" ) + end if + call assert_msg( 467090427, size( this%AA_ ) == size( this%BB_ ) .and. & + size( this%AA_ ) == size( this%lp_ ), & + "Arrays AA, BB, and lp must be the same size for "// & + "temperature-based cross sections." ) + call config%get( "minimum wavelength", this%min_wavelength_, my_name, & + default = 0.0_dk ) + call config%get( "maximum wavelength", this%max_wavelength_, my_name, & + default = huge(1.0_dk) ) + this%min_wavelength_index_ = 1 + do while( wavelengths%mid_( this%min_wavelength_index_ ) & + < this%min_wavelength_ & + .and. this%min_wavelength_index_ <= wavelengths%ncells_ ) + this%min_wavelength_index_ = this%min_wavelength_index_ + 1 + end do + call assert_msg( 286143383, & + wavelengths%mid_( this%min_wavelength_index_ ) & + >= this%min_wavelength_, & + "Minimum wavelength for temperature-based cross section is "// & + "outside the bounds of the wavelength grid." ) + this%max_wavelength_index_ = wavelengths%ncells_ + do while( wavelengths%mid_( this%max_wavelength_index_ ) & + > this%max_wavelength_ & + .and. this%max_wavelength_index_ >= 1 ) + this%max_wavelength_index_ = this%max_wavelength_index_ - 1 + end do + call assert_msg( 490175140, & + wavelengths%mid_( this%max_wavelength_index_ ) & + <= this%max_wavelength_, & + "Maximum wavelength for temperature-based cross section is "// & + "outside the bounds of the wavelength grid." ) + call config%get( "temperature ranges", temp_ranges, my_name, & + found = found ) + if( .not. found ) then + allocate( this%ranges_( 1 ) ) + return + end if + allocate( this%ranges_( temp_ranges%number_of_children( ) ) ) + iter => temp_ranges%get_iterator( ) + i_range = 0 + do while( iter%next( ) ) + i_range = i_range + 1 + call temp_ranges%get( iter, temp_range, my_name ) + this%ranges_( i_range ) = temperature_range_t( temp_range ) + end do + deallocate( iter ) + + end function temperature_parameterization_constructor + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function merge_wavelength_grids( this, input_grid, tuv_grid ) & + result( merged_grid ) + ! Merges wavelength grid from NetCDF input data with parameterization + ! grid (same as the TUV-x grid). + ! Where they overlap, the parameterization is used. + ! Updates the parameterization wavelength indices for new grid. + ! Returns merged wavelength grid. + ! + ! NOTE: Uses mid-points on the TUV-x wavelength grid + + use musica_assert, only : assert + use tuvx_grid, only : grid_t + + class(temperature_parameterization_t), intent(inout) :: this + real(kind=dk), intent(in) :: input_grid(:) + class(grid_t), intent(in) :: tuv_grid + real(kind=dk), allocatable :: merged_grid(:) + + logical :: found_min + integer :: i_wl, n_wl, i_input_wl, i_tuv_wl, n_tuv_wl + + if( size( input_grid ) == 0 ) then + merged_grid = tuv_grid%mid_ + return + end if + + associate( wl_min_index => this%min_wavelength_index_, & + wl_max_index => this%max_wavelength_index_, & + min_wl => this%min_wavelength_, & + max_wl => this%max_wavelength_ ) + n_wl = 0 + do i_input_wl = 1, size( input_grid(:) ) + if( min_wl > input_grid( i_input_wl ) .or. & + max_wl < input_grid( i_input_wl ) ) n_wl = n_wl + 1 + end do + i_tuv_wl = wl_min_index + n_tuv_wl = wl_max_index + n_wl = n_wl + ( n_tuv_wl - i_tuv_wl + 1 ) + allocate( merged_grid( n_wl ) ) + i_input_wl = 1 + i_wl = 1 + found_min = .false. + do + if( i_wl > n_wl ) then + ! end of merged grid + exit + else if( i_tuv_wl > n_tuv_wl .and. & + input_grid( i_input_wl ) <= max_wl ) then + ! skipping input data wavelengths in parameterization range + i_input_wl = i_input_wl + 1 + else if( .not. ( min_wl <= input_grid( i_input_wl ) .and. & + max_wl >= input_grid( i_input_wl ) ) ) then + ! adding input data wavelengths outside of parameterization range + merged_grid( i_wl ) = input_grid( i_input_wl ) + i_input_wl = i_input_wl + 1 + i_wl = i_wl + 1 + else if( i_tuv_wl <= n_tuv_wl ) then + ! adding TUV-x wavelengths in parameterization range + ! + ! TODO This follows logic from original TUV, but perhaps should + ! be modified to assign TUV-x wavelength edges + merged_grid( i_wl ) = tuv_grid%mid_( i_tuv_wl ) + if( .not. found_min ) then + found_min = .true. + wl_min_index = i_wl + end if + wl_max_index = i_wl + i_tuv_wl = i_tuv_wl + 1 + i_wl = i_wl + 1 + end if + end do + call assert( 265861594, i_tuv_wl == n_tuv_wl + 1 ) + call assert( 537808229, i_input_wl <= size( input_grid ) + 1 ) + call assert( 422870529, i_wl == n_wl + 1 ) + end associate + + end function merge_wavelength_grids + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine temperature_parameterization_calculate( this, temperature, & + wavelengths, cross_section ) + + use tuvx_profile, only : profile_t + + class(temperature_parameterization_t), intent(in) :: this + real(kind=dk), intent(in) :: temperature + real(kind=dk), intent(in) :: wavelengths(:) + real(kind=dk), intent(inout) :: cross_section(:) + + ! local variables + real(kind=dk) :: temp, temp_xs( size( cross_section ) ) + integer :: i_lp, i_range, w_min, w_max + + w_min = this%min_wavelength_index_ + w_max = this%max_wavelength_index_ + do i_range = 1, size( this%ranges_ ) + associate( temp_range => this%ranges_( i_range ) ) + if( temperature < temp_range%min_temperature_ .or. & + temperature > temp_range%max_temperature_ ) cycle + if( temp_range%is_fixed_ ) then + temp = temp_range%fixed_temperature_ + else + temp = temperature + end if + if ( this%is_temperature_inverted_ ) then + temp = this%base_temperature_ - temp + else + temp = temp - this%base_temperature_ + end if + temp_xs(:) = 0.0_dk + do i_lp = 1, size( this%lp_ ) + temp_xs( w_min:w_max ) = temp_xs( w_min:w_max ) + & + ( this%AA_( i_lp ) + temp * this%BB_( i_lp ) ) * & + ( wavelengths( w_min:w_max ) & + - this%base_wavelength_ )**this%lp_( i_lp ) + end do + if (this%is_base_10_) then + cross_section( w_min:w_max ) = cross_section( w_min:w_max ) & + + 10**temp_xs( w_min:w_max ) + else + cross_section( w_min:w_max ) = cross_section( w_min:w_max ) & + + exp( temp_xs( w_min:w_max ) ) + end if + end associate + end do + + end subroutine temperature_parameterization_calculate + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + integer function temperature_parameterization_pack_size( this, comm ) & + result( pack_size ) + ! Returns the size of a character buffer required to pack the + ! parameterization + + use musica_mpi, only : musica_mpi_pack_size + + class(temperature_parameterization_t), intent(in) :: this ! parameterization to be packed + integer, intent(in) :: comm ! MPI communicator + +#ifdef MUSICA_USE_MPI + integer :: i_range + + pack_size = musica_mpi_pack_size( this%AA_, comm ) + & + musica_mpi_pack_size( this%BB_, comm ) + & + musica_mpi_pack_size( this%lp_, comm ) + & + musica_mpi_pack_size( this%base_temperature_, comm ) + & + musica_mpi_pack_size( this%base_wavelength_, comm ) + & + musica_mpi_pack_size( this%is_base_10_, comm ) + & + musica_mpi_pack_size( this%is_temperature_inverted_, comm ) + & + musica_mpi_pack_size( this%min_wavelength_, comm ) + & + musica_mpi_pack_size( this%max_wavelength_, comm ) + & + musica_mpi_pack_size( this%min_wavelength_index_, comm ) + & + musica_mpi_pack_size( this%max_wavelength_index_, comm ) + & + musica_mpi_pack_size( allocated( this%ranges_ ), comm ) + if( allocated( this%ranges_ ) ) then + pack_size = pack_size + & + musica_mpi_pack_size( size( this%ranges_ ), comm ) + do i_range = 1, size( this%ranges_ ) + pack_size = pack_size + this%ranges_( i_range )%pack_size( comm ) + end do + end if +#else + pack_size = 0 +#endif + + end function temperature_parameterization_pack_size + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine temperature_parameterization_mpi_pack( this, buffer, position, & + comm ) + ! Packs the parameterization onto a character buffer + + use musica_assert, only : assert + use musica_mpi, only : musica_mpi_pack + + class(temperature_parameterization_t), intent(in) :: this ! parameterization to be packed + character, intent(inout) :: buffer(:) ! memory buffer + integer, intent(inout) :: position ! current buffer position + integer, intent(in) :: comm ! MPI communicator + +#ifdef MUSICA_USE_MPI + integer :: prev_pos, i_range + + prev_pos = position + call musica_mpi_pack( buffer, position, this%AA_, comm ) + call musica_mpi_pack( buffer, position, this%BB_, comm ) + call musica_mpi_pack( buffer, position, this%lp_, comm ) + call musica_mpi_pack( buffer, position, this%base_temperature_, comm ) + call musica_mpi_pack( buffer, position, this%base_wavelength_, comm ) + call musica_mpi_pack( buffer, position, this%is_base_10_, comm ) + call musica_mpi_pack( buffer, position, this%is_temperature_inverted_, & + comm ) + call musica_mpi_pack( buffer, position, this%min_wavelength_, comm ) + call musica_mpi_pack( buffer, position, this%max_wavelength_, comm ) + call musica_mpi_pack( buffer, position, this%min_wavelength_index_, comm ) + call musica_mpi_pack( buffer, position, this%max_wavelength_index_, comm ) + call musica_mpi_pack( buffer, position, allocated( this%ranges_ ), comm ) + if( allocated( this%ranges_ ) ) then + call musica_mpi_pack( buffer, position, size( this%ranges_ ), comm ) + do i_range = 1, size( this%ranges_ ) + call this%ranges_( i_range )%mpi_pack( buffer, position, comm ) + end do + end if + call assert( 267439201, position - prev_pos <= this%pack_size( comm ) ) +#endif + + end subroutine temperature_parameterization_mpi_pack + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine temperature_parameterization_mpi_unpack( this, buffer, position, & + comm ) + ! Unpacks a parameterization from a character buffer + + use musica_assert, only : assert + use musica_mpi, only : musica_mpi_unpack + + class(temperature_parameterization_t), intent(out) :: this ! parameterization to be unpacked + character, intent(inout) :: buffer(:) ! memory buffer + integer, intent(inout) :: position ! current buffer position + integer, intent(in) :: comm ! MPI communicator + +#ifdef MUSICA_USE_MPI + integer :: prev_pos, i_range, n_ranges + logical :: alloced + + prev_pos = position + call musica_mpi_unpack( buffer, position, this%AA_, comm ) + call musica_mpi_unpack( buffer, position, this%BB_, comm ) + call musica_mpi_unpack( buffer, position, this%lp_, comm ) + call musica_mpi_unpack( buffer, position, this%base_temperature_, comm ) + call musica_mpi_unpack( buffer, position, this%base_wavelength_, comm ) + call musica_mpi_unpack( buffer, position, this%is_base_10_, comm ) + call musica_mpi_unpack( buffer, position, this%is_temperature_inverted_, & + comm ) + call musica_mpi_unpack( buffer, position, this%min_wavelength_, comm ) + call musica_mpi_unpack( buffer, position, this%max_wavelength_, comm ) + call musica_mpi_unpack( buffer, position, this%min_wavelength_index_,comm ) + call musica_mpi_unpack( buffer, position, this%max_wavelength_index_,comm ) + call musica_mpi_unpack( buffer, position, alloced, comm ) + if( alloced ) then + call musica_mpi_unpack( buffer, position, n_ranges, comm ) + allocate( this%ranges_( n_ranges ) ) + do i_range = 1, size( this%ranges_ ) + call this%ranges_( i_range )%mpi_unpack( buffer, position, comm ) + end do + end if + call assert( 483905106, position - prev_pos <= this%pack_size( comm ) ) +#endif + + end subroutine temperature_parameterization_mpi_unpack + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +end module tuvx_temperature_parameterization \ No newline at end of file diff --git a/src/cross_sections/util/temperature_range.F90 b/src/cross_sections/util/temperature_range.F90 new file mode 100644 index 00000000..7406bb5e --- /dev/null +++ b/src/cross_sections/util/temperature_range.F90 @@ -0,0 +1,159 @@ +! Copyright (C) 2020-4 National Center for Atmospheric Research +! SPDX-License-Identifier: Apache-2.0 + +module tuvx_temperature_range +! Defines a temperature range for use in temperature-based cross +! section parameterizations + + ! Including musica_config at the module level to avoid an ICE + ! with Intel 2022.1 compiler + use musica_config, only : config_t + use musica_constants, only : dk => musica_dk + + implicit none + + private + public :: temperature_range_t + + + !> Range for temperature-based calculations + type :: temperature_range_t + !> Minimum temperature [K] for inclusion in range + real(kind=dk) :: min_temperature_ = 0.0_dk + !> Maximum temperature [K] for include in range + real(kind=dk) :: max_temperature_ = huge(1.0_dk) + !> Indicates whether to use a fixed temperature for the + !! parameterization calculation. If FALSE, the actual + !! temperature is used. + logical :: is_fixed_ = .false. + !> Fixed temperature [K] to use in paramterization calculation + !! + !! Is only used if is_fixed == TRUE + real(kind=dk) :: fixed_temperature_ = 0.0_dk + contains + !> Returns the number of bytes required to pack the range onto a + !! character buffer + procedure :: pack_size => temperature_range_pack_size + !> Packs the range onto a character buffer + procedure :: mpi_pack => temperature_range_mpi_pack + !> Unpacks a range from a character buffer + procedure :: mpi_unpack => temperature_range_mpi_unpack + end type temperature_range_t + + !> Constructor for temperature_range_t + interface temperature_range_t + module procedure :: temperature_range_constructor + end interface temperature_range_t + +contains + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function temperature_range_constructor( config ) result( this ) + ! Constructs temperature range objects + + use musica_assert, only : assert_msg + use musica_config, only : config_t + use musica_string, only : string_t + + type(temperature_range_t) :: this + type(config_t), intent(inout) :: config + + character(len=*), parameter :: my_name = "temperature range constructor" + type(string_t) :: required_keys(0), optional_keys(3) + logical :: found + + optional_keys(1) = "minimum" + optional_keys(2) = "maximum" + optional_keys(3) = "fixed value" + call assert_msg( 355912601, & + config%validate( required_keys, optional_keys ), & + "Bad configuration for temperature range" ) + + call config%get( "minimum", this%min_temperature_, my_name, & + default = 0.0_dk ) + call config%get( "maximum", this%max_temperature_, my_name, & + default = huge(1.0_dk) ) + call config%get( "fixed value", this%fixed_temperature_, my_name, & + found = found ) + this%is_fixed_ = found + + end function temperature_range_constructor + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + integer function temperature_range_pack_size( this, comm ) & + result( pack_size ) + ! Returns the size of a character buffer required to pack the range + + use musica_mpi, only : musica_mpi_pack_size + + class(temperature_range_t), intent(in) :: this ! temperature range to be packed + integer, intent(in) :: comm ! MPI communicator + +#ifdef MUSICA_USE_MPI + pack_size = musica_mpi_pack_size( this%min_temperature_, comm ) + & + musica_mpi_pack_size( this%max_temperature_, comm ) + & + musica_mpi_pack_size( this%is_fixed_, comm ) + & + musica_mpi_pack_size( this%fixed_temperature_, comm ) +#else + pack_size = 0 +#endif + + end function temperature_range_pack_size + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine temperature_range_mpi_pack( this, buffer, position, comm ) + ! Packs the temperature range onto a character buffer + + use musica_assert, only : assert + use musica_mpi, only : musica_mpi_pack + + class(temperature_range_t), intent(in) :: this ! temperature range to be packed + character, intent(inout) :: buffer(:) ! memory buffer + integer, intent(inout) :: position ! current buffer position + integer, intent(in) :: comm ! MPI communicator + +#ifdef MUSICA_USE_MPI + integer :: prev_pos + + prev_pos = position + call musica_mpi_pack( buffer, position, this%min_temperature_, comm ) + call musica_mpi_pack( buffer, position, this%max_temperature_, comm ) + call musica_mpi_pack( buffer, position, this%is_fixed_, comm ) + call musica_mpi_pack( buffer, position, this%fixed_temperature_, comm ) + call assert( 409699380, position - prev_pos <= this%pack_size( comm ) ) +#endif + + end subroutine temperature_range_mpi_pack + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine temperature_range_mpi_unpack( this, buffer, position, comm ) + ! Unpacks a temperature range from a character buffer + + use musica_assert, only : assert + use musica_mpi, only : musica_mpi_unpack + + class(temperature_range_t), intent(out) :: this ! temperature range to be unpacked + character, intent(inout) :: buffer(:) ! memory buffer + integer, intent(inout) :: position ! current buffer position + integer, intent(in) :: comm ! MPI communicator + +#ifdef MUSICA_USE_MPI + integer :: prev_pos + + prev_pos = position + call musica_mpi_unpack( buffer, position, this%min_temperature_, comm ) + call musica_mpi_unpack( buffer, position, this%max_temperature_, comm ) + call musica_mpi_unpack( buffer, position, this%is_fixed_, comm ) + call musica_mpi_unpack( buffer, position, this%fixed_temperature_, comm ) + call assert( 164457757, position - prev_pos <= this%pack_size( comm ) ) +#endif + + end subroutine temperature_range_mpi_unpack + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +end module tuvx_temperature_range \ No newline at end of file diff --git a/test/unit/tuv_doug/JCALC/CMakeLists.txt b/test/unit/tuv_doug/JCALC/CMakeLists.txt index 634894a2..8f251b71 100644 --- a/test/unit/tuv_doug/JCALC/CMakeLists.txt +++ b/test/unit/tuv_doug/JCALC/CMakeLists.txt @@ -4,6 +4,7 @@ target_sources(tuv_doug PRIVATE XSQY_BRO.f + XSQY_BRONO2.f XSQY_CF2CL2.f XSQY_CFC113.f XSQY_CFC114.f @@ -21,6 +22,7 @@ target_sources(tuv_doug XSQY_HCFC141b.f XSQY_HCFC142b.f XSQY_HNO3.f + XSQY_HO2NO2.f XSQY_N2O5.f ) diff --git a/test/unit/tuv_doug/JCALC/XSQY_BRONO2.f b/test/unit/tuv_doug/JCALC/XSQY_BRONO2.f new file mode 100644 index 00000000..5083671a --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_BRONO2.f @@ -0,0 +1,176 @@ + subroutine XSQY_BRONO2(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product (cross section) x (quantum yield) for brono2 photolysis: ! +! BrONO2 + hv -> products ! +! ! +! cross section: jpl 06 recommendation ! +! quantum yield: jpl 06 recommendation ! +!-----------------------------------------------------------------------------! +! parameters: ! +! nw - integer, number of specified intervals + 1 in working (i) ! +! wavelength grid ! +! wl - real, vector of lower limits of wavelength intervals in (i) ! +! working wavelength grid ! +! wc - real, vector of center points of wavelength intervals in (i) ! +! working wavelength grid ! +! nz - integer, number of altitude levels in working altitude grid (i) ! +! tlev - real, temperature (k) at each specified altitude level (i) ! +! airlev - real, air density (molec/cc) at each altitude level (i) ! +! j - integer, counter for number of weighting functions defined (io) ! +! sq - real, cross section x quantum yield (cm^2) for each (o) ! +! photolysis reaction defined, at each defined wavelength and ! +! at each defined altitude level ! +! jlabel - character*60, string identifier for each photolysis reaction (o) ! +! defined ! +!-----------------------------------------------------------------------------! +! edit history: ! +! 07/27/07: Doug Kinnison ! +!-----------------------------------------------------------------------------! + implicit none + include 'params' + +!-----------------------------------------------------------------------------! +! ... input ! +!-----------------------------------------------------------------------------! + real, intent(in) :: wl(kw) + real, intent(in) :: wc(kw) + real, intent(in) :: tlev(kz) + real, intent(in) :: airlev(kz) + + integer, intent(in) :: nz + integer, intent(in) :: nw + + character*80, intent(in) :: pn + character*60, intent(out) :: jlabel(kj) + + real, intent(out) :: sq(kj,kz,kw) + +!-----------------------------------------------------------------------------! +! ... input/output ! +!-----------------------------------------------------------------------------! + integer, intent(inout) :: j + +!-----------------------------------------------------------------------------! +! ... local ! +!-----------------------------------------------------------------------------! + integer kdata + parameter(kdata=200) + integer i, iw, n, n1, idum, ierr, iz + real x1 (kdata), y1(kdata) + real xin (kdata) + real a1 (kdata), a2(kdata) + real ytmp(nz,kdata) + real ytd (nz,kw) + real yg1 (kw) + real tin (nz) + real qy1 + real qy2 +!----------------------------------------------- +! ... tin set to tlev +!----------------------------------------------- + tin(:) = tlev(:) + +!----------------------------------------------- +! ... jlabel(j) = 'BrONO2 + hv -> Br + NO3' +! ... jlabel(j) = 'BrONO2 + hv -> BrO + NO2' +!----------------------------------------------- + j = j+1 + jlabel(j) = 'BrONO2 + hv -> Br + NO3' + +!----------------------------------------------- +! ... cross sections from JPL06 +!----------------------------------------------- + open(kin, + & file=TRIM(pn)//'XS_BRONO2_JPL06.txt',status='old') + + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + + do iw = 1, n + read(kin,*) xin(iw), y1(iw) + enddo + close(kin) + +!----------------------------------------------- +! ... Read in temperature dep coeff +!----------------------------------------------- + open(kin, + & file=TRIM(pn)//'XS_BRONO2_td_JPL06.txt',status='old') + + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + + do iw = 1, n + read(kin,*) xin(iw), a1(iw), a2(iw) + enddo + close(kin) + +!----------------------------------------------- +! ... derive T-dep (200-296K) +!----------------------------------------------- + do iz = 1, nz + do iw = 1 , n + if ((tin(iz) .GE. 200.) .AND. (tin(iz) .LE. 296.)) Then + ytmp(iz,iw) = y1(iw) * + & ( 1. + + & A1(iw)* (tin(iz)-296.) + + & A2(iw)*((tin(iz)-296.)**2) + & ) + endif + if (tin(iz) .LT. 200.) then + ytmp(iz,iw) = y1(iw) * + & ( 1. + + & A1(iw)* (200.-296.) + + & A2(iw)*((200.-296.)**2) + & ) + endif + if (tin(iz) .GT. 296.) then + ytmp(iz,iw) = y1(iw) + endif + enddo + enddo + +!----------------------------------------------- +! ... Interpolate +!----------------------------------------------- + do iz = 1, nz + n1 = n + y1 = ytmp(iz,:) + x1 = xin + + call addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n1, 0.,0.) + call addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n1, 1e38,0.) + call inter2(nw,wl,yg1,n1,x1,y1,ierr) + ytd(iz,:) = yg1(:) + + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + enddo + +!---------------------------------------------- +! ...Quantum yields JPL06 +! ...This recommendation is only for >300nm +! However, it is used at all wavelengths +!---------------------------------------------- + qy1 = 0.85 + qy2 = 0.15 + do iw = 1, nw-1 + do iz = 1, nz + sq(j,iz,iw) = qy1 * ytd(iz,iw) + sq(j+1,iz,iw) = qy2 * ytd(iz,iw) + enddo + enddo + + j = j+1 + jlabel(j) = 'BrONO2 + hv -> BrO + NO2' + + end subroutine XSQY_BRONO2 diff --git a/test/unit/tuv_doug/JCALC/XSQY_HO2NO2.f b/test/unit/tuv_doug/JCALC/XSQY_HO2NO2.f new file mode 100644 index 00000000..5865bf03 --- /dev/null +++ b/test/unit/tuv_doug/JCALC/XSQY_HO2NO2.f @@ -0,0 +1,210 @@ + subroutine XSQY_HO2NO2(nw,wl,wc,nz,tlev,airlev,j,sq,jlabel,pn) +!-----------------------------------------------------------------------------! +! purpose: ! +! provide product of (cross section) x (quantum yield) for hno4 photolysis ! +! 1) HO2NO2 + hv -> HO2 + NO2 ! +! 2) HO2NO2 + hv -> OH + NO3 ! +! cross sections and QY from JPL06 ! +!-----------------------------------------------------------------------------! +! parameters: ! +! nw - integer, number of specified intervals + 1 in working (i) ! +! wavelength grid ! +! wl - real, vector of lower limits of wavelength intervals in (i) ! +! working wavelength grid ! +! wc - real, vector of center points of wavelength intervals in (i) ! +! working wavelength grid ! +! nz - integer, number of altitude levels in working altitude grid (i) ! +! tlev - real, temperature (k) at each specified altitude level (i) ! +! airlev - real, air density (molec/cc) at each altitude level (i) ! +! j - integer, counter for number of weighting functions defined (io) ! +! sq - real, cross section x quantum yield (cm^2) for each (o) ! +! photolysis reaction defined, at each defined wavelength and ! +! at each defined altitude level ! +! jlabel - character*60, string identifier for each photolysis reaction (o) ! +! defined ! +!-----------------------------------------------------------------------------! +! edit history: ! +! 05/98 original, adapted from former jspec1 subroutine ! +! 06/01 modified by doug kinnison ! +! 01/08 modified by Doug Kinnison ! +!-----------------------------------------------------------------------------! + implicit none + include 'params' + +!-----------------------------------------------------------------------------! +! ... input ! +!-----------------------------------------------------------------------------! + real, intent(in) :: wl(kw) + real, intent(in) :: wc(kw) + real, intent(in) :: tlev(kz) + real, intent(in) :: airlev(kz) + + integer, intent(in) :: nz + integer, intent(in) :: nw + + character*80, intent(in) :: pn + character*60, intent(out) :: jlabel(kj) + real, intent(out) :: sq(kj,kz,kw) + +!-----------------------------------------------------------------------------! +! ... input/output ! +!-----------------------------------------------------------------------------! + integer, intent(inout) :: j + +!-----------------------------------------------------------------------------! +! ... local ! +!-----------------------------------------------------------------------------! + integer kdata + parameter(kdata=100) + integer i, iw, iz, n, n1, idum, ierr, icnt + real x1 (kdata), x2(kdata), wcb(kdata) + real y1 (kdata), aa(kdata), bb (kdata) + real ytmp (nz,kdata), ycomb(nz,kdata) + real ytd (nz,kw), yg(kw) + real Q(nz), tin(nz), t + +!---------------------------------------------- +! ... tin set to tlev +!---------------------------------------------- + tin(:) = tlev(:) + +!---------------------------------------------- +! ... jlabel(j) = 'HO2NO2 -> HO2 + NO2 +! jlabel(j) = 'HO2NO2 -> OH + NO3 +!---------------------------------------------- + j = j + 1 + jlabel(j) = 'HO2NO2 + hv -> OH + NO3' + +!---------------------------------------------- +! ...ho2no2 cross sections plus T-dep. +! (Burkholder et al., 2002.) +!---------------------------------------------- + open(kin,file=TRIM(pn)//'XS_HO2NO2_JPL06.txt',status='old') + +!... read lambda and cross sections + read(kin,*) idum, n + do i = 1, idum-2 + read(kin,*) + enddo + do i = 1, n + read(kin,*) x1(i), y1(i) + enddo + +!... read lambda and T-dep coeff. + read(kin,*) + read(kin,*) idum, n1 + do i = 1, n1 + read(kin,*) x2(i), aa(i), bb(i) + enddo + close(kin) + +!---------------------------------------------- +! ...Derive T-dep Burkholder et al., 2002.) +!---------------------------------------------- + do iz = 1, nz + do iw = 1, n1 + t = MAX(280.,MIN(tin(iz),350.)) + Q(iz) = 1 + exp(-988./(0.69*t)) + ytmp(iz,iw) = ( aa(iw)/Q(iz) + bb(iw)*(1-1/Q(iz)))*1e-20 + enddo + enddo +!---------------------------------------------- +! ... Check routine +! iz = 1 +! do iw = 1, n1 +! print*, iw, x2(iw), ytmp(iz,iw) +! enddo +! stop +!---------------------------------------------- +! ... Combine cross sections + do iz = 1, nz + icnt = 1 + +! ... < 280 nm +! ... x1(iw) goes from 190-350nm + do iw = 1, n + IF (x1(iw) .LT. 280.) THEN + ycomb(iz,icnt) = y1(iw) + wcb (icnt) = x1(iw) + icnt = icnt + 1 + ENDIF + enddo +! ... 280-350 nm + do iw = 1, n1 + ycomb(iz,icnt) = ytmp(iz,iw) + wcb (icnt) = x2 (iw) + icnt = icnt+1 + enddo + enddo + +!... Test No TD +! do iz = 1, nz +! icnt = 1 +! do iw = 1, n +! ycomb(iz,icnt) = y1(iw) +! wcb (icnt) = x1(iw) +! icnt = icnt + 1 +! enddo +! enddo +!---------------------------------------------- +! ... Check routine +! iz = 1 +! print*,"tin=", tin(iz) +! do iw = 1, icnt-1 +! print*, iw, wcb(iw), ycomb(iz,iw) +! enddo +! stop +!---------------------------------------------- +! ... Interpolate Combine cross sections + do iz = 1, nz + n = icnt-1 + y1 = ycomb(iz,:) + x1 = wcb + + call addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) + call addpnt(x1,y1,kdata,n, 0.,0.) + call addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) + call addpnt(x1,y1,kdata,n, 1.e+38,0.) + call inter2(nw,wl,yg,n,x1,y1,ierr) + ytd(iz,:) = yg(:) + + if (ierr .ne. 0) then + write(*,*) ierr, jlabel(j) + stop + endif + enddo + +!------------------------------------------------- +! iz = 1 +! do iw = 24, 80 +! print*, iw, wc(iw), ytd(iz,iw), tin(iz) +! enddo +! stop +!------------------------------------------------- + do iw = 1, nw - 1 + IF (wc(iw) .LT. 200.0) THEN + do iz = 1, nz + sq(j, iz,iw) = 0.30 * ytd(iz,iw) + sq(j+1,iz,iw) = 0.70 * ytd(iz,iw) + enddo + ENDIF + IF (wc(iw) .GE. 200.0) THEN + do iz = 1, nz + sq(j, iz,iw) = 0.20 * ytd(iz,iw) + sq(j+1,iz,iw) = 0.80 * ytd(iz,iw) + enddo + ENDIF + enddo + +!-------------------------------------------------- +! iz = 1 +! do iw = 24, 80 +! print*, wc(iw), sq(j,iz,iw), sq(j+1,iz,iw) +! print*, sq(j,iz,iw)+sq(j+1,iz,iw) +! enddo +! stop +!------------------------------------------------- + j = j + 1 + jlabel(j) = 'HO2NO2 + hv -> HO2 + NO2' + + end subroutine XSQY_HO2NO2 diff --git a/test/unit/tuv_doug/driver.F90 b/test/unit/tuv_doug/driver.F90 index 761a9730..6dbdb8c7 100644 --- a/test/unit/tuv_doug/driver.F90 +++ b/test/unit/tuv_doug/driver.F90 @@ -189,6 +189,12 @@ subroutine calculate( label, temperature, air_density, xsqy ) case( "HCFC142b + hv -> Cl" ) call XSQY_HCFC142b(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) xsqy(:,:) = l_xsqy(1,:nz,:nw) + case( "BrONO2 + hv -> Br + NO3" ) + call XSQY_BRONO2(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) + case( "HO2NO2 + hv -> OH + NO3" ) + call XSQY_HO2NO2(nw,wl,wc,nz,temperature,air_density,j,l_xsqy,all_labels,pn) + xsqy(:,:) = l_xsqy(1,:nz,:nw) case default call die( 946669022 ) end select