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