From 862ac09096bda64fd3617364bc326203daa437b7 Mon Sep 17 00:00:00 2001 From: Jiwon Gim <55209567+boulderdaze@users.noreply.github.com> Date: Wed, 30 Oct 2024 10:56:59 -0600 Subject: [PATCH 1/5] Update surface albedo in TUV-x prior to calculating rate constants (#141) Originator(s): @mattldawson, @boulderdaze Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): - Updates TUV-x wrapper, tests, and configuration data to allow setting the wavelength grid from the host model. This implementation assumes that there will be a standard-named variable that defines the edges of the wavelength grid to use for photolysis rate constant calculations that is defined by the host model. - Updates the profile of surface albedos in TUV-x for each column prior to calculating photolysis rate constants. - Closes #96 and #123 Describe any changes made to the namelist: N/A List all files eliminated and why: N/A List all files added and what they do: - A schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 - A schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 - A schemes/musica/tuvx/musica_ccpp_tuvx_wavelength_grid.F90 - A test/musica/tuvx/test_tuvx_surface_albedo.F90 - A test/musica/tuvx/test_tuvx_wavelength_grid.F90 List all existing files that have been modified, and describe the changes: - M schemes/musica/musica_ccpp.F90 - M schemes/musica/musica_ccpp.meta - M schemes/musica/tuvx/musica_ccpp_tuvx.F90 - M test/musica/micm/test_micm_util.F90 - M test/musica/test_musica_api.F90 - M test/musica/tuvx/CMakeLists.txt - M test/musica/tuvx/configs/ts1_tsmlt.json - M test/musica/tuvx/test_tuvx_height_grid.F90 - M test/musica/tuvx/test_tuvx_temperature.F90 List any test failures: N/A Is this a science-changing update? New physics package, algorithm change, tuning changes, etc? No --------- Co-authored-by: Matt Dawson --- schemes/musica/micm/musica_ccpp_micm.F90 | 8 +- schemes/musica/musica_ccpp.F90 | 25 +-- schemes/musica/musica_ccpp.meta | 12 ++ schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 170 +++++++++++------- .../tuvx/musica_ccpp_tuvx_surface_albedo.F90 | 76 ++++++++ schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 | 55 ++++++ .../tuvx/musica_ccpp_tuvx_wavelength_grid.F90 | 59 ++++++ test/musica/micm/test_micm_util.F90 | 50 +++--- test/musica/test_musica_api.F90 | 155 +++++++++++++--- test/musica/tuvx/CMakeLists.txt | 59 +++++- test/musica/tuvx/configs/ts1_tsmlt.json | 16 -- test/musica/tuvx/test_tuvx_height_grid.F90 | 80 ++++----- test/musica/tuvx/test_tuvx_surface_albedo.F90 | 57 ++++++ test/musica/tuvx/test_tuvx_temperature.F90 | 42 ++--- .../musica/tuvx/test_tuvx_wavelength_grid.F90 | 53 ++++++ 15 files changed, 712 insertions(+), 205 deletions(-) create mode 100644 schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 create mode 100644 schemes/musica/tuvx/musica_ccpp_tuvx_util.F90 create mode 100644 schemes/musica/tuvx/musica_ccpp_tuvx_wavelength_grid.F90 create mode 100644 test/musica/tuvx/test_tuvx_surface_albedo.F90 create mode 100644 test/musica/tuvx/test_tuvx_wavelength_grid.F90 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 From de3c61b3eb77802e382b309cc2b52108b981be1a Mon Sep 17 00:00:00 2001 From: Courtney Peverley Date: Tue, 5 Nov 2024 10:56:19 -0700 Subject: [PATCH 2/5] Replace dynamic constituent routine with register phase (#147) Originator(s): peverwhee Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): Modifies existing dynamic constituent handling for the MUSICA interface to use the new register phase Describe any changes made to the namelist: N/A List all files eliminated and why: N/A List all files added and what they do: N/A List all existing files that have been modified, and describe the changes: (Helpful git command: git diff --name-status development...) M schemes/musica/musica_ccpp.F90 M schemes/musica/muscia_ccpp.meta - remove dynamic_constituent_routine from metadata - add necessary comment header to Fortran - add metadata for register routine List any test failures: Is this a science-changing update? New physics package, algorithm change, tuning changes, etc? No --------- Co-authored-by: Courtney Peverley --- schemes/musica/musica_ccpp.F90 | 2 ++ schemes/musica/musica_ccpp.meta | 38 +++++++++++++++++++++++++++++++-- 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index 09364397..1ef2ea63 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -10,6 +10,8 @@ module musica_ccpp contains + !> \section arg_table_musica_ccpp_register Argument Table + !! \htmlinclude musica_ccpp_register.html subroutine musica_ccpp_register(solver_type, num_grid_cells, constituent_props, errmsg, errcode) use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t diff --git a/schemes/musica/musica_ccpp.meta b/schemes/musica/musica_ccpp.meta index 30d19bad..85e01f3e 100644 --- a/schemes/musica/musica_ccpp.meta +++ b/schemes/musica/musica_ccpp.meta @@ -2,7 +2,41 @@ name = musica_ccpp type = scheme dependencies = micm/musica_ccpp_micm.F90,micm/musica_ccpp_micm_util.F90,tuvx/musica_ccpp_tuvx.F90,tuvx/musica_ccpp_tuvx_height_grid.F90,musica_ccpp_util.F90 - dynamic_constituent_routine = musica_ccpp_register + +[ccpp-arg-table] + name = musica_ccpp_register + type = scheme +[ solver_type ] + standard_name = micm_solver_type + units = none + type = integer + dimensions = () + intent = in +[ num_grid_cells ] + standard_name = number_of_grid_cells + units = count + type = integer + dimensions = () + intent = in +[ constituent_props ] + standard_name = dynamic_constituents_for_musica_ccpp + units = none + dimensions = (:) + allocatable = True + type = ccpp_constituent_properties_t + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errcode ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out [ccpp-arg-table] name = musica_ccpp_init @@ -140,4 +174,4 @@ units = 1 type = integer dimensions = () - intent = out \ No newline at end of file + intent = out From 026c4c5e01afb15b26d3d754997bf8a7d0ab4378 Mon Sep 17 00:00:00 2001 From: mwaxmonsky <137746677+mwaxmonsky@users.noreply.github.com> Date: Thu, 14 Nov 2024 13:53:43 -0500 Subject: [PATCH 3/5] Unit test infrastructure setup (#144) Originator(s): @mwaxmonsky Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): Sets up the basic infrastructure to start enabling unit testing. Addresses #82 Describe any changes made to the namelist: N/A List all files eliminated and why: N/A List all files added and what they do: A .github/workflows/code-coverage.yaml A test/unit-test/CMakeLists.txt A test/unit-test/include/ccpp_kinds.F90 A test/unit-test/tests/CMakeLists.txt A test/unit-test/tests/utilities/CMakeLists.txt A test/unit-test/tests/utilities/test_state_converters.pf - Initial round of infrastructure to start enabling unit tests. List all existing files that have been modified, and describe the changes: N/A (Helpful git command: `git diff --name-status development...`) List any test failures: N/A Is this a science-changing update? New physics package, algorithm change, tuning changes, etc? No --- .github/workflows/unit-tests.yaml | 72 +++++++++++++++++++ test/unit-test/CMakeLists.txt | 39 ++++++++++ test/unit-test/README.md | 26 +++++++ test/unit-test/include/ccpp_kinds.F90 | 11 +++ test/unit-test/tests/CMakeLists.txt | 2 + test/unit-test/tests/utilities/CMakeLists.txt | 5 ++ .../tests/utilities/test_state_converters.pf | 27 +++++++ 7 files changed, 182 insertions(+) create mode 100644 .github/workflows/unit-tests.yaml create mode 100644 test/unit-test/CMakeLists.txt create mode 100644 test/unit-test/README.md create mode 100644 test/unit-test/include/ccpp_kinds.F90 create mode 100644 test/unit-test/tests/CMakeLists.txt create mode 100644 test/unit-test/tests/utilities/CMakeLists.txt create mode 100644 test/unit-test/tests/utilities/test_state_converters.pf diff --git a/.github/workflows/unit-tests.yaml b/.github/workflows/unit-tests.yaml new file mode 100644 index 00000000..d9003cee --- /dev/null +++ b/.github/workflows/unit-tests.yaml @@ -0,0 +1,72 @@ +name: unit-test-code-coverage + +on: + push: + branches: + - development + - main + workflow_dispatch: + pull_request: + +concurrency: + group: ${{ github.workflow }}-${{ github.ref || github.run_id }} + cancel-in-progress: true + +jobs: + gcc-toolchain: + runs-on: ubuntu-latest + steps: + - name: Checkout atmospheric_physics + uses: actions/checkout@v4 + + - name: Install dependencies + run: | + sudo apt update && sudo apt -y install libopenmpi-dev openmpi-bin + + - name: Build pFUnit + run: | + git clone --depth 1 --branch v4.10.0 https://github.com/Goddard-Fortran-Ecosystem/pFUnit.git + cd pFUnit + pwd + cmake -B./build -S. + cd build + make install + + - name: Build atmospheric_physics + run: | + cmake \ + -DCMAKE_PREFIX_PATH=/home/runner/work/atmospheric_physics/atmospheric_physics/pFUnit/build/installed \ + -DATMOSPHERIC_PHYSICS_ENABLE_CODE_COVERAGE=ON \ + -B./build \ + -S./test/unit-test + cd build + make + + - name: Run tests + run: | + cd build && ctest -V --output-on-failure --output-junit test_results.xml + + - name: Upload unit test results + uses: actions/upload-artifact@v4 + with: + name: unit-test-results + path: build/test_results.xml + + - name: Setup GCov + run: | + python3 -m venv venv + source venv/bin/activate + pip3 install gcovr + + - name: Run Gcov + run: | + source venv/bin/activate + cd build + gcovr -r .. --filter '\.\./schemes' --html atmospheric_physics_code_coverage.html --txt + + - name: Upload code coverage results + uses: actions/upload-artifact@v4 + with: + name: code-coverage-results + path: build/atmospheric_physics_code_coverage.html + diff --git a/test/unit-test/CMakeLists.txt b/test/unit-test/CMakeLists.txt new file mode 100644 index 00000000..d307934e --- /dev/null +++ b/test/unit-test/CMakeLists.txt @@ -0,0 +1,39 @@ +cmake_minimum_required(VERSION 3.17) + +project(atmospheric_physics VERSION 0.0.1 LANGUAGES Fortran) + +find_package(PFUNIT REQUIRED) + +if(NOT ATMOSPHERIC_PHYSICS_IS_TOP_LEVEL) + message(WARNING "atmospheric-physics is not integrated into the CMake build of any top level " + "project yet and this CMake is for testing purposes only. " + "Making a change to this project's CMake will not impact the build of " + "a parent project at this time.") +endif() + +option(ATMOSPHERIC_PHYSICS_ENABLE_TESTS "Run pFUnit unit tests" OFF) +option(ATMOSPHERIC_PHYSICS_ENABLE_CODE_COVERAGE "Run code coverage tool" OFF) + +if(ATMOSPHERIC_PHYSICS_ENABLE_CODE_COVERAGE) + add_compile_options(-O0 --coverage) + add_link_options(--coverage) +endif() + +set(CMAKE_BUILD_TYPE Debug) + +set(UTILITIES_SRC + ../../schemes/utilities/state_converters.F90 + ../../schemes/utilities/static_energy.F90 + ../../schemes/utilities/physics_tendency_updaters.F90 + include/ccpp_kinds.F90 +) + +add_library(utilities ${UTILITIES_SRC}) +target_compile_options(utilities PRIVATE -ffree-line-length-none) +target_include_directories(utilities PUBLIC ${CMAKE_CURRENT_BINARY_DIR}) + +if(ATMOSPHERIC_PHYSICS_ENABLE_TESTS OR ATMOSPHERIC_PHYSICS_ENABLE_CODE_COVERAGE) + enable_testing() + add_subdirectory(tests) +endif() + diff --git a/test/unit-test/README.md b/test/unit-test/README.md new file mode 100644 index 00000000..d46c7afd --- /dev/null +++ b/test/unit-test/README.md @@ -0,0 +1,26 @@ +# Unit Tests + +To add or update unit tests, follow the instructions from the [development guide](https://escomp.github.io/CAM-SIMA-docs/atmospheric_physics/development_workflow/#5-unit-testing). Also, make sure [pFUnit](https://github.com/Goddard-Fortran-Ecosystem/pFUnit) is built and installed following the [build directions](https://github.com/Goddard-Fortran-Ecosystem/pFUnit?tab=readme-ov-file#building-and-installing-pfunit) (see the atmospheric_physics [github workflow file](../../.github/workflows/unit-tests.yaml) for a more detailed example of how to build pFUnit): + +```bash +$ git clone --depth 1 --branch v4.10.0 https://github.com/Goddard-Fortran-Ecosystem/pFUnit.git +$ cd pFUnit +$ cmake -B./build -S. +$ cd build +$ make install +``` + +To run the tests, from the root directory of your clone, run: + +```bash +$ cmake \ + -DCMAKE_PREFIX_PATH=/build/installed \ + -DATMOSPHERIC_PHYSICS_ENABLE_CODE_COVERAGE=ON \ + -B./build \ + -S./test/unit-test +$ cd build +$ make +$ ctest -V --output-on-failure +``` + +Where `` is the path to your pfunit repository. The install path of pFUnit may be different depending on how you've built your local verison (see the example in the [workflow file](../../.github/workflows/unit-tests.yaml)). diff --git a/test/unit-test/include/ccpp_kinds.F90 b/test/unit-test/include/ccpp_kinds.F90 new file mode 100644 index 00000000..d56949b3 --- /dev/null +++ b/test/unit-test/include/ccpp_kinds.F90 @@ -0,0 +1,11 @@ +module ccpp_kinds + + use ISO_FORTRAN_ENV, only: kind_phys => REAL64 + + implicit none + private + + public :: kind_phys + + end module ccpp_kinds + diff --git a/test/unit-test/tests/CMakeLists.txt b/test/unit-test/tests/CMakeLists.txt new file mode 100644 index 00000000..b4c4714c --- /dev/null +++ b/test/unit-test/tests/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory(utilities) + diff --git a/test/unit-test/tests/utilities/CMakeLists.txt b/test/unit-test/tests/utilities/CMakeLists.txt new file mode 100644 index 00000000..bca98a41 --- /dev/null +++ b/test/unit-test/tests/utilities/CMakeLists.txt @@ -0,0 +1,5 @@ +add_pfunit_ctest(utilities_tests + TEST_SOURCES test_state_converters.pf + LINK_LIBRARIES utilities +) + diff --git a/test/unit-test/tests/utilities/test_state_converters.pf b/test/unit-test/tests/utilities/test_state_converters.pf new file mode 100644 index 00000000..bc56f038 --- /dev/null +++ b/test/unit-test/tests/utilities/test_state_converters.pf @@ -0,0 +1,27 @@ +@test +subroutine test_temp_to_potential_temp() + use funit + use state_converters, only : temp_to_potential_temp_run + use ccpp_kinds, only: kind_phys + + integer, parameter :: ncol = 5 + integer, parameter :: nz = 5 + + real(kind_phys) :: temp(ncol, nz) + real(kind_phys) :: exner(ncol, nz) + real(kind_phys) :: theta(ncol, nz) + character(len=512) :: errmsg + integer :: errflg + + temp = 1 + exner = 1 + theta = 1 + + errmsg = "" + errflg = 0 + + call temp_to_potential_temp_run(ncol, nz, temp, exner, theta, errmsg, errflg) + + @assertEqual(0, errflg) + +end subroutine test_temp_to_potential_temp From 0ddf4d25b51d22767e58f786c855640f92c87e5d Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Thu, 14 Nov 2024 11:23:12 -0800 Subject: [PATCH 4/5] Add Chapman and Terminator configurations for chemistry (#151) Adds configurations for Terminator and Chapman chemistry and basic integration tests. This includes some data files needed for TUV-x to calculate photolysis rate constants for these mechanisms, but I wasn't sure if this is the best place for them. I'm happy to remove the configuration data from this PR if there is a better place to keep them. closes #150 Also: - Adds calculation of photolysis rate constants from TUV-x - Maps TUV-x rate constants to MICM rate parameters - Removes some array copying using updated MUSICA functions that accept pointers to Fortran arrays - Sets up a gitignore file for the repo --- .gitignore | 7 + schemes/musica/micm/musica_ccpp_micm.F90 | 47 ++- schemes/musica/micm/musica_ccpp_micm_util.F90 | 72 +--- schemes/musica/musica_ccpp.F90 | 82 ++-- schemes/musica/musica_ccpp_namelist.xml | 12 + schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 130 +++++-- test/docker/Dockerfile.musica | 14 +- test/docker/Dockerfile.musica.no_install | 14 +- test/musica/CMakeLists.txt | 8 +- test/musica/micm/configs/chapman/config.json | 6 - .../micm/configs/chapman/reactions.json | 97 ----- test/musica/micm/configs/chapman/species.json | 39 -- test/musica/micm/test_micm_util.F90 | 61 --- test/musica/musica_ccpp_namelist.F90 | 8 +- test/musica/test_musica_api.F90 | 366 +++++++++++++++--- 15 files changed, 504 insertions(+), 459 deletions(-) create mode 100644 .gitignore delete mode 100644 test/musica/micm/configs/chapman/config.json delete mode 100644 test/musica/micm/configs/chapman/reactions.json delete mode 100644 test/musica/micm/configs/chapman/species.json diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..83c5d6be --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +**/*.DS_Store +**/*.pyc +**/build/ +**/CMakeCache.txt +**/CMakeFiles/ +.vscode +xcode/ \ No newline at end of file diff --git a/schemes/musica/micm/musica_ccpp_micm.F90 b/schemes/musica/micm/musica_ccpp_micm.F90 index 1ef50d5a..7b9fbae9 100644 --- a/schemes/musica/micm/musica_ccpp_micm.F90 +++ b/schemes/musica/micm/musica_ccpp_micm.F90 @@ -2,17 +2,18 @@ module musica_ccpp_micm ! Note: "micm_t" is included in an external pre-built MICM library that the host ! model is responsible for linking to during compilation - use musica_micm, only: micm_t + use ccpp_kinds, only: kind_phys use musica_ccpp_util, only: has_error_occurred use musica_ccpp_namelist, only: filename_of_micm_configuration - use ccpp_kinds, only: kind_phys + use musica_micm, only: micm_t implicit none private - public :: micm_register, micm_init, micm_run, micm_final + public :: micm_register, micm_init, micm_run, micm_final, micm, number_of_rate_parameters type(micm_t), pointer :: micm => null( ) + integer :: number_of_rate_parameters = 0 contains @@ -36,7 +37,7 @@ subroutine micm_register(solver_type, num_grid_cells, constituent_props, errmsg, logical :: is_advected integer :: i, species_index - micm => micm_t(filename_of_micm_configuration, solver_type, num_grid_cells, error) + micm => micm_t(trim(filename_of_micm_configuration), solver_type, num_grid_cells, error) if (has_error_occurred(error, errmsg, errcode)) return allocate(constituent_props(micm%species_ordering%size()), stat=errcode) @@ -73,6 +74,7 @@ subroutine micm_register(solver_type, num_grid_cells, constituent_props, errmsg, if (errcode /= 0) return end associate ! map end do + number_of_rate_parameters = micm%user_defined_reaction_rates%size() end subroutine micm_register @@ -91,34 +93,31 @@ subroutine micm_run(time_step, temperature, pressure, dry_air_density, & user_defined_rate_parameters, constituents, errmsg, errcode) use musica_micm, only: solver_stats_t use musica_util, only: string_t, error_t - use iso_c_binding, only: c_double + use iso_c_binding, only: c_double, c_loc - real(kind_phys), intent(in) :: time_step ! s - real(c_double), intent(in) :: temperature(:) ! K - real(c_double), intent(in) :: pressure(:) ! Pa - real(c_double), intent(in) :: dry_air_density(:) ! kg m-3 - real(c_double), intent(in) :: user_defined_rate_parameters(:) ! various units - real(c_double), intent(inout) :: constituents(:) ! mol m-3 - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode + real(kind_phys), intent(in) :: time_step ! s + real(kind_phys), target, intent(in) :: temperature(:,:) ! K + real(kind_phys), target, intent(in) :: pressure(:,:) ! Pa + real(kind_phys), target, intent(in) :: dry_air_density(:,:) ! kg m-3 + real(kind_phys), target, intent(in) :: user_defined_rate_parameters(:,:,:) ! various units + real(kind_phys), target, intent(inout) :: constituents(:,:,:) ! mol m-3 + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode ! local variables type(string_t) :: solver_state type(solver_stats_t) :: solver_stats type(error_t) :: error - real(c_double) :: c_time_step integer :: i_elem - c_time_step = real(time_step, c_double) - - call micm%solve(c_time_step, & - temperature, & - pressure, & - dry_air_density, & - constituents, & - user_defined_rate_parameters, & - solver_state, & - solver_stats, & + call micm%solve(real(time_step, kind=c_double), & + c_loc(temperature), & + c_loc(pressure), & + c_loc(dry_air_density), & + c_loc(constituents), & + c_loc(user_defined_rate_parameters), & + solver_state, & + solver_stats, & error) if (has_error_occurred(error, errmsg, errcode)) return diff --git a/schemes/musica/micm/musica_ccpp_micm_util.F90 b/schemes/musica/micm/musica_ccpp_micm_util.F90 index 47b92c81..2e12de18 100644 --- a/schemes/musica/micm/musica_ccpp_micm_util.F90 +++ b/schemes/musica/micm/musica_ccpp_micm_util.F90 @@ -2,80 +2,10 @@ module musica_ccpp_micm_util implicit none private - public :: reshape_into_micm_arr, reshape_into_ccpp_arr, convert_to_mol_per_cubic_meter, & - convert_to_mass_mixing_ratio + public :: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio contains - !> Reshape array (2D/3D -> 1D) and convert type (kind_phys -> c_double) - subroutine reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, & - micm_temperature, micm_pressure, micm_dry_air_density, & - micm_constituents) - use iso_c_binding, only: c_double - use ccpp_kinds, only: kind_phys - - real(kind_phys), intent(in) :: temperature(:,:) ! K - real(kind_phys), intent(in) :: pressure(:,:) ! Pa - real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 - real(kind_phys), intent(in) :: constituents(:,:,:) ! kg kg-1 - real(c_double), intent(out) :: micm_temperature(:) ! K - real(c_double), intent(out) :: micm_pressure(:) ! Pa - real(c_double), intent(out) :: micm_dry_air_density(:) ! kg m-3 - real(c_double), intent(out) :: micm_constituents(:) ! kg kg-1 - - ! local variables - integer :: num_columns, num_layers, num_constituents - integer :: i_column, i_layer, i_elem, i_constituents - - num_columns = size(constituents, dim=1) - num_layers = size(constituents, dim=2) - num_constituents = size(constituents, dim=3) - - ! Reshape into 1-D arry in species-column first order, referring to - ! state.variables_[i_cell][i_species] = concentrations[i_species_elem++] - i_elem = 1 - i_constituents = 1 - do i_layer = 1, num_layers - do i_column = 1, num_columns - micm_temperature(i_elem) = real(temperature(i_column, i_layer), c_double) - micm_pressure(i_elem) = real(pressure(i_column, i_layer), c_double) - micm_dry_air_density(i_elem) = real(dry_air_density(i_column, i_layer), c_double) - micm_constituents(i_constituents : i_constituents + num_constituents - 1) & - = real(constituents(i_column, i_layer, :), c_double) - i_elem = i_elem + 1 - i_constituents = i_constituents + num_constituents - end do - end do - - end subroutine reshape_into_micm_arr - - !> Reshape array (1D -> 3D) and convert type (c_double -> kind_phys) - subroutine reshape_into_ccpp_arr(micm_constituents, constituents) - use iso_c_binding, only: c_double - use ccpp_kinds, only: kind_phys - - real(c_double), intent(in) :: micm_constituents(:) ! kg kg-1 - real(kind_phys), intent(inout) :: constituents(:,:,:) ! kg kg-1 - - ! local variables - integer :: num_columns, num_layers, num_constituents - integer :: i_column, i_layer, i_constituents - - num_columns = size(constituents, dim=1) - num_layers = size(constituents, dim=2) - num_constituents = size(constituents, dim=3) - - i_constituents = 1 - do i_layer = 1, num_layers - do i_column = 1, num_columns - constituents(i_column, i_layer, :) & - = real(micm_constituents(i_constituents : i_constituents + num_constituents - 1), kind_phys) - i_constituents = i_constituents + num_constituents - end do - end do - - end subroutine reshape_into_ccpp_arr - ! Convert CAM-SIMA unit to MICM unit (kg kg-1 -> mol m-3) subroutine convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents) use ccpp_kinds, only: kind_phys diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index 1ef2ea63..d3512f4a 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -1,7 +1,9 @@ !> Top-level wrapper for MUSICA chemistry components module musica_ccpp - use musica_ccpp_micm, only: micm_register, micm_init, micm_run, micm_final - use musica_ccpp_tuvx, only: tuvx_init, tuvx_run, tuvx_final + use musica_ccpp_micm, only: micm_register, micm_init, micm_run, micm_final + use musica_ccpp_namelist, only: filename_of_tuvx_micm_mapping_configuration + use musica_ccpp_tuvx, only: tuvx_init, tuvx_run, tuvx_final + use musica_util, only: index_mappings_t implicit none private @@ -30,18 +32,20 @@ end subroutine musica_ccpp_register subroutine musica_ccpp_init(vertical_layer_dimension, vertical_interface_dimension, & photolysis_wavelength_grid_interfaces, errmsg, errcode) use ccpp_kinds, only : kind_phys - + use musica_ccpp_micm, only: micm + use musica_ccpp_util, only: has_error_occurred 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, & - photolysis_wavelength_grid_interfaces, errmsg, errcode) - if (errcode /= 0) return call micm_init(errmsg, errcode) if (errcode /= 0) return + call tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & + photolysis_wavelength_grid_interfaces, & + micm%user_defined_reaction_rates, errmsg, errcode) + if (errcode /= 0) return end subroutine musica_ccpp_init @@ -56,53 +60,42 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co geopotential_height_wrt_surface_at_interface, surface_temperature, & 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 use ccpp_kinds, only: kind_phys - use iso_c_binding, only: c_double + use musica_ccpp_micm, only: number_of_rate_parameters + use musica_ccpp_micm_util, only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio - real(kind_phys), intent(in) :: time_step ! s - real(kind_phys), intent(in) :: temperature(:,:) ! K - real(kind_phys), intent(in) :: pressure(:,:) ! Pa - real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 + real(kind_phys), intent(in) :: time_step ! s + real(kind_phys), target, intent(in) :: temperature(:,:) ! K + real(kind_phys), target, intent(in) :: pressure(:,:) ! Pa + real(kind_phys), target, intent(in) :: dry_air_density(:,:) ! kg m-3 type(ccpp_constituent_prop_ptr_t), & - intent(in) :: constituent_props(:) - real(kind_phys), intent(inout) :: constituents(:,:,:) ! kg kg-1 - real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_midpoint(:,:) ! m (column, layer) - real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface) - real(kind_phys), intent(in) :: surface_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 + intent(in) :: constituent_props(:) + real(kind_phys), target, intent(inout) :: constituents(:,:,:) ! kg kg-1 + real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_midpoint(:,:) ! m (column, layer) + real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface) + real(kind_phys), intent(in) :: surface_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 ! local variables - real(c_double), dimension(size(temperature, dim=1) & - * size(temperature, dim=2)) :: micm_temperature - real(c_double), dimension(size(pressure, dim=1) & - * size(pressure, dim=2)) :: micm_pressure - real(c_double), dimension(size(dry_air_density, dim=1) & - * size(dry_air_density, dim=2)) :: micm_dry_air_density - real(c_double), dimension(size(constituents, dim=1) & - * size(constituents, dim=2) & - * size(constituents, dim=3)) :: micm_constituents ! mol m-3 real(kind_phys), dimension(size(constituents, dim=3)) :: molar_mass_arr ! kg mol-1 - - ! temporarily dimensioned to Chapman mechanism until mapping between MICM and TUV-x is implemented - real(c_double), dimension(size(constituents, dim=1) & - * size(constituents, dim=2) & - * 3) :: photolysis_rate_constants ! s-1 + real(kind_phys), dimension(size(constituents, dim=1), & + size(constituents, dim=2), & + number_of_rate_parameters) :: rate_parameters ! various units integer :: i_elem + ! Calculate photolysis rate constants using TUV-x call 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, & + rate_parameters, & errmsg, errcode) ! Get the molar mass that is set in the call to instantiate() @@ -127,16 +120,9 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co ! Convert CAM-SIMA unit to MICM unit (kg kg-1 -> mol m-3) call convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents) - ! Reshape array (2D/3D -> 1D) and convert type (kind_phys -> c_double) - call reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, & - micm_temperature, micm_pressure, micm_dry_air_density, micm_constituents) - - ! temporarily pass in unmapped photolysis rate constants until mapping between MICM and TUV-x is implemented - call micm_run(time_step, micm_temperature, micm_pressure, micm_dry_air_density, & - photolysis_rate_constants, micm_constituents, errmsg, errcode) - - ! Reshape array (1D -> 3D) and convert type (c_double -> kind_phys) - call reshape_into_ccpp_arr(micm_constituents, constituents) + ! Solve chemistry at the current time step + call micm_run(time_step, temperature, pressure, dry_air_density, rate_parameters, & + constituents, errmsg, errcode) ! Convert MICM unit back to CAM-SIMA unit (mol m-3 -> kg kg-1) call convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constituents) diff --git a/schemes/musica/musica_ccpp_namelist.xml b/schemes/musica/musica_ccpp_namelist.xml index a8417876..5015f2d3 100644 --- a/schemes/musica/musica_ccpp_namelist.xml +++ b/schemes/musica/musica_ccpp_namelist.xml @@ -103,5 +103,17 @@ UNSET_PATH + + char*512 + musica_ccpp + musica_ccpp + filename_of_tuvx_micm_mapping_configuration + none + + A configuration file for the mapping of TUV-x photolysis rates to MICM custom rate parameters + + + UNSET_PATH + diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index bcc7be09..5a12a3ee 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -2,10 +2,11 @@ module musica_ccpp_tuvx ! 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 use ccpp_kinds, only: kind_phys use musica_ccpp_namelist, only: filename_of_tuvx_configuration + use musica_ccpp_util, only: has_error_occurred + use musica_tuvx, only: tuvx_t, grid_t, profile_t + use musica_util, only: mappings_t, index_mappings_t implicit none private @@ -18,13 +19,18 @@ module musica_ccpp_tuvx type(profile_t), pointer :: temperature_profile => null() type(profile_t), pointer :: surface_albedo_profile => null() + type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( ) + integer :: number_of_photolysis_rate_constants = 0 + contains !> Initializes TUV-x subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & - wavelength_grid_interfaces, errmsg, errcode) + wavelength_grid_interfaces, micm_rate_parameter_ordering, & + errmsg, errcode) use musica_tuvx, only: grid_map_t, profile_map_t, radiator_map_t - use musica_util, only: error_t + use musica_util, only: error_t, configuration_t + use musica_ccpp_namelist, only: filename_of_tuvx_micm_mapping_configuration use musica_ccpp_tuvx_util, only: tuvx_deallocate use musica_ccpp_tuvx_height_grid, & only: create_height_grid, height_grid_label, height_grid_unit @@ -38,14 +44,17 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & integer, intent(in) :: vertical_layer_dimension ! (count) integer, intent(in) :: vertical_interface_dimension ! (count) real(kind_phys), intent(in) :: wavelength_grid_interfaces(:) ! m + type(mappings_t), intent(in) :: micm_rate_parameter_ordering ! index mappings for MICM rate parameters character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode ! local variables - type(grid_map_t), pointer :: grids - type(profile_map_t), pointer :: profiles - type(radiator_map_t), pointer :: radiators - type(error_t) :: error + type(grid_map_t), pointer :: grids + type(profile_map_t), pointer :: profiles + type(radiator_map_t), pointer :: radiators + type(configuration_t) :: config + type(mappings_t), pointer :: photolysis_rate_constants_ordering + type(error_t) :: error grids => grid_map_t( error ) if (has_error_occurred( error, errmsg, errcode )) return @@ -122,7 +131,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & return end if - tuvx => tuvx_t( filename_of_tuvx_configuration, grids, profiles, & + tuvx => tuvx_t( trim(filename_of_tuvx_configuration), grids, profiles, & radiators, error ) if (has_error_occurred( error, errmsg, errcode )) then call tuvx_deallocate( grids, profiles, radiators, null(), height_grid, & @@ -134,11 +143,7 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & wavelength_grid, temperature_profile, surface_albedo_profile ) grids => tuvx%get_grids( error ) - if (has_error_occurred( error, errmsg, errcode )) then - deallocate( tuvx ) - tuvx => null() - return - end if + if (has_error_occurred( error, errmsg, errcode )) return height_grid => grids%get( height_grid_label, height_grid_unit, error ) if (has_error_occurred( error, errmsg, errcode )) then @@ -179,6 +184,25 @@ subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & call tuvx_deallocate( grids, profiles, null(), null(), null(), null(), & null(), null() ) + photolysis_rate_constants_ordering => & + tuvx%get_photolysis_rate_constants_ordering( error ) + if (has_error_occurred( error, errmsg, errcode )) return + number_of_photolysis_rate_constants = photolysis_rate_constants_ordering%size() + + call config%load_from_file( trim(filename_of_tuvx_micm_mapping_configuration), error ) + if (has_error_occurred( error, errmsg, errcode )) then + deallocate( photolysis_rate_constants_ordering ) + photolysis_rate_constants_ordering => null() + return + end if + + photolysis_rate_constants_mapping => & + index_mappings_t( config, photolysis_rate_constants_ordering, & + micm_rate_parameter_ordering, error ) + deallocate( photolysis_rate_constants_ordering ) + photolysis_rate_constants_ordering => null() + if (has_error_occurred( error, errmsg, errcode )) return + end subroutine tuvx_init !> Calculates photolysis rate constants for the current model conditions @@ -188,34 +212,43 @@ subroutine tuvx_run(temperature, dry_air_density, & 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 + rate_parameters, 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_ccpp_util, only: has_error_occurred 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) - real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_midpoint(:,:) ! m (column, layer) - real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface) - real(kind_phys), intent(in) :: surface_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 - real(kind_phys), intent(out) :: photolysis_rate_constants(:) ! s-1 (column, reaction) - character(len=512), intent(out) :: errmsg - integer, intent(out) :: errcode + real(kind_phys), intent(in) :: temperature(:,:) ! K (column, layer) + real(kind_phys), intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer) + real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_midpoint(:,:) ! m (column, layer) + real(kind_phys), intent(in) :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface) + real(kind_phys), intent(in) :: surface_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 + real(kind_phys), intent(inout) :: rate_parameters(:,:,:) ! various units (column, layer, reaction) + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode ! local variables real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_midpoint, dim = 2)) :: height_midpoints real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_interface, dim = 2)) :: height_interfaces real(kind_phys) :: reciprocal_of_gravitational_acceleration ! s2 m-1 - integer :: i_col + real(kind_phys), dimension(size(rate_parameters, dim=2)+2, & + number_of_photolysis_rate_constants) :: photolysis_rate_constants, & ! s-1 + heating_rates ! K s-1 (TODO: check units) + real(kind_phys) :: solar_zenith_angle ! degrees + real(kind_phys) :: earth_sun_distance ! AU + type(error_t) :: error + integer :: i_col, i_level reciprocal_of_gravitational_acceleration = 1.0_kind_phys / standard_gravitational_acceleration + ! 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 + do i_col = 1, size(temperature, dim=1) call calculate_heights( geopotential_height_wrt_surface_at_midpoint(i_col,:), & geopotential_height_wrt_surface_at_interface(i_col,:), & @@ -229,14 +262,28 @@ subroutine tuvx_run(temperature, dry_air_density, & call set_temperature_values( temperature_profile, temperature(i_col,:), & surface_temperature(i_col), errmsg, errcode ) 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 + ! temporary values until these are available from the host model + solar_zenith_angle = 0.0_kind_phys + earth_sun_distance = 1.0_kind_phys + + ! calculate photolysis rate constants and heating rates + call tuvx%run( solar_zenith_angle, earth_sun_distance, & + photolysis_rate_constants(:,:), heating_rates(:,:), & + error ) + if (has_error_occurred( error, errmsg, errcode )) return + + ! filter out negative photolysis rate constants + photolysis_rate_constants(:,:) = & + max( photolysis_rate_constants(:,:), 0.0_kind_phys ) + + ! map photolysis rate constants to the host model's rate parameters and vertical grid + do i_level = 1, size(rate_parameters, dim=2) + call photolysis_rate_constants_mapping%copy_data( & + photolysis_rate_constants(size(rate_parameters, dim=2)-i_level+2,:), & + rate_parameters(i_col,i_level,:) ) + end do + end do end subroutine tuvx_run @@ -273,6 +320,11 @@ subroutine tuvx_final(errmsg, errcode) tuvx => null() end if + if (associated( photolysis_rate_constants_mapping )) then + deallocate( photolysis_rate_constants_mapping ) + photolysis_rate_constants_mapping => null() + end if + end subroutine tuvx_final end module musica_ccpp_tuvx \ No newline at end of file diff --git a/test/docker/Dockerfile.musica b/test/docker/Dockerfile.musica index 65940f1a..83dea485 100644 --- a/test/docker/Dockerfile.musica +++ b/test/docker/Dockerfile.musica @@ -5,7 +5,7 @@ FROM ubuntu:22.04 -ARG MUSICA_GIT_TAG=dbbdb22f5f2807e27c2695db85291951ba178634 +ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8 RUN apt update \ && apt install -y sudo \ @@ -68,6 +68,12 @@ RUN cd musica \ COPY . atmospheric_physics RUN sudo chown -R test_user:test_user atmospheric_physics +# Clone the MUSICA chemistry data set repository +RUN git clone https://github.com/NCAR/cam-sima-chemistry-data.git \ + && cd cam-sima-chemistry-data \ + && git checkout 144d982715f5bd6fa61e76f255a52db1843c55f2 \ + && mv mechanisms /home/test_user/atmospheric_physics/schemes/musica/configurations + # Must make ccpp-framework available before building test RUN cd atmospheric_physics/test \ && mkdir lib \ @@ -84,10 +90,4 @@ RUN cd atmospheric_physics/test \ -D CCPP_ENABLE_MEMCHECK=ON \ && cmake --build ./build -RUN cd atmospheric_physics/test \ - && mkdir third_party \ - && cd third_party \ - && git clone --depth 1 https://github.com/NCAR/tuv-x.git \ - && cp -r tuv-x/data ../build/data - WORKDIR /home/test_user/atmospheric_physics/test/build \ No newline at end of file diff --git a/test/docker/Dockerfile.musica.no_install b/test/docker/Dockerfile.musica.no_install index b59489af..e8f70fb0 100644 --- a/test/docker/Dockerfile.musica.no_install +++ b/test/docker/Dockerfile.musica.no_install @@ -8,7 +8,7 @@ FROM ubuntu:22.04 -ARG MUSICA_GIT_TAG=dbbdb22f5f2807e27c2695db85291951ba178634 +ARG MUSICA_GIT_TAG=326b5119768d5be9654baf96ae3bd6a1b757fdc8 RUN apt update \ && apt install -y sudo \ @@ -55,6 +55,12 @@ ENV MUSICA_GIT_TAG=${MUSICA_GIT_TAG} COPY . atmospheric_physics RUN sudo chown -R test_user:test_user atmospheric_physics +# Clone the MUSICA chemistry data set repository +RUN git clone https://github.com/NCAR/cam-sima-chemistry-data.git \ + && cd cam-sima-chemistry-data \ + && git checkout 144d982715f5bd6fa61e76f255a52db1843c55f2 \ + && mv mechanisms /home/test_user/atmospheric_physics/schemes/musica/configurations + # Must make ccpp-framework available before building test RUN cd atmospheric_physics/test \ && mkdir lib \ @@ -71,10 +77,4 @@ RUN cd atmospheric_physics/test \ -D CCPP_ENABLE_MEMCHECK=ON \ && cmake --build ./build -RUN cd atmospheric_physics/test \ - && mkdir third_party \ - && cd third_party \ - && git clone --depth 1 https://github.com/NCAR/tuv-x.git \ - && cp -r tuv-x/data ../build/data - WORKDIR /home/test_user/atmospheric_physics/test/build \ No newline at end of file diff --git a/test/musica/CMakeLists.txt b/test/musica/CMakeLists.txt index e41405c1..74b29693 100644 --- a/test/musica/CMakeLists.txt +++ b/test/musica/CMakeLists.txt @@ -56,14 +56,10 @@ add_test( add_memory_check_test(test_musica_api $ "" ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}) +# copy MUSICA configuration data into the build folder add_custom_target( copy_micm_configs ALL ${CMAKE_COMMAND} -E copy_directory - ${CMAKE_CURRENT_SOURCE_DIR}/micm/configs/chapman ${CMAKE_BINARY_DIR}/chapman -) - -add_custom_target( - copy_tuvx_configs ALL ${CMAKE_COMMAND} -E copy_directory - ${CMAKE_CURRENT_SOURCE_DIR}/tuvx/configs ${CMAKE_BINARY_DIR}/configs + ${CMAKE_CURRENT_SOURCE_DIR}/../../schemes/musica/configurations ${CMAKE_BINARY_DIR}/musica_configurations ) add_subdirectory(micm) diff --git a/test/musica/micm/configs/chapman/config.json b/test/musica/micm/configs/chapman/config.json deleted file mode 100644 index 04d0ef28..00000000 --- a/test/musica/micm/configs/chapman/config.json +++ /dev/null @@ -1,6 +0,0 @@ -{ - "camp-files": [ - "species.json", - "reactions.json" - ] -} diff --git a/test/musica/micm/configs/chapman/reactions.json b/test/musica/micm/configs/chapman/reactions.json deleted file mode 100644 index 32af8982..00000000 --- a/test/musica/micm/configs/chapman/reactions.json +++ /dev/null @@ -1,97 +0,0 @@ -{ - "camp-data": [ - { - "name": "Chapman", - "type": "MECHANISM", - "reactions": [ - { - "type": "PHOTOLYSIS", - "reactants": { - "O2": {} - }, - "products": { - "O": { - "yield": 2.0 - } - }, - "MUSICA name": "jO2" - }, - { - "type": "ARRHENIUS", - "A": 8.018e-17, - "reactants": { - "O": {}, - "O2": {} - }, - "products": { - "O3": {} - }, - "MUSICA name": "R2" - }, - { - "type": "PHOTOLYSIS", - "reactants": { - "O3": {} - }, - "products": { - "O": {}, - "O2": {} - }, - "MUSICA name": "jO3->O" - }, - { - "type": "ARRHENIUS", - "A": 1.576e-15, - "reactants": { - "O": {}, - "O3": {} - }, - "products": { - "O2": { - "yield": 2.0 - } - }, - "MUSICA name": "R4" - }, - { - "type": "PHOTOLYSIS", - "reactants": { - "O3": {} - }, - "products": { - "O1D": {}, - "O2": {} - }, - "MUSICA name": "jO3->O1D" - }, - { - "type": "ARRHENIUS", - "A": 7.11e-11, - "reactants": { - "O1D": {}, - "M": {} - }, - "products": { - "O": {}, - "M": {} - }, - "MUSICA name": "R6" - }, - { - "type": "ARRHENIUS", - "A": 1.2e-10, - "reactants": { - "O1D": {}, - "O3": {} - }, - "products": { - "O2": { - "yield": 2.0 - } - }, - "MUSICA name": "R7" - } - ] - } - ] -} \ No newline at end of file diff --git a/test/musica/micm/configs/chapman/species.json b/test/musica/micm/configs/chapman/species.json deleted file mode 100644 index aa2c2dcf..00000000 --- a/test/musica/micm/configs/chapman/species.json +++ /dev/null @@ -1,39 +0,0 @@ -{ - "camp-data": [ - { - "name": "M", - "type": "CHEM_SPEC", - "tracer type": "THIRD_BODY", - "molecular weight [kg mol-1]": 0.029, - "__is advected": false - }, - { - "name": "O2", - "type": "CHEM_SPEC", - "tracer type": "CONSTANT", - "molecular weight [kg mol-1]": 0.032, - "__is advected": false - }, - { - "name": "O", - "type": "CHEM_SPEC", - "absolute tolerance": 1e-12, - "molecular weight [kg mol-1]": 0.016, - "__is advected": false - }, - { - "name": "O1D", - "type": "CHEM_SPEC", - "absolute tolerance": 1e-12, - "molecular weight [kg mol-1]": 0.016, - "__is advected": false - }, - { - "name": "O3", - "type": "CHEM_SPEC", - "absolute tolerance": 1e-12, - "molecular weight [kg mol-1]": 0.048, - "__is advected": true - } - ] -} diff --git a/test/musica/micm/test_micm_util.F90 b/test/musica/micm/test_micm_util.F90 index 798f1408..3b5cf4fa 100644 --- a/test/musica/micm/test_micm_util.F90 +++ b/test/musica/micm/test_micm_util.F90 @@ -7,71 +7,10 @@ program test_micm_util #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_reshape() call test_unit_conversion() contains - subroutine test_reshape() - use iso_c_binding, only: c_double - use ccpp_kinds, only: kind_phys - - integer, parameter :: NUM_SPECIES = 4 - integer, parameter :: NUM_COLUMNS = 2 - integer, parameter :: NUM_LAYERS = 2 - 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) - real(kind_phys) :: constituents(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) - real(c_double) :: micm_temperature(NUM_GRID_CELLS) - 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) - 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 /) - pressure(:,1) = (/ 100._kind_phys, 200._kind_phys /) - pressure(:,2) = (/ 300._kind_phys, 400._kind_phys /) - dry_air_density(:,1) = (/ 100._kind_phys, 200._kind_phys /) - dry_air_density(:,2) = (/ 300._kind_phys, 400._kind_phys /) - constituents(1,1,:) = (/ 0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys, 0.4_kind_phys /) - 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 /) - - 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) == 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), expected_constituents(i_elem), ABS_ERROR) - end do - - call reshape_into_ccpp_arr(micm_constituents, constituents) - - i_arr = 1 - 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), expected_constituents(i_arr), ABS_ERROR) - i_arr = i_arr + 1 - end do - end do - end do - - end subroutine test_reshape - subroutine test_unit_conversion() use ccpp_kinds, only: kind_phys diff --git a/test/musica/musica_ccpp_namelist.F90 b/test/musica/musica_ccpp_namelist.F90 index 69c803b3..d27cc754 100644 --- a/test/musica/musica_ccpp_namelist.F90 +++ b/test/musica/musica_ccpp_namelist.F90 @@ -4,9 +4,9 @@ module musica_ccpp_namelist implicit none private - public :: filename_of_micm_configuration, filename_of_tuvx_configuration - - character(len=*), parameter :: filename_of_micm_configuration = 'chapman' - character(len=*), parameter :: filename_of_tuvx_configuration = 'configs/ts1_tsmlt.json' + + character(len=250), public :: filename_of_micm_configuration = 'musica_configurations/chapman/micm/config.json' + character(len=250), public :: filename_of_tuvx_configuration = 'musica_configurations/chapman/tuvx/config.json' + character(len=250), public :: filename_of_tuvx_micm_mapping_configuration = 'musica_configurations/chapman/tuvx_micm_mapping.json' end module musica_ccpp_namelist diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index 6b4d4154..61cca2f7 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -5,51 +5,25 @@ program run_test_musica_ccpp 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_musica_ccpp_api() + call test_chapman() + call test_terminator() contains - subroutine test_musica_ccpp_api() - use musica_micm, only: Rosenbrock, RosenbrockStandardOrder - use ccpp_kinds, only: kind_phys - use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t - use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + !> Returns the wavelength edges used in the CAM-Chem photolysis rate constant lookup table + !! 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. + subroutine get_wavelength_edges(edges) + use ccpp_kinds, only: kind_phys - integer, parameter :: NUM_SPECIES = 4 - integer, parameter :: NUM_COLUMNS = 2 - integer, parameter :: NUM_LAYERS = 2 - 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 = 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(:) - 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 = (/ & + real(kind_phys), dimension(:), intent(out) :: edges + + edges = (/ & 120.0e-9_kind_phys, & 121.4e-9_kind_phys, & 121.9e-9_kind_phys, & @@ -154,6 +128,60 @@ subroutine test_musica_ccpp_api() 700.0e-9_kind_phys, & 750.0e-9_kind_phys & /) + + end subroutine get_wavelength_edges + + !> Tests the Chapman chemistry scheme + subroutine test_chapman() + use musica_micm, only: Rosenbrock, RosenbrockStandardOrder + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_micm, only: micm + use musica_ccpp_namelist, only: filename_of_micm_configuration, & + filename_of_tuvx_configuration, & + filename_of_tuvx_micm_mapping_configuration + + implicit none + + integer, parameter :: NUM_SPECIES = 5 + ! This test requires that the number of grid cells = 4, which is the default + ! vector dimension for MICM. This restriction will be removed once + ! https://github.com/NCAR/musica/issues/217 is finished. + integer, parameter :: NUM_COLUMNS = 2 + integer, parameter :: NUM_LAYERS = 2 + 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 = 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 + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: initial_constituents ! kg kg-1 + type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:) + type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:) + type(ccpp_constituent_properties_t), pointer :: const_prop + real(kind_phys) :: molar_mass, base_conc + character(len=512) :: species_name, units + character(len=:), allocatable :: micm_species_name + logical :: tmp_bool, is_advected + integer :: i, j, k + integer :: N2_index, O2_index, O_index, O1D_index, O3_index + real(kind_phys) :: total_O, total_O_init + + call get_wavelength_edges(photolysis_wavelength_grid_interfaces) + solver_type = Rosenbrock + time_step = 60._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 /) @@ -168,12 +196,16 @@ subroutine test_musica_ccpp_api() pressure(:,2) = (/ 8000.04_kind_phys, 9000.04_kind_phys /) dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) - constituents(1,1,:) = (/ 0.1_kind_phys, 0.2_kind_phys, 0.3_kind_phys, 0.4_kind_phys /) - 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 /) - call musica_ccpp_register(solver_type, num_grid_cells, constituent_props, errmsg, errcode) + filename_of_micm_configuration = 'musica_configurations/chapman/micm/config.json' + filename_of_tuvx_configuration = 'musica_configurations/chapman/tuvx/config.json' + filename_of_tuvx_micm_mapping_configuration = 'musica_configurations/chapman/tuvx_micm_mapping.json' + + call musica_ccpp_register(solver_type, NUM_GRID_CELLS, constituent_props, errmsg, errcode) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif ASSERT(allocated(constituent_props)) ASSERT(size(constituent_props) == NUM_SPECIES) do i = 1, size(constituent_props) @@ -185,13 +217,213 @@ subroutine test_musica_ccpp_api() ASSERT(errcode == 0) call constituent_props(i)%is_advected(is_advected, errcode, errmsg) ASSERT(errcode == 0) - tmp_bool = (trim(species_name) == "O2" .and. molar_mass == 0.032_kind_phys .and. .not. is_advected) .or. & - (trim(species_name) == "O" .and. molar_mass == 0.016_kind_phys .and. .not. is_advected) .or. & - (trim(species_name) == "O1D" .and. molar_mass == 0.016_kind_phys .and. .not. is_advected) .or. & - (trim(species_name) == "O3" .and. molar_mass == 0.048_kind_phys .and. is_advected) + tmp_bool = (trim(species_name) == "O2" .and. molar_mass == 0.0319988_kind_phys .and. is_advected) .or. & + (trim(species_name) == "O" .and. molar_mass == 0.0159994_kind_phys .and. .not. is_advected) .or. & + (trim(species_name) == "O1D" .and. molar_mass == 0.0159994_kind_phys .and. .not. is_advected) .or. & + (trim(species_name) == "O3" .and. molar_mass == 0.0479982_kind_phys .and. is_advected) .or. & + (trim(species_name) == "N2" .and. molar_mass == 0.0280134_kind_phys .and. is_advected) ASSERT(tmp_bool) call constituent_props(i)%units(units, errcode, errmsg) + if (errcode /= 0) then + write(*,*) errcode, trim(errmsg) + stop 3 + endif + ASSERT(trim(units) == 'kg kg-1') + end do + if (errcode /= 0) then + write(*,*) errcode, trim(errmsg) + stop 3 + end if + + allocate(constituent_props_ptr(size(constituent_props))) + do i = 1, size(constituent_props) + const_prop => constituent_props(i) + call constituent_props_ptr(i)%set(const_prop, errcode, errmsg) + end do + + call musica_ccpp_init(NUM_LAYERS, NUM_LAYERS+1, photolysis_wavelength_grid_interfaces, & + errmsg, errcode) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + + do i = 1, micm%species_ordering%size() + micm_species_name = micm%species_ordering%name(i) + if (micm_species_name == "O2") then + O2_index = i + base_conc = 0.21_kind_phys + else if (micm_species_name == "O") then + O_index = i + base_conc = 1.0e-9_kind_phys + else if (micm_species_name == "O1D") then + O1D_index = i + base_conc = 1.0e-9_kind_phys + else if (micm_species_name == "O3") then + O3_index = i + base_conc = 1.0e-4_kind_phys + else if (micm_species_name == "N2") then + N2_index = i + base_conc = 0.79_kind_phys + else + write(*,*) "Unknown species: ", micm_species_name + stop 3 + endif + do j = 1, NUM_COLUMNS + do k = 1, NUM_LAYERS + constituents(j,k,i) = base_conc * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) + end do + end do + end do + initial_constituents(:,:,:) = constituents(:,:,:) + + write(*,*) "[MUSICA INFO] Initial Time Step" + write(*,fmt="(1x,f10.2)") time_step + write(*,*) "[MUSICA INFO] Initial Temperature" + write(*,fmt="(4(1x,f10.4))") temperature + write(*,*) "[MUSICA INFO] Initial Pressure" + write(*,fmt="(4(1x,f10.4))") pressure + 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, surface_albedo, standard_gravitational_acceleration, & + errmsg, errcode) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + + write(*,*) "[MUSICA INFO] Solved Concentrations" + write(*,fmt="(4(3x,e13.6))") constituents + + call musica_ccpp_final(errmsg, errcode) + + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + + ! Check the results + do i = 1, NUM_COLUMNS + do j = 1, NUM_LAYERS + ! N2 should be unchanged + ASSERT_NEAR(constituents(i,j,N2_index), initial_constituents(i,j,N2_index), 1.0e-13) + ! O3 and O2 should be relatively unchanged + ASSERT_NEAR(constituents(i,j,O3_index), initial_constituents(i,j,O3_index), 1.0e-4) + ASSERT_NEAR(constituents(i,j,O2_index), initial_constituents(i,j,O2_index), 1.0e-4) + ! O and O1D should be < 1e-10 kg kg-1 + ASSERT(constituents(i,j,O_index) < 1.0e-10) + ASSERT(constituents(i,j,O1D_index) < 1.0e-10) + ! total O mass should be conserved + total_O = constituents(i,j,O_index) + constituents(i,j,O1D_index) + & + constituents(i,j,O2_index) + constituents(i,j,O3_index) + total_O_init = initial_constituents(i,j,O_index) + initial_constituents(i,j,O1D_index) + & + initial_constituents(i,j,O2_index) + initial_constituents(i,j,O3_index) + ASSERT_NEAR(total_O, total_O_init, 1.0e-13) + end do + end do + + deallocate(constituent_props_ptr) + + end subroutine test_chapman + + !> Tests the simple Terminator chemistry scheme + subroutine test_terminator() + use musica_micm, only: Rosenbrock, RosenbrockStandardOrder + use ccpp_kinds, only: kind_phys + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t + use musica_ccpp_micm, only: micm + use musica_ccpp_namelist, only: filename_of_micm_configuration, & + filename_of_tuvx_configuration, & + filename_of_tuvx_micm_mapping_configuration + + implicit none + + integer, parameter :: NUM_SPECIES = 2 + ! This test requires that the number of grid cells = 4, which is the default + ! vector dimension for MICM. This restriction will be removed once + ! https://github.com/NCAR/musica/issues/217 is finished. + integer, parameter :: NUM_COLUMNS = 2 + integer, parameter :: NUM_LAYERS = 2 + 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 = 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 + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS,NUM_SPECIES) :: initial_constituents ! kg kg-1 + type(ccpp_constituent_prop_ptr_t), allocatable :: constituent_props_ptr(:) + type(ccpp_constituent_properties_t), allocatable, target :: constituent_props(:) + type(ccpp_constituent_properties_t), pointer :: const_prop + real(kind_phys) :: molar_mass, base_conc + character(len=512) :: species_name, units + character(len=:), allocatable :: micm_species_name + logical :: tmp_bool, is_advected + integer :: i, j, k + integer :: Cl_index, Cl2_index + real(kind_phys) :: total_Cl, total_Cl_init + + call get_wavelength_edges(photolysis_wavelength_grid_interfaces) + solver_type = Rosenbrock + time_step = 60._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 /) + pressure(:,1) = (/ 6000.04_kind_phys, 7000.04_kind_phys /) + pressure(:,2) = (/ 8000.04_kind_phys, 9000.04_kind_phys /) + dry_air_density(:,1) = (/ 3.5_kind_phys, 4.5_kind_phys /) + dry_air_density(:,2) = (/ 5.5_kind_phys, 6.5_kind_phys /) + + filename_of_micm_configuration = 'musica_configurations/terminator/micm/config.json' + filename_of_tuvx_configuration = 'musica_configurations/terminator/tuvx/config.json' + filename_of_tuvx_micm_mapping_configuration = 'musica_configurations/terminator/tuvx_micm_mapping.json' + + call musica_ccpp_register(solver_type, NUM_GRID_CELLS, constituent_props, errmsg, errcode) + if (errcode /= 0) then + write(*,*) trim(errmsg) + stop 3 + endif + ASSERT(allocated(constituent_props)) + ASSERT(size(constituent_props) == NUM_SPECIES) + do i = 1, size(constituent_props) + ASSERT(constituent_props(i)%is_instantiated(errcode, errmsg)) + ASSERT(errcode == 0) + call constituent_props(i)%standard_name(species_name, errcode, errmsg) + ASSERT(errcode == 0) + call constituent_props(i)%molar_mass(molar_mass, errcode, errmsg) ASSERT(errcode == 0) + call constituent_props(i)%is_advected(is_advected, errcode, errmsg) + ASSERT(errcode == 0) + tmp_bool = (trim(species_name) == "Cl" .and. molar_mass == 0.035453_kind_phys .and. is_advected) .or. & + (trim(species_name) == "Cl2" .and. molar_mass == 0.070906_kind_phys .and. is_advected) + ASSERT(tmp_bool) + call constituent_props(i)%units(units, errcode, errmsg) + if (errcode /= 0) then + write(*,*) errcode, trim(errmsg) + stop 3 + endif ASSERT(trim(units) == 'kg kg-1') end do if (errcode /= 0) then @@ -212,6 +444,26 @@ subroutine test_musica_ccpp_api() stop 3 endif + do i = 1, micm%species_ordering%size() + micm_species_name = micm%species_ordering%name(i) + if (micm_species_name == "Cl") then + Cl_index = i + base_conc = 1.0e-10_kind_phys + else if (micm_species_name == "Cl2") then + Cl2_index = i + base_conc = 1.0e-6_kind_phys + else + write(*,*) "Unknown species: ", micm_species_name + stop 3 + endif + do j = 1, NUM_COLUMNS + do k = 1, NUM_LAYERS + constituents(j,k,i) = base_conc * (1.0 + 0.1 * (j-1) + 0.01 * (k-1)) + end do + end do + end do + initial_constituents(:,:,:) = constituents(:,:,:) + write(*,*) "[MUSICA INFO] Initial Time Step" write(*,fmt="(1x,f10.2)") time_step write(*,*) "[MUSICA INFO] Initial Temperature" @@ -241,8 +493,22 @@ subroutine test_musica_ccpp_api() stop 3 endif + ! Check the results + do i = 1, NUM_COLUMNS + do j = 1, NUM_LAYERS + ! Cl2 should not be unchanged + ASSERT(constituents(i,j,Cl2_index) .ne. initial_constituents(i,j,Cl2_index)) + ! Cl should not be unchanged + ASSERT(constituents(i,j,Cl_index) .ne. initial_constituents(i,j,Cl_index)) + ! total Cl mass should be conserved + total_Cl = constituents(i,j,Cl_index) + constituents(i,j,Cl2_index) + total_Cl_init = initial_constituents(i,j,Cl_index) + initial_constituents(i,j,Cl2_index) + ASSERT_NEAR(total_Cl, total_Cl_init, 1.0e-13) + end do + end do + deallocate(constituent_props_ptr) - end subroutine test_musica_ccpp_api + end subroutine test_terminator end program run_test_musica_ccpp \ No newline at end of file From 79df94abe0283f6717defaf9368bfb7e8eaff4f2 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 18 Nov 2024 15:51:47 -0500 Subject: [PATCH 5/5] Implement CCPPized check_energy_chng and check_energy_fix (#145) Fixes #114 * Companion PR in CAM: https://github.com/ESCOMP/CAM/pull/1180 * Companion PR in CAM-SIMA: https://github.com/ESCOMP/CAM-SIMA/pull/316 This PR fixes the following NCAR/atmospheric_physics Github issues: #114 - Implements check_energy_chng. The routine computes total energies using physics and dycore formula and takes in boundary fluxes of vapor, liquid+ice, ice, sensible heat, and writes to log significant energy conservation errors. In order to take in the scheme name and fluxes from the last calling physics scheme (usually zero) a check_energy_zero_fluxes has to be ran before the scheme, then the physics scheme to be checked, then check_energy_chng. - Implements check_energy_fix. The routine computes the heating rate required for global mean total energy conservation which is later applied by apply_heating_rate. Supporting routines include the global mean calculator for total energy (check_energy_gmean) - stored in separate folder to prevent CAM from building it. The global means are very useful for diagnosing model state (as the computed energy numbers include all model state incl. constituents) through a simple one-line printout. check_energy_gmean_diagnostics implements this. At the end of the timestep check_energy_save_teout has to be ran in order to save total energies at the end of the timestep for use in the next timestep. - Diagnostics of intermediate quantities used in check_energy_diagnostics. List all existing files that have been added (A), modified (M), or deleted (D), and describe the changes: ``` - Docs: M doc/ChangeLog - check_energy_chng and supporting routines: A schemes/check_energy/check_energy_chng.F90 A schemes/check_energy/check_energy_chng.meta A schemes/check_energy/check_energy_chng_namelist.xml A schemes/check_energy/check_energy_scaling.F90 A schemes/check_energy/check_energy_scaling.meta A schemes/check_energy/check_energy_zero_fluxes.F90 A schemes/check_energy/check_energy_zero_fluxes.meta - check_energy_fix and supporting routines: A schemes/check_energy/check_energy_fix.F90 A schemes/check_energy/check_energy_fix.meta A schemes/check_energy/check_energy_gmean/check_energy_gmean.F90 A schemes/check_energy/check_energy_gmean/check_energy_gmean.meta A schemes/check_energy/check_energy_save_teout.F90 A schemes/check_energy/check_energy_save_teout.meta - check_energy related diagnostics: A schemes/sima_diagnostics/check_energy_diagnostics.F90 A schemes/sima_diagnostics/check_energy_diagnostics.meta - check_energy_gmean "nstep, te" atm.log output: A schemes/sima_diagnostics/check_energy_gmean_diagnostics.F90 A schemes/sima_diagnostics/check_energy_gmean_diagnostics.meta - test SDF including all functionality: A test/test_suites/suite_check_energy.xml ``` List and Describe any test failures: N/A Summarize any changes to answers: none --- doc/ChangeLog | 72 +++ doc/NamesNotInDictionary.txt | 289 ++++++++---- schemes/check_energy/check_energy_chng.F90 | 414 ++++++++++++++++++ schemes/check_energy/check_energy_chng.meta | 413 +++++++++++++++++ .../check_energy_chng_namelist.xml | 19 + schemes/check_energy/check_energy_fix.F90 | 42 ++ schemes/check_energy/check_energy_fix.meta | 55 +++ .../check_energy_gmean/check_energy_gmean.F90 | 70 +++ .../check_energy_gmean.meta | 86 ++++ .../check_energy/check_energy_save_teout.F90 | 32 ++ .../check_energy/check_energy_save_teout.meta | 25 ++ schemes/check_energy/check_energy_scaling.F90 | 42 ++ .../check_energy/check_energy_scaling.meta | 43 ++ .../check_energy/check_energy_zero_fluxes.F90 | 34 ++ .../check_energy_zero_fluxes.meta | 43 ++ .../dycore_energy_consistency_adjust.F90 | 46 ++ .../dycore_energy_consistency_adjust.meta | 43 ++ .../check_energy_diagnostics.F90 | 86 ++++ .../check_energy_diagnostics.meta | 83 ++++ .../check_energy_gmean_diagnostics.F90 | 61 +++ .../check_energy_gmean_diagnostics.meta | 59 +++ ...heck_energy_gmean_diagnostics_namelist.xml | 19 + suites/suite_adiabatic.xml | 30 ++ suites/suite_cam7.xml | 31 +- suites/suite_kessler.xml | 21 + suites/suite_tj2016.xml | 22 + 26 files changed, 2091 insertions(+), 89 deletions(-) create mode 100644 schemes/check_energy/check_energy_chng.F90 create mode 100644 schemes/check_energy/check_energy_chng.meta create mode 100644 schemes/check_energy/check_energy_chng_namelist.xml create mode 100644 schemes/check_energy/check_energy_fix.F90 create mode 100644 schemes/check_energy/check_energy_fix.meta create mode 100644 schemes/check_energy/check_energy_gmean/check_energy_gmean.F90 create mode 100644 schemes/check_energy/check_energy_gmean/check_energy_gmean.meta create mode 100644 schemes/check_energy/check_energy_save_teout.F90 create mode 100644 schemes/check_energy/check_energy_save_teout.meta create mode 100644 schemes/check_energy/check_energy_scaling.F90 create mode 100644 schemes/check_energy/check_energy_scaling.meta create mode 100644 schemes/check_energy/check_energy_zero_fluxes.F90 create mode 100644 schemes/check_energy/check_energy_zero_fluxes.meta create mode 100644 schemes/check_energy/dycore_energy_consistency_adjust.F90 create mode 100644 schemes/check_energy/dycore_energy_consistency_adjust.meta create mode 100644 schemes/sima_diagnostics/check_energy_diagnostics.F90 create mode 100644 schemes/sima_diagnostics/check_energy_diagnostics.meta create mode 100644 schemes/sima_diagnostics/check_energy_gmean_diagnostics.F90 create mode 100644 schemes/sima_diagnostics/check_energy_gmean_diagnostics.meta create mode 100644 schemes/sima_diagnostics/check_energy_gmean_diagnostics_namelist.xml create mode 100644 suites/suite_adiabatic.xml diff --git a/doc/ChangeLog b/doc/ChangeLog index 721b9f0d..b08a0987 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,5 +1,77 @@ =============================================================== +Tag name: atmos_phys0_07_000 +Originator(s): jimmielin +Date: November 18, 2024 +One-line Summary: Implement CCPPized check_energy_chng and check_energy_fix +Github PR URL: + +This PR fixes the following NCAR/atmospheric_physics Github issues: #114 + +- Implements check_energy_chng. The routine computes total energies using physics and dycore formula and takes in boundary fluxes of vapor, liquid+ice, ice, sensible heat, and writes to log significant energy conservation errors. +In order to take in the scheme name and fluxes from the last calling physics scheme (usually zero) a check_energy_zero_fluxes has to be ran before the scheme, then the physics scheme to be checked, then check_energy_chng. + +- Implements check_energy_fix. The routine computes the heating rate required for global mean total energy conservation which is later applied by apply_heating_rate. +Supporting routines include the global mean calculator for total energy (check_energy_gmean) - stored in separate folder to prevent CAM from building it. +The global means are very useful for diagnosing model state (as the computed energy numbers include all model state incl. constituents) through a simple one-line printout. check_energy_gmean_diagnostics implements this. +At the end of the timestep check_energy_save_teout has to be ran in order to save total energies at the end of the timestep for use in the next timestep. + +- Diagnostics of intermediate quantities used in check_energy_diagnostics. + +- Implements dycore_energy_consistency_adjust which adjusts temperature and temperature tendencies for the MPAS and SE dynamical cores. + +Code reviewed by: + +List all existing files that have been added (A), modified (M), or deleted (D), +and describe the changes: + +- Docs: +M doc/ChangeLog + +- check_energy_chng and supporting routines: +A schemes/check_energy/check_energy_chng.F90 +A schemes/check_energy/check_energy_chng.meta +A schemes/check_energy/check_energy_chng_namelist.xml +A schemes/check_energy/check_energy_scaling.F90 +A schemes/check_energy/check_energy_scaling.meta +A schemes/check_energy/check_energy_zero_fluxes.F90 +A schemes/check_energy/check_energy_zero_fluxes.meta + +- check_energy_fix and supporting routines: +A schemes/check_energy/check_energy_fix.F90 +A schemes/check_energy/check_energy_fix.meta +A schemes/check_energy/check_energy_gmean/check_energy_gmean.F90 +A schemes/check_energy/check_energy_gmean/check_energy_gmean.meta +A schemes/check_energy/check_energy_save_teout.F90 +A schemes/check_energy/check_energy_save_teout.meta + +- check_energy related diagnostics: +A schemes/sima_diagnostics/check_energy_diagnostics.F90 +A schemes/sima_diagnostics/check_energy_diagnostics.meta + +- check_energy_gmean "nstep, te" atm.log output: +A schemes/sima_diagnostics/check_energy_gmean_diagnostics.F90 +A schemes/sima_diagnostics/check_energy_gmean_diagnostics.meta + +- dycore_energy_consistency_adjust and applications in simple physics: +A schemes/check_energy/dycore_energy_consistency_adjust.F90 +A schemes/check_energy/dycore_energy_consistency_adjust.meta +M suites/suite_held_suarez_1994.xml +M suites/suite_kessler.xml + +- add check_energy to SDFs: +M suites/suite_cam7.xml + +- adiabatic SDF: +A suites/suite_adiabatic.xml + + +List and Describe any test failures: N/A + +Summarize any changes to answers: none + +=============================================================== + Tag name: Originator(s): jimmielin Date: October 17, 2024 diff --git a/doc/NamesNotInDictionary.txt b/doc/NamesNotInDictionary.txt index c934e247..57566b4d 100644 --- a/doc/NamesNotInDictionary.txt +++ b/doc/NamesNotInDictionary.txt @@ -1,19 +1,42 @@ ####################### Date/time of when script was run: -2024-10-10 14:17:36.423628 +2024-11-18 15:44:54.993446 ####################### Non-dictionary standard names found in the following metadata files: -------------------------- +atmospheric_physics/schemes/sima_diagnostics/check_energy_gmean_diagnostics.meta + + - flag_for_energy_global_means_output + - global_mean_heating_rate_correction_for_energy_conservation + - global_mean_vertically_integrated_total_energy_at_end_of_physics_timestep + - global_mean_vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep + +-------------------------- + +atmospheric_physics/schemes/sima_diagnostics/check_energy_diagnostics.meta + + - cumulative_total_energy_boundary_flux_using_physics_energy_formula + - cumulative_total_water_boundary_flux + - ratio_of_specific_heat_of_air_used_in_physics_energy_formula_to_specific_heat_of_air_used_in_dycore_energy_formula + - specific_heat_of_air_used_in_dycore + - vertically_integrated_total_energy_at_end_of_physics_timestep + - vertically_integrated_total_energy_using_dycore_energy_formula + - vertically_integrated_total_energy_using_physics_energy_formula + - vertically_integrated_total_water + +-------------------------- + atmospheric_physics/schemes/sima_diagnostics/sima_state_diagnostics.meta - air_pressure_at_interface - air_pressure_of_dry_air_at_interface - ln_air_pressure_at_interface - ln_air_pressure_of_dry_air_at_interface + - surface_air_pressure -------------------------- @@ -45,65 +68,49 @@ atmospheric_physics/schemes/sima_diagnostics/tropopause_diagnostics.meta -------------------------- -atmospheric_physics/schemes/dry_adiabatic_adjust/dadadj.meta +atmospheric_physics/schemes/tj2016/tj2016_precip.meta - - air_pressure_at_interface - - binary_indicator_for_dry_adiabatic_adjusted_grid_cell - - number_of_iterations_for_dry_adiabatic_adjustment_algorithm_convergence - - number_of_vertical_levels_from_model_top_where_dry_adiabatic_adjustment_occurs - - tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + - gas_constant_of_water_vapor + - lwe_large_scale_precipitation_rate_at_surface + - ratio_of_water_vapor_to_dry_air_molecular_weights + - sum_of_sigma_pressure_hybrid_coordinate_a_coefficient_and_sigma_pressure_hybrid_coordinate_b_coefficient -------------------------- -atmospheric_physics/schemes/dry_adiabatic_adjust/dadadj_apply_qv_tendency.meta +atmospheric_physics/schemes/tj2016/tj2016_sfc_pbl_hs.meta - - tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water + - air_pressure_at_interface + - eddy_heat_diffusivity + - eddy_momentum_diffusivity + - gas_constant_of_water_vapor + - ln_air_pressure_at_interface + - pi_constant + - ratio_of_water_vapor_to_dry_air_molecular_weights + - sum_of_sigma_pressure_hybrid_coordinate_a_coefficient_and_sigma_pressure_hybrid_coordinate_b_coefficient + - surface_air_pressure + - surface_eastward_wind_stress + - surface_evaporation_rate + - surface_northward_wind_stress + - surface_upward_sensible_heat_flux + - tendency_of_air_temperature_due_to_diabatic_heating + - tendency_of_air_temperature_due_to_vertical_diffusion + - tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_due_to_vertical_diffusion -------------------------- -atmospheric_physics/schemes/utilities/geopotential_temp.meta +atmospheric_physics/schemes/dry_adiabatic_adjust/dadadj.meta - air_pressure_at_interface - - ln_air_pressure_at_interface + - binary_indicator_for_dry_adiabatic_adjusted_grid_cell + - number_of_iterations_for_dry_adiabatic_adjustment_algorithm_convergence + - number_of_vertical_levels_from_model_top_where_dry_adiabatic_adjustment_occurs + - tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water -------------------------- -atmospheric_physics/schemes/tropopause_find/tropopause_find.meta +atmospheric_physics/schemes/dry_adiabatic_adjust/dadadj_apply_qv_tendency.meta - - air_pressure_at_interface - - fill_value_for_diagnostic_output - - fractional_calendar_days_on_end_of_current_timestep - - pi_constant - - ratio_of_dry_air_gas_constant_to_specific_heat_of_dry_air_at_constant_pressure - - tropopause_air_pressure - - tropopause_air_pressure_from_chemical_method - - tropopause_air_pressure_from_climatological_method - - tropopause_air_pressure_from_climatology_dataset - - tropopause_air_pressure_from_cold_point_method - - tropopause_air_pressure_from_hybrid_stobie_linoz_with_climatological_backup_method - - tropopause_air_pressure_from_lapse_rate_method - - tropopause_air_temperature - - tropopause_air_temperature_from_chemical_method - - tropopause_air_temperature_from_climatological_method - - tropopause_air_temperature_from_cold_point_method - - tropopause_air_temperature_from_hybrid_stobie_linoz_with_climatological_backup_method - - tropopause_air_temperature_from_lapse_rate_method - - tropopause_calendar_days_from_climatology - - tropopause_geopotential_height_wrt_surface - - tropopause_geopotential_height_wrt_surface_from_chemical_method - - tropopause_geopotential_height_wrt_surface_from_climatological_method - - tropopause_geopotential_height_wrt_surface_from_cold_point_method - - tropopause_geopotential_height_wrt_surface_from_hybrid_stobie_linoz_with_climatological_backup_method - - tropopause_geopotential_height_wrt_surface_from_lapse_rate_method - - tropopause_vertical_layer_index - - tropopause_vertical_layer_index_from_chemical_method - - tropopause_vertical_layer_index_from_climatological_method - - tropopause_vertical_layer_index_from_cold_point_method - - tropopause_vertical_layer_index_from_hybrid_stobie_linoz_with_climatological_backup_method - - tropopause_vertical_layer_index_from_hybrid_stobie_linoz_with_climatological_backup_method_for_chemistry - - tropopause_vertical_layer_index_from_lapse_rate_method - - vertical_layer_index_lower_bound_from_hybrid_stobie_linoz_with_climatological_backup_method_for_linearized_ozone_chemistry - - vertical_layer_index_lower_bound_from_hybrid_stobie_linoz_with_climatological_backup_method_for_stratospheric_chemistry + - tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water -------------------------- @@ -134,6 +141,32 @@ atmospheric_physics/schemes/zhang_mcfarlane/zm_conv_momtran.meta -------------------------- +atmospheric_physics/schemes/zhang_mcfarlane/zm_conv_evap.meta + + - + - cloud_area_fraction + - flag_for_zhang_mcfarlane_convective_organization_parameterization? + - freezing_point_of_water? + - frozen_precipitation_mass_flux_at_interface_due_to_deep_convection? + - heating_rate + - latent_heat_of_fusion_of_water_at_0c? + - latent_heat_of_vaporization_of_water_at_0c? + - lwe_frozen_precipitation_rate_at_surface_due_to_deep_convection + - lwe_precipitation_rate_at_surface_due_to_deep_convection + - precipitation_mass_flux_at_interface_due_to_deep_convection? + - pressure_thickness + - specific_heat_of_dry_air_at_constant_pressure? + - tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_frozen_precipitation_melt? + - tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_frozen_precipitation_production_in_deep_convection? + - tendency_of_frozen_precipitation_wrt_moist_air_and_condensed_water_due_to_deep_convection? + - tendency_of_precipitation_wrt_moist_air_and_condensed_water_due_to_deep_convection? + - tendency_of_precipitation_wrt_moist_air_and_condensed_water_due_to_deep_convection_excluding_subcloud_evaporation + - tendency_of_water_vapor_mixing_ratio_wrt_moist_air and_condensed_water? + - tunable_evaporation_efficiency_for_land_in_zhang_mcfarlane_deep_convection_scheme? + - tunable_evaporation_efficiency_in_zhang_mcfarlane_deep_convection_scheme? + +-------------------------- + atmospheric_physics/schemes/zhang_mcfarlane/zm_convr.meta - air_pressure_at_interface @@ -212,57 +245,139 @@ atmospheric_physics/schemes/zhang_mcfarlane/zm_conv_convtran.meta -------------------------- -atmospheric_physics/schemes/zhang_mcfarlane/zm_conv_evap.meta +atmospheric_physics/schemes/utilities/geopotential_temp.meta - - - - cloud_area_fraction - - flag_for_zhang_mcfarlane_convective_organization_parameterization? - - freezing_point_of_water? - - frozen_precipitation_mass_flux_at_interface_due_to_deep_convection? - - heating_rate - - latent_heat_of_fusion_of_water_at_0c? - - latent_heat_of_vaporization_of_water_at_0c? - - lwe_frozen_precipitation_rate_at_surface_due_to_deep_convection - - lwe_precipitation_rate_at_surface_due_to_deep_convection - - precipitation_mass_flux_at_interface_due_to_deep_convection? - - pressure_thickness - - specific_heat_of_dry_air_at_constant_pressure? - - tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_frozen_precipitation_melt? - - tendency_of_dry_air_enthalpy_at_constant_pressure_due_to_frozen_precipitation_production_in_deep_convection? - - tendency_of_frozen_precipitation_wrt_moist_air_and_condensed_water_due_to_deep_convection? - - tendency_of_precipitation_wrt_moist_air_and_condensed_water_due_to_deep_convection? - - tendency_of_precipitation_wrt_moist_air_and_condensed_water_due_to_deep_convection_excluding_subcloud_evaporation - - tendency_of_water_vapor_mixing_ratio_wrt_moist_air and_condensed_water? - - tunable_evaporation_efficiency_for_land_in_zhang_mcfarlane_deep_convection_scheme? - - tunable_evaporation_efficiency_in_zhang_mcfarlane_deep_convection_scheme? + - air_pressure_at_interface + - ln_air_pressure_at_interface -------------------------- -atmospheric_physics/schemes/tj2016/tj2016_precip.meta +atmospheric_physics/schemes/check_energy/check_energy_save_teout.meta - - gas_constant_of_water_vapor - - lwe_large_scale_precipitation_rate_at_surface - - ratio_of_water_vapor_to_dry_air_molecular_weights - - sum_of_sigma_pressure_hybrid_coordinate_a_coefficient_and_sigma_pressure_hybrid_coordinate_b_coefficient + - vertically_integrated_total_energy_at_end_of_physics_timestep + - vertically_integrated_total_energy_using_dycore_energy_formula -------------------------- -atmospheric_physics/schemes/tj2016/tj2016_sfc_pbl_hs.meta +atmospheric_physics/schemes/check_energy/check_energy_chng.meta + + - air_pressure_of_dry_air_at_interface + - air_temperature_at_start_of_physics_timestep + - cumulative_total_energy_boundary_flux_using_physics_energy_formula + - cumulative_total_water_boundary_flux + - flag_for_energy_conservation_warning + - geopotential_height_wrt_surface_at_start_of_physics_timestep + - latent_heat_of_fusion_of_water_at_0c + - net_liquid_and_lwe_ice_fluxes_through_top_and_bottom_of_atmosphere_column + - net_lwe_ice_fluxes_through_top_and_bottom_of_atmosphere_column + - net_sensible_heat_flux_through_top_and_bottom_of_atmosphere_column + - net_water_vapor_fluxes_through_top_and_bottom_of_atmosphere_column + - number_of_atmosphere_columns_with_significant_energy_or_water_imbalances + - ratio_of_specific_heat_of_air_used_in_physics_energy_formula_to_specific_heat_of_air_used_in_dycore_energy_formula + - specific_heat_of_air_used_in_dycore + - total_energy_formula_for_dycore + - total_energy_formula_for_physics + - vertically_integrated_total_energy_at_end_of_physics_timestep + - vertically_integrated_total_energy_using_dycore_energy_formula + - vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep + - vertically_integrated_total_energy_using_physics_energy_formula + - vertically_integrated_total_energy_using_physics_energy_formula_at_start_of_physics_timestep + - vertically_integrated_total_water + - vertically_integrated_total_water_at_start_of_physics_timestep + +-------------------------- + +atmospheric_physics/schemes/check_energy/check_energy_zero_fluxes.meta + + - net_liquid_and_lwe_ice_fluxes_through_top_and_bottom_of_atmosphere_column + - net_lwe_ice_fluxes_through_top_and_bottom_of_atmosphere_column + - net_sensible_heat_flux_through_top_and_bottom_of_atmosphere_column + - net_water_vapor_fluxes_through_top_and_bottom_of_atmosphere_column + +-------------------------- + +atmospheric_physics/schemes/check_energy/dycore_energy_consistency_adjust.meta + + - flag_for_dycore_energy_consistency_adjustment + - ratio_of_specific_heat_of_air_used_in_physics_energy_formula_to_specific_heat_of_air_used_in_dycore_energy_formula + +-------------------------- + +atmospheric_physics/schemes/check_energy/check_energy_scaling.meta + + - ratio_of_specific_heat_of_air_used_in_physics_energy_formula_to_specific_heat_of_air_used_in_dycore_energy_formula + - specific_heat_of_air_used_in_dycore + +-------------------------- + +atmospheric_physics/schemes/check_energy/check_energy_fix.meta - air_pressure_at_interface - - eddy_heat_diffusivity - - eddy_momentum_diffusivity - - gas_constant_of_water_vapor - - ln_air_pressure_at_interface + - global_mean_heating_rate_correction_for_energy_conservation + - net_sensible_heat_flux_through_top_and_bottom_of_atmosphere_column + +-------------------------- + +atmospheric_physics/schemes/check_energy/check_energy_gmean/check_energy_gmean.meta + + - air_pressure_at_interface + - global_mean_air_pressure_at_top_of_atmosphere_model + - global_mean_heating_rate_correction_for_energy_conservation + - global_mean_surface_air_pressure + - global_mean_total_energy_correction_for_energy_conservation + - global_mean_vertically_integrated_total_energy_at_end_of_physics_timestep + - global_mean_vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep + - vertically_integrated_total_energy_at_end_of_physics_timestep + - vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep + +-------------------------- + +atmospheric_physics/schemes/tropopause_find/tropopause_find.meta + + - air_pressure_at_interface + - fill_value_for_diagnostic_output + - fractional_calendar_days_on_end_of_current_timestep - pi_constant - - ratio_of_water_vapor_to_dry_air_molecular_weights - - sum_of_sigma_pressure_hybrid_coordinate_a_coefficient_and_sigma_pressure_hybrid_coordinate_b_coefficient - - surface_eastward_wind_stress - - surface_evaporation_rate - - surface_northward_wind_stress - - surface_upward_sensible_heat_flux - - tendency_of_air_temperature_due_to_diabatic_heating - - tendency_of_air_temperature_due_to_vertical_diffusion - - tendency_of_water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_due_to_vertical_diffusion + - ratio_of_dry_air_gas_constant_to_specific_heat_of_dry_air_at_constant_pressure + - tropopause_air_pressure + - tropopause_air_pressure_from_chemical_method + - tropopause_air_pressure_from_climatological_method + - tropopause_air_pressure_from_climatology_dataset + - tropopause_air_pressure_from_cold_point_method + - tropopause_air_pressure_from_hybrid_stobie_linoz_with_climatological_backup_method + - tropopause_air_pressure_from_lapse_rate_method + - tropopause_air_temperature + - tropopause_air_temperature_from_chemical_method + - tropopause_air_temperature_from_climatological_method + - tropopause_air_temperature_from_cold_point_method + - tropopause_air_temperature_from_hybrid_stobie_linoz_with_climatological_backup_method + - tropopause_air_temperature_from_lapse_rate_method + - tropopause_calendar_days_from_climatology + - tropopause_geopotential_height_wrt_surface + - tropopause_geopotential_height_wrt_surface_from_chemical_method + - tropopause_geopotential_height_wrt_surface_from_climatological_method + - tropopause_geopotential_height_wrt_surface_from_cold_point_method + - tropopause_geopotential_height_wrt_surface_from_hybrid_stobie_linoz_with_climatological_backup_method + - tropopause_geopotential_height_wrt_surface_from_lapse_rate_method + - tropopause_vertical_layer_index + - tropopause_vertical_layer_index_from_chemical_method + - tropopause_vertical_layer_index_from_climatological_method + - tropopause_vertical_layer_index_from_cold_point_method + - tropopause_vertical_layer_index_from_hybrid_stobie_linoz_with_climatological_backup_method + - tropopause_vertical_layer_index_from_hybrid_stobie_linoz_with_climatological_backup_method_for_chemistry + - tropopause_vertical_layer_index_from_lapse_rate_method + - vertical_layer_index_lower_bound_from_hybrid_stobie_linoz_with_climatological_backup_method_for_linearized_ozone_chemistry + - vertical_layer_index_lower_bound_from_hybrid_stobie_linoz_with_climatological_backup_method_for_stratospheric_chemistry + +-------------------------- + +atmospheric_physics/schemes/musica/musica_ccpp.meta + + - blackbody_temperature_at_surface + - dynamic_constituents_for_musica_ccpp + - micm_solver_type + - number_of_grid_cells + - photolysis_wavelength_grid_interfaces + - surface_albedo_due_to_UV_and_VIS_direct ####################### diff --git a/schemes/check_energy/check_energy_chng.F90 b/schemes/check_energy/check_energy_chng.F90 new file mode 100644 index 00000000..91af5151 --- /dev/null +++ b/schemes/check_energy/check_energy_chng.F90 @@ -0,0 +1,414 @@ +module check_energy_chng + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: check_energy_chng_init + public :: check_energy_chng_timestep_init + public :: check_energy_chng_run + + ! Private module options. + logical :: print_energy_errors = .false. ! Turn on verbose output identifying columns that fail + ! energy/water checks? + +contains + +!> \section arg_table_check_energy_chng_init Argument Table +!! \htmlinclude arg_table_check_energy_chng_init.html + subroutine check_energy_chng_init(print_energy_errors_in) + ! Input arguments + logical, intent(in) :: print_energy_errors_in + + print_energy_errors = print_energy_errors_in + end subroutine check_energy_chng_init + + ! Compute initial values of energy and water integrals, + ! and zero out cumulative boundary tendencies. +!> \section arg_table_check_energy_chng_timestep_init Argument Table +!! \htmlinclude arg_table_check_energy_chng_timestep_init.html + subroutine check_energy_chng_timestep_init( & + ncol, pver, pcnst, & + is_first_timestep, & + q, pdel, & + u, v, T, & + pintdry, phis, zm, & + cp_phys, & ! cpairv generally, cpair fixed size for subcolumns code + cp_or_cv_dycore, & + te_ini_phys, te_ini_dyn, & + tw_ini, & + te_cur_phys, te_cur_dyn, & + tw_cur, & + tend_te_tnd, tend_tw_tnd, & + temp_ini, z_ini, & + count, & + teout, & + energy_formula_physics, energy_formula_dycore, & + errmsg, errflg) + + ! This scheme is non-portable due to dependencies on cam_thermo + ! for hydrostatic energy calculation (physics and dycore formulas) + use cam_thermo, only: get_hydrostatic_energy + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS + + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver ! number of vertical layers + integer, intent(in) :: pcnst ! number of ccpp constituents + logical, intent(in) :: is_first_timestep ! is first step of initial run? + real(kind_phys), intent(in) :: q(:,:,:) ! constituent mass mixing ratios [kg kg-1] + real(kind_phys), intent(in) :: pdel(:,:) ! layer thickness [Pa] + real(kind_phys), intent(in) :: u(:,:) ! zonal wind [m s-1] + real(kind_phys), intent(in) :: v(:,:) ! meridional wind [m s-1] + real(kind_phys), intent(in) :: T(:,:) ! temperature [K] + real(kind_phys), intent(in) :: pintdry(:,:) ! interface pressure dry [Pa] + real(kind_phys), intent(in) :: phis(:) ! surface geopotential [m2 s-2] + real(kind_phys), intent(in) :: zm(:,:) ! geopotential height at layer midpoints [m] + real(kind_phys), intent(in) :: cp_phys(:,:) ! enthalpy (cpairv generally) [J kg-1 K-1] + real(kind_phys), intent(in) :: cp_or_cv_dycore(:,:) ! enthalpy or heat capacity, dycore dependent [J K-1 kg-1] + integer, intent(in) :: energy_formula_physics! total energy formulation physics + integer, intent(in) :: energy_formula_dycore ! total energy formulation dycore + + ! Output arguments + real(kind_phys), intent(out) :: temp_ini(:,:) ! initial temperature [K] + real(kind_phys), intent(out) :: z_ini(:,:) ! initial geopotential height [m] + integer, intent(out) :: count ! count of values with significant energy or water imbalances [1] + real(kind_phys), intent(out) :: teout(:) ! total energy for global fixer in next timestep [J m-2] + real(kind_phys), intent(out) :: tend_te_tnd(:) ! total energy tendency [J m-2 s-1] + real(kind_phys), intent(out) :: tend_tw_tnd(:) ! total water tendency [kg m-2 s-1] + + ! Input/Output arguments + real(kind_phys), intent(inout) :: te_ini_phys(:) ! physics formula: initial total energy [J m-2] + real(kind_phys), intent(inout) :: te_ini_dyn (:) ! dycore formula: initial total energy [J m-2] + real(kind_phys), intent(inout) :: tw_ini (:) ! initial total water [kg m-2] + real(kind_phys), intent(inout) :: te_cur_phys(:) ! physics formula: current total energy [J m-2] + real(kind_phys), intent(inout) :: te_cur_dyn (:) ! dycore formula: current total energy [J m-2] + real(kind_phys), intent(inout) :: tw_cur (:) ! current total water [kg m-2] + + ! Output arguments + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + !------------------------------------------------ + ! Physics total energy. + !------------------------------------------------ + call get_hydrostatic_energy( & + tracer = q(1:ncol,1:pver,1:pcnst), & ! moist mixing ratios + moist_mixing_ratio = .true., & + pdel_in = pdel (1:ncol,1:pver), & + cp_or_cv = cp_phys (1:ncol,1:pver), & + U = u (1:ncol,1:pver), & + V = v (1:ncol,1:pver), & + T = T (1:ncol,1:pver), & + vcoord = energy_formula_physics, & ! energy formula for physics + ptop = pintdry (1:ncol,1), & + phis = phis (1:ncol), & + te = te_ini_phys(1:ncol), & ! vertically integrated total energy + H2O = tw_ini (1:ncol) & ! v.i. total water + ) + + ! Save initial state temperature and geopotential height for use in run phase + temp_ini(:ncol,:) = T (:ncol, :) + z_ini (:ncol,:) = zm(:ncol, :) + + !------------------------------------------------ + ! Dynamical core total energy. + !------------------------------------------------ + if (energy_formula_dycore == ENERGY_FORMULA_DYCORE_SE) then + ! SE dycore specific hydrostatic energy (enthalpy) + call get_hydrostatic_energy( & + tracer = q(1:ncol,1:pver,1:pcnst), & ! moist mixing ratios + moist_mixing_ratio = .true., & + pdel_in = pdel (1:ncol,1:pver), & + cp_or_cv = cp_or_cv_dycore(1:ncol,1:pver), & + U = u (1:ncol,1:pver), & + V = v (1:ncol,1:pver), & + T = T (1:ncol,1:pver), & + vcoord = energy_formula_dycore, & ! energy formula for dycore + ptop = pintdry (1:ncol,1), & + phis = phis (1:ncol), & + te = te_ini_dyn (1:ncol) & ! WRITE OPERATION - vertically integrated total energy + ) + + else if (energy_formula_dycore == ENERGY_FORMULA_DYCORE_MPAS) then + ! MPAS dycore: compute cv if vertical coordinate is height: cv = cp - R (internal energy) + call get_hydrostatic_energy( & + tracer = q(1:ncol,1:pver,1:pcnst), & ! moist mixing ratios + moist_mixing_ratio = .true., & + pdel_in = pdel (1:ncol,1:pver), & + cp_or_cv = cp_or_cv_dycore(1:ncol,1:pver), & + U = u (1:ncol,1:pver), & + V = v (1:ncol,1:pver), & + T = T (1:ncol,1:pver), & ! enthalpy-scaled temperature for energy consistency + vcoord = energy_formula_dycore, & ! energy formula for dycore + ptop = pintdry (1:ncol,1), & + phis = phis (1:ncol), & + z_mid = z_ini (1:ncol,:), & ! unique for MPAS + te = te_ini_dyn (1:ncol) & ! WRITE OPERATION - vertically integrated total energy + ) + else + ! FV dycore: dycore energy is the same as physics + te_ini_dyn(:ncol) = te_ini_phys(:ncol) + endif + + ! Set current state to be the same as initial + te_cur_phys(:ncol) = te_ini_phys(:ncol) + te_cur_dyn (:ncol) = te_ini_dyn (:ncol) + tw_cur (:ncol) = tw_ini (:ncol) + + ! Zero out current energy unbalances count + count = 0 + + ! Zero out cumulative boundary fluxes + tend_te_tnd(:ncol) = 0._kind_phys + tend_tw_tnd(:ncol) = 0._kind_phys + + ! If first timestep, initialize value of teout + if(is_first_timestep) then + teout(:ncol) = te_ini_dyn(:ncol) + endif + + end subroutine check_energy_chng_timestep_init + + + ! Check that energy and water change matches the boundary fluxes. +!> \section arg_table_check_energy_chng_run Argument Table +!! \htmlinclude arg_table_check_energy_chng_run.html + subroutine check_energy_chng_run( & + ncol, pver, pcnst, & + iulog, & + q, pdel, & + u, v, T, & + pintdry, phis, zm, & + cp_phys, & ! cpairv generally, cpair fixed size for subcolumns code + cp_or_cv_dycore, & + scaling_dycore, & ! From check_energy_scaling to work around subcolumns code + te_cur_phys, te_cur_dyn, & + tw_cur, & + tend_te_tnd, tend_tw_tnd, & + temp_ini, z_ini, & + count, ztodt, & + latice, latvap, & + energy_formula_physics, energy_formula_dycore, & + name, flx_vap, flx_cnd, flx_ice, flx_sen, & + errmsg, errflg) + + ! This scheme is non-portable due to dependencies on cam_thermo + ! for hydrostatic energy calculation (physics and dycore formulas) + use cam_thermo, only: get_hydrostatic_energy + + ! Dependency for energy formula used by physics and dynamical cores + use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_FV, ENERGY_FORMULA_DYCORE_SE, ENERGY_FORMULA_DYCORE_MPAS + + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver ! number of vertical layers + integer, intent(in) :: pcnst ! number of ccpp constituents + integer, intent(in) :: iulog ! log output unit + real(kind_phys), intent(in) :: q(:,:,:) ! constituent mass mixing ratios [kg kg-1] + real(kind_phys), intent(in) :: pdel(:,:) ! layer thickness [Pa] + real(kind_phys), intent(in) :: u(:,:) ! zonal wind [m s-1] + real(kind_phys), intent(in) :: v(:,:) ! meridional wind [m s-1] + real(kind_phys), intent(in) :: T(:,:) ! temperature [K] + real(kind_phys), intent(in) :: pintdry(:,:) ! interface pressure dry [Pa] + real(kind_phys), intent(in) :: phis(:) ! surface geopotential [m2 s-2] + real(kind_phys), intent(in) :: zm(:,:) ! geopotential height at layer midpoints [m] + real(kind_phys), intent(in) :: temp_ini(:,:) ! initial temperature [K] + real(kind_phys), intent(in) :: z_ini(:,:) ! initial geopotential height [m] + real(kind_phys), intent(in) :: cp_phys(:,:) ! enthalpy (cpairv generally) [J kg-1 K-1] + real(kind_phys), intent(in) :: cp_or_cv_dycore(:,:) ! enthalpy or heat capacity, dycore dependent [J K-1 kg-1] + real(kind_phys), intent(in) :: scaling_dycore(:,:) ! scaling for conversion of temperature increment [1] + real(kind_phys), intent(in) :: ztodt ! physics timestep [s] + real(kind_phys), intent(in) :: latice ! constant, latent heat of fusion of water at 0 C [J kg-1] + real(kind_phys), intent(in) :: latvap ! constant, latent heat of vaporization of water at 0 C [J kg-1] + integer, intent(in) :: energy_formula_physics! total energy formulation physics + integer, intent(in) :: energy_formula_dycore ! total energy formulation dycore + + ! Input from CCPP-scheme being checked: + ! parameterization name; surface fluxes of (1) vapor, (2) liquid+ice, (3) ice, (4) sensible heat + ! to pass in the values to be checked, call check_energy_zero_input_fluxes to reset these values + ! before a parameterization that is checked, then update these values as-needed + ! (can be all zero; in fact, most parameterizations calling _chng call with zero arguments) + ! + ! Original comment from BAB: + ! Note that the precip and ice fluxes are in precip units (m/s). + ! I would prefer to have kg/m2/s. + ! I would also prefer liquid (not total) and ice fluxes + character(len=*), intent(in) :: name ! parameterization name for fluxes + real(kind_phys), intent(in) :: flx_vap(:) ! boundary flux of vapor [kg m-2 s-1] + real(kind_phys), intent(in) :: flx_cnd(:) ! boundary flux of liquid+ice (precip?) [m s-1] + real(kind_phys), intent(in) :: flx_ice(:) ! boundary flux of ice [m s-1] + real(kind_phys), intent(in) :: flx_sen(:) ! boundary flux of sensible heat [W m-2] + + ! Input/Output arguments + real(kind_phys), intent(inout) :: te_cur_phys(:) ! physics formula: current total energy [J m-2] + real(kind_phys), intent(inout) :: te_cur_dyn (:) ! dycore formula: current total energy [J m-2] + real(kind_phys), intent(inout) :: tw_cur (:) ! current total water [kg m-2] + integer, intent(inout) :: count ! count of columns with significant energy or water imbalances [1] + real(kind_phys), intent(inout) :: tend_te_tnd(:) ! total energy tendency [J m-2 s-1] + real(kind_phys), intent(inout) :: tend_tw_tnd(:) ! total water tendency [kg m-2 s-1] + + ! Output arguments + character(len=512), intent(out) :: errmsg ! error message + integer, intent(out) :: errflg ! error flag + + ! Local variables + real(kind_phys) :: te_xpd(ncol) ! expected value (f0 + dt*boundary_flux) + real(kind_phys) :: te_dif(ncol) ! energy of input state - original energy + real(kind_phys) :: te_tnd(ncol) ! tendency from last process + real(kind_phys) :: te_rer(ncol) ! relative error in energy column + + real(kind_phys) :: tw_xpd(ncol) ! expected value (w0 + dt*boundary_flux) + real(kind_phys) :: tw_dif(ncol) ! tw_inp - original water + real(kind_phys) :: tw_tnd(ncol) ! tendency from last process + real(kind_phys) :: tw_rer(ncol) ! relative error in water column + + real(kind_phys) :: te(ncol) ! vertical integral of total energy + real(kind_phys) :: tw(ncol) ! vertical integral of total water + real(kind_phys) :: temp(ncol,pver) ! temperature + + real(kind_phys) :: se(ncol) ! enthalpy or internal energy (J/m2) + real(kind_phys) :: po(ncol) ! surface potential or potential energy (J/m2) + real(kind_phys) :: ke(ncol) ! kinetic energy (J/m2) + real(kind_phys) :: wv(ncol) ! column integrated vapor (kg/m2) + real(kind_phys) :: liq(ncol) ! column integrated liquid (kg/m2) + real(kind_phys) :: ice(ncol) ! column integrated ice (kg/m2) + + integer :: i + + !------------------------------------------------ + ! Physics total energy. + !------------------------------------------------ + call get_hydrostatic_energy( & + tracer = q(1:ncol,1:pver,1:pcnst), & ! moist mixing ratios + moist_mixing_ratio = .true., & + pdel_in = pdel (1:ncol,1:pver), & + cp_or_cv = cp_phys(1:ncol,1:pver), & + U = u (1:ncol,1:pver), & + V = v (1:ncol,1:pver), & + T = T (1:ncol,1:pver), & + vcoord = energy_formula_physics, & ! energy formula for physics + ptop = pintdry(1:ncol,1), & + phis = phis (1:ncol), & + te = te (1:ncol), & ! vertically integrated total energy + H2O = tw (1:ncol), & ! v.i. total water + se = se (1:ncol), & ! v.i. enthalpy + po = po (1:ncol), & ! v.i. PHIS term + ke = ke (1:ncol), & ! v.i. kinetic energy + wv = wv (1:ncol), & ! v.i. water vapor + liq = liq (1:ncol), & ! v.i. liquid + ice = ice (1:ncol) & ! v.i. ice + ) + + ! compute expected values and tendencies + do i = 1, ncol + ! change in static energy and total water + te_dif(i) = te(i) - te_cur_phys(i) + tw_dif(i) = tw(i) - tw_cur (i) + + ! expected tendencies from boundary fluxes for last process + te_tnd(i) = flx_vap(i)*(latvap+latice) - (flx_cnd(i) - flx_ice(i))*1000._kind_phys*latice + flx_sen(i) + tw_tnd(i) = flx_vap(i) - flx_cnd(i) *1000._kind_phys + + ! cummulative tendencies from boundary fluxes + tend_te_tnd(i) = tend_te_tnd(i) + te_tnd(i) + tend_tw_tnd(i) = tend_tw_tnd(i) + tw_tnd(i) + + ! expected new values from previous state plus boundary fluxes + te_xpd(i) = te_cur_phys(i) + te_tnd(i)*ztodt + tw_xpd(i) = tw_cur (i) + tw_tnd(i)*ztodt + + ! relative error, expected value - input state / previous state + te_rer(i) = (te_xpd(i) - te(i)) / te_cur_phys(i) + end do + + ! relative error for total water (allow for dry atmosphere) + tw_rer = 0._kind_phys + where (tw_cur(:ncol) > 0._kind_phys) + tw_rer(:ncol) = (tw_xpd(:ncol) - tw(:ncol)) / tw_cur(:ncol) + end where + + ! error checking + if (print_energy_errors) then + if (any(abs(te_rer(1:ncol)) > 1.E-14_kind_phys .or. abs(tw_rer(1:ncol)) > 1.E-10_kind_phys)) then + do i = 1, ncol + ! the relative error threshold for the water budget has been reduced to 1.e-10 + ! to avoid messages generated by QNEG3 calls + ! PJR- change to identify if error in energy or water + if (abs(te_rer(i)) > 1.E-14_kind_phys ) then + count = count + 1 + write(iulog,*) "significant energy conservation error after ", name, & + " count", count, "col", i + write(iulog,*) te(i),te_xpd(i),te_dif(i),tend_te_tnd(i)*ztodt, & + te_tnd(i)*ztodt,te_rer(i) + endif + if ( abs(tw_rer(i)) > 1.E-10_kind_phys) then + count = count + 1 + write(iulog,*) "significant water conservation error after ", name, & + " count", count, "col", i + write(iulog,*) tw(i),tw_xpd(i),tw_dif(i),tend_tw_tnd(i)*ztodt, & + tw_tnd(i)*ztodt,tw_rer(i) + end if + end do + end if + end if + + ! WRITE OPERATION - copy new value to state, including total water. + ! the total water operations are consistent regardless of energy formula, + ! so it only has to be written once. + do i = 1, ncol + te_cur_phys(i) = te(i) + tw_cur(i) = tw(i) + end do + + !------------------------------------------------ + ! Dynamical core total energy. + !------------------------------------------------ + if (energy_formula_dycore == ENERGY_FORMULA_DYCORE_SE) then + ! SE dycore specific hydrostatic energy + + ! enthalpy scaling for energy consistency + temp(1:ncol,:) = temp_ini(1:ncol,:)+scaling_dycore(1:ncol,:)*(T(1:ncol,:)-temp_ini(1:ncol,:)) + + call get_hydrostatic_energy( & + tracer = q(1:ncol,1:pver,1:pcnst), & ! moist mixing ratios + moist_mixing_ratio = .true., & + pdel_in = pdel (1:ncol,1:pver), & + cp_or_cv = cp_or_cv_dycore(1:ncol,1:pver), & + U = u (1:ncol,1:pver), & + V = v (1:ncol,1:pver), & + T = temp (1:ncol,1:pver), & ! enthalpy-scaled temperature for energy consistency + vcoord = energy_formula_dycore, & ! energy formula for dycore + ptop = pintdry (1:ncol,1), & + phis = phis (1:ncol), & + te = te_cur_dyn (1:ncol) & ! WRITE OPERATION - vertically integrated total energy + ) + + else if (energy_formula_dycore == ENERGY_FORMULA_DYCORE_MPAS) then + ! MPAS dycore: compute cv if vertical coordinate is height: cv = cp - R + + ! REMOVECAM: note this scaling is different with subcols off/on which is why it was put into separate scheme (hplin, 9/5/24) + temp(1:ncol,:) = temp_ini(1:ncol,:)+scaling_dycore(1:ncol,:)*(T(1:ncol,:)-temp_ini(1:ncol,:)) + + call get_hydrostatic_energy( & + tracer = q(1:ncol,1:pver,1:pcnst), & ! moist mixing ratios + moist_mixing_ratio = .true., & + pdel_in = pdel (1:ncol,1:pver), & + cp_or_cv = cp_or_cv_dycore(1:ncol,1:pver), & + U = u (1:ncol,1:pver), & + V = v (1:ncol,1:pver), & + T = temp (1:ncol,1:pver), & ! enthalpy-scaled temperature for energy consistency + vcoord = energy_formula_dycore, & ! energy formula for dycore + ptop = pintdry (1:ncol,1), & + phis = phis (1:ncol), & + z_mid = z_ini (1:ncol,:), & ! unique for MPAS + te = te_cur_dyn (1:ncol) & ! WRITE OPERATION - vertically integrated total energy + ) + + else + ! FV dycore + te_cur_dyn(1:ncol) = te(1:ncol) + end if + end subroutine check_energy_chng_run + +end module check_energy_chng diff --git a/schemes/check_energy/check_energy_chng.meta b/schemes/check_energy/check_energy_chng.meta new file mode 100644 index 00000000..f69792df --- /dev/null +++ b/schemes/check_energy/check_energy_chng.meta @@ -0,0 +1,413 @@ +[ccpp-table-properties] + name = check_energy_chng + type = scheme + dependencies = ../../../../data/cam_thermo.F90,../../../../data/cam_thermo_formula.F90 + + +[ccpp-arg-table] + name = check_energy_chng_init + type = scheme +[ print_energy_errors_in ] + standard_name = flag_for_energy_conservation_warning + units = flag + type = logical + dimensions = () + intent = in + +[ccpp-arg-table] + name = check_energy_chng_timestep_init + type = scheme +[ ncol ] + standard_name = horizontal_dimension + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ pcnst ] + standard_name = number_of_ccpp_constituents + units = count + type = integer + dimensions = () + intent = in +[ is_first_timestep ] + standard_name = is_first_timestep + units = flag + type = logical + dimensions = () + intent = in +[ q ] + standard_name = ccpp_constituents + units = none + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension, number_of_ccpp_constituents) + intent = in +[ pdel ] + standard_name = air_pressure_thickness + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension) + intent = in +[ u ] + standard_name = eastward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension) + intent = in +[ v ] + standard_name = northward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension) + intent = in +[ T ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension) + intent = in +[ pintdry ] + standard_name = air_pressure_of_dry_air_at_interface + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_interface_dimension) + intent = in +[ phis ] + standard_name = surface_geopotential + units = m2 s-2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = in +[ zm ] + standard_name = geopotential_height_wrt_surface + units = m + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension) + intent = in +[ cp_phys ] + standard_name = composition_dependent_specific_heat_of_dry_air_at_constant_pressure + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension) + intent = in +[ cp_or_cv_dycore ] + standard_name = specific_heat_of_air_used_in_dycore + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension) + intent = in +[ te_ini_phys ] + standard_name = vertically_integrated_total_energy_using_physics_energy_formula_at_start_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = inout +[ te_ini_dyn ] + standard_name = vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = inout +[ tw_ini ] + standard_name = vertically_integrated_total_water_at_start_of_physics_timestep + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = inout +[ te_cur_phys ] + standard_name = vertically_integrated_total_energy_using_physics_energy_formula + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = inout +[ te_cur_dyn ] + standard_name = vertically_integrated_total_energy_using_dycore_energy_formula + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = inout +[ tw_cur ] + standard_name = vertically_integrated_total_water + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = inout +[ tend_te_tnd ] + standard_name = cumulative_total_energy_boundary_flux_using_physics_energy_formula + units = J m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = out +[ tend_tw_tnd ] + standard_name = cumulative_total_water_boundary_flux + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = out +[ temp_ini ] + standard_name = air_temperature_at_start_of_physics_timestep + units = K + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension) + intent = out +[ z_ini ] + standard_name = geopotential_height_wrt_surface_at_start_of_physics_timestep + units = m + type = real | kind = kind_phys + dimensions = (horizontal_dimension, vertical_layer_dimension) + intent = out +[ count ] + standard_name = number_of_atmosphere_columns_with_significant_energy_or_water_imbalances + units = count + type = integer + dimensions = () + intent = out +[ teout ] + standard_name = vertically_integrated_total_energy_at_end_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = out +[ energy_formula_physics ] + standard_name = total_energy_formula_for_physics + units = 1 + type = integer + dimensions = () + intent = in +[ energy_formula_dycore ] + standard_name = total_energy_formula_for_dycore + units = 1 + type = integer + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = check_energy_chng_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ pcnst ] + standard_name = number_of_ccpp_constituents + units = count + type = integer + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ q ] + standard_name = ccpp_constituents + units = none + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents) + intent = in +[ pdel ] + standard_name = air_pressure_thickness + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ u ] + standard_name = eastward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ v ] + standard_name = northward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ T ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ pintdry ] + standard_name = air_pressure_of_dry_air_at_interface + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ phis ] + standard_name = surface_geopotential + units = m2 s-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ zm ] + standard_name = geopotential_height_wrt_surface + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cp_phys ] + standard_name = composition_dependent_specific_heat_of_dry_air_at_constant_pressure + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cp_or_cv_dycore ] + standard_name = specific_heat_of_air_used_in_dycore + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ scaling_dycore ] + standard_name = ratio_of_specific_heat_of_air_used_in_physics_energy_formula_to_specific_heat_of_air_used_in_dycore_energy_formula + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ te_cur_phys ] + standard_name = vertically_integrated_total_energy_using_physics_energy_formula + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ te_cur_dyn ] + standard_name = vertically_integrated_total_energy_using_dycore_energy_formula + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ tw_cur ] + standard_name = vertically_integrated_total_water + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ tend_te_tnd ] + standard_name = cumulative_total_energy_boundary_flux_using_physics_energy_formula + units = J m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ tend_tw_tnd ] + standard_name = cumulative_total_water_boundary_flux + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ temp_ini ] + standard_name = air_temperature_at_start_of_physics_timestep + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ z_ini ] + standard_name = geopotential_height_wrt_surface_at_start_of_physics_timestep + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ count ] + standard_name = number_of_atmosphere_columns_with_significant_energy_or_water_imbalances + units = count + type = integer + dimensions = () + intent = inout +[ ztodt ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ latice ] + standard_name = latent_heat_of_fusion_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ latvap ] + standard_name = latent_heat_of_vaporization_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ energy_formula_physics ] + standard_name = total_energy_formula_for_physics + units = 1 + type = integer + dimensions = () + intent = in +[ energy_formula_dycore ] + standard_name = total_energy_formula_for_dycore + units = 1 + type = integer + dimensions = () + intent = in +[ name ] + standard_name = scheme_name + units = none + type = character | kind = len=* + dimensions = () + intent = in +[ flx_vap ] + standard_name = net_water_vapor_fluxes_through_top_and_bottom_of_atmosphere_column + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ flx_cnd ] + standard_name = net_liquid_and_lwe_ice_fluxes_through_top_and_bottom_of_atmosphere_column + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ flx_ice ] + standard_name = net_lwe_ice_fluxes_through_top_and_bottom_of_atmosphere_column + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ flx_sen ] + standard_name = net_sensible_heat_flux_through_top_and_bottom_of_atmosphere_column + units = W m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/check_energy/check_energy_chng_namelist.xml b/schemes/check_energy/check_energy_chng_namelist.xml new file mode 100644 index 00000000..e6ff1bbc --- /dev/null +++ b/schemes/check_energy/check_energy_chng_namelist.xml @@ -0,0 +1,19 @@ + + + + + + + logical + diagnostics + check_energy_nl + flag_for_energy_conservation_warning + flag + + Turn on verbose output identifying columns that fail energy/water conservation checks. Default: FALSE + + + false + + + diff --git a/schemes/check_energy/check_energy_fix.F90 b/schemes/check_energy/check_energy_fix.F90 new file mode 100644 index 00000000..1bdc7ae8 --- /dev/null +++ b/schemes/check_energy/check_energy_fix.F90 @@ -0,0 +1,42 @@ +module check_energy_fix + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: check_energy_fix_run + +contains + + ! Add heating rate required for global mean total energy conservation +!> \section arg_table_check_energy_fix_run Argument Table +!! \htmlinclude arg_table_check_energy_fix_run.html + subroutine check_energy_fix_run(ncol, pver, pint, gravit, heat_glob, ptend_s, eshflx, scheme_name) + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver ! number of vertical layers + real(kind_phys), intent(in) :: pint(:,:) ! interface pressure [Pa] + real(kind_phys), intent(in) :: gravit ! gravitational acceleration [m s-2] + real(kind_phys), intent(in) :: heat_glob ! global mean heating rate [J kg-1 s-1] + real(kind_phys), intent(out) :: ptend_s(:,:) ! physics tendency heating rate [J kg-1 s-1] + real(kind_phys), intent(out) :: eshflx(:) ! effective sensible heat flux [W m-2] + ! for check_energy_chng + + character(len=64), intent(out) :: scheme_name ! scheme name + + ! Local variables + integer :: i + + ! Set scheme name for check_energy_chng + scheme_name = "check_energy_fix" + + ! add (-) global mean total energy difference as heating + ptend_s(:ncol, :pver) = heat_glob + + ! compute effective sensible heat flux + do i = 1, ncol + eshflx(i) = heat_glob * (pint(i,pver+1) - pint(i,1)) / gravit + end do + end subroutine check_energy_fix_run + +end module check_energy_fix diff --git a/schemes/check_energy/check_energy_fix.meta b/schemes/check_energy/check_energy_fix.meta new file mode 100644 index 00000000..b5f935fb --- /dev/null +++ b/schemes/check_energy/check_energy_fix.meta @@ -0,0 +1,55 @@ +[ccpp-table-properties] + name = check_energy_fix + type = scheme + +[ccpp-arg-table] + name = check_energy_fix_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ pint ] + standard_name = air_pressure_at_interface + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ heat_glob ] + standard_name = global_mean_heating_rate_correction_for_energy_conservation + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ ptend_s ] + standard_name = tendency_of_dry_air_enthalpy_at_constant_pressure + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ eshflx ] + standard_name = net_sensible_heat_flux_through_top_and_bottom_of_atmosphere_column + units = W m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ scheme_name ] + standard_name = scheme_name + units = none + type = character | kind = len=64 + dimensions = () + intent = out diff --git a/schemes/check_energy/check_energy_gmean/check_energy_gmean.F90 b/schemes/check_energy/check_energy_gmean/check_energy_gmean.F90 new file mode 100644 index 00000000..4d50c428 --- /dev/null +++ b/schemes/check_energy/check_energy_gmean/check_energy_gmean.F90 @@ -0,0 +1,70 @@ +module check_energy_gmean + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: check_energy_gmean_run + +contains + + ! Compute global mean total energy of physics input and output states + ! computed consistently with dynamical core vertical coordinate + ! (under hydrostatic assumption) +!> \section arg_table_check_energy_gmean_run Argument Table +!! \htmlinclude arg_table_check_energy_gmean_run.html + subroutine check_energy_gmean_run( & + ncol, pver, dtime, & + gravit, & + pint, & + te_ini_dyn, teout, & + tedif_glob, heat_glob, & + teinp_glob, teout_glob, psurf_glob, ptopb_glob) + + ! This scheme is non-portable due to dependency: Global mean module gmean from src/utils + use gmean_mod, only: gmean + + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver ! number of vertical layers + real(kind_phys), intent(in) :: dtime ! physics time step [s] + real(kind_phys), intent(in) :: gravit ! gravitational acceleration [m s-2] + real(kind_phys), intent(in) :: pint(:,:) ! interface pressure [Pa] + real(kind_phys), intent(in) :: te_ini_dyn(:) ! dycore formula: initial total energy [J m-2] + real(kind_phys), intent(in) :: teout(:) ! total energy for global fixer in next timestep [J m-2] + + ! Output arguments + real(kind_phys), intent(out) :: tedif_glob ! global mean energy difference [J m-2] + real(kind_phys), intent(out) :: heat_glob ! global mean heating rate [J kg-1 s-1] + + ! Output for check_energy_gmean_diagnostics + real(kind_phys), intent(out) :: teinp_glob ! global mean energy of input state [J m-2] + real(kind_phys), intent(out) :: teout_glob ! global mean energy of output state [J m-2] + real(kind_phys), intent(out) :: psurf_glob ! global mean surface pressure [Pa] + real(kind_phys), intent(out) :: ptopb_glob ! global mean top boundary pressure [Pa] + + ! Local variables + real(kind_phys) :: te(ncol, 4) ! total energy of input/output states (copy) + real(kind_phys) :: te_glob(4) ! global means of total energy + + ! Copy total energy out of input and output states. + ! These four fields will have their global means calculated respectively + te(:ncol, 1) = te_ini_dyn(:ncol) ! Input energy using dycore energy formula [J m-2] + te(:ncol, 2) = teout(:ncol) ! Total energy from end of physics timestep [J m-2] + te(:ncol, 3) = pint(:ncol, pver+1) ! Surface pressure for heating rate [Pa] + te(:ncol, 4) = pint(:ncol, 1) ! Model top pressure for heating rate [Pa] + ! not constant for z-based vertical coordinate + + ! Compute global means of input and output energies and of surface pressure for heating rate. + call gmean(te, te_glob, 4) + teinp_glob = te_glob(1) + teout_glob = te_glob(2) + psurf_glob = te_glob(3) + ptopb_glob = te_glob(4) + + ! Compute global mean total energy difference for check_energy_fix + tedif_glob = teinp_glob - teout_glob + heat_glob = -tedif_glob/dtime * gravit / (psurf_glob - ptopb_glob) ! [J kg-1 s-1] + end subroutine check_energy_gmean_run + +end module check_energy_gmean diff --git a/schemes/check_energy/check_energy_gmean/check_energy_gmean.meta b/schemes/check_energy/check_energy_gmean/check_energy_gmean.meta new file mode 100644 index 00000000..8ebf6f83 --- /dev/null +++ b/schemes/check_energy/check_energy_gmean/check_energy_gmean.meta @@ -0,0 +1,86 @@ +[ccpp-table-properties] + name = check_energy_gmean + type = scheme + dependencies = ../../../../../utils/gmean_mod.F90 + +[ccpp-arg-table] + name = check_energy_gmean_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ dtime ] + standard_name = timestep_for_physics + units = s + type = real | kind = kind_phys + dimensions = () + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ pint ] + standard_name = air_pressure_at_interface + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ te_ini_dyn ] + standard_name = vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ teout ] + standard_name = vertically_integrated_total_energy_at_end_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ tedif_glob ] + standard_name = global_mean_total_energy_correction_for_energy_conservation + units = J m-2 + type = real | kind = kind_phys + dimensions = () + intent = out +[ heat_glob ] + standard_name = global_mean_heating_rate_correction_for_energy_conservation + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = () + intent = out +[ teinp_glob ] + standard_name = global_mean_vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = () + intent = out +[ teout_glob ] + standard_name = global_mean_vertically_integrated_total_energy_at_end_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = () + intent = out +[ psurf_glob ] + standard_name = global_mean_surface_air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = () + intent = out +[ ptopb_glob ] + standard_name = global_mean_air_pressure_at_top_of_atmosphere_model + units = Pa + type = real | kind = kind_phys + dimensions = () + intent = out diff --git a/schemes/check_energy/check_energy_save_teout.F90 b/schemes/check_energy/check_energy_save_teout.F90 new file mode 100644 index 00000000..6ccdc15c --- /dev/null +++ b/schemes/check_energy/check_energy_save_teout.F90 @@ -0,0 +1,32 @@ +! save total energy for global fixer in next timestep +! this must be called after the last parameterization and physics_update, +! and after a final check_energy_chng to compute te_cur. +module check_energy_save_teout + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: check_energy_save_teout_run + +contains + +!> \section arg_table_check_energy_save_teout_run Argument Table +!! \htmlinclude arg_table_check_energy_save_teout_run.html + subroutine check_energy_save_teout_run(ncol, te_cur_dyn, teout) + + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + real(kind_phys), intent(in) :: te_cur_dyn (:) ! dycore formula: current total energy [J m-2] + + ! Output arguments + real(kind_phys), intent(out) :: teout(:) ! total energy for global fixer in next timestep [J m-2] + + ! nb hplin: note that in physpkg.F90, the pbuf is updated to the previous dyn timestep + ! through itim_old. Need to check if we need to replicate such pbuf functionality + ! in the CAM-SIMA/CCPP infrastructure. + teout(:ncol) = te_cur_dyn(:ncol) + + end subroutine check_energy_save_teout_run + +end module check_energy_save_teout diff --git a/schemes/check_energy/check_energy_save_teout.meta b/schemes/check_energy/check_energy_save_teout.meta new file mode 100644 index 00000000..57a712a0 --- /dev/null +++ b/schemes/check_energy/check_energy_save_teout.meta @@ -0,0 +1,25 @@ +[ccpp-table-properties] + name = check_energy_save_teout + type = scheme + +[ccpp-arg-table] + name = check_energy_save_teout_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ te_cur_dyn ] + standard_name = vertically_integrated_total_energy_using_dycore_energy_formula + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ teout ] + standard_name = vertically_integrated_total_energy_at_end_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out diff --git a/schemes/check_energy/check_energy_scaling.F90 b/schemes/check_energy/check_energy_scaling.F90 new file mode 100644 index 00000000..c10818c1 --- /dev/null +++ b/schemes/check_energy/check_energy_scaling.F90 @@ -0,0 +1,42 @@ +module check_energy_scaling + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: check_energy_scaling_run + +contains + + ! CCPP routine to get scaling factor for conversion of temperature increment. + ! This is extracted to a separate subroutine so that scaling_dycore can be passed + ! directly to the CCPP-ized check_energy_chng_run subroutine from CAM with subcolumns. + ! + ! When subcolumns are removed from CAM, this dummy scheme can be removed, and + ! scaling_dycore can just be calculated in check_energy_chng. (hplin, 9/5/24) +!> \section arg_table_check_energy_scaling_run Argument Table +!! \htmlinclude arg_table_check_energy_scaling_run.html + subroutine check_energy_scaling_run( & + ncol, & + cp_or_cv_dycore, cpairv, & + scaling_dycore, & + errmsg, errflg) + + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + real(kind_phys), intent(in) :: cp_or_cv_dycore(:,:) ! cp or cv from dycore [J kg-1 K-1] + real(kind_phys), intent(in) :: cpairv(:,:) ! specific heat of dry air at constant pressure [J kg-1 K-1] + + ! Output arguments + real(kind_phys), intent(out) :: scaling_dycore(:,:) ! scaling for conversion of temperature increment [1] + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + scaling_dycore(:ncol,:) = cpairv(:ncol,:) / cp_or_cv_dycore(:ncol,:) + + end subroutine check_energy_scaling_run + +end module check_energy_scaling diff --git a/schemes/check_energy/check_energy_scaling.meta b/schemes/check_energy/check_energy_scaling.meta new file mode 100644 index 00000000..400e97e5 --- /dev/null +++ b/schemes/check_energy/check_energy_scaling.meta @@ -0,0 +1,43 @@ +[ccpp-table-properties] + name = check_energy_scaling + type = scheme + +[ccpp-arg-table] + name = check_energy_scaling_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ cp_or_cv_dycore ] + standard_name = specific_heat_of_air_used_in_dycore + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cpairv ] + standard_name = composition_dependent_specific_heat_of_dry_air_at_constant_pressure + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ scaling_dycore ] + standard_name = ratio_of_specific_heat_of_air_used_in_physics_energy_formula_to_specific_heat_of_air_used_in_dycore_energy_formula + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/check_energy/check_energy_zero_fluxes.F90 b/schemes/check_energy/check_energy_zero_fluxes.F90 new file mode 100644 index 00000000..bda0d74c --- /dev/null +++ b/schemes/check_energy/check_energy_zero_fluxes.F90 @@ -0,0 +1,34 @@ +! zeros input fluxes to check_energy +! before running other schemes +module check_energy_zero_fluxes + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: check_energy_zero_fluxes_run + +contains + +!> \section arg_table_check_energy_zero_fluxes_run Argument Table +!! \htmlinclude arg_table_check_energy_zero_fluxes_run.html + subroutine check_energy_zero_fluxes_run(ncol, name, flx_vap, flx_cnd, flx_ice, flx_sen) + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + + ! Output arguments + character(len=64), intent(out) :: name ! parameterization name for fluxes + real(kind_phys), intent(out) :: flx_vap(:) ! boundary flux of vapor [kg m-2 s-1] + real(kind_phys), intent(out) :: flx_cnd(:) ! boundary flux of liquid+ice (precip?) [m s-1] + real(kind_phys), intent(out) :: flx_ice(:) ! boundary flux of ice (snow?) [m s-1] + real(kind_phys), intent(out) :: flx_sen(:) ! boundary flux of sensible heat [W m-2] + + ! reset values to zero + name = '' + flx_vap(:) = 0._kind_phys + flx_cnd(:) = 0._kind_phys + flx_ice(:) = 0._kind_phys + flx_sen(:) = 0._kind_phys + end subroutine check_energy_zero_fluxes_run + +end module check_energy_zero_fluxes diff --git a/schemes/check_energy/check_energy_zero_fluxes.meta b/schemes/check_energy/check_energy_zero_fluxes.meta new file mode 100644 index 00000000..65782599 --- /dev/null +++ b/schemes/check_energy/check_energy_zero_fluxes.meta @@ -0,0 +1,43 @@ +[ccpp-table-properties] + name = check_energy_zero_fluxes + type = scheme + +[ccpp-arg-table] + name = check_energy_zero_fluxes_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ name ] + standard_name = scheme_name + units = none + type = character | kind = len=64 + dimensions = () + intent = out +[ flx_vap ] + standard_name = net_water_vapor_fluxes_through_top_and_bottom_of_atmosphere_column + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ flx_cnd ] + standard_name = net_liquid_and_lwe_ice_fluxes_through_top_and_bottom_of_atmosphere_column + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ flx_ice ] + standard_name = net_lwe_ice_fluxes_through_top_and_bottom_of_atmosphere_column + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ flx_sen ] + standard_name = net_sensible_heat_flux_through_top_and_bottom_of_atmosphere_column + units = W m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out diff --git a/schemes/check_energy/dycore_energy_consistency_adjust.F90 b/schemes/check_energy/dycore_energy_consistency_adjust.F90 new file mode 100644 index 00000000..14fce382 --- /dev/null +++ b/schemes/check_energy/dycore_energy_consistency_adjust.F90 @@ -0,0 +1,46 @@ +! MPAS and SE dynamical core specific +! 1) scaling of temperature for enforcing energy consistency +! 2) and to ensure correct computation of temperature dependent diagnostic tendencies, e.g., dtcore +module dycore_energy_consistency_adjust + use ccpp_kinds, only: kind_phys + implicit none + private + + public :: dycore_energy_consistency_adjust_run + +contains +!> \section arg_table_dycore_energy_consistency_adjust_run Argument Table +!! \htmlinclude arg_table_dycore_energy_consistency_adjust_run.html + subroutine dycore_energy_consistency_adjust_run( & + ncol, pver, & + do_consistency_adjust, & + scaling_dycore, & + tend_dTdt, & + tend_dTdt_local) + + ! Input arguments + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver ! number of vertical layers + logical, intent(in) :: do_consistency_adjust ! do energy consistency adjustment? + real(kind_phys), intent(in) :: scaling_dycore(:,:) ! scaling for conversion of temperature increment [1] + real(kind_phys), intent(in) :: tend_dTdt(:,:) ! model physics temperature tendency [K s-1] + + ! Output arguments + real(kind_phys), intent(out) :: tend_dTdt_local(:,:) ! (scheme) temperature tendency [K s-1] + + if (do_consistency_adjust) then + ! original formula for scaling of temperature: + ! T(:ncol,:) = temp_ini(:ncol,:) + & + ! scaling_dycore(:ncol,:) * (T(:ncol,:) - temp_ini(:ncol,:)) + ! and temperature tendency due to model physics: + ! tend_dTdt(:ncol,:) = scaling_dycore(:ncol,:) * tend_dTdt(:ncol,:) + ! + ! the terms can be arranged for this scaling to be applied through scheme tendencies + ! at the cost of a round-off level difference + tend_dTdt_local(:ncol,:) = (scaling_dycore(:ncol,:) - 1._kind_phys) * tend_dTdt(:ncol,:) + endif + ! do nothing for dynamical cores with energy consistent with CAM physics + + end subroutine dycore_energy_consistency_adjust_run + +end module dycore_energy_consistency_adjust diff --git a/schemes/check_energy/dycore_energy_consistency_adjust.meta b/schemes/check_energy/dycore_energy_consistency_adjust.meta new file mode 100644 index 00000000..e1afb294 --- /dev/null +++ b/schemes/check_energy/dycore_energy_consistency_adjust.meta @@ -0,0 +1,43 @@ +[ccpp-table-properties] + name = dycore_energy_consistency_adjust + type = scheme + +[ccpp-arg-table] + name = dycore_energy_consistency_adjust_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ do_consistency_adjust ] + standard_name = flag_for_dycore_energy_consistency_adjustment + units = flag + type = logical + dimensions = () + intent = in +[ scaling_dycore ] + standard_name = ratio_of_specific_heat_of_air_used_in_physics_energy_formula_to_specific_heat_of_air_used_in_dycore_energy_formula + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ tend_dTdt ] + standard_name = tendency_of_air_temperature_due_to_model_physics + units = K s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ tend_dTdt_local ] + standard_name = tendency_of_air_temperature + units = K s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out diff --git a/schemes/sima_diagnostics/check_energy_diagnostics.F90 b/schemes/sima_diagnostics/check_energy_diagnostics.F90 new file mode 100644 index 00000000..be52eb8d --- /dev/null +++ b/schemes/sima_diagnostics/check_energy_diagnostics.F90 @@ -0,0 +1,86 @@ +! Diagnostic scheme for check_energy +! Not all quantities are needed as diagnostics; this module is designed to ease development +module check_energy_diagnostics + use ccpp_kinds, only: kind_phys + + implicit none + private + save + + public :: check_energy_diagnostics_init ! init routine + public :: check_energy_diagnostics_run ! main routine + +CONTAINS + + !> \section arg_table_check_energy_diagnostics_init Argument Table + !! \htmlinclude check_energy_diagnostics_init.html + subroutine check_energy_diagnostics_init(errmsg, errflg) + use cam_history, only: history_add_field + use cam_history_support, only: horiz_only + + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Local variables: + + errmsg = '' + errflg = 0 + + ! History add field calls + call history_add_field('cp_or_cv_dycore', 'specific_heat_of_air_used_in_dycore', 'lev', 'inst', 'J kg-1 K-1') + call history_add_field('scaling_dycore', 'ratio_of_specific_heat_of_air_used_in_physics_energy_formula_to_specific_heat_of_air_used_in_dycore_energy_formula', 'lev', 'inst', '1') + + call history_add_field('te_cur_phys', 'vertically_integrated_total_energy_using_physics_energy_formula', horiz_only, 'inst', 'J m-2') + call history_add_field('te_cur_dyn', 'vertically_integrated_total_energy_using_dycore_energy_formula', horiz_only, 'inst', 'J m-2') + call history_add_field('tw_cur', 'vertically_integrated_total_water', horiz_only, 'inst', 'kg m-2') + + call history_add_field('tend_te_tnd', 'cumulative_total_energy_boundary_flux_using_physics_energy_formula', horiz_only, 'inst', 'J m-2 s-1') + call history_add_field('tend_tw_tnd', 'cumulative_total_water_boundary_flux', horiz_only, 'inst', 'kg m-2 s-1') + + call history_add_field('teout', 'vertically_integrated_total_energy_at_end_of_physics_timestep', horiz_only, 'inst', 'J m-2') + + end subroutine check_energy_diagnostics_init + + !> \section arg_table_check_energy_diagnostics_run Argument Table + !! \htmlinclude check_energy_diagnostics_run.html + subroutine check_energy_diagnostics_run( & + cp_or_cv_dycore, scaling_dycore, & + te_cur_phys, te_cur_dyn, tw_cur, & + tend_te_tnd, tend_tw_tnd, teout, & + errmsg, errflg) + + use cam_history, only: history_out_field + !------------------------------------------------ + ! Input / output parameters + !------------------------------------------------ + ! State variables + real(kind_phys), intent(in) :: cp_or_cv_dycore(:,:) + real(kind_phys), intent(in) :: scaling_dycore(:,:) + real(kind_phys), intent(in) :: te_cur_phys(:) + real(kind_phys), intent(in) :: te_cur_dyn(:) + real(kind_phys), intent(in) :: tw_cur(:) + real(kind_phys), intent(in) :: tend_te_tnd(:) + real(kind_phys), intent(in) :: tend_tw_tnd(:) + real(kind_phys), intent(in) :: teout(:) + + + ! CCPP error handling variables + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + ! History out field calls + call history_out_field('cp_or_cv_dycore', cp_or_cv_dycore) + call history_out_field('scaling_dycore', scaling_dycore) + call history_out_field('te_cur_phys', te_cur_phys) + call history_out_field('te_cur_dyn', te_cur_dyn) + call history_out_field('tw_cur', tw_cur) + call history_out_field('tend_te_tnd', tend_te_tnd) + call history_out_field('tend_tw_tnd', tend_tw_tnd) + call history_out_field('teout', teout) + + end subroutine check_energy_diagnostics_run + +end module check_energy_diagnostics diff --git a/schemes/sima_diagnostics/check_energy_diagnostics.meta b/schemes/sima_diagnostics/check_energy_diagnostics.meta new file mode 100644 index 00000000..cd5e1532 --- /dev/null +++ b/schemes/sima_diagnostics/check_energy_diagnostics.meta @@ -0,0 +1,83 @@ +[ccpp-table-properties] + name = check_energy_diagnostics + type = scheme + +[ccpp-arg-table] + name = check_energy_diagnostics_init + type = scheme +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = check_energy_diagnostics_run + type = scheme +[ cp_or_cv_dycore ] + standard_name = specific_heat_of_air_used_in_dycore + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ scaling_dycore ] + standard_name = ratio_of_specific_heat_of_air_used_in_physics_energy_formula_to_specific_heat_of_air_used_in_dycore_energy_formula + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ te_cur_phys ] + standard_name = vertically_integrated_total_energy_using_physics_energy_formula + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ te_cur_dyn ] + standard_name = vertically_integrated_total_energy_using_dycore_energy_formula + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ tw_cur ] + standard_name = vertically_integrated_total_water + units = kg m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ tend_te_tnd ] + standard_name = cumulative_total_energy_boundary_flux_using_physics_energy_formula + units = J m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ tend_tw_tnd ] + standard_name = cumulative_total_water_boundary_flux + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ teout ] + standard_name = vertically_integrated_total_energy_at_end_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/sima_diagnostics/check_energy_gmean_diagnostics.F90 b/schemes/sima_diagnostics/check_energy_gmean_diagnostics.F90 new file mode 100644 index 00000000..cd582b2e --- /dev/null +++ b/schemes/sima_diagnostics/check_energy_gmean_diagnostics.F90 @@ -0,0 +1,61 @@ +! Diagnostic scheme for check_energy_gmean +! This creates a debug printout with: +! - global mean input energy using dycore energy formula +! - global mean output energy at end of physics timestep using dycore energy formula +! - global mean heating rate for energy conservation +! These numbers are very useful for matching models bit-for-bit because they include "everything" in the model. +module check_energy_gmean_diagnostics + use ccpp_kinds, only: kind_phys + + implicit none + private + save + + public :: check_energy_gmean_diagnostics_init + public :: check_energy_gmean_diagnostics_run ! main routine + + ! Private module options. + logical :: print_global_means = .false. + +contains + +!> \section arg_table_check_energy_gmean_diagnostics_init Argument Table +!! \htmlinclude arg_table_check_energy_gmean_diagnostics_init.html + subroutine check_energy_gmean_diagnostics_init(print_global_means_in) + ! Input arguments + logical, intent(in) :: print_global_means_in + + print_global_means = print_global_means_in + end subroutine check_energy_gmean_diagnostics_init + + !> \section arg_table_check_energy_gmean_diagnostics_run Argument Table + !! \htmlinclude check_energy_gmean_diagnostics_run.html + subroutine check_energy_gmean_diagnostics_run( & + amIRoot, & + iulog, & + teinp_glob, & + teout_glob, & + heat_glob, & + errmsg, errflg) + + logical, intent(in) :: amIRoot ! are we on the MPI root task? + integer, intent(in) :: iulog ! log output unit + + real(kind_phys), intent(in) :: teinp_glob ! global mean energy of input state [J m-2] + real(kind_phys), intent(in) :: teout_glob ! global mean energy of output state [J m-2] + real(kind_phys), intent(in) :: heat_glob ! global mean heating rate [J kg-1 s-1] + + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + if (print_global_means .and. amIRoot) then + write(iulog,'(1x,a26,1x,e25.17,1x,a15,1x,e25.17)') 'global mean input energy =', teinp_glob, 'output energy =', teout_glob + write(iulog,'(1x,a38,1x,e25.17)') 'heating rate for energy conservation =', heat_glob + endif + + end subroutine check_energy_gmean_diagnostics_run + +end module check_energy_gmean_diagnostics diff --git a/schemes/sima_diagnostics/check_energy_gmean_diagnostics.meta b/schemes/sima_diagnostics/check_energy_gmean_diagnostics.meta new file mode 100644 index 00000000..7eff0709 --- /dev/null +++ b/schemes/sima_diagnostics/check_energy_gmean_diagnostics.meta @@ -0,0 +1,59 @@ +[ccpp-table-properties] + name = check_energy_gmean_diagnostics + type = scheme + +[ccpp-arg-table] + name = check_energy_gmean_diagnostics_init + type = scheme +[ print_global_means_in ] + standard_name = flag_for_energy_global_means_output + units = flag + type = logical + dimensions = () + intent = in + +[ccpp-arg-table] + name = check_energy_gmean_diagnostics_run + type = scheme +[ amIRoot ] + standard_name = flag_for_mpi_root + units = flag + type = logical + dimensions = () + intent = in +[ iulog ] + standard_name = log_output_unit + units = 1 + type = integer + dimensions = () + intent = in +[ teinp_glob ] + standard_name = global_mean_vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ teout_glob ] + standard_name = global_mean_vertically_integrated_total_energy_at_end_of_physics_timestep + units = J m-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ heat_glob ] + standard_name = global_mean_heating_rate_correction_for_energy_conservation + units = J kg-1 s-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/sima_diagnostics/check_energy_gmean_diagnostics_namelist.xml b/schemes/sima_diagnostics/check_energy_gmean_diagnostics_namelist.xml new file mode 100644 index 00000000..cd5896f6 --- /dev/null +++ b/schemes/sima_diagnostics/check_energy_gmean_diagnostics_namelist.xml @@ -0,0 +1,19 @@ + + + + + + + logical + diagnostics + check_energy_gmean_nl + flag_for_energy_global_means_output + flag + + Turn on output of global means of total energy and heating rate for energy conservation. Default: TRUE + + + true + + + diff --git a/suites/suite_adiabatic.xml b/suites/suite_adiabatic.xml new file mode 100644 index 00000000..eda9d892 --- /dev/null +++ b/suites/suite_adiabatic.xml @@ -0,0 +1,30 @@ + + + + + + check_energy_gmean + + check_energy_gmean_diagnostics + + + check_energy_zero_fluxes + check_energy_fix + apply_heating_rate + geopotential_temp + + + check_energy_scaling + check_energy_chng + + + check_energy_save_teout + + + sima_state_diagnostics + + + + sima_tend_diagnostics + + diff --git a/suites/suite_cam7.xml b/suites/suite_cam7.xml index 1bacbe94..7351223d 100644 --- a/suites/suite_cam7.xml +++ b/suites/suite_cam7.xml @@ -2,21 +2,48 @@ + + check_energy_gmean + + check_energy_gmean_diagnostics + + + check_energy_zero_fluxes + check_energy_fix + apply_heating_rate + geopotential_temp + + + check_energy_scaling + check_energy_chng + dadadj dadadj_apply_qv_tendency apply_heating_rate qneg geopotential_temp - - sima_state_diagnostics + + sima_state_diagnostics + tropopause_find tropopause_diagnostics + + + check_energy_save_teout + + + check_energy_scaling + dycore_energy_consistency_adjust + apply_tendency_of_air_temperature + sima_tend_diagnostics diff --git a/suites/suite_kessler.xml b/suites/suite_kessler.xml index d8d58adb..792c6ee3 100644 --- a/suites/suite_kessler.xml +++ b/suites/suite_kessler.xml @@ -16,10 +16,31 @@ kessler_update qneg geopotential_temp + + + check_energy_zero_fluxes + check_energy_scaling + check_energy_chng + + sima_state_diagnostics kessler_diagnostics + + + + check_energy_scaling + dycore_energy_consistency_adjust + apply_tendency_of_air_temperature + + sima_tend_diagnostics diff --git a/suites/suite_tj2016.xml b/suites/suite_tj2016.xml index 8aa50be6..a99bca79 100644 --- a/suites/suite_tj2016.xml +++ b/suites/suite_tj2016.xml @@ -5,6 +5,16 @@ tj2016_precip apply_heating_rate qneg + + + check_energy_zero_fluxes + check_energy_scaling + check_energy_chng + + sima_state_diagnostics @@ -13,6 +23,18 @@ apply_tendency_of_eastward_wind apply_tendency_of_northward_wind qneg + + + + + check_energy_scaling + dycore_energy_consistency_adjust + apply_tendency_of_air_temperature + + sima_tend_diagnostics