diff --git a/schemes/musica/micm/musica_ccpp_micm.F90 b/schemes/musica/micm/musica_ccpp_micm.F90 index 2419d40e..1ef50d5a 100644 --- a/schemes/musica/micm/musica_ccpp_micm.F90 +++ b/schemes/musica/micm/musica_ccpp_micm.F90 @@ -16,7 +16,7 @@ module musica_ccpp_micm contains - !> Register MICM constituent properties with the CCPP + !> Registers MICM constituent properties with the CCPP subroutine micm_register(solver_type, num_grid_cells, constituent_props, errmsg, errcode) use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t use musica_micm, only: Rosenbrock, RosenbrockStandardOrder @@ -76,7 +76,7 @@ subroutine micm_register(solver_type, num_grid_cells, constituent_props, errmsg, end subroutine micm_register - !> Intitialize MICM + !> Initializes MICM subroutine micm_init(errmsg, errcode) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -86,7 +86,7 @@ subroutine micm_init(errmsg, errcode) end subroutine micm_init - !> Solve chemistry at the current time step + !> Solves chemistry at the current time step subroutine micm_run(time_step, temperature, pressure, dry_air_density, & user_defined_rate_parameters, constituents, errmsg, errcode) use musica_micm, only: solver_stats_t @@ -124,7 +124,7 @@ subroutine micm_run(time_step, temperature, pressure, dry_air_density, & end subroutine micm_run - !> Finalize MICM + !> Finalizes MICM subroutine micm_final(errmsg, errcode) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index 7b76ee55..09364397 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -26,14 +26,17 @@ end subroutine musica_ccpp_register !> \section arg_table_musica_ccpp_init Argument Table !! \htmlinclude musica_ccpp_init.html subroutine musica_ccpp_init(vertical_layer_dimension, vertical_interface_dimension, & - errmsg, errcode) - integer, intent(in) :: vertical_layer_dimension ! (count) - integer, intent(in) :: vertical_interface_dimension ! (count) + photolysis_wavelength_grid_interfaces, errmsg, errcode) + use ccpp_kinds, only : kind_phys + + integer, intent(in) :: vertical_layer_dimension ! (count) + integer, intent(in) :: vertical_interface_dimension ! (count) + real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! m character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & - errmsg, errcode) + photolysis_wavelength_grid_interfaces, errmsg, errcode) if (errcode /= 0) return call micm_init(errmsg, errcode) if (errcode /= 0) return @@ -43,13 +46,14 @@ end subroutine musica_ccpp_init !> \section arg_table_musica_ccpp_run Argument Table !! \htmlinclude musica_ccpp_run.html !! - !! The standard name for the variable 'surface_tempearture' is - !! `blackbody_temperature_at_surface` because this is what we have as - !! the standard name for `cam_in%ts`, whcih represents the same quantity. + !! The standard name for the variable 'surface_temperature' is + !! 'blackbody_temperature_at_surface' because this is what we have as + !! the standard name for 'cam_in%ts', whcih represents the same quantity. subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props, & constituents, geopotential_height_wrt_surface_at_midpoint, & geopotential_height_wrt_surface_at_interface, surface_temperature, & - surface_geopotential, standard_gravitational_acceleration, errmsg, errcode) + surface_geopotential, surface_albedo, & + standard_gravitational_acceleration, errmsg, errcode) use musica_ccpp_micm_util, only: reshape_into_micm_arr, reshape_into_ccpp_arr use musica_ccpp_micm_util, only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t @@ -67,6 +71,7 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface) real(kind_phys), intent(in) :: surface_temperature(:) ! K real(kind_phys), intent(in) :: surface_geopotential(:) ! m2 s-2 + real(kind_phys), intent(in) :: surface_albedo ! unitless real(kind_phys), intent(in) :: standard_gravitational_acceleration ! m s-2 character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -92,8 +97,8 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co call tuvx_run(temperature, dry_air_density, & geopotential_height_wrt_surface_at_midpoint, & geopotential_height_wrt_surface_at_interface, & - surface_temperature, & - surface_geopotential, & + surface_temperature, surface_geopotential, & + surface_albedo, & standard_gravitational_acceleration, & photolysis_rate_constants, & errmsg, errcode) diff --git a/schemes/musica/musica_ccpp.meta b/schemes/musica/musica_ccpp.meta index e4a38460..30d19bad 100644 --- a/schemes/musica/musica_ccpp.meta +++ b/schemes/musica/musica_ccpp.meta @@ -19,6 +19,12 @@ type = integer dimensions = () intent = in +[ photolysis_wavelength_grid_interfaces ] + standard_name = photolysis_wavelength_grid_interfaces + units = m + type = real | kind = kind_phys + dimensions = (photolysis_wavelength_grid_interface_dimension) + intent = in [ errmsg ] standard_name = ccpp_error_message units = none @@ -95,6 +101,12 @@ units = m2 s-2 dimensions = (horizontal_loop_extent) intent = in +[ surface_albedo ] + standard_name = surface_albedo_due_to_UV_and_VIS_direct + type = real | kind = kind_phys + units = None + dimensions = () + intent = in [ standard_gravitational_acceleration ] standard_name = standard_gravitational_acceleration units = m s-2 diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index 9ac1c553..bcc7be09 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -1,6 +1,6 @@ module musica_ccpp_tuvx - ! Note: "tuvx_t" is included in an external pre-built tuvx library that the host + ! Note: "tuvx_t" is included in an external pre-built TUV-x library that the host ! model is responsible for linking to during compilation use musica_tuvx, only: tuvx_t, grid_t, profile_t use musica_ccpp_util, only: has_error_occurred @@ -12,23 +12,32 @@ module musica_ccpp_tuvx public :: tuvx_init, tuvx_run, tuvx_final - type(tuvx_t), pointer :: tuvx => null( ) - type(grid_t), pointer :: height_grid => null( ) - type(profile_t), pointer :: temperature_profile => null( ) + type(tuvx_t), pointer :: tuvx => null() + type(grid_t), pointer :: height_grid => null() + type(grid_t), pointer :: wavelength_grid => null() + type(profile_t), pointer :: temperature_profile => null() + type(profile_t), pointer :: surface_albedo_profile => null() contains - !> Intitialize TUV-x + !> Initializes TUV-x subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & - errmsg, errcode) + wavelength_grid_interfaces, errmsg, errcode) use musica_tuvx, only: grid_map_t, profile_map_t, radiator_map_t use musica_util, only: error_t - use musica_ccpp_tuvx_height_grid, only: create_height_grid, & - height_grid_label, height_grid_unit - use musica_ccpp_tuvx_temperature, only: create_temperature_profile, & - temperature_label, temperature_unit - integer, intent(in) :: vertical_layer_dimension ! (count) - integer, intent(in) :: vertical_interface_dimension ! (count) + use musica_ccpp_tuvx_util, only: tuvx_deallocate + use musica_ccpp_tuvx_height_grid, & + only: create_height_grid, height_grid_label, height_grid_unit + use musica_ccpp_tuvx_wavelength_grid, & + only: create_wavelength_grid, wavelength_grid_label, wavelength_grid_unit + use musica_ccpp_tuvx_temperature, & + only: create_temperature_profile, temperature_label, temperature_unit + use musica_ccpp_tuvx_surface_albedo, & + only: create_surface_albedo_profile, surface_albedo_label, surface_albedo_unit + + integer, intent(in) :: vertical_layer_dimension ! (count) + integer, intent(in) :: vertical_interface_dimension ! (count) + real(kind_phys), intent(in) :: wavelength_grid_interfaces(:) ! m character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -50,71 +59,79 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & call grids%add( height_grid, error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( grids ) - deallocate( height_grid ) - height_grid => null() + call tuvx_deallocate( grids, null(), null(), null(), height_grid, null(), & + null(), null() ) + return + end if + + wavelength_grid => create_wavelength_grid( wavelength_grid_interfaces, & + errmsg, errcode ) + if (errcode /= 0) then + call tuvx_deallocate( grids, null(), null(), null(), height_grid, null(), & + null(), null() ) + return + endif + + call grids%add( wavelength_grid, error ) + if (has_error_occurred( error, errmsg, errcode )) then + call tuvx_deallocate( grids, null(), null(), null(), height_grid, & + wavelength_grid, null(), null() ) return end if profiles => profile_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( grids ) - deallocate( height_grid ) - height_grid => null() + call tuvx_deallocate( grids, null(), null(), null(), height_grid, & + wavelength_grid, null(), null() ) return end if temperature_profile => create_temperature_profile( height_grid, errmsg, errcode ) if (errcode /= 0) then - deallocate( grids ) - deallocate( height_grid ) - height_grid => null() - deallocate( profiles ) + call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & + wavelength_grid, null(), null() ) return endif call profiles%add( temperature_profile, error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( grids ) - deallocate( height_grid ) - height_grid => null() - deallocate( profiles ) - deallocate( temperature_profile ) - temperature_profile => null() + call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & + wavelength_grid, temperature_profile, null() ) + return + end if + + surface_albedo_profile => create_surface_albedo_profile( wavelength_grid, & + errmsg, errcode ) + if (errcode /= 0) then + call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & + wavelength_grid, temperature_profile, null() ) + return + endif + + call profiles%add( surface_albedo_profile, error ) + if (has_error_occurred( error, errmsg, errcode )) then + call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & + wavelength_grid, temperature_profile, surface_albedo_profile ) return end if radiators => radiator_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( grids ) - deallocate( height_grid ) - height_grid => null() - deallocate( profiles ) - deallocate( temperature_profile ) - temperature_profile => null() + call tuvx_deallocate( grids, profiles, null(), null(), height_grid, & + wavelength_grid, temperature_profile, surface_albedo_profile ) return end if tuvx => tuvx_t( filename_of_tuvx_configuration, grids, profiles, & radiators, error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( grids ) - deallocate( height_grid ) - height_grid => null() - deallocate( profiles ) - deallocate( temperature_profile ) - temperature_profile => null() - deallocate( radiators ) + call tuvx_deallocate( grids, profiles, radiators, null(), height_grid, & + wavelength_grid, temperature_profile, surface_albedo_profile ) return end if - deallocate( grids ) - deallocate( height_grid ) - height_grid => null() - deallocate( profiles ) - deallocate( temperature_profile ) - temperature_profile => null() - deallocate( radiators ) + call tuvx_deallocate( grids, profiles, radiators, null(), height_grid, & + wavelength_grid, temperature_profile, surface_albedo_profile ) grids => tuvx%get_grids( error ) if (has_error_occurred( error, errmsg, errcode )) then @@ -125,30 +142,42 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & height_grid => grids%get( height_grid_label, height_grid_unit, error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( tuvx ) - tuvx => null() - deallocate( grids ) + call tuvx_deallocate( grids, null(), null(), tuvx, null(), null(), & + null(), null() ) return end if - deallocate( grids ) + wavelength_grid => grids%get( wavelength_grid_label, wavelength_grid_unit, & + error ) + if (has_error_occurred( error, errmsg, errcode )) then + call tuvx_deallocate( grids, null(), null(), tuvx, height_grid, null(), & + null(), null() ) + return + end if profiles => tuvx%get_profiles( error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( tuvx ) - tuvx => null() + call tuvx_deallocate( grids, null(), null(), tuvx, height_grid, & + wavelength_grid, null(), null() ) return end if temperature_profile => profiles%get( temperature_label, temperature_unit, error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( tuvx ) - tuvx => null() - deallocate( profiles ) + call tuvx_deallocate( grids, profiles, null(), tuvx, height_grid, & + wavelength_grid, null(), null() ) + return + end if + + surface_albedo_profile => profiles%get( surface_albedo_label, surface_albedo_unit, error ) + if (has_error_occurred( error, errmsg, errcode )) then + call tuvx_deallocate( grids, profiles, null(), tuvx, height_grid, & + wavelength_grid, temperature_profile, null() ) return end if - deallocate( profiles ) + call tuvx_deallocate( grids, profiles, null(), null(), null(), null(), & + null(), null() ) end subroutine tuvx_init @@ -157,11 +186,13 @@ subroutine tuvx_run(temperature, dry_air_density, & geopotential_height_wrt_surface_at_midpoint, & geopotential_height_wrt_surface_at_interface, & surface_temperature, surface_geopotential, & + surface_albedo, & standard_gravitational_acceleration, & photolysis_rate_constants, errmsg, errcode) - use musica_util, only: error_t - use musica_ccpp_tuvx_height_grid, only: set_height_grid_values, calculate_heights - use musica_ccpp_tuvx_temperature, only: set_temperature_values + use musica_util, only: error_t + use musica_ccpp_tuvx_height_grid, only: set_height_grid_values, calculate_heights + use musica_ccpp_tuvx_temperature, only: set_temperature_values + use musica_ccpp_tuvx_surface_albedo, only: set_surface_albedo_values real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer) real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer) @@ -169,6 +200,7 @@ subroutine tuvx_run(temperature, dry_air_density, & real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface) real(kind_phys), intent(in) :: surface_temperature(:) ! K real(kind_phys), intent(in) :: surface_geopotential(:) ! m2 s-2 + real(kind_phys), intent(in) :: surface_albedo ! unitless real(kind_phys), intent(in) :: standard_gravitational_acceleration ! m s-2 ! temporarily set to Chapman mechanism and 1 dimension ! until mapping between MICM and TUV-x is implemented @@ -199,12 +231,16 @@ subroutine tuvx_run(temperature, dry_air_density, & if (errcode /= 0) return end do + ! surface albedo with respect to direct UV/visible radiation + call set_surface_albedo_values( surface_albedo_profile, surface_albedo, errmsg, errcode ) + if (errcode /= 0) return + ! stand-in until actual photolysis rate constants are calculated photolysis_rate_constants(:) = 1.0e-6_kind_phys end subroutine tuvx_run - !> Finalize tuvx + !> Finalizes TUV-x subroutine tuvx_final(errmsg, errcode) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -217,11 +253,21 @@ subroutine tuvx_final(errmsg, errcode) height_grid => null() end if + if (associated( wavelength_grid )) then + deallocate( wavelength_grid ) + wavelength_grid => null() + end if + if (associated( temperature_profile )) then deallocate( temperature_profile ) temperature_profile => null() end if + if (associated( surface_albedo_profile )) then + deallocate( surface_albedo_profile ) + surface_albedo_profile => null() + end if + if (associated( tuvx )) then deallocate( tuvx ) tuvx => null() diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 new file mode 100644 index 00000000..751d18c6 --- /dev/null +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 @@ -0,0 +1,76 @@ +module musica_ccpp_tuvx_surface_albedo + implicit none + + private + public :: create_surface_albedo_profile, set_surface_albedo_values, & + surface_albedo_label, surface_albedo_unit + + !> Label for surface albedo in TUV-x + character(len=*), parameter :: surface_albedo_label = "surface albedo" + !> Unit for surface albedo in TUV-x + character(len=*), parameter :: surface_albedo_unit = "none" + !> Default value of number of wavelength bins + integer, parameter :: DEFAULT_NUM_WAVELENGTH_BINS = 0 + !> Number of wavelength bins + integer, protected :: num_wavelength_bins = DEFAULT_NUM_WAVELENGTH_BINS + +contains + + !> Creates a TUV-x surface albedo profile from the host-model wavelength grid + function create_surface_albedo_profile( wavelength_grid, errmsg, errcode ) & + result( profile ) + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_tuvx_profile, only: profile_t + use musica_util, only: error_t + + type(grid_t), intent(inout) :: wavelength_grid + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + type(profile_t), pointer :: profile + + ! local variables + type(error_t) :: error + + num_wavelength_bins = wavelength_grid%number_of_sections( error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + profile => profile_t( surface_albedo_label, surface_albedo_unit, & + wavelength_grid, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end function create_surface_albedo_profile + + !> Sets TUV-x surface albedo values + !! + !! CAM uses a single value for surface albedo at all wavelengths + subroutine set_surface_albedo_values( profile, host_surface_albedo, & + errmsg, errcode ) + use ccpp_kinds, only: kind_phys + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_profile, only: profile_t + use musica_util, only: error_t + + type(profile_t), intent(inout) :: profile + real(kind_phys), intent(in) :: host_surface_albedo ! unitless + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! local variables + type(error_t) :: error + real(kind_phys) :: surface_albedo_interfaces(num_wavelength_bins + 1) + + if (size(surface_albedo_interfaces) <= DEFAULT_NUM_WAVELENGTH_BINS + 1) then + errmsg = "[MUSICA Error] Invalid size of TUV-x wavelength interfaces." + errcode = 1 + return + end if + + surface_albedo_interfaces(:) = host_surface_albedo + + call profile%set_edge_values( surface_albedo_interfaces, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end subroutine set_surface_albedo_values + +end module musica_ccpp_tuvx_surface_albedo \ No newline at end of file diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 new file mode 100644 index 00000000..25ac1a28 --- /dev/null +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 @@ -0,0 +1,55 @@ +module musica_ccpp_tuvx_util + implicit none + + private + public :: tuvx_deallocate + +contains + + !> This is a helper subroutine created to deallocate objects associated with TUV-x + subroutine tuvx_deallocate(grids, profiles, radiators, tuvx, height_grid, & + wavelength_grid, temperature_profile, surface_albedo_profile) + use musica_tuvx, only: tuvx_t, grid_map_t, profile_map_t, radiator_map_t, & + grid_t, profile_t + + type(grid_map_t), pointer :: grids + type(profile_map_t), pointer :: profiles + type(radiator_map_t), pointer :: radiators + type(tuvx_t), pointer :: tuvx + type(grid_t), pointer :: height_grid + type(grid_t), pointer :: wavelength_grid + type(profile_t), pointer :: temperature_profile + type(profile_t), pointer :: surface_albedo_profile + + if (associated( grids )) deallocate( grids ) + if (associated( profiles )) deallocate( profiles ) + if (associated( radiators )) deallocate( radiators ) + + if (associated( tuvx )) then + deallocate( tuvx ) + tuvx => null() + end if + + if (associated( height_grid )) then + deallocate( height_grid ) + height_grid => null() + end if + + if (associated( wavelength_grid )) then + deallocate( wavelength_grid ) + wavelength_grid => null() + end if + + if (associated( temperature_profile )) then + deallocate( temperature_profile ) + temperature_profile => null() + end if + + if (associated( surface_albedo_profile )) then + deallocate( surface_albedo_profile ) + surface_albedo_profile => null() + end if + + end subroutine tuvx_deallocate + +end module musica_ccpp_tuvx_util \ No newline at end of file diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_wavelength_grid.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_wavelength_grid.F90 new file mode 100644 index 00000000..96528d5d --- /dev/null +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_wavelength_grid.F90 @@ -0,0 +1,59 @@ +module musica_ccpp_tuvx_wavelength_grid + + implicit none + + private + public :: create_wavelength_grid + + ! TUV-x Wavelegnth grid notes + ! + !----------------------------------------------------------------------- + ! The wavelength grid used with TUV-x is based on the grid used in the + ! CAM-Chem photolysis rate constant lookup tables. Slight modifications + ! were made to the grid in the Shumann-Runge and Lyman-alpha regions to + ! work with the expectations of the TUV-x code. + ! + ! The wavelength grid is defined by the host model. Any wavelength- + ! resolved quantities passed to TUV-x must be on this grid. + + !> Label for wavelength grid in TUV-x + character(len=*), parameter, public :: wavelength_grid_label = "wavelength" + !> Unit for wavelength grid in TUV-x + character(len=*), parameter, public :: wavelength_grid_unit = "nm" + +contains + + !> Creates a TUV-x wavelength grid + function create_wavelength_grid( wavelength_grid_interfaces, errmsg, errcode ) & + result( wavelength_grid ) + + use ccpp_kinds, only: kind_phys + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx_grid, only: grid_t + use musica_util, only: error_t + + real(kind_phys), intent(in) :: wavelength_grid_interfaces(:) ! m + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errcode + type(grid_t), pointer :: wavelength_grid + + ! local variables + real(kind_phys) :: interfaces( size( wavelength_grid_interfaces ) ) ! nm + reaL(kind_phys) :: midpoints( size( wavelength_grid_interfaces ) - 1 ) ! nm + type(error_t) :: error + + interfaces(:) = wavelength_grid_interfaces(:) * 1.0e9 + midpoints(:) = & + 0.5 * ( interfaces( 1: size( interfaces ) - 1 ) & + + interfaces( 2: size( interfaces ) ) ) + wavelength_grid => grid_t( wavelength_grid_label, wavelength_grid_unit, & + size( midpoints ), error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + call wavelength_grid%set_edges( interfaces, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + call wavelength_grid%set_midpoints( midpoints, error ) + if ( has_error_occurred( error, errmsg, errcode ) ) return + + end function create_wavelength_grid + +end module musica_ccpp_tuvx_wavelength_grid \ No newline at end of file diff --git a/test/musica/micm/test_micm_util.F90 b/test/musica/micm/test_micm_util.F90 index 5581d5d8..798f1408 100644 --- a/test/musica/micm/test_micm_util.F90 +++ b/test/musica/micm/test_micm_util.F90 @@ -19,7 +19,11 @@ subroutine test_reshape() integer, parameter :: NUM_SPECIES = 4 integer, parameter :: NUM_COLUMNS = 2 integer, parameter :: NUM_LAYERS = 2 - integer, parameter :: NUM_GRID_CELLS = 4 + integer, parameter :: NUM_GRID_CELLS = NUM_COLUMNS * NUM_LAYERS + real, parameter :: ABS_ERROR = 1e-7 + real(c_double) :: expected_conditions(NUM_GRID_CELLS) = [100.0, 200.0, 300.0, 400.0] + real(c_double) :: expected_constituents(NUM_GRID_CELLS*NUM_SPECIES) = & + [0.1, 0.2, 0.3, 0.4, 0.21, 0.22, 0.23, 0.24, 0.41, 0.42, 0.43, 0.44, 0.31, 0.32, 0.33, 0.34] real(kind_phys) :: temperature(NUM_COLUMNS,NUM_LAYERS) real(kind_phys) :: pressure(NUM_COLUMNS,NUM_LAYERS) real(kind_phys) :: dry_air_density(NUM_COLUMNS,NUM_LAYERS) @@ -28,12 +32,7 @@ subroutine test_reshape() real(c_double) :: micm_pressure(NUM_GRID_CELLS) real(c_double) :: micm_dry_air_density(NUM_GRID_CELLS) real(c_double) :: micm_constituents(NUM_GRID_CELLS*NUM_SPECIES) - - ! local variables - real(c_double), dimension(NUM_GRID_CELLS) :: arr_conditions - real(c_double), dimension(NUM_GRID_CELLS*NUM_SPECIES) :: arr_constituents - integer :: i_column, i_layer, i_elem, i_arr - real :: abs_error = 1e-7 + integer :: i_column, i_layer, i_elem, i_arr temperature(:,1) = (/ 100._kind_phys, 200._kind_phys /) temperature(:,2) = (/ 300._kind_phys, 400._kind_phys /) @@ -45,20 +44,18 @@ subroutine test_reshape() constituents(1,2,:) = (/ 0.41_kind_phys, 0.42_kind_phys, 0.43_kind_phys, 0.44_kind_phys /) constituents(2,1,:) = (/ 0.21_kind_phys, 0.22_kind_phys, 0.23_kind_phys, 0.24_kind_phys /) constituents(2,2,:) = (/ 0.31_kind_phys, 0.32_kind_phys, 0.33_kind_phys, 0.34_kind_phys /) - arr_conditions = (/ 100.0, 200.0, 300.0, 400.0 /) - arr_constituents = (/ 0.1, 0.2, 0.3, 0.4, 0.21, 0.22, 0.23, 0.24, 0.41, 0.42, 0.43, 0.44, 0.31, 0.32, 0.33, 0.34 /) call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, & micm_temperature, micm_pressure, micm_dry_air_density, micm_constituents) do i_elem = 1, NUM_GRID_CELLS - ASSERT(micm_temperature(i_elem) == arr_conditions(i_elem)) - ASSERT(micm_pressure(i_elem) == arr_conditions(i_elem)) - ASSERT(micm_dry_air_density(i_elem) == arr_conditions(i_elem)) + ASSERT(micm_temperature(i_elem) == expected_conditions(i_elem)) + ASSERT(micm_pressure(i_elem) == expected_conditions(i_elem)) + ASSERT(micm_dry_air_density(i_elem) == expected_conditions(i_elem)) end do do i_elem = 1, size(micm_constituents) - ASSERT_NEAR(micm_constituents(i_elem), arr_constituents(i_elem), abs_error) + ASSERT_NEAR(micm_constituents(i_elem), expected_constituents(i_elem), ABS_ERROR) end do call reshape_into_ccpp_arr(micm_constituents, constituents) @@ -67,7 +64,7 @@ subroutine test_reshape() do i_layer = 1, NUM_LAYERS do i_column = 1, NUM_COLUMNS do i_elem = 1, NUM_SPECIES - ASSERT_NEAR(constituents(i_column, i_layer, i_elem), arr_constituents(i_arr), abs_error) + ASSERT_NEAR(constituents(i_column, i_layer, i_elem), expected_constituents(i_arr), ABS_ERROR) i_arr = i_arr + 1 end do end do @@ -78,16 +75,16 @@ end subroutine test_reshape subroutine test_unit_conversion() use ccpp_kinds, only: kind_phys - integer, parameter :: NUM_COLUMNS = 2 - integer, parameter :: NUM_LAYERS = 2 - integer, parameter :: NUM_SPECIES = 4 - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: dry_air_density ! kg m-3 - real(kind_phys), dimension(NUM_SPECIES) :: molar_mass_arr - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: constituents - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: ccpp_constituents ! kg kg-1 - real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: micm_constituents ! mol m-3 - integer :: i_column, i_layer, i_elem - real :: abs_error = 1e-3 + integer, parameter :: NUM_COLUMNS = 2 + integer, parameter :: NUM_LAYERS = 2 + integer, parameter :: NUM_SPECIES = NUM_COLUMNS * NUM_LAYERS + real, parameter :: ABS_ERROR = 1e-3 + real(kind_phys) :: dry_air_density(NUM_COLUMNS,NUM_LAYERS) + real(kind_phys) :: molar_mass_arr(NUM_SPECIES) + real(kind_phys) :: constituents(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) + real(kind_phys) :: ccpp_constituents(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) ! kg kg-1 + real(kind_phys) :: micm_constituents(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) ! mol m-3 + integer :: i_column, i_layer, i_elem dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) @@ -112,7 +109,7 @@ subroutine test_unit_conversion() do i_column = 1, NUM_COLUMNS do i_layer = 1, NUM_LAYERS do i_elem = 1, NUM_SPECIES - ASSERT_NEAR(constituents(i_column, i_layer, i_elem), micm_constituents(i_column, i_layer, i_elem), abs_error) + ASSERT_NEAR(constituents(i_column, i_layer, i_elem), micm_constituents(i_column, i_layer, i_elem), ABS_ERROR) end do end do end do @@ -122,10 +119,11 @@ subroutine test_unit_conversion() do i_column = 1, NUM_COLUMNS do i_layer = 1, NUM_LAYERS do i_elem = 1, NUM_SPECIES - ASSERT_NEAR(constituents(i_column, i_layer, i_elem), ccpp_constituents(i_column, i_layer, i_elem), abs_error) + ASSERT_NEAR(constituents(i_column, i_layer, i_elem), ccpp_constituents(i_column, i_layer, i_elem), ABS_ERROR) end do end do end do + end subroutine test_unit_conversion end program test_micm_util \ No newline at end of file diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index ac37e30c..6b4d4154 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -16,44 +16,151 @@ subroutine test_musica_ccpp_api() use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t - implicit none - integer, parameter :: NUM_SPECIES = 4 integer, parameter :: NUM_COLUMNS = 2 integer, parameter :: NUM_LAYERS = 2 - integer :: solver_type + integer, parameter :: NUM_WAVELENGTH_BINS = 102 + integer :: num_grid_cells = NUM_COLUMNS * NUM_LAYERS + integer :: solver_type = Rosenbrock integer :: errcode character(len=512) :: errmsg - real(kind_phys) :: time_step ! s + real(kind_phys) :: time_step = 60._kind_phys ! s + real(kind_phys), dimension(NUM_WAVELENGTH_BINS+1) :: photolysis_wavelength_grid_interfaces ! m real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: geopotential_height_wrt_surface_at_midpoint ! m real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m real(kind_phys), dimension(NUM_COLUMNS) :: surface_temperature ! K real(kind_phys), dimension(NUM_COLUMNS) :: surface_geopotential ! m2 s-2 + real(kind_phys) :: surface_albedo ! unitless real(kind_phys) :: standard_gravitational_acceleration ! s2 m-1 real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: temperature ! K real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: pressure ! Pa real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: dry_air_density ! kg m-3 real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: constituents ! kg kg-1 type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:) - - ! local variables - type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:) - type(ccpp_constituent_properties_t), pointer :: const_prop - real(kind_phys) :: molar_mass - character(len=512) :: species_name, units - logical :: tmp_bool, is_advected - integer :: num_grid_cells - integer :: i - - solver_type = Rosenbrock - num_grid_cells = NUM_COLUMNS * NUM_LAYERS - time_step = 60._kind_phys + type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:) + type(ccpp_constituent_properties_t), pointer :: const_prop + real(kind_phys) :: molar_mass + character(len=512) :: species_name, units + logical :: tmp_bool, is_advected + integer :: i + + ! These are the values that will be used in CAM-SIMA and correspond to the wavelength + ! bins used in the CAM-Chem photolysis rate constant lookup table. + ! + ! We're using the actual values here because several of the TS1/TSMLT photolysis + ! rate constant configurations are sensitive to the wavelength grid. + photolysis_wavelength_grid_interfaces = (/ & + 120.0e-9_kind_phys, & + 121.4e-9_kind_phys, & + 121.9e-9_kind_phys, & + 123.5e-9_kind_phys, & + 124.3e-9_kind_phys, & + 125.5e-9_kind_phys, & + 126.3e-9_kind_phys, & + 127.1e-9_kind_phys, & + 130.1e-9_kind_phys, & + 131.1e-9_kind_phys, & + 135.0e-9_kind_phys, & + 140.0e-9_kind_phys, & + 145.0e-9_kind_phys, & + 150.0e-9_kind_phys, & + 155.0e-9_kind_phys, & + 160.0e-9_kind_phys, & + 165.0e-9_kind_phys, & + 168.0e-9_kind_phys, & + 171.0e-9_kind_phys, & + 173.0e-9_kind_phys, & + 174.4e-9_kind_phys, & + 175.4e-9_kind_phys, & + 177.0e-9_kind_phys, & + 178.6e-9_kind_phys, & + 180.2e-9_kind_phys, & + 181.8e-9_kind_phys, & + 183.5e-9_kind_phys, & + 185.2e-9_kind_phys, & + 186.9e-9_kind_phys, & + 188.7e-9_kind_phys, & + 190.5e-9_kind_phys, & + 192.3e-9_kind_phys, & + 194.2e-9_kind_phys, & + 196.1e-9_kind_phys, & + 198.0e-9_kind_phys, & + 200.0e-9_kind_phys, & + 202.0e-9_kind_phys, & + 204.1e-9_kind_phys, & + 206.2e-9_kind_phys, & + 208.0e-9_kind_phys, & + 211.0e-9_kind_phys, & + 214.0e-9_kind_phys, & + 217.0e-9_kind_phys, & + 220.0e-9_kind_phys, & + 223.0e-9_kind_phys, & + 226.0e-9_kind_phys, & + 229.0e-9_kind_phys, & + 232.0e-9_kind_phys, & + 235.0e-9_kind_phys, & + 238.0e-9_kind_phys, & + 241.0e-9_kind_phys, & + 244.0e-9_kind_phys, & + 247.0e-9_kind_phys, & + 250.0e-9_kind_phys, & + 253.0e-9_kind_phys, & + 256.0e-9_kind_phys, & + 259.0e-9_kind_phys, & + 263.0e-9_kind_phys, & + 267.0e-9_kind_phys, & + 271.0e-9_kind_phys, & + 275.0e-9_kind_phys, & + 279.0e-9_kind_phys, & + 283.0e-9_kind_phys, & + 287.0e-9_kind_phys, & + 291.0e-9_kind_phys, & + 295.0e-9_kind_phys, & + 298.5e-9_kind_phys, & + 302.5e-9_kind_phys, & + 305.5e-9_kind_phys, & + 308.5e-9_kind_phys, & + 311.5e-9_kind_phys, & + 314.5e-9_kind_phys, & + 317.5e-9_kind_phys, & + 322.5e-9_kind_phys, & + 327.5e-9_kind_phys, & + 332.5e-9_kind_phys, & + 337.5e-9_kind_phys, & + 342.5e-9_kind_phys, & + 347.5e-9_kind_phys, & + 350.0e-9_kind_phys, & + 355.0e-9_kind_phys, & + 360.0e-9_kind_phys, & + 365.0e-9_kind_phys, & + 370.0e-9_kind_phys, & + 375.0e-9_kind_phys, & + 380.0e-9_kind_phys, & + 385.0e-9_kind_phys, & + 390.0e-9_kind_phys, & + 395.0e-9_kind_phys, & + 400.0e-9_kind_phys, & + 405.0e-9_kind_phys, & + 410.0e-9_kind_phys, & + 415.0e-9_kind_phys, & + 420.0e-9_kind_phys, & + 430.0e-9_kind_phys, & + 440.0e-9_kind_phys, & + 450.0e-9_kind_phys, & + 500.0e-9_kind_phys, & + 550.0e-9_kind_phys, & + 600.0e-9_kind_phys, & + 650.0e-9_kind_phys, & + 700.0e-9_kind_phys, & + 750.0e-9_kind_phys & + /) geopotential_height_wrt_surface_at_midpoint(1,:) = (/ 2000.0_kind_phys, 500.0_kind_phys /) geopotential_height_wrt_surface_at_midpoint(2,:) = (/ 2000.0_kind_phys, -500.0_kind_phys /) geopotential_height_wrt_surface_at_interface(1,:) = (/ 3000.0_kind_phys, 1000.0_kind_phys, 0.0_kind_phys /) geopotential_height_wrt_surface_at_interface(2,:) = (/ 3000.0_kind_phys, 500.0_kind_phys, -1500.0_kind_phys /) surface_temperature = (/ 300.0_kind_phys, 300.0_kind_phys /) surface_geopotential = (/ 100.0_kind_phys, 200.0_kind_phys /) + surface_albedo = 0.10_kind_phys standard_gravitational_acceleration = 10.0_kind_phys temperature(:,1) = (/ 100._kind_phys, 200._kind_phys /) temperature(:,2) = (/ 300._kind_phys, 400._kind_phys /) @@ -98,7 +205,8 @@ subroutine test_musica_ccpp_api() call constituent_props_ptr(i)%set(const_prop, errcode, errmsg) end do - call musica_ccpp_init(NUM_LAYERS, NUM_LAYERS+1, errmsg, errcode) + call musica_ccpp_init(NUM_LAYERS, NUM_LAYERS+1, photolysis_wavelength_grid_interfaces, & + errmsg, errcode) if (errcode /= 0) then write(*,*) trim(errmsg) stop 3 @@ -113,10 +221,11 @@ subroutine test_musica_ccpp_api() write(*,*) "[MUSICA INFO] Initial Concentrations" write(*,fmt="(4(3x,e13.6))") constituents - call musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props_ptr, & - constituents, geopotential_height_wrt_surface_at_midpoint, & - geopotential_height_wrt_surface_at_interface, surface_temperature, & - surface_geopotential, standard_gravitational_acceleration, errmsg, errcode) + call musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props_ptr, & + constituents, geopotential_height_wrt_surface_at_midpoint, & + geopotential_height_wrt_surface_at_interface, surface_temperature, & + surface_geopotential, surface_albedo, standard_gravitational_acceleration, & + errmsg, errcode) if (errcode /= 0) then write(*,*) trim(errmsg) stop 3 @@ -132,6 +241,8 @@ subroutine test_musica_ccpp_api() stop 3 endif + deallocate(constituent_props_ptr) + end subroutine test_musica_ccpp_api end program run_test_musica_ccpp \ No newline at end of file diff --git a/test/musica/tuvx/CMakeLists.txt b/test/musica/tuvx/CMakeLists.txt index b863be59..2c8623fd 100644 --- a/test/musica/tuvx/CMakeLists.txt +++ b/test/musica/tuvx/CMakeLists.txt @@ -26,6 +26,34 @@ add_test( add_memory_check_test(test_tuvx_height_grid $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) +# Wavelength grid +add_executable(test_tuvx_wavelength_grid test_tuvx_wavelength_grid.F90) + +target_sources(test_tuvx_wavelength_grid + PUBLIC + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_wavelength_grid.F90 + ${MUSICA_SRC_PATH}/musica_ccpp_util.F90 + ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 +) + +target_link_libraries(test_tuvx_wavelength_grid + PRIVATE + musica::musica-fortran +) + +set_target_properties(test_tuvx_wavelength_grid + PROPERTIES + LINKER_LANGUAGE Fortran +) + +add_test( + NAME test_tuvx_wavelength_grid + COMMAND $ + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} +) + +add_memory_check_test(test_tuvx_wavelength_grid $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) + # Temperature add_executable(test_tuvx_temperature test_tuvx_temperature.F90) @@ -53,4 +81,33 @@ add_test( WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} ) -add_memory_check_test(test_tuvx_temperature $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) \ No newline at end of file +add_memory_check_test(test_tuvx_temperature $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) + +# Surface albedo +add_executable(test_tuvx_surface_albedo test_tuvx_surface_albedo.F90) + +target_sources(test_tuvx_surface_albedo + PUBLIC + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_wavelength_grid.F90 + ${MUSICA_SRC_PATH}/tuvx/musica_ccpp_tuvx_surface_albedo.F90 + ${MUSICA_SRC_PATH}/musica_ccpp_util.F90 + ${CCPP_TEST_SRC_PATH}/ccpp_kinds.F90 +) + +target_link_libraries(test_tuvx_surface_albedo + PRIVATE + musica::musica-fortran +) + +set_target_properties(test_tuvx_surface_albedo + PROPERTIES + LINKER_LANGUAGE Fortran +) + +add_test( + NAME test_tuvx_surface_albedo + COMMAND $ + WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} +) + +add_memory_check_test(test_tuvx_surface_albedo $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) \ No newline at end of file diff --git a/test/musica/tuvx/configs/ts1_tsmlt.json b/test/musica/tuvx/configs/ts1_tsmlt.json index 378a38fe..b2a2f986 100644 --- a/test/musica/tuvx/configs/ts1_tsmlt.json +++ b/test/musica/tuvx/configs/ts1_tsmlt.json @@ -4,12 +4,6 @@ "cross section parameters file": "data/cross_sections/O2_parameters.txt" }, "grids": [ - { - "name": "wavelength", - "type": "from csv file", - "units": "nm", - "file path": "data/grids/wavelength/cam.csv" - }, { "name": "time", "type": "from config file", @@ -54,16 +48,6 @@ "month": 3, "day": 21 }, - { - "name": "surface albedo", - "type": "from config file", - "units": "none", - "uniform value": 0.10, - "grid": { - "name": "wavelength", - "units": "nm" - } - }, { "name": "extraterrestrial flux", "enable diagnostics" : true, diff --git a/test/musica/tuvx/test_tuvx_height_grid.F90 b/test/musica/tuvx/test_tuvx_height_grid.F90 index c0f5c6ed..6da26960 100644 --- a/test/musica/tuvx/test_tuvx_height_grid.F90 +++ b/test/musica/tuvx/test_tuvx_height_grid.F90 @@ -7,6 +7,8 @@ program test_tuvx_height_grid #define ASSERT(x) if (.not.(x)) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: x"; stop 1; endif #define ASSERT_NEAR( a, b, abs_error ) if( (abs(a - b) >= abs_error) .and. (abs(a - b) /= 0.0) ) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: a, b"; stop 1; endif + real, parameter :: ABS_ERROR = 1e-5 + call test_create_height_grid() call test_calculate_height_grid_values() @@ -17,20 +19,19 @@ subroutine test_create_height_grid() use musica_tuvx_grid, only: grid_t use ccpp_kinds, only: kind_phys - integer, parameter :: NUM_HOST_MIDPOINTS = 2 - integer, parameter :: NUM_HOST_INTERFACES = 3 - real(kind_phys), target :: host_midpoints(NUM_HOST_MIDPOINTS) = [200.8_kind_phys, 100.6_kind_phys] - real(kind_phys), target :: host_interfaces(NUM_HOST_INTERFACES) = [250.3_kind_phys, 150.2_kind_phys, 50.1_kind_phys] - type(grid_t), pointer :: height_grid - character(len=512) :: errmsg - integer :: errcode - real(kind_phys) :: abs_error = 1e-5 - integer :: i - - ! local variables - real(kind_phys), dimension(NUM_HOST_MIDPOINTS+1) :: midpoints - real(kind_phys), dimension(NUM_HOST_INTERFACES+1) :: interfaces - type(error_t) :: error + integer, parameter :: NUM_HOST_MIDPOINTS = 2 + integer, parameter :: NUM_HOST_INTERFACES = 3 + real(kind_phys) :: host_midpoints(NUM_HOST_MIDPOINTS) = [200.8_kind_phys, 100.6_kind_phys] + real(kind_phys) :: host_interfaces(NUM_HOST_INTERFACES) = [250.3_kind_phys, 150.2_kind_phys, 50.1_kind_phys] + real(kind_phys) :: expected_midpoints(NUM_HOST_MIDPOINTS+1) = [(100.6 + 50.1) * 0.5, 150.2, (250.3 + 200.8) * 0.5] + real(kind_phys) :: expected_interfaces(NUM_HOST_INTERFACES+1) = [50.1_kind_phys, 100.6_kind_phys, 200.8_kind_phys, 250.3_kind_phys] + real(kind_phys) :: midpoints(NUM_HOST_MIDPOINTS+1) + real(kind_phys) :: interfaces(NUM_HOST_INTERFACES+1) + type(grid_t), pointer :: height_grid => null() + type(error_t) :: error + character(len=512) :: errmsg + integer :: errcode + integer :: i height_grid => create_height_grid(-1, 0, errmsg, errcode) ASSERT(errcode == 1) @@ -54,50 +55,47 @@ subroutine test_create_height_grid() call height_grid%get_midpoints(midpoints, error) ASSERT(error%is_success()) + do i = 1, size(midpoints) + ASSERT_NEAR(midpoints(i), expected_midpoints(i), ABS_ERROR) + end do call height_grid%get_edges(interfaces, error) ASSERT(error%is_success()) - - ASSERT_NEAR(midpoints(1), (100.6 + 50.1) * 0.5, abs_error) - ASSERT_NEAR(midpoints(2), 150.2, abs_error) - ASSERT_NEAR(midpoints(3), (250.3 + 200.8) * 0.5, abs_error) - ASSERT_NEAR(interfaces(1), 50.1, abs_error) - ASSERT_NEAR(interfaces(2), 100.6, abs_error) - ASSERT_NEAR(interfaces(3), 200.8, abs_error) - ASSERT_NEAR(interfaces(4), 250.3, abs_error) + do i = 1, size(interfaces) + ASSERT_NEAR(interfaces(i), expected_interfaces(i), ABS_ERROR) + end do deallocate( height_grid ) end subroutine test_create_height_grid subroutine test_calculate_height_grid_values() - use ccpp_kinds, only: kind_phys - integer, parameter :: NUM_LAYERS = 2 - real(kind_phys), dimension(NUM_LAYERS) :: geopotential_height_wrt_surface_at_midpoint ! m - real(kind_phys), dimension(NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m - real(kind_phys) :: surface_geopotential ! m2 s-2 - real(kind_phys) :: reciprocal_of_gravitational_acceleration ! s2 m-1 - real(kind_phys), dimension(NUM_LAYERS) :: height_midpoints ! km - real(kind_phys), dimension(NUM_LAYERS+1) :: height_interfaces ! km - - geopotential_height_wrt_surface_at_midpoint(:) = (/ 2000.0_kind_phys, 500.0_kind_phys /) - geopotential_height_wrt_surface_at_interface(:) = (/ 3000.0_kind_phys, 1000.0_kind_phys, 0.0_kind_phys /) - surface_geopotential = 100.0_kind_phys - reciprocal_of_gravitational_acceleration = 10.0_kind_phys + integer, parameter :: NUM_LAYERS = 2 + real(kind_phys) :: geopotential_height_wrt_surface_at_midpoint(NUM_LAYERS) = [2000.0_kind_phys, 500.0_kind_phys] + real(kind_phys) :: geopotential_height_wrt_surface_at_interface(NUM_LAYERS+1) = [3000.0_kind_phys, 1000.0_kind_phys, 0.0_kind_phys] + real(kind_phys) :: surface_geopotential = 100.0_kind_phys + real(kind_phys) :: reciprocal_of_gravitational_acceleration = 10.0_kind_phys + real(kind_phys) :: expected_height_midpoints(NUM_LAYERS) = [3.0_kind_phys, 1.5_kind_phys] + real(kind_phys) :: expected_height_interfaces(NUM_LAYERS+1) = [4.0_kind_phys, 2.0_kind_phys, 1.0_kind_phys] + real(kind_phys) :: height_midpoints(NUM_LAYERS) + real(kind_phys) :: height_interfaces(NUM_LAYERS+1) + integer :: i call calculate_heights(geopotential_height_wrt_surface_at_midpoint, & geopotential_height_wrt_surface_at_interface, & surface_geopotential, reciprocal_of_gravitational_acceleration, & height_midpoints, height_interfaces) - ASSERT_NEAR(height_midpoints(1), 3.0, 1e-5) - ASSERT_NEAR(height_midpoints(2), 1.5, 1e-5) - ASSERT_NEAR(height_interfaces(1), 4.0, 1e-5) - ASSERT_NEAR(height_interfaces(2), 2.0, 1e-5) - ASSERT_NEAR(height_interfaces(3), 1.0, 1e-5) + do i = 1, size(height_midpoints) + ASSERT_NEAR(height_midpoints(i), expected_height_midpoints(i), ABS_ERROR) + end do + + do i = 1, size(height_interfaces) + ASSERT_NEAR(height_interfaces(i), expected_height_interfaces(i), ABS_ERROR) + end do end subroutine test_calculate_height_grid_values -end program test_tuvx_height_grid +end program test_tuvx_height_grid \ No newline at end of file diff --git a/test/musica/tuvx/test_tuvx_surface_albedo.F90 b/test/musica/tuvx/test_tuvx_surface_albedo.F90 new file mode 100644 index 00000000..c47ebc63 --- /dev/null +++ b/test/musica/tuvx/test_tuvx_surface_albedo.F90 @@ -0,0 +1,57 @@ +program test_tuvx_surface_albedo + + use musica_ccpp_tuvx_surface_albedo + + implicit none + +#define ASSERT(x) if (.not.(x)) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: x"; stop 1; endif +#define ASSERT_NEAR( a, b, abs_error ) if( (abs(a - b) >= abs_error) .and. (abs(a - b) /= 0.0) ) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: a, b"; stop 1; endif + + call test_update_surface_albedo() + +contains + + subroutine test_update_surface_albedo() + use musica_ccpp_tuvx_wavelength_grid, only: create_wavelength_grid + use musica_util, only: error_t + use musica_tuvx_grid, only: grid_t + use musica_tuvx_profile, only: profile_t + use ccpp_kinds, only: kind_phys + + integer, parameter :: NUM_WAVELENGTH_BINS = 4 + real, parameter :: ABS_ERROR = 1e-6 + real(kind_phys) :: wavelength_grid_interfaces(NUM_WAVELENGTH_BINS + 1) = & + [200.0e-9_kind_phys, 210.0e-9_kind_phys, 240.0e-9_kind_phys, 300.0e-9_kind_phys, 400.0e-9_kind_phys] + real(kind_phys) :: host_surface_albedo = 0.3_kind_phys + real(kind_phys) :: expected_surface_albedo_interfaces(NUM_WAVELENGTH_BINS + 1) = 0.3_kind_phys + real(kind_phys) :: surface_albedo_interfaces(NUM_WAVELENGTH_BINS + 1) + type(grid_t), pointer :: wavelength_grid + type(profile_t), pointer :: profile + type(error_t) :: error + character(len=512) :: errmsg + integer :: errcode + integer :: i + + wavelength_grid => create_wavelength_grid(wavelength_grid_interfaces, errmsg, errcode) + ASSERT(errcode == 0) + ASSERT(associated(wavelength_grid)) + + profile => create_surface_albedo_profile( wavelength_grid, errmsg, errcode ) + ASSERT(errcode == 0) + ASSERT(associated(profile)) + + call set_surface_albedo_values( profile, host_surface_albedo, errmsg, errcode ) + ASSERT(errcode == 0) + + call profile%get_edge_values( surface_albedo_interfaces, error) + ASSERT(error%is_success()) + do i = 1, size(surface_albedo_interfaces) + ASSERT_NEAR(surface_albedo_interfaces(i), expected_surface_albedo_interfaces(i), ABS_ERROR) + end do + + deallocate( profile ) + deallocate( wavelength_grid ) + + end subroutine test_update_surface_albedo + +end program test_tuvx_surface_albedo \ No newline at end of file diff --git a/test/musica/tuvx/test_tuvx_temperature.F90 b/test/musica/tuvx/test_tuvx_temperature.F90 index 2f78eb56..cecde8b4 100644 --- a/test/musica/tuvx/test_tuvx_temperature.F90 +++ b/test/musica/tuvx/test_tuvx_temperature.F90 @@ -18,20 +18,21 @@ subroutine test_update_temperature() use musica_tuvx_profile, only: profile_t use ccpp_kinds, only: kind_phys - integer, parameter :: NUM_HOST_MIDPOINTS = 5 - integer, parameter :: NUM_HOST_INTERFACES = 6 - real(kind_phys) :: host_midpoint_temperature(NUM_HOST_MIDPOINTS) & - = [800.8_kind_phys, 700.7_kind_phys, 600.6_kind_phys, 500.5_kind_phys, 400.4_kind_phys] - real(kind_phys) :: host_surface_temperature = 300.3_kind_phys - real(kind_phys) :: interface_temperatures(NUM_HOST_MIDPOINTS+2) - real(kind_phys) :: expected_interface_temperatures(NUM_HOST_MIDPOINTS+2) & - = [300.3_kind_phys, 400.4_kind_phys, 500.5_kind_phys, 600.6_kind_phys, 700.7_kind_phys, 800.8_kind_phys, 800.8_kind_phys] - type(grid_t), pointer :: height_grid - type(profile_t), pointer :: profile - character(len=512) :: errmsg - integer :: errcode - type(error_t) :: error - real(kind_phys) :: abs_error = 1e-4 + integer, parameter :: NUM_HOST_MIDPOINTS = 5 + integer, parameter :: NUM_HOST_INTERFACES = 6 + real, parameter :: ABS_ERROR = 1e-4 + real(kind_phys) :: host_midpoint_temperature(NUM_HOST_MIDPOINTS) = & + [800.8_kind_phys, 700.7_kind_phys, 600.6_kind_phys, 500.5_kind_phys, 400.4_kind_phys] + real(kind_phys) :: host_surface_temperature = 300.3_kind_phys + real(kind_phys) :: expected_temperature_interfaces(NUM_HOST_MIDPOINTS+2) = & + [300.3_kind_phys, 400.4_kind_phys, 500.5_kind_phys, 600.6_kind_phys, 700.7_kind_phys, 800.8_kind_phys, 800.8_kind_phys] + real(kind_phys) :: temperature_interfaces(NUM_HOST_MIDPOINTS+2) + type(grid_t), pointer :: height_grid + type(profile_t), pointer :: profile + character(len=512) :: errmsg + integer :: errcode + type(error_t) :: error + integer :: i height_grid => create_height_grid(NUM_HOST_MIDPOINTS, NUM_HOST_INTERFACES, errmsg, errcode) ASSERT(errcode == 0) @@ -45,16 +46,11 @@ subroutine test_update_temperature() host_surface_temperature, errmsg, errcode ) ASSERT(errcode == 0) - call profile%get_edge_values( interface_temperatures, error) + call profile%get_edge_values( temperature_interfaces, error) ASSERT(error%is_success()) - - ASSERT_NEAR(interface_temperatures(1), expected_interface_temperatures(1), abs_error) - ASSERT_NEAR(interface_temperatures(2), expected_interface_temperatures(2), abs_error) - ASSERT_NEAR(interface_temperatures(3), expected_interface_temperatures(3), abs_error) - ASSERT_NEAR(interface_temperatures(4), expected_interface_temperatures(4), abs_error) - ASSERT_NEAR(interface_temperatures(5), expected_interface_temperatures(5), abs_error) - ASSERT_NEAR(interface_temperatures(6), expected_interface_temperatures(6), abs_error) - ASSERT_NEAR(interface_temperatures(7), expected_interface_temperatures(7), abs_error) + do i = 1, size(temperature_interfaces) + ASSERT_NEAR(temperature_interfaces(i), expected_temperature_interfaces(i), ABS_ERROR) + end do deallocate( profile ) deallocate( height_grid ) diff --git a/test/musica/tuvx/test_tuvx_wavelength_grid.F90 b/test/musica/tuvx/test_tuvx_wavelength_grid.F90 new file mode 100644 index 00000000..e90d49c8 --- /dev/null +++ b/test/musica/tuvx/test_tuvx_wavelength_grid.F90 @@ -0,0 +1,53 @@ +program test_tuvx_wavelength_grid + + use musica_ccpp_tuvx_wavelength_grid + + implicit none + +#define ASSERT(x) if (.not.(x)) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: x"; stop 1; endif +#define ASSERT_NEAR( a, b, abs_error ) if( (abs(a - b) >= abs_error) .and. (abs(a - b) /= 0.0) ) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: a, b"; stop 1; endif + + call test_create_wavelength_grid() + +contains + + subroutine test_create_wavelength_grid() + use musica_util, only: error_t + use musica_tuvx_grid, only: grid_t + use ccpp_kinds, only: kind_phys + + integer, parameter :: NUM_WAVELENGTH_GRID_MIDPOINTS = 2 + integer, parameter :: NUM_WAVELENGTH_GRID_INTERFACES = 3 + real, parameter :: ABS_ERROR = 1e-5 + real(kind_phys) :: host_interfaces(NUM_WAVELENGTH_GRID_INTERFACES) = [180.0e-9_kind_phys, 200.0e-9_kind_phys, 240.0e-9_kind_phys] + real(kind_phys) :: expected_interfaces(NUM_WAVELENGTH_GRID_INTERFACES) = [180.0_kind_phys, 200.0_kind_phys, 240.0_kind_phys] + real(kind_phys) :: expected_midpoints(NUM_WAVELENGTH_GRID_MIDPOINTS) = [190.0_kind_phys, 220.0_kind_phys] + real(kind_phys) :: interfaces(NUM_WAVELENGTH_GRID_INTERFACES) + real(kind_phys) :: midpoints(NUM_WAVELENGTH_GRID_MIDPOINTS) + type(grid_t), pointer :: wavelength_grid => null() + character(len=512) :: errmsg + integer :: errcode + type(error_t) :: error + integer :: i + + wavelength_grid => create_wavelength_grid(host_interfaces, errmsg, errcode) + ASSERT(errcode == 0) + ASSERT(associated(wavelength_grid)) + + call wavelength_grid%get_edges(interfaces, error) + ASSERT(error%is_success()) + do i = 1, NUM_WAVELENGTH_GRID_INTERFACES + ASSERT_NEAR(interfaces(i), expected_interfaces(i), ABS_ERROR) + end do + + call wavelength_grid%get_midpoints(midpoints, error) + ASSERT(error%is_success()) + do i = 1, NUM_WAVELENGTH_GRID_MIDPOINTS + ASSERT_NEAR(midpoints(i), expected_midpoints(i), ABS_ERROR) + end do + + deallocate(wavelength_grid) + + end subroutine test_create_wavelength_grid + +end program test_tuvx_wavelength_grid \ No newline at end of file