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/.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/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/musica/micm/musica_ccpp_micm.F90 b/schemes/musica/micm/musica_ccpp_micm.F90 index 2419d40e..7b9fbae9 100644 --- a/schemes/musica/micm/musica_ccpp_micm.F90 +++ b/schemes/musica/micm/musica_ccpp_micm.F90 @@ -2,21 +2,22 @@ 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 - !> 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 @@ -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,10 +74,11 @@ 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 - !> Intitialize MICM + !> Initializes MICM subroutine micm_init(errmsg, errcode) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -86,45 +88,42 @@ 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 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 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/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 7b76ee55..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 @@ -10,6 +12,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 @@ -26,76 +30,72 @@ 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 + 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, & - 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 !> \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) - 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 + surface_geopotential, surface_albedo, & + standard_gravitational_acceleration, errmsg, errcode) 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) :: 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_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() @@ -120,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.meta b/schemes/musica/musica_ccpp.meta index e4a38460..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 @@ -19,6 +53,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 +135,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 @@ -128,4 +174,4 @@ units = 1 type = integer dimensions = () - intent = out \ No newline at end of file + intent = out 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 9ac1c553..5a12a3ee 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -1,42 +1,60 @@ 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 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 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() + + type(index_mappings_t), pointer :: photolysis_rate_constants_mapping => null( ) + integer :: number_of_photolysis_rate_constants = 0 contains - !> Intitialize TUV-x + !> Initializes TUV-x subroutine tuvx_init(vertical_layer_dimension, vertical_interface_dimension, & + 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_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_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 + 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 + 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 @@ -50,105 +68,140 @@ 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, & + tuvx => tuvx_t( trim(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 )) return + + height_grid => grids%get( height_grid_label, height_grid_unit, error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( tuvx ) - tuvx => null() + call tuvx_deallocate( grids, null(), null(), tuvx, null(), null(), & + null(), null() ) return end if - height_grid => grids%get( height_grid_label, height_grid_unit, error ) + wavelength_grid => grids%get( wavelength_grid_label, wavelength_grid_unit, & + error ) if (has_error_occurred( error, errmsg, errcode )) then - deallocate( tuvx ) - tuvx => null() - deallocate( grids ) + call tuvx_deallocate( grids, null(), null(), tuvx, height_grid, null(), & + null(), null() ) return end if - deallocate( grids ) - 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() ) + + 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 @@ -157,33 +210,45 @@ 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) + 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 - - 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) :: 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 + 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 + 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,:), & @@ -197,14 +262,32 @@ 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 - ! 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 - !> Finalize tuvx + !> Finalizes TUV-x subroutine tuvx_final(errmsg, errcode) character(len=512), intent(out) :: errmsg integer, intent(out) :: errcode @@ -217,16 +300,31 @@ 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() 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/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/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 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 5581d5d8..3b5cf4fa 100644 --- a/test/musica/micm/test_micm_util.F90 +++ b/test/musica/micm/test_micm_util.F90 @@ -7,87 +7,23 @@ 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 + subroutine test_unit_conversion() + 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 = 4 - real(kind_phys) :: temperature(NUM_COLUMNS,NUM_LAYERS) - real(kind_phys) :: pressure(NUM_COLUMNS,NUM_LAYERS) + 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(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) - - ! 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 - - 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 /) - 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)) - end do - - do i_elem = 1, size(micm_constituents) - ASSERT_NEAR(micm_constituents(i_elem), arr_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), arr_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 - - 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 + 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 +48,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 +58,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/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 ac37e30c..61cca2f7 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -5,48 +5,182 @@ 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() + !> 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 + + real(kind_phys), dimension(:), intent(out) :: edges + + edges = (/ & + 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 & + /) + + 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 = 4 + 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 :: 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 + 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 - num_grid_cells = NUM_COLUMNS * NUM_LAYERS 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 /) @@ -54,6 +188,7 @@ subroutine test_musica_ccpp_api() 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 /) @@ -61,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) @@ -78,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 @@ -98,12 +437,33 @@ 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 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" @@ -113,10 +473,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 +493,22 @@ subroutine test_musica_ccpp_api() stop 3 endif - end subroutine test_musica_ccpp_api + ! 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_terminator 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 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