From 0e417de23ae73a1f3dac7aa8de07ad839a568eb1 Mon Sep 17 00:00:00 2001 From: Courtney Peverley Date: Tue, 4 Feb 2025 11:23:14 -0700 Subject: [PATCH 1/6] Add diagnostic output for ALL constituents (#199) Originator(s): peverwhee Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): The sima_state_diagnostics scheme only output the "known" constituents. This PR adds history add and out field calls for ALL other constituents. I'm open to discussion about the diagnostic names for these! Describe any changes made to the namelist: n/a List all files eliminated and why: n/a List all files added and what they do: n/a List all existing files that have been modified, and describe the changes: (Helpful git command: `git diff --name-status development...`) M schemes/sima_diagnostics/sima_state_diagnostics.F90 - add "history_add_field" and "history_out_field" calls for constituents outside of the core expected ones. List any test failures: all CAM-SIMA tests pass Is this a science-changing update? New physics package, algorithm change, tuning changes, etc? no - diagnostics update --- .../sima_state_diagnostics.F90 | 68 ++++++++++++++++++- 1 file changed, 66 insertions(+), 2 deletions(-) diff --git a/schemes/sima_diagnostics/sima_state_diagnostics.F90 b/schemes/sima_diagnostics/sima_state_diagnostics.F90 index 42786f07..1d1b5500 100644 --- a/schemes/sima_diagnostics/sima_state_diagnostics.F90 +++ b/schemes/sima_diagnostics/sima_state_diagnostics.F90 @@ -2,6 +2,7 @@ module sima_state_diagnostics use ccpp_kinds, only: kind_phys use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use cam_history_support, only: fieldname_len implicit none private @@ -41,7 +42,13 @@ subroutine sima_state_diagnostics_init(const_props, errmsg, errflg) integer :: const_idx, name_idx integer :: const_num_found + character(len=64) :: units + logical :: is_wet character(len=512) :: standard_name + character(len=fieldname_len) :: diagnostic_name + character(len=512) :: long_name + character(len=3) :: mixing_ratio_type + logical :: const_found(size(const_props)) errmsg = '' errflg = 0 @@ -70,14 +77,20 @@ subroutine sima_state_diagnostics_init(const_props, errmsg, errflg) call history_add_field('LNPINT', 'ln_air_pressure_at_interfaces', 'ilev', 'avg', '1') call history_add_field('LNPINTDRY', 'ln_air_pressure_of_dry_air_at_interfaces', 'ilev', 'avg', '1') call history_add_field('ZI', 'geopotential_height_wrt_surface_at_interfaces', 'ilev', 'avg', 'm') - ! Add constituent fields + + ! Add expected constituent fields const_num_found = 0 + const_found = .false. do const_idx = 1, size(const_props) call const_props(const_idx)%standard_name(standard_name, errflg, errmsg) + if (errflg /= 0) then + return + end if do name_idx = 1, size(const_std_names) if (trim(standard_name) == trim(const_std_names(name_idx))) then call history_add_field(trim(const_diag_names(name_idx)), trim(const_std_names(name_idx)), 'lev', 'avg', 'kg kg-1', mixing_ratio='wet') const_num_found = const_num_found + 1 + const_found(const_idx) = .true. end if end do if (const_num_found == size(const_std_names)) then @@ -85,6 +98,36 @@ subroutine sima_state_diagnostics_init(const_props, errmsg, errflg) end if end do + ! Add fields for all other constituents + do const_idx = 1, size(const_props) + if (.not. const_found(const_idx)) then + call const_props(const_idx)%standard_name(standard_name, errflg, errmsg) + if (errflg /= 0) then + return + end if + ! truncate the standard name if necessary + diagnostic_name = standard_name + call const_props(const_idx)%units(units, errflg, errmsg) + if (errflg /= 0) then + return + end if + call const_props(const_idx)%is_wet(is_wet, errflg, errmsg) + if (errflg /= 0) then + return + end if + call const_props(const_idx)%long_name(long_name, errflg, errmsg) + if (errflg /= 0) then + return + end if + if (is_wet) then + mixing_ratio_type = 'wet' + else + mixing_ratio_type = 'dry' + end if + call history_add_field(trim(diagnostic_name), trim(long_name), 'lev', 'avg', trim(units), mixing_ratio=mixing_ratio_type) + end if + end do + end subroutine sima_state_diagnostics_init !> \section arg_table_sima_state_diagnostics_run Argument Table @@ -130,8 +173,10 @@ subroutine sima_state_diagnostics_run(ps, psdry, phis, T, u, v, dse, omega, & integer, intent(out) :: errflg character(len=512) :: standard_name + character(len=fieldname_len) :: diagnostic_name integer :: const_idx, name_idx integer :: const_num_found + logical :: const_found(size(const_props)) errmsg = '' errflg = 0 @@ -161,14 +206,21 @@ subroutine sima_state_diagnostics_run(ps, psdry, phis, T, u, v, dse, omega, & call history_out_field('LNPINTDRY', lnpintdry) call history_out_field('ZI' , zi) - ! Capture constituent fields + ! Capture expected constituent fields const_num_found = 0 + const_found = .false. do const_idx = 1, size(const_props) call const_props(const_idx)%standard_name(standard_name, errflg, errmsg) + if (errflg /= 0) then + return + end if do name_idx = 1, size(const_std_names) if (trim(standard_name) == trim(const_std_names(name_idx))) then call history_out_field(trim(const_diag_names(name_idx)), const_array(:,:,const_idx)) const_num_found = const_num_found + 1 + const_found(const_idx) = .true. + else + call history_out_field(trim(standard_name), const_array(:,:,const_idx)) end if end do if (const_num_found == size(const_std_names)) then @@ -176,5 +228,17 @@ subroutine sima_state_diagnostics_run(ps, psdry, phis, T, u, v, dse, omega, & end if end do + ! Capture all other constituent fields + do const_idx = 1, size(const_props) + if (.not. const_found(const_idx)) then + call const_props(const_idx)%standard_name(standard_name, errflg, errmsg) + if (errflg /= 0) then + return + end if + diagnostic_name = standard_name + call history_out_field(trim(diagnostic_name), const_array(:,:,const_idx)) + end if + end do + end subroutine sima_state_diagnostics_run end module sima_state_diagnostics From 60b71f3fd9f10b362caf64afffce12de10854d4d Mon Sep 17 00:00:00 2001 From: Jiwon Gim <55209567+boulderdaze@users.noreply.github.com> Date: Tue, 4 Feb 2025 13:21:17 -0700 Subject: [PATCH 2/6] Build MUSICA in CAM-SIMA (#198) Originator(s): @boulderdaze Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): - Closes #196 - Fixed bugs for the CAM-SIMA build Describe any changes made to the namelist: ``` M schemes/musica/musica_ccpp_namelist.xml -> Specified the path to directory for the configuration ``` List all files eliminated and why: N/A List all files added and what they do: N/A List all existing files that have been modified, and describe the changes: ``` M schemes/musica/musica_ccpp.F90 M schemes/musica/musica_ccpp.meta M schemes/musica/musica_ccpp_namelist.xml M schemes/musica/tuvx/musica_ccpp_tuvx.F90 M schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 M test/musica/test_musica_api.F90 ``` List any test failures: N/A Is this a science-changing update? New physics package, algorithm change, tuning changes, etc? No --- schemes/musica/musica_ccpp.F90 | 24 +++++++++---------- schemes/musica/musica_ccpp.meta | 9 +++++-- schemes/musica/musica_ccpp_namelist.xml | 6 ++--- schemes/musica/tuvx/musica_ccpp_tuvx.F90 | 23 +++++++++--------- .../tuvx/musica_ccpp_tuvx_surface_albedo.F90 | 6 ++--- test/musica/test_musica_api.F90 | 12 +++++----- 6 files changed, 43 insertions(+), 37 deletions(-) diff --git a/schemes/musica/musica_ccpp.F90 b/schemes/musica/musica_ccpp.F90 index 70d9791d..27a8897a 100644 --- a/schemes/musica/musica_ccpp.F90 +++ b/schemes/musica/musica_ccpp.F90 @@ -124,22 +124,22 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co extract_subset_constituents, update_constituents 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) :: temperature(:,:) ! K (column, layer) + real(kind_phys), target, intent(in) :: pressure(:,:) ! Pa (column, layer) + real(kind_phys), target, intent(in) :: dry_air_density(:,:) ! kg m-3 (column, layer) type(ccpp_constituent_prop_ptr_t), & - intent(in) :: constituent_props(:) - real(kind_phys), target, intent(inout) :: constituents(:,:,:) ! kg kg-1 + intent(in) :: constituent_props(:) ! (constituent) + real(kind_phys), target, intent(inout) :: constituents(:,:,:) ! kg kg-1 (column, layer, constituent) 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_geopotential(:) ! m2 s-2 - real(kind_phys), intent(in) :: surface_temperature(:) ! K - real(kind_phys), intent(in) :: surface_albedo ! unitless - real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! nm - real(kind_phys), intent(in) :: extraterrestrial_flux(:) ! photons cm-2 s-1 nm-1 + real(kind_phys), intent(in) :: surface_geopotential(:) ! m2 s-2 (column) + real(kind_phys), intent(in) :: surface_temperature(:) ! K (column) + real(kind_phys), intent(in) :: surface_albedo(:) ! fraction (column) + real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! nm (wavelength interface) + real(kind_phys), intent(in) :: extraterrestrial_flux(:) ! photons cm-2 s-1 nm-1 (wavelength interface) real(kind_phys), intent(in) :: standard_gravitational_acceleration ! m s-2 - real(kind_phys), intent(in) :: cloud_area_fraction(:,:) ! unitless (column, level) - real(kind_phys), intent(in) :: air_pressure_thickness(:,:) ! Pa (column, level) + real(kind_phys), intent(in) :: cloud_area_fraction(:,:) ! fraction (column, layer) + real(kind_phys), intent(in) :: air_pressure_thickness(:,:) ! Pa (column, layer) real(kind_phys), intent(in) :: solar_zenith_angle(:) ! radians (column) real(kind_phys), intent(in) :: earth_sun_distance ! AU character(len=512), intent(out) :: errmsg diff --git a/schemes/musica/musica_ccpp.meta b/schemes/musica/musica_ccpp.meta index bad6f271..bb251532 100644 --- a/schemes/musica/musica_ccpp.meta +++ b/schemes/musica/musica_ccpp.meta @@ -1,7 +1,12 @@ [ccpp-table-properties] 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 + dependencies = musica_ccpp.F90,musica_ccpp_util.F90,musica_ccpp_species.F90 + dependencies = micm/musica_ccpp_micm.F90,micm/musica_ccpp_micm_util.F90 + dependencies = tuvx/musica_ccpp_tuvx_aerosol_optics.F90,tuvx/musica_ccpp_tuvx_cloud_optics.F90,tuvx/musica_ccpp_tuvx_extraterrestrial_flux.F90 + dependencies = tuvx/musica_ccpp_tuvx_gas_species.F90,tuvx/musica_ccpp_tuvx_height_grid.F90,tuvx/musica_ccpp_tuvx_load_species.F90 + dependencies = tuvx/musica_ccpp_tuvx_surface_albedo.F90,tuvx/musica_ccpp_tuvx_temperature.F90,tuvx/musica_ccpp_tuvx_wavelength_grid.F90 + dependencies = tuvx/musica_ccpp_tuvx.F90 [ccpp-arg-table] name = musica_ccpp_register @@ -139,7 +144,7 @@ standard_name = surface_albedo_due_to_UV_and_VIS_direct type = real | kind = kind_phys units = fraction - dimensions = () + dimensions = (horizontal_loop_extent) intent = in [ photolysis_wavelength_grid_interfaces ] standard_name = photolysis_wavelength_grid_interfaces diff --git a/schemes/musica/musica_ccpp_namelist.xml b/schemes/musica/musica_ccpp_namelist.xml index c36d1b81..88ac00e6 100644 --- a/schemes/musica/musica_ccpp_namelist.xml +++ b/schemes/musica/musica_ccpp_namelist.xml @@ -99,7 +99,7 @@ A configuration file for the MICM chemical solver. - UNSET_PATH + ${CASEROOT} @@ -112,7 +112,7 @@ A configuration file for the TUV-x photolysis rate calculator - UNSET_PATH + ${CASEROOT} @@ -125,7 +125,7 @@ A configuration file for the mapping of TUV-x photolysis rates to MICM custom rate parameters - UNSET_PATH + ${CASEROOT} \ No newline at end of file diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 index 1732ffab..293a2905 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx.F90 @@ -533,15 +533,15 @@ subroutine tuvx_run(temperature, dry_air_density, & real(kind_phys), intent(in) :: constituents(:,:,:) ! kg kg-1 (column, layer, constituent) 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_geopotential(:) ! m2 s-2 - real(kind_phys), intent(in) :: surface_temperature(:) ! K - real(kind_phys), intent(in) :: surface_albedo ! unitless - real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! nm - real(kind_phys), intent(in) :: extraterrestrial_flux(:) ! photons cm-2 s-1 nm-1 + real(kind_phys), intent(in) :: surface_geopotential(:) ! m2 s-2 (column) + real(kind_phys), intent(in) :: surface_temperature(:) ! K (column) + real(kind_phys), intent(in) :: surface_albedo(:) ! fraction (column) + real(kind_phys), intent(in) :: photolysis_wavelength_grid_interfaces(:) ! nm (wavelength interface) + real(kind_phys), intent(in) :: extraterrestrial_flux(:) ! photons cm-2 s-1 nm-1 (wavelength interface) real(kind_phys), intent(in) :: standard_gravitational_acceleration ! m s-2 - real(kind_phys), intent(in) :: cloud_area_fraction(:,:) ! unitless (column, layer) + real(kind_phys), intent(in) :: cloud_area_fraction(:,:) ! fraction (column, layer) real(kind_phys), intent(in) :: air_pressure_thickness(:,:) ! Pa (column, layer) - real(kind_phys), intent(in) :: solar_zenith_angle(:) ! radians + real(kind_phys), intent(in) :: solar_zenith_angle(:) ! radians (column) real(kind_phys), intent(in) :: earth_sun_distance ! m real(kind_phys), intent(inout) :: rate_parameters(:,:,:) ! various units (column, layer, reaction) character(len=512), intent(out) :: errmsg @@ -561,10 +561,6 @@ subroutine tuvx_run(temperature, dry_air_density, & 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 - call set_extraterrestrial_flux_values( extraterrestrial_flux_profile, & photolysis_wavelength_grid_interfaces, & extraterrestrial_flux, errmsg, errcode ) @@ -587,6 +583,11 @@ subroutine tuvx_run(temperature, dry_air_density, & height_deltas, errmsg, errcode ) if (errcode /= 0) return + ! surface albedo with respect to direct UV/visible radiation + call set_surface_albedo_values( surface_albedo_profile, surface_albedo(i_col), & + errmsg, errcode ) + if (errcode /= 0) return + call set_temperature_values( temperature_profile, temperature(i_col,:), & surface_temperature(i_col), errmsg, errcode ) if (errcode /= 0) return diff --git a/schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 b/schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 index 8608a12d..a52894c3 100644 --- a/schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 +++ b/schemes/musica/tuvx/musica_ccpp_tuvx_surface_albedo.F90 @@ -43,7 +43,7 @@ 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, & + subroutine set_surface_albedo_values( profile, surface_albedo, & errmsg, errcode ) use ccpp_kinds, only: kind_phys use musica_ccpp_util, only: has_error_occurred @@ -51,7 +51,7 @@ subroutine set_surface_albedo_values( profile, host_surface_albedo, & use musica_util, only: error_t type(profile_t), intent(inout) :: profile - real(kind_phys), intent(in) :: host_surface_albedo ! unitless + real(kind_phys), intent(in) :: surface_albedo ! fraction character(len=*), intent(out) :: errmsg integer, intent(out) :: errcode @@ -65,7 +65,7 @@ subroutine set_surface_albedo_values( profile, host_surface_albedo, & return end if - surface_albedo_interfaces(:) = host_surface_albedo + surface_albedo_interfaces(:) = surface_albedo call profile%set_edge_values( surface_albedo_interfaces, error ) if ( has_error_occurred( error, errmsg, errcode ) ) return diff --git a/test/musica/test_musica_api.F90 b/test/musica/test_musica_api.F90 index 76af8099..4ffae4a8 100644 --- a/test/musica/test_musica_api.F90 +++ b/test/musica/test_musica_api.F90 @@ -169,7 +169,7 @@ subroutine test_chapman() real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m real(kind_phys), dimension(NUM_COLUMNS) :: surface_geopotential ! m2 s-2 real(kind_phys), dimension(NUM_COLUMNS) :: surface_temperature ! K - real(kind_phys) :: surface_albedo ! unitless + real(kind_phys), dimension(NUM_COLUMNS) :: surface_albedo ! fraction integer, parameter :: num_photolysis_wavelength_grid_sections = 8 ! (count) real(kind_phys), dimension(num_photolysis_wavelength_grid_sections+1) :: flux_data_photolysis_wavelength_interfaces ! nm real(kind_phys), dimension(num_photolysis_wavelength_grid_sections) :: extraterrestrial_flux ! photons cm-2 s-1 nm-1 @@ -177,7 +177,7 @@ subroutine test_chapman() 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) :: cloud_area_fraction ! unitless + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: cloud_area_fraction ! fraction real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: air_pressure_thickness ! Pa real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & NUM_SPECIES+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) :: constituents ! kg kg-1 @@ -204,7 +204,7 @@ subroutine test_chapman() 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 + 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 /) @@ -411,7 +411,7 @@ subroutine test_terminator() real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m real(kind_phys), dimension(NUM_COLUMNS) :: surface_geopotential ! m2 s-2 real(kind_phys), dimension(NUM_COLUMNS) :: surface_temperature ! K - real(kind_phys) :: surface_albedo ! unitless + real(kind_phys), dimension(NUM_COLUMNS) :: surface_albedo ! fraction integer, parameter :: num_photolysis_wavelength_grid_sections = 8 ! (count) real(kind_phys), dimension(num_photolysis_wavelength_grid_sections+1) :: flux_data_photolysis_wavelength_interfaces ! nm real(kind_phys), dimension(num_photolysis_wavelength_grid_sections) :: extraterrestrial_flux ! photons cm-2 s-1 nm-1 @@ -419,7 +419,7 @@ subroutine test_terminator() 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) :: cloud_area_fraction ! unitless + real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: cloud_area_fraction ! fraction real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS) :: air_pressure_thickness ! Pa real(kind_phys), dimension(NUM_COLUMNS,NUM_LAYERS, & NUM_SPECIES+NUM_TUVX_CONSTITUENTS+NUM_TUVX_ONLY_GAS_SPECIES) :: constituents ! kg kg-1 @@ -446,7 +446,7 @@ subroutine test_terminator() 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 + 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 /) From 0c545134ae003e214a0dd152a91ae163b60b893c Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Thu, 13 Feb 2025 21:29:10 -0500 Subject: [PATCH 3/6] Implement thermo_water_update scheme (non-portable wrapper to cam_thermo_water_update) (#178) Originator(s): @jimmielin Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): Closes #177, implementing thermo_water_update scheme to call cam_thermo_water_update. Describe any changes made to the namelist: N/A List all files eliminated and why: N/A List all files added and what they do: Implement thermo_water_update scheme A schemes/thermo_water_update/thermo_water_update.F90 A schemes/thermo_water_update/thermo_water_update.meta List all existing files that have been modified, and describe the changes: (Helpful git command: `git diff --name-status development...`) Include thermo_water_update scheme in CAM7 suite and simple physics suites M suites/suite_cam7.xml M suites/suite_kessler.xml M suites/suite_tj2016.xml List any test failures: Is this a science-changing update? New physics package, algorithm change, tuning changes, etc? Yes --------- Co-authored-by: Jesse Nusbaumer --- .../thermo_water_update.F90 | 50 +++++++++++++++++++ .../thermo_water_update.meta | 50 +++++++++++++++++++ suites/suite_cam7.xml | 3 ++ suites/suite_kessler.xml | 4 ++ suites/suite_tj2016.xml | 3 ++ 5 files changed, 110 insertions(+) create mode 100644 schemes/thermo_water_update/thermo_water_update.F90 create mode 100644 schemes/thermo_water_update/thermo_water_update.meta diff --git a/schemes/thermo_water_update/thermo_water_update.F90 b/schemes/thermo_water_update/thermo_water_update.F90 new file mode 100644 index 00000000..d2623d6e --- /dev/null +++ b/schemes/thermo_water_update/thermo_water_update.F90 @@ -0,0 +1,50 @@ +! This is a non-portable wrapper subroutine for cam_thermo_water_update +! in the cam_thermo module. +module thermo_water_update + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: thermo_water_update_run + +contains + + ! Update water dependent properties +!> \section arg_table_thermo_water_update_run Argument Table +!! \htmlinclude arg_table_thermo_water_update_run.html + subroutine thermo_water_update_run( & + mmr, & + ncol, pver, & + energy_formula_dycore, & + pdel, pdeldry, & + cp_or_cv_dycore) + + ! This scheme is non-portable due to dependencies on cam_thermo + ! for the actual logic of cam_thermo_water_update, which depends on air_composition + ! and a series of other subroutines/module properties + use cam_thermo, only: cam_thermo_water_update + + ! Input arguments + real(kind_phys), intent(in) :: mmr(:,:,:) ! constituent mass mixing ratios [kg kg-1] + integer, intent(in) :: ncol ! number of atmospheric columns + integer, intent(in) :: pver ! number of vertical layers + integer, intent(in) :: energy_formula_dycore ! total energy formulation used by dycore + real(kind_phys), intent(in) :: pdel(:,:) ! layer thickness [Pa] + real(kind_phys), intent(in) :: pdeldry(:,:) ! dry layer thickness [Pa] + + ! Output arguments + real(kind_phys), intent(out) :: cp_or_cv_dycore(:,:) ! enthalpy or heat capacity, dycore dependent [J K-1 kg-1] + + call cam_thermo_water_update( & + mmr = mmr, & ! mmr*factor is a dry mixing ratio + ncol = ncol, & + pver = pver, & + energy_formula = energy_formula_dycore, & + cp_or_cv_dycore = cp_or_cv_dycore(:ncol,:), & + to_dry_factor = pdel(:ncol,:)/pdeldry(:ncol,:) & ! factor to convert to dry + ) + + end subroutine thermo_water_update_run + +end module thermo_water_update diff --git a/schemes/thermo_water_update/thermo_water_update.meta b/schemes/thermo_water_update/thermo_water_update.meta new file mode 100644 index 00000000..395b8ae5 --- /dev/null +++ b/schemes/thermo_water_update/thermo_water_update.meta @@ -0,0 +1,50 @@ +[ccpp-table-properties] + name = thermo_water_update + type = scheme + dependencies = ../../../../data/cam_thermo.F90 + +[ccpp-arg-table] + name = thermo_water_update_run + type = scheme +[ mmr ] + standard_name = ccpp_constituents + units = none + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents) + intent = in +[ 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 +[ energy_formula_dycore ] + standard_name = total_energy_formula_for_dycore + units = 1 + type = integer + dimensions = () + intent = in +[ pdel ] + standard_name = air_pressure_thickness + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ pdeldry ] + standard_name = air_pressure_thickness_of_dry_air + units = Pa + 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 = out diff --git a/suites/suite_cam7.xml b/suites/suite_cam7.xml index ae0daa42..b949654e 100644 --- a/suites/suite_cam7.xml +++ b/suites/suite_cam7.xml @@ -81,6 +81,9 @@ check_energy_save_teout + + thermo_water_update + diff --git a/suites/suite_kessler.xml b/suites/suite_kessler.xml index 792c6ee3..787e543d 100644 --- a/suites/suite_kessler.xml +++ b/suites/suite_kessler.xml @@ -30,6 +30,10 @@ kessler_diagnostics + + + thermo_water_update + + thermo_water_update + check_energy_gmean_diagnostics - + check_energy_zero_fluxes check_energy_fix apply_heating_rate geopotential_temp - - + check_energy_scaling check_energy_chng + + check_energy_fix_diagnostics check_energy_save_teout diff --git a/suites/suite_cam4.xml b/suites/suite_cam4.xml deleted file mode 100644 index c7832514..00000000 --- a/suites/suite_cam4.xml +++ /dev/null @@ -1,202 +0,0 @@ - - - - - - - physics_state_check - - - qneg3 - - physics_state_check - diag_state_b4_phys_write - - - calc_te_and_aam_budgets - check_energy_fix - check_energy_cam_update_pre_chng - check_energy_chng - check_energy_cam_update_post_chng - calc_te_and_aam_budgets - - diag_conv_tend_ini - calc_dtcore - - - - dadadj_calc_update - - - - - - zm_convr - zm_convr_cam_update - zm_conv_evap - zm_conv_evap_cam_update - momtran - momtran_cam_update - convtran - convtran_cam_update - - check_energy_chng - - - cmfmca - cmfmca_cam_update - zm_conv_evap - zm_conv_evap_cam_update - - check_energy_chng - check_tracers_chng - sslt_rebin_adv - - - tropopause_find - cld_sediment_vel - cld_sediment_tend - rk_sed_cam_update - add_conv_detrain_rliq - cldfrc - cldfrc_rh2 - rhdfda_calc - cldfrc_fice - strat_cond_repartition - pcond - rk_pcond_cam_update - cldfrc - cldfrc_diag_output - cldefr - - convtran - check_tracers_chng - diag_phys_writeout - diag_conv - cloud_diagnostics_calc - - - calc_solar_zenith_ang - group_day_night - calc_col_mean_co2 - eccf_calc - tropopause_find - calc_col_mean_o2 - aer_rad_props_sw - radcswmx - rad_sw_diag_output - aer_rad_props_lw - radclwmx - rad_lw_diag_output - cloud_cover_diags_output - rad_qdp_q_calc - rad_data_write - radheat_tend - radheat_diag_output - rad_q_qdp_calc - set_srf_net_sw - - check_energy_chng - tropopause_output - cam_export - diag_export - - - - - - flux_avg_run - physics_state_check - calc_flx_net - - - qneg4 - aoa_tracers_timestep_tend - check_tracers_chng - co2_cycle_set_ptend - - - - set_dry_to_wet - tint_calc - ubc_get_vals - set_top_tint_val - rhoi_calc - compute_tms - tms_diag_output - compute_blj - blj_diag_output - temp_to_potential_temp - virtem - calc_ustar - calc_obklen - - - trbintd - pblintd - austausch_atm - austausch_pbl - - pbl_diag_calc - dse_top_calc - compute_molec_diff_wet - compute_vdiff_wet - compute_molec_diff_dry - compute_vdiff_dry - diff_flux_diag - diff_flux_tend_dry - set_wet_to_dry - pbl_diag_calc - diff_mass_check - pbl_output - - rayleigh_friction_tend - check_energy_chng - check_tracers_chng - charge_balance - - - - set_dry_to_wet - gw_prof - - - effgw_calc - gw_src_calc - gw_drag_prof - calc_taucd - egwdffi_tot_calc - momentum_flux_calc - energy_change - energy_fixer - gw_cam_update - - - check_energy_chng - lunar_tides_tend - check_energy_chng - calc_te_and_aam_budgets - nudging_timestep_tend - check_energy_chng - set_dry_to_wet - - - qfac_state_adjust - geopotential_dse - - calc_te_and_aam_budgets - dtcore_reset - diag_phys_tend_writeout - - - diff --git a/suites/suite_cam6.xml b/suites/suite_cam6.xml deleted file mode 100644 index d995b70b..00000000 --- a/suites/suite_cam6.xml +++ /dev/null @@ -1,364 +0,0 @@ - - - - - - - physics_state_check - clybry_fam_adj - - - qneg - - physics_state_check - diag_state_b4_phys_write - - - calc_te_and_aam_budgets - check_energy_fix - check_energy_cam_update_pre_chng - check_energy_chng - check_energy_cam_update_post_chng - calc_te_and_aam_budgets - - diag_conv_tend_ini - calc_dtcore - - - - dadadj_calc_update - - - - zm_convr - zm_convr_cam_update - zm_conv_evap - zm_conv_evap_cam_update - momtran - momtran_cam_update - convtran - convtran_cam_update - - check_energy_chng - check_tracers_chng - sslt_rebin_adv - - - - - - clubb_input_prepare - clubb_var_vert_inv - - stats_begin_timestep_api - advance_clubb_core_api - update_xp2_mc_intr - calculate_thlp2_rad_intr - stats_end_timestep_clubb - - - - clubb_output_prepare - clubb_var_vert_rev - clubb_upper_diss - clubb_cam_update - liquid_macro_tend - liquid_macro_CAM_update - conv_cond_detrain_calc - conv_cond_detrain_cam_update - set_wet_to_dry - clubb_diag_output - - check_energy_chng - - - - - hetfrz_classnuc_cam_save_cbaero - aero_get_num_mmr - aero_calc_wsub - nucleate_ice_cam_calc - lcldm_min_check - aero_cam_drop_activate - aero_cam_contact_freezing - ndrop_bam_ccn - hetfrz_classnuc_cam_calc - microp_aero_diag_output - - - - micro_mg_get_cols3_0 - calc_incloud_LWP - micro_calc_tropopause - micro_mg_tend3_0 - mg_calc_outputs - calc_atm_density - - - size_dist_param_liq - micro_eff_radius_liq - size_dist_param_liq_const - calc_ncic_grid - size_dist_param_liq - micro_eff_radius_liq - micro_eff_radius_rain - micro_eff_radius_snow - micro_eff_radius_graupel - calc_niic_grid - size_dist_param_basic - micro_eff_radius_ice - - calc_micro_column_vars - calc_prec_efficiency - micro_diag_output - - massless_droplet_destroyer - - check_energy_chng - - diag_clip_tend_writeout - check_energy_chng - calc_prec_sum - - calc_prec_avg - modal_aero_calcsize_diag - - - - - modal_aero_wateruptake_dr - - - - modal_aero_calcsize_sub - modal_aero_wateruptake_dr - wetdep_prec_calc - coarse_fact_calc - - - aero_fact_calc - wetdepa_v2 - wetdep_diag_output - - - ma_convproc_intr - set_srf_wetdep - - convtran - check_tracers_chng - diag_phys_writeout - diag_conv - cloud_diagnostics_calc - - - calc_solar_zenith_ang - group_day_night - tropopause_find - - - rrtmg_state_prepare - - rrtmg_state_update - aer_rad_props_sw - rad_rrtmg_sw - rad_sw_diag_output - - rad_cnst_out - - rrtmg_state_update - aer_rad_props_lw - rad_rrtmg_lw - rad_lw_diag_output - - cosp_input_prepare - cosp_simulator - rad_qdp_q_calc - rad_data_write - radheat_tend - radheat_diag_output - rad_q_qdp_calc - set_srf_net_sw - - check_energy_chng - tropopause_output - cam_export - diag_export - - - - - - flux_avg_run - physics_state_check - calc_flx_net - - - aero_model_emissions - calc_MEGAN_fluxes - set_srf_emissions - srf_emis_diag_output - fire_emissions_srf - - - - qneg_surface - - aoa_tracers_timestep_tend - check_tracers_chng - co2_cycle_set_ptend - - - short_lived_species_writeic - get_curr_calday - chem_tropopause_find - neu_wetdep_tend - chem_calc_cldw - gas_phase_chemdr - nitro_srf_flx - chem_diag_output - - check_energy_chng - check_tracers_chng - - - set_dry_to_wet - tint_calc - ubc_get_vals - set_top_tint_val - rhoi_calc - compute_tms - tms_diag_output - compute_blj - blj_diag_output - theta_calc - virtem - calc_ustar - calc_obklen - pbl_diag_calc - dse_top_calc - compute_molec_diff_wet - compute_vdiff_wet - compute_molec_diff_dry - compute_vdiff_dry - aero_srf_flx_add - diff_flux_diag - diff_flux_tend_dry - set_wet_to_dry - pbl_diag_calc - diff_mass_check - pbl_output - - rayleigh_friction_tend - check_energy_chng - check_tracers_chng - - - calcram - calcram_diag_output - calc_atm_density - - modal_aero_depvel_part - - - - calc_aero_vars - - modal_aero_depvel_part - - aero_tracer_indx - depvel_m_to_pa - dust_sediment_tend - drydep_diag_output - - - set_srf_drydep - - charge_balance - - - - set_dry_to_wet - gw_prof - - - effgw_calc - gw_src_calc - gw_drag_prof - calc_taucd - egwdffi_tot_calc - momentum_flux_calc - energy_change - energy_fixer - gw_cam_update - - - check_energy_chng - lunar_tides_tend - check_energy_chng - calc_te_and_aam_budgets - nudging_timestep_tend - check_energy_chng - set_dry_to_wet - - - qfac_state_adjust - geopotential_dse - - calc_te_and_aam_budgets - dtcore_reset - diag_phys_tend_writeout - clybry_fam_set - - - diff --git a/suites/suite_cam6_silhs.xml b/suites/suite_cam6_silhs.xml deleted file mode 100644 index 9bd45239..00000000 --- a/suites/suite_cam6_silhs.xml +++ /dev/null @@ -1,379 +0,0 @@ - - - - - - - physics_state_check - clybry_fam_adj - - - qneg3 - - physics_state_check - diag_state_b4_phys_write - - - calc_te_and_aam_budgets - check_energy_fix - check_energy_cam_update_pre_chng - check_energy_chng - check_energy_cam_update_post_chng - calc_te_and_aam_budgets - - diag_conv_tend_ini - calc_dtcore - - - - dadadj_calc_update - - - - zm_convr - zm_convr_cam_update - zm_conv_evap - zm_conv_evap_cam_update - momtran - momtran_cam_update - convtran - convtran_cam_update - - check_energy_chng - check_tracers_chng - sslt_rebin_adv - - - - - - clubb_input_prepare - clubb_var_vert_inv - - stats_begin_timestep_api - advance_clubb_core_api - update_xp2_mc_intr - calculate_thlp2_rad_intr - stats_end_timestep_clubb - - - - clubb_output_prepare - clubb_var_vert_rev - clubb_upper_diss - clubb_cam_update - liquid_macro_tend - liquid_macro_CAM_update - conv_cond_detrain_calc - conv_cond_detrain_cam_update - set_wet_to_dry - clubb_diag_output - - check_energy_chng - - - - - hetfrz_classnuc_cam_save_cbaero - aero_get_num_mmr - aero_calc_wsub - nucleate_ice_cam_calc - lcldm_min_check - aero_cam_drop_activate - aero_cam_contact_freezing - ndrop_bam_ccn - hetfrz_classnuc_cam_calc - microp_aero_diag_output - - - - micro_mg_get_cols3_0 - calc_incloud_LWP - micro_calc_tropopause - micro_mg_tend3_0 - mg_calc_outputs - calc_atm_density - - - size_dist_param_liq - micro_eff_radius_liq - size_dist_param_liq_const - calc_ncic_grid - size_dist_param_liq - micro_eff_radius_liq - micro_eff_radius_rain - micro_eff_radius_snow - micro_eff_radius_graupel - calc_niic_grid - size_dist_param_basic - micro_eff_radius_ice - - calc_micro_column_vars - calc_prec_efficiency - micro_diag_output - - - - calc_total_moisture - calc_dry_static_density - convert_omega_to_w - convert_dse_to_temp - exner_silhs - thl_calc - silhs_var_vert_inv - silhs_get_subcol_wgts - lh_microphys_var_covar_driver_api - zero_upper_level - - subcol_SILHS_fill_holes_conserv_calc - massless_droplet_destroyer - subcol_SILHS_hydromet_conc_tend_lim - - check_energy_chng - - diag_clip_tend_writeout - check_energy_chng - calc_prec_sum - - calc_prec_avg - modal_aero_calcsize_diag - - - - - modal_aero_wateruptake_dr - - - - modal_aero_calcsize_sub - modal_aero_wateruptake_dr - wetdep_prec_calc - coarse_fact_calc - - - aero_fact_calc - wetdepa_v2 - wetdep_diag_output - - - ma_convproc_intr - set_srf_wetdep - - convtran - check_tracers_chng - diag_phys_writeout - diag_conv - cloud_diagnostics_calc - - - calc_solar_zenith_ang - group_day_night - tropopause_find - - - rrtmg_state_prepare - - rrtmg_state_update - aer_rad_props_sw - rad_rrtmg_sw - rad_sw_diag_output - - rad_cnst_out - - rrtmg_state_update - aer_rad_props_lw - rad_rrtmg_lw - rad_lw_diag_output - - cosp_input_prepare - cosp_simulator - rad_qdp_q_calc - rad_data_write - radheat_tend - radheat_diag_output - rad_q_qdp_calc - set_srf_net_sw - - check_energy_chng - tropopause_output - cam_export - diag_export - - - - - - flux_avg_run - physics_state_check - calc_flx_net - - - aero_model_emissions - calc_MEGAN_fluxes - set_srf_emissions - srf_emis_diag_output - fire_emissions_srf - - - - qneg4 - - aoa_tracers_timestep_tend - check_tracers_chng - co2_cycle_set_ptend - - - short_lived_species_writeic - get_curr_calday - chem_tropopause_find - neu_wetdep_tend - chem_calc_cldw - gas_phase_chemdr - nitro_srf_flx - chem_diag_output - - check_energy_chng - check_tracers_chng - - - set_dry_to_wet - tint_calc - ubc_get_vals - set_top_tint_val - rhoi_calc - compute_tms - tms_diag_output - compute_blj - blj_diag_output - theta_calc - virtem - calc_ustar - calc_obklen - pbl_diag_calc - dse_top_calc - compute_molec_diff_wet - compute_vdiff_wet - compute_molec_diff_dry - compute_vdiff_dry - aero_srf_flx_add - diff_flux_diag - diff_flux_tend_dry - set_wet_to_dry - pbl_diag_calc - diff_mass_check - pbl_output - - rayleigh_friction_tend - check_energy_chng - check_tracers_chng - - - calcram - calcram_diag_output - calc_atm_density - - modal_aero_depvel_part - - - - calc_aero_vars - - modal_aero_depvel_part - - aero_tracer_indx - depvel_m_to_pa - dust_sediment_tend - drydep_diag_output - - - set_srf_drydep - - charge_balance - - - - set_dry_to_wet - gw_prof - - - effgw_calc - gw_src_calc - gw_drag_prof - calc_taucd - egwdffi_tot_calc - momentum_flux_calc - energy_change - energy_fixer - gw_cam_update - - - check_energy_chng - lunar_tides_tend - check_energy_chng - calc_te_and_aam_budgets - nudging_timestep_tend - check_energy_chng - set_dry_to_wet - - - qfac_state_adjust - geopotential_dse - - calc_te_and_aam_budgets - dtcore_reset - diag_phys_tend_writeout - clybry_fam_set - - - diff --git a/suites/suite_cam7.xml b/suites/suite_cam7.xml index b949654e..49c345cb 100644 --- a/suites/suite_cam7.xml +++ b/suites/suite_cam7.xml @@ -10,15 +10,19 @@ check_energy_gmean_diagnostics - + check_energy_zero_fluxes check_energy_fix apply_heating_rate geopotential_temp - - + check_energy_scaling check_energy_chng + + check_energy_fix_diagnostics dadadj From 252b500a93c89f36ece7d8ba08fd8eb025279eaa Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 14 Feb 2025 11:50:33 -0500 Subject: [PATCH 5/6] Add kessler_test Test SDF for FPHYStest to split off from FKESSLER test (#204) Originator(s): jimmielin Description (include issue title and the keyword ['closes', 'fixes', 'resolves'] and issue number): The FPHYStest test for Kessler in SIMA right now uses the same `suite_kessler.xml` as FKESSLER, which means it includes `physics_after_coupler` and other `check_energy` schemes we do not want to include in the FPHYStest test which should only be testing `kessler_tend`. This PR includes a new `suite_kessler_test.xml` that removes all the other schemes from the SDF and keeps only what we want to test, and consistent with how the snapshots are taken in CAM: ```fortran if (trim(cam_take_snapshot_before) == "kessler_tend") then call cam_snapshot_all_outfld(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf) end if call kessler_tend(state, ptend, ztodt, pbuf) if ( (trim(cam_take_snapshot_after) == "kessler_tend") .and. & (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then call cam_snapshot_ptend_outfld(ptend, lchnk) end if call physics_update(state, ptend, ztodt, tend) if (trim(cam_take_snapshot_after) == "kessler_tend") then call cam_snapshot_all_outfld(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf) end if ``` List all namelist files that were added or changed: N/A List all files eliminated and why: N/A List all files added and what they do: ``` A test/test_suites/suite_kessler_test.xml - new kessler_test suite ``` List all existing files that have been modified, and describe the changes: (Helpful git command: `git diff --name-status development...`) N/A List all automated tests that failed, as well as an explanation for why they weren't fixed: N/A Is this an answer-changing PR? If so, is it a new physics package, algorithm change, tuning change, etc? Yes for FPHYStest kessler test If yes to the above question, describe how this code was validated with the new/modified features: BFB when using this suite --- test/test_suites/suite_kessler_test.xml | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 test/test_suites/suite_kessler_test.xml diff --git a/test/test_suites/suite_kessler_test.xml b/test/test_suites/suite_kessler_test.xml new file mode 100644 index 00000000..9feab1bf --- /dev/null +++ b/test/test_suites/suite_kessler_test.xml @@ -0,0 +1,24 @@ + + + + + calc_exner + temp_to_potential_temp + calc_dry_air_ideal_gas_density + wet_to_dry_water_vapor + wet_to_dry_cloud_liquid_water + wet_to_dry_rain + kessler + potential_temp_to_temp + dry_to_wet_water_vapor + dry_to_wet_cloud_liquid_water + dry_to_wet_rain + kessler_update + qneg + geopotential_temp + + + sima_state_diagnostics + kessler_diagnostics + + From 603de701e21c052e42d31975d09ef140b5bd03cf Mon Sep 17 00:00:00 2001 From: mwaxmonsky <137746677+mwaxmonsky@users.noreply.github.com> Date: Wed, 19 Feb 2025 16:08:50 -0500 Subject: [PATCH 6/6] PBL utils transfer (#193) Originator(s): @mwaxmonsky Summary (include the keyword ['closes', 'fixes', 'resolves'] and issue number): Adds pbl_utils capability from CAM. Describe any changes made to the namelist: N/A List all files eliminated and why: N/A List all files added and what they do: M phys_utils/atmos_phys_pbl_utils.F90 - Transferred all functionality (except for `compute_radf(...)`) from ESCOMP/cam/src/physics/cam/pbl_utils.F90 M .github/workflows/unit-tests.yaml - Added phys_utils to code coverage. M test/unit-test/CMakeLists.txt M test/unit-test/tests/CMakeLists.txt M test/unit-test/tests/phys_utils/CMakeLists.txt M test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf - Creating CMake target library for `phys_utils` and adding sample unit tests. List all existing files that have been modified, and describe the changes: N/A (Helpful git command: `git diff --name-status development...`) List any test failures: N/A Is this a science-changing update? New physics package, algorithm change, tuning changes, etc? No --------- Co-authored-by: Jesse Nusbaumer --- .github/workflows/unit-tests.yaml | 2 +- phys_utils/atmos_phys_pbl_utils.F90 | 216 ++++++++++++++++++ test/unit-test/CMakeLists.txt | 9 + test/unit-test/tests/CMakeLists.txt | 2 +- .../unit-test/tests/phys_utils/CMakeLists.txt | 4 + .../tests/phys_utils/test_atmos_pbl_utils.pf | 25 ++ 6 files changed, 256 insertions(+), 2 deletions(-) create mode 100644 phys_utils/atmos_phys_pbl_utils.F90 create mode 100644 test/unit-test/tests/phys_utils/CMakeLists.txt create mode 100644 test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf diff --git a/.github/workflows/unit-tests.yaml b/.github/workflows/unit-tests.yaml index d9003cee..45ef3a57 100644 --- a/.github/workflows/unit-tests.yaml +++ b/.github/workflows/unit-tests.yaml @@ -62,7 +62,7 @@ jobs: run: | source venv/bin/activate cd build - gcovr -r .. --filter '\.\./schemes' --html atmospheric_physics_code_coverage.html --txt + gcovr -r .. --filter '\.\./schemes' --filter '\.\./phys_utils' --html atmospheric_physics_code_coverage.html --txt - name: Upload code coverage results uses: actions/upload-artifact@v4 diff --git a/phys_utils/atmos_phys_pbl_utils.F90 b/phys_utils/atmos_phys_pbl_utils.F90 new file mode 100644 index 00000000..c4499493 --- /dev/null +++ b/phys_utils/atmos_phys_pbl_utils.F90 @@ -0,0 +1,216 @@ +module atmos_phys_pbl_utils + ! Planetary boundary layer related functions used for vertical diffusion schemes. + + use ccpp_kinds, only: kind_phys + + implicit none + private + + public :: calc_ideal_gas_rrho + public :: calc_friction_velocity + public :: calc_kinematic_heat_flux + public :: calc_kinematic_water_vapor_flux + public :: calc_kinematic_buoyancy_flux + public :: calc_obukhov_length + public :: calc_virtual_temperature + public :: calc_eddy_flux_coefficient + public :: calc_free_atm_eddy_flux_coefficient + + real(kind_phys), parameter :: minimum_friction_velocity = 0.01_kind_phys ! Assuming minimum for coarse grids + real(kind_phys), parameter :: minimum_eddy_flux_coefficient = 0.01_kind_phys ! CCM1 2.f.14 + +contains + + pure elemental function calc_ideal_gas_rrho(rair, surface_temperature, pmid) result(rrho) + ! air density reciprocal + ! Taken from https://glossary.ametsoc.org/wiki/Equation_of_state + ! where \alpha = rrho + real(kind_phys), intent(in) :: rair ! gas constant for dry air [ J kg-1 K-1 ] + real(kind_phys), intent(in) :: surface_temperature ! [ K ] + real(kind_phys), intent(in) :: pmid ! midpoint pressure (bottom level) [ Pa ] + + real(kind_phys) :: rrho ! 1./bottom level density [ m3 kg-1 ] + + rrho = rair * surface_temperature / pmid + end function calc_ideal_gas_rrho + + pure elemental function calc_friction_velocity(taux, tauy, rrho) result(friction_velocity) + ! https://glossary.ametsoc.org/wiki/Friction_velocity + ! NOTE: taux,tauy come form the expansion of the Reynolds stress + ! + ! Also found in: + ! Stull, Roland B. An Introduction to Boundary Layer Meteorology. Springer Kluwer Academic Publishers, 1988. Print. + ! DOI: https://doi.org/10.1007/978-94-009-3027-8 + ! Equation 2.10b, page 67 + + real(kind_phys), intent(in) :: taux ! surface u stress [ N m-2 ] + real(kind_phys), intent(in) :: tauy ! surface v stress [ N m-2 ] + real(kind_phys), intent(in) :: rrho ! 1./bottom level density [ m3 kg-1 ] + + real(kind_phys) :: friction_velocity ! surface friction velocity [ m s-1 ] + + friction_velocity = max( sqrt( sqrt(taux**2 + tauy**2)*rrho ), minimum_friction_velocity ) + end function calc_friction_velocity + + pure elemental function calc_kinematic_heat_flux(shflx, rrho, cpair) result(khfs) + real(kind_phys), intent(in) :: shflx ! surface heat flux [ W m-2 ] + real(kind_phys), intent(in) :: rrho ! 1./bottom level density [ m3 kg-1 ] + real(kind_phys), intent(in) :: cpair ! specific heat of dry air [ J kg-1 K-1 ] + + real(kind_phys) :: khfs ! surface kinematic heat flux [ m K s-1 ] + + khfs = shflx*rrho/cpair + end function calc_kinematic_heat_flux + + pure elemental function calc_kinematic_water_vapor_flux(qflx, rrho) result(kqfs) + real(kind_phys), intent(in) :: qflx ! water vapor flux [ kg m-2 s-1 ] + real(kind_phys), intent(in) :: rrho ! 1./bottom level density [ m3 kg-1 ] + + real(kind_phys) :: kqfs ! surface kinematic water vapor flux [ kg kg-1 m s-1 ] + + kqfs = qflx*rrho + end function calc_kinematic_water_vapor_flux + + pure elemental function calc_kinematic_buoyancy_flux(khfs, zvir, ths, kqfs) result(kbfs) + real(kind_phys), intent(in) :: khfs ! surface kinematic heat flux [ m K s-1 ] + real(kind_phys), intent(in) :: zvir ! rh2o/rair - 1 + real(kind_phys), intent(in) :: ths ! potential temperature at surface [ K ] + real(kind_phys), intent(in) :: kqfs ! surface kinematic water vapor flux [ kg kg-1 m s-1 ] + + ! (`kbfs = \overline{(w' \theta'_v)}_s`) + real(kind_phys) :: kbfs ! surface kinematic buoyancy flux [ m K s-1 ] + + kbfs = khfs + zvir*ths*kqfs + end function calc_kinematic_buoyancy_flux + + pure elemental function calc_obukhov_length(thvs, ustar, g, karman, kbfs) result(obukhov_length) + ! Stull, Roland B. An Introduction to Boundary Layer Meteorology. Springer Kluwer Academic Publishers, 1988. Print. + ! DOI: https://doi.org/10.1007/978-94-009-3027-8 + ! Equation 5.7c, page 181 + ! \frac{-\theta*u_*^3}{g*k*\overline{(w' \theta_v')}_s} = frac{-\theta*u_*^3}{g*k*kbfs} + + real(kind_phys), intent(in) :: thvs ! virtual potential temperature at surface [ K ] + real(kind_phys), intent(in) :: ustar ! Surface friction velocity [ m s-1 ] + real(kind_phys), intent(in) :: g ! acceleration of gravity [ m s-2 ] + real(kind_phys), intent(in) :: karman ! Von Karman's constant (unitless) + real(kind_phys), intent(in) :: kbfs ! surface kinematic buoyancy flux [ m K s-1 ] + + real(kind_phys) :: obukhov_length ! Obukhov length [ m ] + + ! Added sign(...) term to prevent division by 0 and using the fact that + ! `kbfs = \overline{(w' \theta_v')}_s` + obukhov_length = -thvs * ustar**3 / & + (g*karman*(kbfs + sign(1.e-10_kind_phys,kbfs))) + end function calc_obukhov_length + + pure elemental function calc_virtual_temperature(temperature, specific_humidity, zvir) result(virtual_temperature) + ! Williamson, D., Kiehl, J., Ramanathan, V., Dickinson, R., & Hack, J. (1987). + ! Description of the NCAR Community Climate Model (CCM1). + ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987) + ! Equation 2.a.7 + + real(kind_phys), intent(in) :: temperature + real(kind_phys), intent(in) :: specific_humidity + real(kind_phys), intent(in) :: zvir ! rh2o/rair - 1 + + real(kind_phys) :: virtual_temperature + + virtual_temperature = temperature * (1.0_kind_phys + zvir*specific_humidity) + end function calc_virtual_temperature + + pure elemental function calc_eddy_flux_coefficient(mixing_length_squared, & + richardson_number, & + shear_squared) & + result(kvf) + ! Computes exchange coefficients for turbulent flows. + ! + ! The stable case (Richardson number, Ri>0) is taken from Holtslag and + ! Beljaars (1989), ECMWF proceedings. + ! The unstable case (Richardson number, Ri<0) is taken from CCM1. + + real(kind_phys), intent(in) :: mixing_length_squared ! [ m2 ] + real(kind_phys), intent(in) :: richardson_number ! [ 1 ] + real(kind_phys), intent(in) :: shear_squared ! [ s-2 ] + + real(kind_phys) :: kvf ! Eddy diffusivity ! [ m2 s-1 ] + + real(kind_phys) :: fofri ! f(ri) + real(kind_phys) :: kvn ! Neutral Kv + + if( richardson_number < 0.0_kind_phys ) then + fofri = unstable_gradient_richardson_stability_parameter(richardson_number) + else + fofri = stable_gradient_richardson_stability_parameter(richardson_number) + end if + kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared) + kvf = max( minimum_eddy_flux_coefficient, kvn * fofri ) + end function calc_eddy_flux_coefficient + + pure elemental function calc_free_atm_eddy_flux_coefficient(mixing_length_squared, & + richardson_number, & + shear_squared) & + result(kvf) + ! same as austausch_atm but only mixing for Ri<0 + ! i.e. no background mixing and mixing for Ri>0 + + real(kind_phys), intent(in) :: mixing_length_squared ! [ m2 ] + real(kind_phys), intent(in) :: richardson_number ! [ 1 ] + real(kind_phys), intent(in) :: shear_squared ! [ s-2 ] + + real(kind_phys) :: kvf ! [ m2 s-1 ] + + real(kind_phys) :: fofri ! f(ri) + real(kind_phys) :: kvn ! Neutral Kv + + kvf = 0.0_kind_phys + if( richardson_number < 0.0_kind_phys ) then + fofri = unstable_gradient_richardson_stability_parameter(richardson_number) + kvn = neutral_exchange_coefficient(mixing_length_squared, shear_squared) + kvf = kvn * fofri + end if + end function calc_free_atm_eddy_flux_coefficient + + pure elemental function unstable_gradient_richardson_stability_parameter(richardson_number) result(modifier) + ! Williamson, D., Kiehl, J., Ramanathan, V., Dickinson, R., & Hack, J. (1987). + ! Description of the NCAR Community Climate Model (CCM1). + ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987) + ! Equation 2.f.13 + ! \sqrt{ 1-18*Ri } + + real(kind_phys), intent(in) :: richardson_number + + real(kind_phys) :: modifier + + modifier = sqrt( 1._kind_phys - 18._kind_phys * richardson_number ) + end function unstable_gradient_richardson_stability_parameter + + pure elemental function stable_gradient_richardson_stability_parameter(richardson_number) result(modifier) + ! Holtslag, A. A. M., and Beljaars A. C. M. , 1989: Surface flux parameterization schemes: Developments and experiences at KNMI. + ! ECMWF Workshop on Parameterization of Fluxes and Land Surface, Reading, United Kingdom, ECMWF, 121–147. + ! equation 20, page 140 + ! Originally used published equation from CCM1, 2.f.12, page 11 + ! \frac{1}{1+10*Ri*(1+8*Ri)} + + real(kind_phys), intent(in) :: richardson_number + + real(kind_phys) :: modifier + + modifier = 1.0_kind_phys / & + ( 1.0_kind_phys + 10.0_kind_phys * richardson_number * ( 1.0_kind_phys + 8.0_kind_phys * richardson_number ) ) + end function stable_gradient_richardson_stability_parameter + + pure elemental function neutral_exchange_coefficient(mixing_length_squared, shear_squared) result(neutral_k) + ! Williamson, D., Kiehl, J., Ramanathan, V., Dickinson, R., & Hack, J. (1987). + ! Description of the NCAR Community Climate Model (CCM1). + ! University Corporation for Atmospheric Research. https://doi.org/10.5065/D6TB14WH (Original work published 1987) + ! Equation 2.f.15, page 12 + ! NOTE: shear_squared variable currently (01/2025) computed in hb_diff.F90 (s2 in trbintd(...)) matches referenced equation. + + real(kind_phys), intent(in) :: mixing_length_squared ! [ m2 ] + real(kind_phys), intent(in) :: shear_squared ! [ s-2 ] + + real(kind_phys) :: neutral_k + + neutral_k = mixing_length_squared * sqrt(shear_squared) + end function neutral_exchange_coefficient +end module atmos_phys_pbl_utils diff --git a/test/unit-test/CMakeLists.txt b/test/unit-test/CMakeLists.txt index d307934e..49d17506 100644 --- a/test/unit-test/CMakeLists.txt +++ b/test/unit-test/CMakeLists.txt @@ -32,6 +32,15 @@ add_library(utilities ${UTILITIES_SRC}) target_compile_options(utilities PRIVATE -ffree-line-length-none) target_include_directories(utilities PUBLIC ${CMAKE_CURRENT_BINARY_DIR}) +set(PHYS_UTILS_SRC + ../../phys_utils/atmos_phys_pbl_utils.F90 + include/ccpp_kinds.F90 +) + +add_library(phys_utils ${PHYS_UTILS_SRC}) +target_compile_options(phys_utils PRIVATE -ffree-line-length-none) +target_include_directories(phys_utils PUBLIC ${CMAKE_CURRENT_BINARY_DIR}) + if(ATMOSPHERIC_PHYSICS_ENABLE_TESTS OR ATMOSPHERIC_PHYSICS_ENABLE_CODE_COVERAGE) enable_testing() add_subdirectory(tests) diff --git a/test/unit-test/tests/CMakeLists.txt b/test/unit-test/tests/CMakeLists.txt index b4c4714c..705189f1 100644 --- a/test/unit-test/tests/CMakeLists.txt +++ b/test/unit-test/tests/CMakeLists.txt @@ -1,2 +1,2 @@ add_subdirectory(utilities) - +add_subdirectory(phys_utils) diff --git a/test/unit-test/tests/phys_utils/CMakeLists.txt b/test/unit-test/tests/phys_utils/CMakeLists.txt new file mode 100644 index 00000000..5e3020cc --- /dev/null +++ b/test/unit-test/tests/phys_utils/CMakeLists.txt @@ -0,0 +1,4 @@ +add_pfunit_ctest(phys_utils_tests + TEST_SOURCES test_atmos_pbl_utils.pf + LINK_LIBRARIES phys_utils +) diff --git a/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf new file mode 100644 index 00000000..d7bb2683 --- /dev/null +++ b/test/unit-test/tests/phys_utils/test_atmos_pbl_utils.pf @@ -0,0 +1,25 @@ +@test +subroutine test_free_eddy_coef_is_zero_when_ri_equals_zero() + use funit + use atmos_phys_pbl_utils, only : calc_free_atm_eddy_flux_coefficient + use ccpp_kinds, only: kind_phys + + real(kind_phys) :: kvf + + kvf = calc_free_atm_eddy_flux_coefficient(30.0_kind_phys, 0.0_kind_phys, 0.01_kind_phys) + + @assertEqual(0.0_kind_phys, kvf) +end subroutine test_free_eddy_coef_is_zero_when_ri_equals_zero + +@test +subroutine test_free_eddy_atm_coef_is_zero_when_ri_stable_near_zero() + use funit + use atmos_phys_pbl_utils, only : calc_free_atm_eddy_flux_coefficient + use ccpp_kinds, only: kind_phys + + real(kind_phys) :: kvf + + kvf = calc_free_atm_eddy_flux_coefficient(30.0_kind_phys, nearest(0.0_kind_phys, 1.0_kind_phys), 0.01_kind_phys) + + @assertEqual(0.0_kind_phys, kvf) +end subroutine test_free_eddy_atm_coef_is_zero_when_ri_stable_near_zero