From ff30858f55edcc7c70b6aea706ffb1d0bbdee190 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 16 Dec 2024 14:46:55 -0700 Subject: [PATCH 01/24] Tidy up MPAS-related changes introduced in PR #316 * Adjust wording and keep code comments up-to-date * Concentrate all calls to `mark_as_initialized` in one place * Fix up code style inconsistencies --- src/dynamics/mpas/dyn_comp.F90 | 37 +++++++++++++++++------------- src/dynamics/mpas/dyn_coupling.F90 | 17 +++++--------- 2 files changed, 27 insertions(+), 27 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 3911af4f..a911531a 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -187,9 +187,8 @@ end subroutine dyn_readnl ! ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) - use cam_thermo_formula, only: energy_formula_dycore, ENERGY_FORMULA_DYCORE_MPAS - use phys_vars_init_check, only: mark_as_initialized - use physics_types, only: dycore_energy_consistency_adjust + use cam_thermo_formula, only: energy_formula_dycore, energy_formula_dycore_mpas + use physics_types, only: dycore_energy_consistency_adjust type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in @@ -207,14 +206,13 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) nullify(pio_init_file) nullify(pio_topo_file) - ! Set dynamical core energy formula for use in cam_thermo. - energy_formula_dycore = ENERGY_FORMULA_DYCORE_MPAS - call mark_as_initialized('total_energy_formula_for_dycore') + ! Set the energy formula of dynamical core to MPAS for use in `cam_thermo`. + energy_formula_dycore = energy_formula_dycore_mpas - ! Dynamical core energy is not consistent with CAM physics and requires - ! temperature and temperature tendency adjustment at end of physics. + ! The total energy of dynamical core, which uses "MPAS formula" as set above, is not consistent with + ! that of CAM physics, which uses "FV formula". Therefore, temperature and temperature tendency adjustments + ! are needed at the end of each physics time step. dycore_energy_consistency_adjust = .true. - call mark_as_initialized('flag_for_dycore_energy_consistency_adjustment') allocate(constituent_name(num_advected), stat=ierr) call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__) @@ -931,15 +929,21 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) end if end subroutine dyn_exchange_constituent_state - !> Mark everything in the `physics_{state,tend}` derived types along with constituents as initialized - !> to prevent physics from attempting to read them from a file. These variables are to be exchanged later - !> during dynamics-physics coupling. + !> Mark everything in the `physics_types` module along with constituents as initialized + !> to prevent physics from attempting to read them from a file. !> (KCW, 2024-05-23) subroutine mark_variable_as_initialized() character(*), parameter :: subname = 'dyn_comp::mark_variable_as_initialized' integer :: i - ! CCPP standard names of physical quantities in the `physics_{state,tend}` derived types. + ! The variables below are managed by dynamics interface. + ! We are responsible for initializing and updating them. + + ! These variables are to be set during dynamics initialization. + call mark_as_initialized('flag_for_dycore_energy_consistency_adjustment') + call mark_as_initialized('total_energy_formula_for_dycore') + + ! These variables are to be set during dynamics-physics coupling. call mark_as_initialized('air_pressure') call mark_as_initialized('air_pressure_at_interface') call mark_as_initialized('air_pressure_of_dry_air') @@ -960,6 +964,7 @@ subroutine mark_variable_as_initialized() call mark_as_initialized('reciprocal_of_air_pressure_thickness') call mark_as_initialized('reciprocal_of_air_pressure_thickness_of_dry_air') call mark_as_initialized('reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure') + call mark_as_initialized('specific_heat_of_air_used_in_dycore') call mark_as_initialized('surface_air_pressure') call mark_as_initialized('surface_geopotential') call mark_as_initialized('surface_pressure_of_dry_air') @@ -972,10 +977,10 @@ subroutine mark_variable_as_initialized() call mark_as_initialized(trim(adjustl(const_name(i)))) end do - call mark_as_initialized('specific_heat_of_air_used_in_dycore') + ! The variables below are not managed by dynamics interface. They are used by external CCPP physics schemes. + ! While we are not responsible for initializing or updating them, we still need to help mark them as initialized. - ! These energy variables are calculated by check_energy_timestep_init - ! but need to be marked here + ! These variables are to be set externally by the `check_energy_chng` CCPP physics scheme. call mark_as_initialized('vertically_integrated_total_energy_at_end_of_physics_timestep') call mark_as_initialized('vertically_integrated_total_energy_using_dycore_energy_formula') call mark_as_initialized('vertically_integrated_total_energy_using_dycore_energy_formula_at_start_of_physics_timestep') diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index a2066e16..14b211f3 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -3,7 +3,7 @@ module dyn_coupling use cam_abortutils, only: check_allocate, endrun use cam_constituents, only: const_is_water_species, const_qmin, num_advected use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update - use cam_thermo_formula, only: ENERGY_FORMULA_DYCORE_MPAS + use cam_thermo_formula, only: energy_formula_dycore_mpas use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_state, reverse, mpas_dynamical_core, & ncells_solve use dynconst, only: constant_cpd => cpair, constant_g => gravit, constant_p0 => pref, & @@ -333,15 +333,10 @@ subroutine set_physics_state_external() call cam_thermo_dry_air_update( & constituents, phys_state % t, ncells_solve, pver, cam_runtime_opts % update_thermodynamic_variables()) - ! update cp_or_cv_dycore in SIMA state. - ! (note: at this point q is dry) + ! Update `cp_or_cv_dycore` by calling `cam_thermo_water_update`. + ! Note that this subroutine expects constituents to be dry. call cam_thermo_water_update( & - mmr = constituents, & ! dry MMR - ncol = ncells_solve, & - pver = pver, & - energy_formula = ENERGY_FORMULA_DYCORE_MPAS, & - cp_or_cv_dycore = cp_or_cv_dycore & - ) + constituents, ncells_solve, pver, energy_formula_dycore_mpas, cp_or_cv_dycore) ! This variable name is really misleading. It actually represents the reciprocal of Exner function ! with respect to surface pressure. This definition is sometimes used for boundary layer work. See @@ -364,7 +359,7 @@ subroutine set_physics_state_external() end if ! Set `zi` (i.e., geopotential height at layer interfaces) and `zm` (i.e., geopotential height at layer midpoints). - ! Note that `rairv` and `zvirv` are updated externally by `cam_thermo_update`. + ! Note that `rairv` and `zvirv` are updated externally by `cam_thermo_dry_air_update`. call geopotential_temp_run( & pver, lagrangian_vertical, pver, 1, pverp, 1, num_advected, & phys_state % lnpint, phys_state % pint, phys_state % pmid, phys_state % pdel, phys_state % rpdel, phys_state % t, & @@ -378,7 +373,7 @@ subroutine set_physics_state_external() end if ! Set `dse` (i.e., dry static energy). - ! Note that `cpairv` is updated externally by `cam_thermo_update`. + ! Note that `cpairv` is updated externally by `cam_thermo_dry_air_update`. call update_dry_static_energy_run( & pver, constant_g, phys_state % t, phys_state % zm, phys_state % phis, phys_state % dse, cpairv, ierr, cerr) From 14f701984173f23e085a3856bc69948c71522c85 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 21 Oct 2024 14:13:27 -0600 Subject: [PATCH 02/24] Add support for MPAS output stream Also fix variable rank information. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 67 +++++++++++++------ 1 file changed, 46 insertions(+), 21 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index ca8e4159..b9188fb8 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -189,6 +189,8 @@ end subroutine model_error_if !> !> Here, it is described as: !> var_info_type(name="xCell", type="real", rank=1) + !> However, note that MPAS treats the "Time" dimension specially. It is implemented as 1-d pointer arrays of + !> custom derived types. For a variable with the "Time" dimension, its rank needs to be subtracted by one. type :: var_info_type private @@ -277,40 +279,54 @@ end subroutine model_error_if ! the `atm_init_coupled_diagnostics` subroutine in MPAS. ! If a variable first appears on the LHS of an equation, it should be in restart. ! If a variable first appears on the RHS of an equation, it should be in input. + ! The remaining ones of interest should be in output. !> This list corresponds to the "input" stream in MPAS registry. !> It consists of variables that are members of the "diag" and "state" struct. !> Only variables that are specific to the "input" stream are included. type(var_info_type), parameter :: input_var_info_list(*) = [ & - var_info_type('Time' , 'real' , 1), & + var_info_type('Time' , 'real' , 0), & var_info_type('initial_time' , 'character' , 0), & - var_info_type('rho' , 'real' , 3), & - var_info_type('rho_base' , 'real' , 3), & + var_info_type('rho' , 'real' , 2), & + var_info_type('rho_base' , 'real' , 2), & var_info_type('scalars' , 'real' , 3), & - var_info_type('theta' , 'real' , 3), & - var_info_type('theta_base' , 'real' , 3), & - var_info_type('u' , 'real' , 3), & - var_info_type('w' , 'real' , 3), & - var_info_type('xtime' , 'character' , 1) & + var_info_type('theta' , 'real' , 2), & + var_info_type('theta_base' , 'real' , 2), & + var_info_type('u' , 'real' , 2), & + var_info_type('w' , 'real' , 2), & + var_info_type('xtime' , 'character' , 0) & ] !> This list corresponds to the "restart" stream in MPAS registry. !> It consists of variables that are members of the "diag" and "state" struct. !> Only variables that are specific to the "restart" stream are included. type(var_info_type), parameter :: restart_var_info_list(*) = [ & - var_info_type('exner' , 'real' , 1), & - var_info_type('exner_base' , 'real' , 1), & - var_info_type('pressure_base' , 'real' , 1), & - var_info_type('pressure_p' , 'real' , 1), & - var_info_type('rho_p' , 'real' , 1), & - var_info_type('rho_zz' , 'real' , 1), & - var_info_type('rtheta_base' , 'real' , 1), & - var_info_type('rtheta_p' , 'real' , 1), & - var_info_type('ru' , 'real' , 1), & - var_info_type('ru_p' , 'real' , 1), & - var_info_type('rw' , 'real' , 1), & - var_info_type('rw_p' , 'real' , 1), & - var_info_type('theta_m' , 'real' , 1) & + var_info_type('exner' , 'real' , 2), & + var_info_type('exner_base' , 'real' , 2), & + var_info_type('pressure_base' , 'real' , 2), & + var_info_type('pressure_p' , 'real' , 2), & + var_info_type('rho_p' , 'real' , 2), & + var_info_type('rho_zz' , 'real' , 2), & + var_info_type('rtheta_base' , 'real' , 2), & + var_info_type('rtheta_p' , 'real' , 2), & + var_info_type('ru' , 'real' , 2), & + var_info_type('ru_p' , 'real' , 2), & + var_info_type('rw' , 'real' , 2), & + var_info_type('rw_p' , 'real' , 2), & + var_info_type('theta_m' , 'real' , 2) & + ] + + !> This list corresponds to the "output" stream in MPAS registry. + !> It consists of variables that are members of the "diag" struct. + !> Only variables that are specific to the "output" stream are included. + type(var_info_type), parameter :: output_var_info_list(*) = [ & + var_info_type('divergence' , 'real' , 2), & + var_info_type('pressure' , 'real' , 2), & + var_info_type('relhum' , 'real' , 2), & + var_info_type('surface_pressure' , 'real' , 1), & + var_info_type('uReconstructMeridional' , 'real' , 2), & + var_info_type('uReconstructZonal' , 'real' , 2), & + var_info_type('vorticity' , 'real' , 2) & ] contains !> Print a debug message with optionally the value(s) of a variable. @@ -1752,6 +1768,8 @@ pure function parse_stream_name_fragment(stream_name_fragment) result(var_info_l allocate(var_info_list, source=input_var_info_list) case ('restart') allocate(var_info_list, source=restart_var_info_list) + case ('output') + allocate(var_info_list, source=output_var_info_list) case default allocate(var_info_list(0)) @@ -1775,6 +1793,13 @@ pure function parse_stream_name_fragment(stream_name_fragment) result(var_info_l var_info_list_buffer = pack(restart_var_info_list, var_name_list == trim(adjustl(stream_name_fragment))) var_info_list = [var_info_list, var_info_list_buffer] end if + + var_name_list = output_var_info_list % name + + if (any(var_name_list == trim(adjustl(stream_name_fragment)))) then + var_info_list_buffer = pack(output_var_info_list, var_name_list == trim(adjustl(stream_name_fragment))) + var_info_list = [var_info_list, var_info_list_buffer] + end if end select end function parse_stream_name_fragment From b7410ce36a17381f4cf71bb923ae5907d4afe9c5 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 23 Oct 2024 14:41:27 -0600 Subject: [PATCH 03/24] Relax restrictions on mixed precision Relax restrictions on mixed precision between MPAS dynamical core and its input data. --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index b9188fb8..d57e7e80 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -2264,11 +2264,15 @@ subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compa cycle end if case ('real') + ! When MPAS dynamical core is compiled at single precision, pairing it with double precision input data + ! is not allowed to prevent loss of precision. if (rkind == r4kind .and. vartype /= pio_real) then cycle end if - if (rkind == r8kind .and. vartype /= pio_double) then + ! When MPAS dynamical core is compiled at double precision, pairing it with single and double precision + ! input data is allowed. + if (rkind == r8kind .and. vartype /= pio_real .and. vartype /= pio_double) then cycle end if case default From 0c0d94d8b3f65ef3e228985e1a79df99176b7327 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Tue, 29 Oct 2024 17:18:07 -0600 Subject: [PATCH 04/24] Use existing functionality to access MPAS configuration `get_variable_pointer` also works for accessing MPAS configuration. It has built-in error checking so no need to do it manually. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 30 ++++--------------- 1 file changed, 5 insertions(+), 25 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index d57e7e80..70ba65d4 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -635,11 +635,7 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, & subname, __LINE__) end select - call mpas_pool_get_config(self % domain_ptr % configs, 'config_calendar_type', config_pointer_c) - - if (.not. associated(config_pointer_c)) then - call self % model_error('Failed to find config "config_calendar_type"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_calendar_type') config_pointer_c = trim(adjustl(mpas_calendar)) call self % debug_print('config_calendar_type = ', [config_pointer_c]) @@ -648,43 +644,27 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, & ! MPAS represents date and time in ISO 8601 format. However, the separator between date and time is `_` ! instead of standard `T`. ! Format in `YYYY-MM-DD_hh:mm:ss` is acceptable. - call mpas_pool_get_config(self % domain_ptr % configs, 'config_start_time', config_pointer_c) - - if (.not. associated(config_pointer_c)) then - call self % model_error('Failed to find config "config_start_time"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_start_time') config_pointer_c = stringify(start_date_time(1:3), '-') // '_' // stringify(start_date_time(4:6), ':') call self % debug_print('config_start_time = ', [config_pointer_c]) nullify(config_pointer_c) - call mpas_pool_get_config(self % domain_ptr % configs, 'config_stop_time', config_pointer_c) - - if (.not. associated(config_pointer_c)) then - call self % model_error('Failed to find config "config_stop_time"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_stop_time') config_pointer_c = stringify(stop_date_time(1:3), '-') // '_' // stringify(stop_date_time(4:6), ':') call self % debug_print('config_stop_time = ', [config_pointer_c]) nullify(config_pointer_c) ! Format in `DD_hh:mm:ss` is acceptable. - call mpas_pool_get_config(self % domain_ptr % configs, 'config_run_duration', config_pointer_c) - - if (.not. associated(config_pointer_c)) then - call self % model_error('Failed to find config "config_run_duration"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_c, 'cfg', 'config_run_duration') config_pointer_c = stringify([run_duration(1)]) // '_' // stringify(run_duration(2:4), ':') call self % debug_print('config_run_duration = ', [config_pointer_c]) nullify(config_pointer_c) ! Reflect current run type to MPAS. - call mpas_pool_get_config(self % domain_ptr % configs, 'config_do_restart', config_pointer_l) - - if (.not. associated(config_pointer_l)) then - call self % model_error('Failed to find config "config_do_restart"', subname, __LINE__) - end if + call self % get_variable_pointer(config_pointer_l, 'cfg', 'config_do_restart') if (initial_run) then ! Run type is initial run. From 2dd1e95f9c1681f160b9568f77a693a9c8dd7caf Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 21 Oct 2024 14:17:42 -0600 Subject: [PATCH 05/24] Implement finalization for MPAS dynamical core --- .../mpas/driver/dyn_mpas_subdriver.F90 | 141 +++++++++++++++++- 1 file changed, 135 insertions(+), 6 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 70ba65d4..788828ea 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -23,13 +23,14 @@ module dyn_mpas_subdriver ! Modules from MPAS. use atm_core, only: atm_compute_output_diagnostics, atm_do_timestep, atm_mpas_init_block use atm_core_interface, only: atm_setup_core, atm_setup_domain - use atm_time_integration, only: mpas_atm_dynamics_init + use atm_time_integration, only: mpas_atm_dynamics_init, mpas_atm_dynamics_finalize use mpas_atm_dimensions, only: mpas_atm_set_dims - use mpas_atm_halos, only: atm_build_halo_groups, exchange_halo_group - use mpas_atm_threading, only: mpas_atm_threading_init + use mpas_atm_halos, only: atm_build_halo_groups, atm_destroy_halo_groups, exchange_halo_group + use mpas_atm_threading, only: mpas_atm_threading_init, mpas_atm_threading_finalize use mpas_attlist, only: mpas_modify_att use mpas_bootstrapping, only: mpas_bootstrap_framework_phase1, mpas_bootstrap_framework_phase2 use mpas_constants, only: mpas_constants_compute_derived + use mpas_decomp, only: mpas_decomp_destroy_decomp_list use mpas_derived_types, only: core_type, domain_type, & mpas_pool_type, mpas_pool_field_info_type, & mpas_pool_character, mpas_pool_real, mpas_pool_integer, & @@ -42,11 +43,12 @@ module dyn_mpas_subdriver use mpas_dmpar, only: mpas_dmpar_exch_halo_field, & mpas_dmpar_max_int, mpas_dmpar_sum_int use mpas_domain_routines, only: mpas_allocate_domain - use mpas_field_routines, only: mpas_allocate_scratch_field - use mpas_framework, only: mpas_framework_init_phase1, mpas_framework_init_phase2 + use mpas_field_routines, only: mpas_allocate_scratch_field, mpas_deallocate_scratch_field + use mpas_framework, only: mpas_framework_init_phase1, mpas_framework_init_phase2, mpas_framework_finalize use mpas_io_streams, only: mpas_createstream, mpas_closestream, mpas_streamaddfield, & mpas_readstream, mpas_writestream, mpas_writestreamatt use mpas_kind_types, only: rkind, r4kind, r8kind, strkind + use mpas_log, only: mpas_log_finalize use mpas_pool_routines, only: mpas_pool_create_pool, mpas_pool_destroy_pool, mpas_pool_get_subpool, & mpas_pool_add_config, mpas_pool_get_config, & mpas_pool_get_array, & @@ -56,9 +58,11 @@ module dyn_mpas_subdriver use mpas_stream_inquiry, only: mpas_stream_inquiry_new_streaminfo use mpas_stream_manager, only: postread_reindex, prewrite_reindex, postwrite_reindex use mpas_string_utils, only: mpas_string_replace - use mpas_timekeeping, only: mpas_advance_clock, mpas_get_clock_time, mpas_get_time, & + use mpas_timekeeping, only: mpas_advance_clock, mpas_destroy_clock, & + mpas_get_clock_time, mpas_get_time, & mpas_set_timeinterval, & operator(+), operator(<) + use mpas_timer, only: mpas_timer_write_header, mpas_timer_write, mpas_timer_finalize use mpas_vector_operations, only: mpas_initialize_vectors implicit none @@ -128,6 +132,7 @@ end subroutine model_error_if procedure, pass, public :: compute_edge_wind => dyn_mpas_compute_edge_wind procedure, pass, public :: init_phase4 => dyn_mpas_init_phase4 procedure, pass, public :: run => dyn_mpas_run + procedure, pass, public :: final => dyn_mpas_final ! Accessor subroutines for users to access internal states of MPAS dynamical core. @@ -3021,6 +3026,130 @@ subroutine dyn_mpas_run(self) call self % debug_print(subname // ' completed') end subroutine dyn_mpas_run + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_final + ! + !> \brief Finalizes MPAS dynamical core as well as its framework + !> \author Michael Duda + !> \date 29 February 2020 + !> \details + !> This subroutine finalizes and cleans up MPAS dynamical core as well as its + !> framework that was set up during initialization. Finalization happens in + !> reverse chronological order. + !> Essentially, it closely follows what is done in `atm_core_finalize` and + !> `mpas_finalize`, except that here, there is no need to call MPAS diagnostics + !> manager or MPAS stream manager. + !> \addenda + !> Ported and refactored for CAM-SIMA. (KCW, 2024-10-10) + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_final(self) + class(mpas_dynamical_core_type), intent(inout) :: self + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_final' + integer :: ierr + type(field2dreal), pointer :: field_2d_real + + call self % debug_print(subname // ' entered') + + nullify(field_2d_real) + + ! First, wind down MPAS dynamical core by calling its own finalization procedures. + + ! Some additional "scratch" fields are needed for interoperability with CAM-SIMA, but they are not finalized by + ! `mpas_atm_dynamics_finalize`. Finalize them below. + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_uzonal', field_2d_real, timelevel=1) + call mpas_deallocate_scratch_field(field_2d_real) + nullify(field_2d_real) + + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, 'tend_umerid', field_2d_real, timelevel=1) + call mpas_deallocate_scratch_field(field_2d_real) + nullify(field_2d_real) + + ! Opposite to `mpas_atm_dynamics_init`. + call mpas_atm_dynamics_finalize(self % domain_ptr) + + ! Opposite to `atm_build_halo_groups`. + call atm_destroy_halo_groups(self % domain_ptr, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to destroy halo exchange groups', subname, __LINE__) + end if + + nullify(exchange_halo_group) + + ! Opposite to `mpas_atm_threading_init`. + call mpas_atm_threading_finalize(self % domain_ptr % blocklist, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to clean up OpenMP threading', subname, __LINE__) + end if + + ! Opposite to `mpas_create_clock`, which was called by `atm_simulation_clock_init`, then `atm_setup_clock`. + call mpas_destroy_clock(self % domain_ptr % clock, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to clean up clock', subname, __LINE__) + end if + + ! Opposite to `mpas_decomp_create_decomp_list`, which was called by `atm_setup_decompositions`. + call mpas_decomp_destroy_decomp_list(self % domain_ptr % decompositions) + + deallocate(self % domain_ptr % streaminfo) + + ! Write timing information to log. + call mpas_timer_write_header() + call mpas_timer_write() + + ! Opposite to `mpas_timer_init`, which was called by `mpas_framework_init_phase2`. + call mpas_timer_finalize(self % domain_ptr) + + ! Opposite to `mpas_log_init`, which was called by `atm_setup_log`. + call mpas_log_finalize(ierr) + + if (ierr /= 0) then + call self % model_error('Failed to clean up log', subname, __LINE__) + end if + + ! Opposite to `mpas_framework_init_phase1` and `mpas_framework_init_phase2`. + call mpas_framework_finalize(self % domain_ptr % dminfo, self % domain_ptr) + + call self % debug_print(subname // ' completed') + + call self % debug_print('Successful finalization of MPAS dynamical core') + + ! Second, clean up this MPAS dynamical core instance. + + ! Initialized by `dyn_mpas_init_phase4`. + self % coupling_time_interval = 0 + self % number_of_time_steps = 0 + + ! Initialized by `dyn_mpas_define_scalar`. + deallocate(self % constituent_name) + deallocate(self % index_constituent_to_mpas_scalar) + deallocate(self % index_mpas_scalar_to_constituent) + deallocate(self % is_water_species) + + ! Initialized by `dyn_mpas_init_phase3`. + self % number_of_constituents = 0 + + ! Initialized by `dyn_mpas_init_phase1`. + self % log_unit = output_unit + self % mpi_comm = mpi_comm_null + self % mpi_rank = 0 + self % mpi_rank_root = .false. + + nullify(self % model_error) + + deallocate(self % corelist % domainlist) + deallocate(self % corelist) + + nullify(self % corelist) + nullify(self % domain_ptr) + + self % debug_output = .false. + end subroutine dyn_mpas_final + !------------------------------------------------------------------------------- ! function dyn_mpas_get_constituent_name ! From f7b459ba1af074b1496fe6a6339d60bda11d2b86 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 21 Oct 2024 14:20:06 -0600 Subject: [PATCH 06/24] Implement `dyn_final` --- src/dynamics/mpas/dyn_comp.F90 | 72 ++++++++++++++++++++++++++++++++-- 1 file changed, 68 insertions(+), 4 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index a911531a..c6eff173 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -51,7 +51,7 @@ module dyn_comp public :: dyn_readnl public :: dyn_init public :: dyn_run - ! public :: dyn_final + public :: dyn_final public :: dyn_debug_print public :: dyn_exchange_constituent_state @@ -999,9 +999,73 @@ subroutine dyn_run() call mpas_dynamical_core % run() end subroutine dyn_run - ! Not used for now. Intended to be called by `stepon_final` in `src/dynamics/mpas/stepon.F90`. - ! subroutine dyn_final() - ! end subroutine dyn_final + !> Finalize MPAS dynamical core as well as its framework. + !> (KCW, 2024-10-04) + subroutine dyn_final() + character(*), parameter :: subname = 'dyn_comp::dyn_final' + + ! Quick hack for dumping variables from MPAS dynamical core. + ! Remove it once history and restart are wired up in CAM-SIMA. + call dyn_variable_dump() + + ! After this point, do not access anything under MPAS dynamical core or runtime errors will ensue. + call mpas_dynamical_core % final() + end subroutine dyn_final + + subroutine dyn_variable_dump() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_instance, only: atm_id + use physics_types, only: phys_state + ! Module(s) from CESM Share. + use shr_pio_mod, only: shr_pio_getioformat, shr_pio_getiosys, shr_pio_getiotype + ! Module(s) from external libraries. + use pio, only: file_desc_t, iosystem_desc_t, pio_createfile, pio_closefile, pio_clobber, pio_noerr + + character(*), parameter :: subname = 'dyn_comp::dyn_variable_dump' + integer :: ierr + integer :: pio_ioformat, pio_iotype + real(kind_r8), pointer :: surface_pressure(:) + type(file_desc_t), pointer :: pio_file + type(iosystem_desc_t), pointer :: pio_iosystem + + nullify(pio_file) + nullify(pio_iosystem) + nullify(surface_pressure) + + call mpas_dynamical_core % get_variable_pointer(surface_pressure, 'diag', 'surface_pressure') + + surface_pressure(1:ncells_solve) = phys_state % ps(:) + + nullify(surface_pressure) + + call mpas_dynamical_core % exchange_halo('surface_pressure') + + allocate(pio_file, stat=ierr) + call check_allocate(ierr, subname, 'pio_file', 'dyn_comp', __LINE__) + + pio_iosystem => shr_pio_getiosys(atm_id) + + pio_ioformat = shr_pio_getioformat(atm_id) + pio_ioformat = ior(pio_ioformat, pio_clobber) + + pio_iotype = shr_pio_getiotype(atm_id) + + ierr = pio_createfile(pio_iosystem, pio_file, pio_iotype, 'dyn_variable_dump.nc', pio_ioformat) + + if (ierr /= pio_noerr) then + call endrun('Failed to create file for variable dumping', subname, __LINE__) + end if + + call mpas_dynamical_core % read_write_stream(pio_file, 'w', 'invariant+input+restart+output') + + call pio_closefile(pio_file) + + deallocate(pio_file) + + nullify(pio_file) + nullify(pio_iosystem) + end subroutine dyn_variable_dump !> Helper function for reversing the order of elements in `array`. !> (KCW, 2024-07-17) From 2e7ac8f34c561f749be6924bccd6f9f7bd3cc70e Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 21 Oct 2024 14:20:32 -0600 Subject: [PATCH 07/24] Wire up `dyn_final` --- src/dynamics/mpas/stepon.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90 index d089636e..d6d86503 100644 --- a/src/dynamics/mpas/stepon.F90 +++ b/src/dynamics/mpas/stepon.F90 @@ -1,7 +1,7 @@ module stepon ! Modules from CAM-SIMA. use camsrfexch, only: cam_out_t - use dyn_comp, only: dyn_import_t, dyn_export_t, dyn_run + use dyn_comp, only: dyn_import_t, dyn_export_t, dyn_run, dyn_final use dyn_coupling, only: dynamics_to_physics_coupling, physics_to_dynamics_coupling use physics_types, only: physics_state, physics_tend use runtime_obj, only: runtime_options @@ -70,5 +70,7 @@ subroutine stepon_final(cam_runtime_opts, dyn_in, dyn_out) type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in type(dyn_export_t), intent(in) :: dyn_out + + call dyn_final() end subroutine stepon_final end module stepon From cbb43cfcaca9deec19d1a883257fbf6a00b228d0 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 11 Nov 2024 17:57:15 -0700 Subject: [PATCH 08/24] Move `use` statements to the smallest scope possible Comply with CAM(-SIMA) coding standards. However, kind and length parameters are exempt from this rule. They are better at module level. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 216 ++++++++++++++---- src/dynamics/mpas/dyn_comp.F90 | 147 ++++++++---- src/dynamics/mpas/dyn_coupling.F90 | 99 +++++--- src/dynamics/mpas/dyn_grid.F90 | 76 ++++-- src/dynamics/mpas/stepon.F90 | 37 ++- 5 files changed, 420 insertions(+), 155 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 788828ea..e6677492 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -8,62 +8,16 @@ module dyn_mpas_subdriver !------------------------------------------------------------------------------- use, intrinsic :: iso_fortran_env, only: output_unit - - ! Modules from external libraries. + ! Module(s) from external libraries. #ifdef MPAS_USE_MPI_F08 use mpi_f08, only: mpi_comm_null, mpi_comm_rank, mpi_success, & mpi_comm_type => mpi_comm, operator(==) #else use mpi, only: mpi_comm_null, mpi_comm_rank, mpi_success #endif - use pio, only: pio_char, pio_int, pio_real, pio_double, & - file_desc_t, iosystem_desc_t, pio_file_is_open, pio_iosystem_is_active, & - pio_inq_varid, pio_inq_varndims, pio_inq_vartype, pio_noerr - - ! Modules from MPAS. - use atm_core, only: atm_compute_output_diagnostics, atm_do_timestep, atm_mpas_init_block - use atm_core_interface, only: atm_setup_core, atm_setup_domain - use atm_time_integration, only: mpas_atm_dynamics_init, mpas_atm_dynamics_finalize - use mpas_atm_dimensions, only: mpas_atm_set_dims - use mpas_atm_halos, only: atm_build_halo_groups, atm_destroy_halo_groups, exchange_halo_group - use mpas_atm_threading, only: mpas_atm_threading_init, mpas_atm_threading_finalize - use mpas_attlist, only: mpas_modify_att - use mpas_bootstrapping, only: mpas_bootstrap_framework_phase1, mpas_bootstrap_framework_phase2 - use mpas_constants, only: mpas_constants_compute_derived - use mpas_decomp, only: mpas_decomp_destroy_decomp_list - use mpas_derived_types, only: core_type, domain_type, & - mpas_pool_type, mpas_pool_field_info_type, & - mpas_pool_character, mpas_pool_real, mpas_pool_integer, & - mpas_stream_type, mpas_stream_noerr, & - mpas_time_type, mpas_timeinterval_type, mpas_now, mpas_start_time, & - mpas_io_native_precision, mpas_io_pnetcdf, mpas_io_read, mpas_io_write, & - field0dchar, field1dchar, & - field0dinteger, field1dinteger, field2dinteger, field3dinteger, & - field0dreal, field1dreal, field2dreal, field3dreal, field4dreal, field5dreal - use mpas_dmpar, only: mpas_dmpar_exch_halo_field, & - mpas_dmpar_max_int, mpas_dmpar_sum_int - use mpas_domain_routines, only: mpas_allocate_domain - use mpas_field_routines, only: mpas_allocate_scratch_field, mpas_deallocate_scratch_field - use mpas_framework, only: mpas_framework_init_phase1, mpas_framework_init_phase2, mpas_framework_finalize - use mpas_io_streams, only: mpas_createstream, mpas_closestream, mpas_streamaddfield, & - mpas_readstream, mpas_writestream, mpas_writestreamatt + ! Module(s) from MPAS. + use mpas_derived_types, only: core_type, domain_type use mpas_kind_types, only: rkind, r4kind, r8kind, strkind - use mpas_log, only: mpas_log_finalize - use mpas_pool_routines, only: mpas_pool_create_pool, mpas_pool_destroy_pool, mpas_pool_get_subpool, & - mpas_pool_add_config, mpas_pool_get_config, & - mpas_pool_get_array, & - mpas_pool_add_dimension, mpas_pool_get_dimension, & - mpas_pool_get_field, mpas_pool_get_field_info, & - mpas_pool_initialize_time_levels, mpas_pool_shift_time_levels - use mpas_stream_inquiry, only: mpas_stream_inquiry_new_streaminfo - use mpas_stream_manager, only: postread_reindex, prewrite_reindex, postwrite_reindex - use mpas_string_utils, only: mpas_string_replace - use mpas_timekeeping, only: mpas_advance_clock, mpas_destroy_clock, & - mpas_get_clock_time, mpas_get_time, & - mpas_set_timeinterval, & - operator(+), operator(<) - use mpas_timer, only: mpas_timer_write_header, mpas_timer_write, mpas_timer_finalize - use mpas_vector_operations, only: mpas_initialize_vectors implicit none @@ -482,6 +436,11 @@ end function stringify ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_phase1(self, mpi_comm, model_error_impl, log_unit, mpas_log_unit) + ! Module(s) from MPAS. + use atm_core_interface, only: atm_setup_core, atm_setup_domain + use mpas_domain_routines, only: mpas_allocate_domain + use mpas_framework, only: mpas_framework_init_phase1 + class(mpas_dynamical_core_type), intent(inout) :: self #ifdef MPAS_USE_MPI_F08 type(mpi_comm_type), intent(in) :: mpi_comm @@ -700,6 +659,12 @@ end subroutine dyn_mpas_read_namelist ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_phase2(self, pio_iosystem) + ! Module(s) from external libraries. + use pio, only: iosystem_desc_t, pio_iosystem_is_active + ! Module(s) from MPAS. + use mpas_framework, only: mpas_framework_init_phase2 + use mpas_stream_inquiry, only: mpas_stream_inquiry_new_streaminfo + class(mpas_dynamical_core_type), intent(in) :: self type(iosystem_desc_t), pointer, intent(in) :: pio_iosystem @@ -789,6 +754,13 @@ end subroutine dyn_mpas_init_phase2 ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file) + ! Module(s) from external libraries. + use pio, only: file_desc_t, pio_file_is_open + ! Module(s) from MPAS. + use mpas_bootstrapping, only: mpas_bootstrap_framework_phase1, mpas_bootstrap_framework_phase2 + use mpas_derived_types, only: mpas_io_pnetcdf, mpas_pool_type + use mpas_pool_routines, only: mpas_pool_add_config, mpas_pool_add_dimension, mpas_pool_get_dimension + class(mpas_dynamical_core_type), intent(inout) :: self integer, intent(in) :: number_of_constituents type(file_desc_t), pointer, intent(in) :: pio_file @@ -885,6 +857,10 @@ end subroutine dyn_mpas_init_phase3 ! !------------------------------------------------------------------------------- subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species) + ! Module(s) from MPAS. + use mpas_derived_types, only: field3dreal, mpas_pool_type + use mpas_pool_routines, only: mpas_pool_add_dimension, mpas_pool_get_field + class(mpas_dynamical_core_type), intent(inout) :: self character(*), intent(in) :: constituent_name(:) logical, intent(in) :: is_water_species(:) @@ -1143,6 +1119,14 @@ end subroutine dyn_mpas_define_scalar ! !------------------------------------------------------------------------------- subroutine dyn_mpas_read_write_stream(self, pio_file, stream_mode, stream_name) + ! Module(s) from external libraries. + use pio, only: file_desc_t + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type, mpas_stream_noerr, mpas_stream_type + use mpas_io_streams, only: mpas_closestream, mpas_readstream, mpas_writestream + use mpas_pool_routines, only: mpas_pool_destroy_pool + use mpas_stream_manager, only: postread_reindex, prewrite_reindex, postwrite_reindex + class(mpas_dynamical_core_type), intent(in) :: self type(file_desc_t), pointer, intent(in) :: pio_file character(*), intent(in) :: stream_mode @@ -1249,6 +1233,17 @@ end subroutine dyn_mpas_read_write_stream ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file, stream_mode, stream_name) + ! Module(s) from external libraries. + use pio, only: file_desc_t, pio_file_is_open + ! Module(s) from MPAS. + use mpas_derived_types, only: field0dchar, field1dchar, & + field0dinteger, field1dinteger, field2dinteger, field3dinteger, & + field0dreal, field1dreal, field2dreal, field3dreal, field4dreal, field5dreal, & + mpas_io_native_precision, mpas_io_pnetcdf, mpas_io_read, mpas_io_write, & + mpas_pool_type, mpas_stream_noerr, mpas_stream_type + use mpas_io_streams, only: mpas_createstream, mpas_streamaddfield + use mpas_pool_routines, only: mpas_pool_add_config, mpas_pool_create_pool, mpas_pool_get_field + class(mpas_dynamical_core_type), intent(in) :: self type(mpas_pool_type), pointer, intent(out) :: mpas_pool type(mpas_stream_type), pointer, intent(out) :: mpas_stream @@ -1584,6 +1579,9 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file !> Helper subroutine for adding a 0-d stream attribute by calling `mpas_writestreamatt` with error checking. !> (KCW, 2024-03-14) subroutine add_stream_attribute_0d(attribute_name, attribute_value) + ! Module(s) from MPAS. + use mpas_io_streams, only: mpas_writestreamatt + character(*), intent(in) :: attribute_name class(*), intent(in) :: attribute_value @@ -1621,6 +1619,9 @@ end subroutine add_stream_attribute_0d !> Helper subroutine for adding a 1-d stream attribute by calling `mpas_writestreamatt` with error checking. !> (KCW, 2024-03-14) subroutine add_stream_attribute_1d(attribute_name, attribute_value) + ! Module(s) from MPAS. + use mpas_io_streams, only: mpas_writestreamatt + character(*), intent(in) :: attribute_name class(*), intent(in) :: attribute_value(:) @@ -1884,6 +1885,16 @@ end function index_unique ! !------------------------------------------------------------------------------- subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compatible, pio_file, var_info) + ! Module(s) from external libraries. + use pio, only: file_desc_t, pio_file_is_open, & + pio_char, pio_int, pio_real, pio_double, & + pio_inq_varid, pio_inq_varndims, pio_inq_vartype, pio_noerr + ! Module(s) from MPAS. + use mpas_derived_types, only: field0dchar, field1dchar, & + field0dinteger, field1dinteger, field2dinteger, field3dinteger, & + field0dreal, field1dreal, field2dreal, field3dreal, field4dreal, field5dreal + use mpas_pool_routines, only: mpas_pool_get_field + class(mpas_dynamical_core_type), intent(in) :: self logical, allocatable, intent(out) :: var_is_present(:) logical, allocatable, intent(out) :: var_is_tkr_compatible(:) @@ -2299,6 +2310,13 @@ end subroutine dyn_mpas_check_variable_status ! !------------------------------------------------------------------------------- subroutine dyn_mpas_exchange_halo(self, field_name) + ! Module(s) from MPAS. + use mpas_derived_types, only: field1dinteger, field2dinteger, field3dinteger, & + field1dreal, field2dreal, field3dreal, field4dreal, field5dreal, & + mpas_pool_field_info_type, mpas_pool_integer, mpas_pool_real + use mpas_dmpar, only: mpas_dmpar_exch_halo_field + use mpas_pool_routines, only: mpas_pool_get_field, mpas_pool_get_field_info + class(mpas_dynamical_core_type), intent(in) :: self character(*), intent(in) :: field_name @@ -2482,6 +2500,10 @@ end subroutine dyn_mpas_exchange_halo ! !------------------------------------------------------------------------------- subroutine dyn_mpas_compute_unit_vector(self) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_vector_operations, only: mpas_initialize_vectors + class(mpas_dynamical_core_type), intent(in) :: self character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_unit_vector' @@ -2666,6 +2688,21 @@ end subroutine dyn_mpas_compute_edge_wind ! !------------------------------------------------------------------------------- subroutine dyn_mpas_init_phase4(self, coupling_time_interval) + ! Module(s) from MPAS. + use atm_core, only: atm_mpas_init_block + use atm_time_integration, only: mpas_atm_dynamics_init + use mpas_atm_dimensions, only: mpas_atm_set_dims + use mpas_atm_halos, only: atm_build_halo_groups, exchange_halo_group + use mpas_atm_threading, only: mpas_atm_threading_init + use mpas_attlist, only: mpas_modify_att + use mpas_constants, only: mpas_constants_compute_derived + use mpas_derived_types, only: field0dreal, field2dreal, mpas_pool_type, mpas_time_type, & + mpas_start_time + use mpas_field_routines, only: mpas_allocate_scratch_field + use mpas_pool_routines, only: mpas_pool_get_field, mpas_pool_initialize_time_levels + use mpas_string_utils, only: mpas_string_replace + use mpas_timekeeping, only: mpas_get_clock_time, mpas_get_time + class(mpas_dynamical_core_type), intent(inout) :: self integer, intent(in) :: coupling_time_interval ! Set the time interval, in seconds, over which MPAS dynamical core ! should integrate each time it is called to run. @@ -2936,6 +2973,15 @@ end subroutine dyn_mpas_init_phase4 ! !------------------------------------------------------------------------------- subroutine dyn_mpas_run(self) + ! Module(s) from MPAS. + use atm_core, only: atm_compute_output_diagnostics, atm_do_timestep + use mpas_derived_types, only: mpas_pool_type, mpas_time_type, mpas_timeinterval_type, & + mpas_now + use mpas_pool_routines, only: mpas_pool_shift_time_levels + use mpas_timekeeping, only: mpas_advance_clock, mpas_get_clock_time, mpas_get_time, & + mpas_set_timeinterval, & + operator(+), operator(<) + class(mpas_dynamical_core_type), intent(inout) :: self character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_run' @@ -3044,6 +3090,19 @@ end subroutine dyn_mpas_run ! !------------------------------------------------------------------------------- subroutine dyn_mpas_final(self) + ! Module(s) from MPAS. + use atm_time_integration, only: mpas_atm_dynamics_finalize + use mpas_atm_halos, only: atm_destroy_halo_groups, exchange_halo_group + use mpas_atm_threading, only: mpas_atm_threading_finalize + use mpas_decomp, only: mpas_decomp_destroy_decomp_list + use mpas_derived_types, only: field2dreal + use mpas_field_routines, only: mpas_deallocate_scratch_field + use mpas_framework, only: mpas_framework_finalize + use mpas_log, only: mpas_log_finalize + use mpas_pool_routines, only: mpas_pool_get_field + use mpas_timekeeping, only: mpas_destroy_clock + use mpas_timer, only: mpas_timer_write_header, mpas_timer_write, mpas_timer_finalize + class(mpas_dynamical_core_type), intent(inout) :: self character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_final' @@ -3369,6 +3428,9 @@ end subroutine dyn_mpas_get_local_mesh_dimension subroutine dyn_mpas_get_global_mesh_dimension(self, & ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max, & sphere_radius) + ! Module(s) from MPAS. + use mpas_dmpar, only: mpas_dmpar_max_int, mpas_dmpar_sum_int + class(mpas_dynamical_core_type), intent(in) :: self integer, intent(out) :: ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max real(rkind), intent(out) :: sphere_radius @@ -3416,6 +3478,10 @@ end subroutine dyn_mpas_get_global_mesh_dimension !> It is used by the `dyn_mpas_get_variable_{pointer,value}_*` subroutines to draw a variable from a pool. !> (KCW, 2024-03-21) subroutine dyn_mpas_get_pool_pointer(self, pool_pointer, pool_name) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_subpool + class(mpas_dynamical_core_type), intent(in) :: self type(mpas_pool_type), pointer, intent(out) :: pool_pointer character(*), intent(in) :: pool_name @@ -3461,6 +3527,10 @@ end subroutine dyn_mpas_get_pool_pointer ! !------------------------------------------------------------------------------- subroutine dyn_mpas_get_variable_pointer_c0(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array, mpas_pool_get_config + class(mpas_dynamical_core_type), intent(in) :: self character(strkind), pointer, intent(out) :: variable_pointer character(*), intent(in) :: pool_name @@ -3489,6 +3559,10 @@ subroutine dyn_mpas_get_variable_pointer_c0(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_c0 subroutine dyn_mpas_get_variable_pointer_c1(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self character(strkind), pointer, intent(out) :: variable_pointer(:) character(*), intent(in) :: pool_name @@ -3511,6 +3585,10 @@ subroutine dyn_mpas_get_variable_pointer_c1(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_c1 subroutine dyn_mpas_get_variable_pointer_i0(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array, mpas_pool_get_config, mpas_pool_get_dimension + class(mpas_dynamical_core_type), intent(in) :: self integer, pointer, intent(out) :: variable_pointer character(*), intent(in) :: pool_name @@ -3542,6 +3620,10 @@ subroutine dyn_mpas_get_variable_pointer_i0(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_i0 subroutine dyn_mpas_get_variable_pointer_i1(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array, mpas_pool_get_dimension + class(mpas_dynamical_core_type), intent(in) :: self integer, pointer, intent(out) :: variable_pointer(:) character(*), intent(in) :: pool_name @@ -3570,6 +3652,10 @@ subroutine dyn_mpas_get_variable_pointer_i1(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_i1 subroutine dyn_mpas_get_variable_pointer_i2(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self integer, pointer, intent(out) :: variable_pointer(:, :) character(*), intent(in) :: pool_name @@ -3592,6 +3678,10 @@ subroutine dyn_mpas_get_variable_pointer_i2(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_i2 subroutine dyn_mpas_get_variable_pointer_i3(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self integer, pointer, intent(out) :: variable_pointer(:, :, :) character(*), intent(in) :: pool_name @@ -3614,6 +3704,10 @@ subroutine dyn_mpas_get_variable_pointer_i3(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_i3 subroutine dyn_mpas_get_variable_pointer_l0(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_config + class(mpas_dynamical_core_type), intent(in) :: self logical, pointer, intent(out) :: variable_pointer character(*), intent(in) :: pool_name @@ -3640,6 +3734,10 @@ subroutine dyn_mpas_get_variable_pointer_l0(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_l0 subroutine dyn_mpas_get_variable_pointer_r0(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array, mpas_pool_get_config + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer character(*), intent(in) :: pool_name @@ -3668,6 +3766,10 @@ subroutine dyn_mpas_get_variable_pointer_r0(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r0 subroutine dyn_mpas_get_variable_pointer_r1(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:) character(*), intent(in) :: pool_name @@ -3690,6 +3792,10 @@ subroutine dyn_mpas_get_variable_pointer_r1(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r1 subroutine dyn_mpas_get_variable_pointer_r2(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:, :) character(*), intent(in) :: pool_name @@ -3712,6 +3818,10 @@ subroutine dyn_mpas_get_variable_pointer_r2(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r2 subroutine dyn_mpas_get_variable_pointer_r3(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:, :, :) character(*), intent(in) :: pool_name @@ -3734,6 +3844,10 @@ subroutine dyn_mpas_get_variable_pointer_r3(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r3 subroutine dyn_mpas_get_variable_pointer_r4(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:, :, :, :) character(*), intent(in) :: pool_name @@ -3756,6 +3870,10 @@ subroutine dyn_mpas_get_variable_pointer_r4(self, variable_pointer, pool_name, v end subroutine dyn_mpas_get_variable_pointer_r4 subroutine dyn_mpas_get_variable_pointer_r5(self, variable_pointer, pool_name, variable_name, time_level) + ! Module(s) from MPAS. + use mpas_derived_types, only: mpas_pool_type + use mpas_pool_routines, only: mpas_pool_get_array + class(mpas_dynamical_core_type), intent(in) :: self real(rkind), pointer, intent(out) :: variable_pointer(:, :, :, :, :) character(*), intent(in) :: pool_name diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index c6eff173..b703f4b5 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -1,47 +1,10 @@ module dyn_comp + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8, & + len_cs => shr_kind_cs + ! Module(s) from MPAS. use dyn_mpas_subdriver, only: mpas_dynamical_core_type - ! Modules from CAM-SIMA. - use air_composition, only: thermodynamic_active_species_num, & - thermodynamic_active_species_liq_num, & - thermodynamic_active_species_ice_num, & - thermodynamic_active_species_idx, thermodynamic_active_species_idx_dycore, & - thermodynamic_active_species_liq_idx, thermodynamic_active_species_liq_idx_dycore, & - thermodynamic_active_species_ice_idx, thermodynamic_active_species_ice_idx_dycore - use cam_abortutils, only: check_allocate, endrun - use cam_constituents, only: const_name, const_is_dry, const_is_water_species, num_advected, readtrace - use cam_control_mod, only: initial_run - use cam_field_read, only: cam_read_field - use cam_grid_support, only: cam_grid_get_latvals, cam_grid_get_lonvals, cam_grid_id - use cam_initfiles, only: initial_file_get_id, topo_file_get_id - use cam_instance, only: atm_id - use cam_logfile, only: debug_output, debugout_none, iulog - use cam_pio_utils, only: clean_iodesc_list - use dyn_tests_utils, only: vc_height - use dynconst, only: constant_cpd => cpair, constant_g => gravit, constant_p0 => pref, constant_pi => pi, & - constant_rd => rair, constant_rv => rh2o, & - deg_to_rad - use inic_analytic, only: analytic_ic_active, dyn_set_inic_col - use runtime_obj, only: runtime_options - use spmd_utils, only: iam, masterproc, mpicom - use string_utils, only: stringify - use time_manager, only: get_start_date, get_stop_date, get_step_size, get_run_duration, timemgr_get_calendar_cf - use vert_coord, only: pver, pverp - - ! Modules from CCPP. - use cam_ccpp_cap, only: cam_constituents_array - use ccpp_kinds, only: kind_phys - use phys_vars_init_check, only: mark_as_initialized, std_name_len - use physics_types, only: phys_state - - ! Modules from CESM Share. - use shr_file_mod, only: shr_file_getunit - use shr_kind_mod, only: kind_cs => shr_kind_cs, kind_r8 => shr_kind_r8 - use shr_pio_mod, only: shr_pio_getiosys - - ! Modules from external libraries. - use pio, only: file_desc_t, iosystem_desc_t, pio_file_is_open - implicit none private @@ -87,6 +50,11 @@ module dyn_comp !> If `printer` is not supplied, the MPI root rank will print. Otherwise, the designated MPI rank will print instead. !> (KCW, 2024-02-03) subroutine dyn_debug_print(message, variable, printer) + ! Module(s) from CAM-SIMA. + use cam_logfile, only: debug_output, debugout_none, iulog + use spmd_utils, only: iam, masterproc + use string_utils, only: stringify + character(*), intent(in) :: message class(*), optional, intent(in) :: variable(:) integer, optional, intent(in) :: printer @@ -121,10 +89,23 @@ end subroutine dyn_debug_print ! ! Called by `read_namelist` in `src/control/runtime_opts.F90`. subroutine dyn_readnl(namelist_path) + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: endrun + use cam_control_mod, only: initial_run + use cam_instance, only: atm_id + use cam_logfile, only: debug_output, debugout_none, iulog + use spmd_utils, only: mpicom + use time_manager, only: get_start_date, get_stop_date, get_run_duration, timemgr_get_calendar_cf + ! Module(s) from CESM Share. + use shr_file_mod, only: shr_file_getunit + use shr_pio_mod, only: shr_pio_getiosys + ! Module(s) from external libraries. + use pio, only: iosystem_desc_t + character(*), intent(in) :: namelist_path character(*), parameter :: subname = 'dyn_comp::dyn_readnl' - character(kind_cs) :: cam_calendar + character(len_cs) :: cam_calendar integer :: log_unit(2) integer :: start_date_time(6), & ! YYYY, MM, DD, hh, mm, ss. stop_date_time(6), & ! YYYY, MM, DD, hh, mm, ss. @@ -187,8 +168,27 @@ end subroutine dyn_readnl ! ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use air_composition, only: thermodynamic_active_species_num, & + thermodynamic_active_species_liq_num, & + thermodynamic_active_species_ice_num, & + thermodynamic_active_species_idx, thermodynamic_active_species_idx_dycore, & + thermodynamic_active_species_liq_idx, thermodynamic_active_species_liq_idx_dycore, & + thermodynamic_active_species_ice_idx, thermodynamic_active_species_ice_idx_dycore + use cam_abortutils, only: check_allocate + use cam_constituents, only: const_name, const_is_water_species, num_advected, readtrace + use cam_control_mod, only: initial_run + use cam_initfiles, only: initial_file_get_id, topo_file_get_id + use cam_pio_utils, only: clean_iodesc_list use cam_thermo_formula, only: energy_formula_dycore, energy_formula_dycore_mpas + use inic_analytic, only: analytic_ic_active use physics_types, only: dycore_energy_consistency_adjust + use runtime_obj, only: runtime_options + use time_manager, only: get_step_size + ! Module(s) from CCPP. + use phys_vars_init_check, only: std_name_len + ! Module(s) from external libraries. + use pio, only: file_desc_t type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in @@ -306,6 +306,14 @@ end subroutine dyn_init !> Otherwise, if topography file is not used, check that the surface geometric height in MPAS is zero. !> (KCW, 2024-05-10) subroutine check_topography_data(pio_file) + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_field_read, only: cam_read_field + use cam_logfile, only: debug_output, debugout_none + use dynconst, only: constant_g => gravit + ! Module(s) from external libraries. + use pio, only: file_desc_t, pio_file_is_open + type(file_desc_t), pointer, intent(in) :: pio_file character(*), parameter :: subname = 'dyn_comp::check_topography_data' @@ -389,9 +397,15 @@ subroutine set_analytic_initial_condition() !> Initialize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines. !> (KCW, 2024-05-13) subroutine init_shared_variables() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_grid_support, only: cam_grid_get_latvals, cam_grid_get_lonvals, cam_grid_id + use dynconst, only: deg_to_rad + use vert_coord, only: pverp + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::init_shared_variables' - integer :: ierr integer :: i + integer :: ierr integer, pointer :: indextocellid(:) real(kind_r8), pointer :: lat_deg(:), lon_deg(:) @@ -459,9 +473,15 @@ end subroutine final_shared_variables !> Set MPAS state `u` (i.e., horizontal velocity at edge interfaces). !> (KCW, 2024-05-13) subroutine set_mpas_state_u() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_tests_utils, only: vc_height + use inic_analytic, only: dyn_set_inic_col + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_u' - integer :: ierr integer :: i + integer :: ierr real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :) call dyn_debug_print('Setting MPAS state "u"') @@ -522,6 +542,13 @@ end subroutine set_mpas_state_w !> Set MPAS state `scalars` (i.e., constituents). !> (KCW, 2024-05-17) subroutine set_mpas_state_scalars() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use cam_constituents, only: num_advected + use dyn_tests_utils, only: vc_height + use inic_analytic, only: dyn_set_inic_col + use vert_coord, only: pver + ! CCPP standard name of `qv`, which denotes water vapor mixing ratio. character(*), parameter :: constituent_qv_standard_name = & 'water_vapor_mixing_ratio_wrt_dry_air' @@ -588,6 +615,13 @@ end subroutine set_mpas_state_scalars !> Set MPAS state `rho` (i.e., dry air density) and `theta` (i.e., potential temperature). !> (KCW, 2024-05-19) subroutine set_mpas_state_rho_theta() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_tests_utils, only: vc_height + use dynconst, only: constant_p0 => pref, constant_rd => rair, constant_rv => rh2o + use inic_analytic, only: dyn_set_inic_col + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_rho_theta' integer :: i, k integer :: ierr @@ -700,6 +734,11 @@ end subroutine set_mpas_state_rho_theta !> Set MPAS state `rho_base` (i.e., base state dry air density) and `theta_base` (i.e., base state potential temperature). !> (KCW, 2024-05-21) subroutine set_mpas_state_rho_base_theta_base() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dynconst, only: constant_p0 => pref, constant_rd => rair + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_rho_base_theta_base' integer :: i, k integer :: ierr @@ -758,6 +797,9 @@ end subroutine set_mpas_state_rho_base_theta_base !> `t_v` is the mean virtual temperature between `z_1` and `z_2`. !> (KCW, 2024-07-02) pure elemental function p_by_hypsometric_equation(p_1, z_1, t_v, z_2) result(p_2) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_g => gravit, constant_rd => rair + real(kind_r8), intent(in) :: p_1, z_1, t_v, z_2 real(kind_r8) :: p_2 @@ -772,6 +814,9 @@ end function p_by_hypsometric_equation !> Poisson equation. !> (KCW, 2024-07-02) pure elemental function theta_by_poisson_equation(p_1, t_1, p_0) result(t_0) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_cpd => cpair, constant_rd => rair + real(kind_r8), intent(in) :: p_1, t_1, p_0 real(kind_r8) :: t_0 @@ -791,6 +836,15 @@ end subroutine set_analytic_initial_condition !> Some procedures in CAM-SIMA expect constituent states to be dry, while the others expect them to be moist. !> (KCW, 2024-09-26) subroutine dyn_exchange_constituent_state(direction, exchange, conversion) + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_constituents, only: const_is_dry, const_is_water_species, num_advected + use physics_types, only: phys_state + use vert_coord, only: pver + ! Module(s) from CCPP. + use cam_ccpp_cap, only: cam_constituents_array + use ccpp_kinds, only: kind_phys + character(*), intent(in) :: direction logical, intent(in) :: exchange logical, intent(in) :: conversion @@ -933,6 +987,11 @@ end subroutine dyn_exchange_constituent_state !> to prevent physics from attempting to read them from a file. !> (KCW, 2024-05-23) subroutine mark_variable_as_initialized() + ! Module(s) from CAM-SIMA. + use cam_constituents, only: const_name, num_advected + ! Module(s) from CCPP. + use phys_vars_init_check, only: mark_as_initialized + character(*), parameter :: subname = 'dyn_comp::mark_variable_as_initialized' integer :: i diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 14b211f3..1b02d7a3 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -1,31 +1,7 @@ module dyn_coupling - ! Modules from CAM-SIMA. - use cam_abortutils, only: check_allocate, endrun - use cam_constituents, only: const_is_water_species, const_qmin, num_advected - use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update - use cam_thermo_formula, only: energy_formula_dycore_mpas - use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_state, reverse, mpas_dynamical_core, & - ncells_solve - use dynconst, only: constant_cpd => cpair, constant_g => gravit, constant_p0 => pref, & - constant_rd => rair, constant_rv => rh2o - use runtime_obj, only: cam_runtime_opts - use vert_coord, only: pver, pverp - - ! Modules from CCPP. - use cam_ccpp_cap, only: cam_constituents_array, cam_model_const_properties - use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t - use ccpp_kinds, only: kind_phys - use geopotential_temp, only: geopotential_temp_run - use physics_types, only: cappav, cpairv, rairv, zvirv, & - dtime_phys, lagrangian_vertical, & - phys_state, phys_tend - use physics_types, only: cp_or_cv_dycore - use qneg, only: qneg_run - use static_energy, only: update_dry_static_energy_run - use string_utils, only: stringify - - ! Modules from CESM Share. - use shr_kind_mod, only: kind_cx => shr_kind_cx, kind_r8 => shr_kind_r8 + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8, & + len_cx => shr_kind_cx implicit none @@ -38,6 +14,9 @@ module dyn_coupling !> The other coupling direction is implemented by its counterpart, `physics_to_dynamics_coupling`. !> (KCW, 2024-07-31) subroutine dynamics_to_physics_coupling() + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_state, ncells_solve + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling' integer :: column_index integer, allocatable :: is_water_species_index(:) @@ -93,6 +72,12 @@ subroutine dynamics_to_physics_coupling() !> `set_physics_state_column` internal subroutines. !> (KCW, 2024-07-20) subroutine init_shared_variables() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use cam_constituents, only: const_is_water_species, num_advected + use dyn_comp, only: mpas_dynamical_core + use vert_coord, only: pver, pverp + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::init_shared_variables' integer :: i integer :: ierr @@ -198,6 +183,10 @@ end subroutine final_shared_variables !> should be called in pairs. !> (KCW, 2024-07-30) subroutine update_shared_variables(i) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_g => gravit, constant_rd => rair, constant_rv => rh2o + use vert_coord, only: pver, pverp + integer, intent(in) :: i character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::update_shared_variables' @@ -259,6 +248,11 @@ end subroutine update_shared_variables !> This subroutine and `update_shared_variables` should be called in pairs. !> (KCW, 2024-07-30) subroutine set_physics_state_column(i) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: reverse + use dynconst, only: constant_g => gravit + use physics_types, only: phys_state + integer, intent(in) :: i character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_column' @@ -294,8 +288,27 @@ end subroutine set_physics_state_column !> Set variables in the `physics_state` derived type by calling external procedures. !> (KCW, 2024-07-30) subroutine set_physics_state_external() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate, endrun + use cam_constituents, only: const_qmin, num_advected + use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update + use cam_thermo_formula, only: energy_formula_dycore_mpas + use dyn_comp, only: mpas_dynamical_core + use dynconst, only: constant_g => gravit + use physics_types, only: cappav, cp_or_cv_dycore, cpairv, lagrangian_vertical, phys_state, rairv, zvirv + use runtime_obj, only: cam_runtime_opts + use string_utils, only: stringify + use vert_coord, only: pver, pverp + ! Module(s) from CCPP. + use cam_ccpp_cap, only: cam_constituents_array, cam_model_const_properties + use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t + use ccpp_kinds, only: kind_phys + use geopotential_temp, only: geopotential_temp_run + use qneg, only: qneg_run + use static_energy, only: update_dry_static_energy_run + character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_external' - character(kind_cx) :: cerr + character(len_cx) :: cerr integer :: i integer :: ierr real(kind_phys), allocatable :: minimum_constituents(:) @@ -394,6 +407,9 @@ end subroutine dynamics_to_physics_coupling !> The other coupling direction is implemented by its counterpart, `dynamics_to_physics_coupling`. !> (KCW, 2024-09-20) subroutine physics_to_dynamics_coupling() + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_exchange_constituent_state + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling' integer, pointer :: index_qv real(kind_r8), allocatable :: qv_prev(:, :) ! Water vapor mixing ratio (kg kg-1) @@ -415,6 +431,11 @@ subroutine physics_to_dynamics_coupling() !> Initialize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines. !> (KCW, 2024-09-13) subroutine init_shared_variables() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, ncells_solve + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::init_shared_variables' integer :: ierr @@ -457,6 +478,10 @@ end subroutine final_shared_variables !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. !> (KCW, 2024-09-11) subroutine set_mpas_physics_tendency_ru() + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_debug_print, reverse, mpas_dynamical_core, ncells_solve + use physics_types, only: phys_tend + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_ru' integer :: i real(kind_r8), pointer :: u_tendency(:, :), v_tendency(:, :) @@ -484,6 +509,9 @@ end subroutine set_mpas_physics_tendency_ru !> In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. !> (KCW, 2024-09-11) subroutine set_mpas_physics_tendency_rho() + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, ncells_solve + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rho' real(kind_r8), pointer :: rho_tendency(:, :) @@ -506,6 +534,13 @@ end subroutine set_mpas_physics_tendency_rho !> due to physics). In MPAS, a "coupled" variable means that it is multiplied by a vertical metric term, `rho_zz`. !> (KCW, 2024-09-19) subroutine set_mpas_physics_tendency_rtheta() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_comp, only: dyn_debug_print, reverse, mpas_dynamical_core, ncells_solve + use dynconst, only: constant_rd => rair, constant_rv => rh2o + use physics_types, only: dtime_phys, phys_tend + use vert_coord, only: pver + character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rtheta' integer :: i integer :: ierr @@ -592,6 +627,10 @@ end subroutine set_mpas_physics_tendency_rtheta !> `t == t_of_theta_rhod_qv(theta_of_t_rhod_qv(t, rhod, qv), rhod, qv)`. !> (KCW, 2024-09-13) pure elemental function t_of_theta_rhod_qv(theta, rhod, qv) result(t) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_cpd => cpair, constant_p0 => pref, & + constant_rd => rair, constant_rv => rh2o + real(kind_r8), intent(in) :: theta, rhod, qv real(kind_r8) :: t @@ -626,6 +665,10 @@ end function t_of_theta_rhod_qv !> `theta == theta_of_t_rhod_qv(t_of_theta_rhod_qv(theta, rhod, qv), rhod, qv)`. !> (KCW, 2024-09-13) pure elemental function theta_of_t_rhod_qv(t, rhod, qv) result(theta) + ! Module(s) from CAM-SIMA. + use dynconst, only: constant_cpd => cpair, constant_p0 => pref, & + constant_rd => rair, constant_rv => rh2o + real(kind_r8), intent(in) :: t, rhod, qv real(kind_r8) :: theta diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 07277b1c..2356fcd4 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -1,32 +1,11 @@ module dyn_grid - ! Modules from CAM-SIMA. - use cam_abortutils, only: check_allocate, endrun - use cam_constituents, only: num_advected - use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register, & - horiz_coord_t, horiz_coord_create, & - max_hcoordname_len - use cam_history_support, only: add_vert_coord - use cam_initfiles, only: initial_file_get_id + ! Module(s) from CAM-SIMA. + use cam_grid_support, only: max_hcoordname_len use cam_map_utils, only: kind_imap => imap - use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, & - ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels, & - ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max, & - sphere_radius - use dynconst, only: constant_p0 => pref, constant_pi => pi, rad_to_deg, dynconst_init - use physics_column_type, only: kind_pcol, physics_column_t - use physics_grid, only: phys_decomp, phys_grid_init - use ref_pres, only: ref_pres_init - use spmd_utils, only: iam - use std_atm_profile, only: std_atm_pres - use string_utils, only: stringify - use vert_coord, only: pver, pverp, vert_coord_init - - ! Modules from CESM Share. + use physics_column_type, only: kind_pcol + ! Module(s) from CESM Share. use shr_kind_mod, only: kind_r8 => shr_kind_r8 - ! Modules from external libraries. - use pio, only: file_desc_t - implicit none private @@ -52,6 +31,20 @@ module dyn_grid ! ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine model_grid_init() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: endrun + use cam_constituents, only: num_advected + use cam_initfiles, only: initial_file_get_id + use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, & + ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels, & + ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max, & + sphere_radius + use dynconst, only: dynconst_init + use string_utils, only: stringify + use vert_coord, only: pver, vert_coord_init + ! Module(s) from external libraries. + use pio, only: file_desc_t + character(*), parameter :: subname = 'dyn_grid::model_grid_init' type(file_desc_t), pointer :: pio_file @@ -127,6 +120,16 @@ end subroutine model_grid_init !> Initialize reference pressure for use by physics. !> (KCW, 2024-03-25) subroutine init_reference_pressure() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use cam_history_support, only: add_vert_coord + use dyn_comp, only: dyn_debug_print, mpas_dynamical_core + use dynconst, only: constant_p0 => pref + use ref_pres, only: ref_pres_init + use std_atm_profile, only: std_atm_pres + use string_utils, only: stringify + use vert_coord, only: pver, pverp + character(*), parameter :: subname = 'dyn_grid::init_reference_pressure' ! Number of pure pressure levels at model top. integer, parameter :: num_pure_p_lev = 0 @@ -228,6 +231,15 @@ end subroutine init_reference_pressure !> Provide grid and mapping information between global and local indexes to physics by calling `phys_grid_init`. !> (KCW, 2024-03-27) subroutine init_physics_grid() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use dyn_comp, only: mpas_dynamical_core, ncells_global, ncells_solve, sphere_radius + use dynconst, only: constant_pi => pi, rad_to_deg + use physics_column_type, only: physics_column_t + use physics_grid, only: phys_grid_init + use spmd_utils, only: iam + use string_utils, only: stringify + character(*), parameter :: subname = 'dyn_grid::init_physics_grid' character(max_hcoordname_len), allocatable :: dyn_attribute_name(:) integer :: hdim1_d, hdim2_d ! First and second horizontal dimensions of physics grid. @@ -308,6 +320,17 @@ end subroutine init_physics_grid !> * "mpas_vertex": Grid that is centered at MPAS "vertex" points. !> (KCW, 2024-03-28) subroutine define_cam_grid() + ! Module(s) from CAM-SIMA. + use cam_abortutils, only: check_allocate + use cam_grid_support, only: cam_grid_attribute_register, cam_grid_register, & + horiz_coord_create, horiz_coord_t + use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, & + ncells_global, nedges_global, nvertices_global, & + ncells_solve, nedges_solve, nvertices_solve, & + sphere_radius + use dynconst, only: constant_pi => pi, rad_to_deg + use string_utils, only: stringify + character(*), parameter :: subname = 'dyn_grid::define_cam_grid' integer :: i integer :: ierr @@ -500,6 +523,9 @@ end subroutine define_cam_grid !> Helper function for returning grid id given its name. !> (KCW, 2024-03-27) pure function dyn_grid_id(name) + ! Module(s) from CAM-SIMA. + use physics_grid, only: phys_decomp + character(*), intent(in) :: name integer :: dyn_grid_id diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90 index d6d86503..dfaf2c18 100644 --- a/src/dynamics/mpas/stepon.F90 +++ b/src/dynamics/mpas/stepon.F90 @@ -1,13 +1,5 @@ module stepon - ! Modules from CAM-SIMA. - use camsrfexch, only: cam_out_t - use dyn_comp, only: dyn_import_t, dyn_export_t, dyn_run, dyn_final - use dyn_coupling, only: dynamics_to_physics_coupling, physics_to_dynamics_coupling - use physics_types, only: physics_state, physics_tend - use runtime_obj, only: runtime_options - use time_manager, only: get_step_size - - ! Modules from CCPP. + ! Module(s) from CCPP. use ccpp_kinds, only: kind_phys implicit none @@ -22,6 +14,10 @@ module stepon contains ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine stepon_init(cam_runtime_opts, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_export_t, dyn_import_t + use runtime_obj, only: runtime_options + type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in type(dyn_export_t), intent(in) :: dyn_out @@ -29,6 +25,13 @@ end subroutine stepon_init ! Called by `cam_timestep_init` in `src/control/cam_comp.F90`. subroutine stepon_timestep_init(dtime_phys, cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_export_t, dyn_import_t + use dyn_coupling, only: dynamics_to_physics_coupling + use physics_types, only: physics_state, physics_tend + use runtime_obj, only: runtime_options + use time_manager, only: get_step_size + real(kind_phys), intent(out) :: dtime_phys type(runtime_options), intent(in) :: cam_runtime_opts type(physics_state), intent(in) :: phys_state @@ -44,6 +47,12 @@ end subroutine stepon_timestep_init ! Called by `cam_run2` in `src/control/cam_comp.F90`. subroutine stepon_run2(cam_runtime_opts, phys_state, phys_tend, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_export_t, dyn_import_t + use dyn_coupling, only: physics_to_dynamics_coupling + use physics_types, only: physics_state, physics_tend + use runtime_obj, only: runtime_options + type(runtime_options), intent(in) :: cam_runtime_opts type(physics_state), intent(in) :: phys_state type(physics_tend), intent(in) :: phys_tend @@ -55,6 +64,12 @@ end subroutine stepon_run2 ! Called by `cam_run3` in `src/control/cam_comp.F90`. subroutine stepon_run3(dtime_phys, cam_runtime_opts, cam_out, phys_state, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use camsrfexch, only: cam_out_t + use dyn_comp, only: dyn_export_t, dyn_import_t, dyn_run + use physics_types, only: physics_state + use runtime_obj, only: runtime_options + real(kind_phys), intent(in) :: dtime_phys type(runtime_options), intent(in) :: cam_runtime_opts type(cam_out_t), intent(in) :: cam_out @@ -67,6 +82,10 @@ end subroutine stepon_run3 ! Called by `cam_final` in `src/control/cam_comp.F90`. subroutine stepon_final(cam_runtime_opts, dyn_in, dyn_out) + ! Module(s) from CAM-SIMA. + use dyn_comp, only: dyn_export_t, dyn_import_t, dyn_final + use runtime_obj, only: runtime_options + type(runtime_options), intent(in) :: cam_runtime_opts type(dyn_import_t), intent(in) :: dyn_in type(dyn_export_t), intent(in) :: dyn_out From 5a4d3bc69bfae94b8c6e1b91daee982f3cfca145 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 13 Nov 2024 17:35:03 -0700 Subject: [PATCH 09/24] Support MPAS dynamical core in single precision mode It is now possible to run CAM-SIMA with MPAS dynamical core in single precision mode. It can be enabled by defining the `SINGLE_PRECISION` macro in `CPPFLAGS`. --- src/dynamics/mpas/Makefile.in.CESM | 4 + .../mpas/driver/dyn_mpas_subdriver.F90 | 4 + src/dynamics/mpas/dyn_comp.F90 | 86 +++++++++---------- src/dynamics/mpas/dyn_coupling.F90 | 61 ++++++------- src/dynamics/mpas/dyn_grid.F90 | 55 ++++++------ 5 files changed, 112 insertions(+), 98 deletions(-) diff --git a/src/dynamics/mpas/Makefile.in.CESM b/src/dynamics/mpas/Makefile.in.CESM index 8e4ba049..b116c79e 100644 --- a/src/dynamics/mpas/Makefile.in.CESM +++ b/src/dynamics/mpas/Makefile.in.CESM @@ -49,6 +49,10 @@ export CPPFLAGS := -D_MPI \ ifneq ($(strip $(PIODEF)),) export CPPFLAGS += $(strip $(PIODEF)) endif +# Uncomment below for MPI Fortran 2008 interface support. There is currently no corresponding option in CIME to enable this. +# export CPPFLAGS += -DMPAS_USE_MPI_F08 +# Uncomment below for single precision mode support. There is currently no corresponding option in CIME to enable this. +# export CPPFLAGS += -DSINGLE_PRECISION export LINKER := $(strip $(FC)) export SCC := $(strip $(CC)) export SCXX := $(strip $(CXX)) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index e6677492..ea793dcf 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -22,6 +22,7 @@ module dyn_mpas_subdriver implicit none private + public :: mpas_dynamical_core_real_kind public :: mpas_dynamical_core_type abstract interface @@ -33,6 +34,9 @@ subroutine model_error_if(message, file, line) end subroutine model_error_if end interface + !> The native floating-point precision of MPAS dynamical core. + integer, parameter :: mpas_dynamical_core_real_kind = rkind + !> The "class" of MPAS dynamical core. !> Important data structures like states of MPAS dynamical core are encapsulated inside this derived type to prevent misuse. !> Type-bound procedures provide well-defined APIs for CAM-SIMA to interact with MPAS dynamical core. diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index b703f4b5..4e28bc8a 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -3,7 +3,7 @@ module dyn_comp use shr_kind_mod, only: kind_r8 => shr_kind_r8, & len_cs => shr_kind_cs ! Module(s) from MPAS. - use dyn_mpas_subdriver, only: mpas_dynamical_core_type + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind, mpas_dynamical_core_type implicit none @@ -44,7 +44,7 @@ module dyn_comp ! Local and global mesh dimensions of MPAS dynamical core. integer :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels integer :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max - real(kind_r8) :: sphere_radius + real(kind_dyn_mpas) :: sphere_radius contains !> Print a debug message with optionally the value(s) of a variable. !> If `printer` is not supplied, the MPI root rank will print. Otherwise, the designated MPI rank will print instead. @@ -322,7 +322,7 @@ subroutine check_topography_data(pio_file) real(kind_r8), parameter :: error_tolerance = 1.0E-3_kind_r8 ! Error tolerance for consistency check. real(kind_r8), allocatable :: surface_geometric_height(:) ! Computed from topography file. real(kind_r8), allocatable :: surface_geopotential(:) ! Read from topography file. - real(kind_r8), pointer :: zgrid(:, :) ! From MPAS. Geometric height (meters) at layer interfaces. + real(kind_dyn_mpas), pointer :: zgrid(:, :) ! From MPAS. Geometric height (meters) at layer interfaces. nullify(zgrid) @@ -354,7 +354,7 @@ subroutine check_topography_data(pio_file) surface_geometric_height(:) = surface_geopotential(:) / constant_g ! Surface geometric height in MPAS should match the values in topography file. - if (any(abs(zgrid(1, 1:ncells_solve) - surface_geometric_height) > error_tolerance)) then + if (any(abs(real(zgrid(1, 1:ncells_solve), kind_r8) - surface_geometric_height(:)) > error_tolerance)) then call endrun('Surface geometric height in MPAS is not consistent with topography data', subname, __LINE__) end if @@ -364,7 +364,7 @@ subroutine check_topography_data(pio_file) call dyn_debug_print('Topography file is not used') ! Surface geometric height in MPAS should be zero. - if (any(abs(zgrid(1, 1:ncells_solve)) > error_tolerance)) then + if (any(abs(real(zgrid(1, 1:ncells_solve), kind_r8)) > error_tolerance)) then call endrun('Surface geometric height in MPAS is not zero', subname, __LINE__) end if end if @@ -379,10 +379,10 @@ subroutine set_analytic_initial_condition() integer, allocatable :: global_grid_index(:) real(kind_r8), allocatable :: buffer_2d_real(:, :), buffer_3d_real(:, :, :) real(kind_r8), allocatable :: lat_rad(:), lon_rad(:) - real(kind_r8), allocatable :: z_int(:, :) ! Geometric height (meters) at layer interfaces. - ! Dimension and vertical index orders follow CAM-SIMA convention. - real(kind_r8), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces. - ! Dimension and vertical index orders follow MPAS convention. + real(kind_r8), allocatable :: z_int(:, :) ! Geometric height (meters) at layer interfaces. + ! Dimension and vertical index orders follow CAM-SIMA convention. + real(kind_dyn_mpas), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces. + ! Dimension and vertical index orders follow MPAS convention. call init_shared_variables() @@ -454,7 +454,7 @@ subroutine init_shared_variables() ! Vertical index order is reversed between CAM-SIMA and MPAS. do i = 1, ncells_solve - z_int(i, :) = reverse(zgrid(:, i)) + z_int(i, :) = reverse(real(zgrid(:, i), kind_r8)) end do end subroutine init_shared_variables @@ -482,7 +482,7 @@ subroutine set_mpas_state_u() character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_u' integer :: i integer :: ierr - real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :) + real(kind_dyn_mpas), pointer :: ucellzonal(:, :), ucellmeridional(:, :) call dyn_debug_print('Setting MPAS state "u"') @@ -500,7 +500,7 @@ subroutine set_mpas_state_u() ! Vertical index order is reversed between CAM-SIMA and MPAS. do i = 1, ncells_solve - ucellzonal(:, i) = reverse(buffer_2d_real(i, :)) + ucellzonal(:, i) = real(reverse(buffer_2d_real(i, :)), kind_dyn_mpas) end do buffer_2d_real(:, :) = 0.0_kind_r8 @@ -509,7 +509,7 @@ subroutine set_mpas_state_u() ! Vertical index order is reversed between CAM-SIMA and MPAS. do i = 1, ncells_solve - ucellmeridional(:, i) = reverse(buffer_2d_real(i, :)) + ucellmeridional(:, i) = real(reverse(buffer_2d_real(i, :)), kind_dyn_mpas) end do deallocate(buffer_2d_real) @@ -523,7 +523,7 @@ end subroutine set_mpas_state_u !> (KCW, 2024-05-13) subroutine set_mpas_state_w() character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_w' - real(kind_r8), pointer :: w(:, :) + real(kind_dyn_mpas), pointer :: w(:, :) call dyn_debug_print('Setting MPAS state "w"') @@ -531,7 +531,7 @@ subroutine set_mpas_state_w() call mpas_dynamical_core % get_variable_pointer(w, 'state', 'w', time_level=1) - w(:, 1:ncells_solve) = 0.0_kind_r8 + w(:, 1:ncells_solve) = 0.0_kind_dyn_mpas nullify(w) @@ -558,7 +558,7 @@ subroutine set_mpas_state_scalars() integer :: ierr integer, allocatable :: constituent_index(:) integer, pointer :: index_qv - real(kind_r8), pointer :: scalars(:, :, :) + real(kind_dyn_mpas), pointer :: scalars(:, :, :) call dyn_debug_print('Setting MPAS state "scalars"') @@ -585,7 +585,7 @@ subroutine set_mpas_state_scalars() do j = 1, num_advected ! Vertical index order is reversed between CAM-SIMA and MPAS. scalars(j, :, i) = & - reverse(buffer_3d_real(i, :, mpas_dynamical_core % map_constituent_index(j))) + real(reverse(buffer_3d_real(i, :, mpas_dynamical_core % map_constituent_index(j))), kind_dyn_mpas) end do end do @@ -599,7 +599,7 @@ subroutine set_mpas_state_scalars() ! Convert specific humidity to water vapor mixing ratio. scalars(index_qv, :, 1:ncells_solve) = & - scalars(index_qv, :, 1:ncells_solve) / (1.0_kind_r8 - scalars(index_qv, :, 1:ncells_solve)) + scalars(index_qv, :, 1:ncells_solve) / (1.0_kind_dyn_mpas - scalars(index_qv, :, 1:ncells_solve)) end if deallocate(buffer_3d_real) @@ -635,9 +635,9 @@ subroutine set_mpas_state_rho_theta() real(kind_r8), allocatable :: tm_mid_col(:) ! Modified "moist" temperature (K) at layer midpoints of each column. ! Be advised that it is not virtual temperature. ! See doi:10.5065/1DFH-6P97 and doi:10.1175/MWR-D-11-00215.1 for details. - real(kind_r8), pointer :: rho(:, :) - real(kind_r8), pointer :: theta(:, :) - real(kind_r8), pointer :: scalars(:, :, :) + real(kind_dyn_mpas), pointer :: rho(:, :) + real(kind_dyn_mpas), pointer :: theta(:, :) + real(kind_dyn_mpas), pointer :: scalars(:, :, :) call dyn_debug_print('Setting MPAS state "rho" and "theta"') @@ -686,16 +686,16 @@ subroutine set_mpas_state_rho_theta() ! Set `rho` and `theta` column by column. This way, peak memory usage can be reduced. do i = 1, ncells_solve - qv_mid_col(:) = scalars(index_qv, :, i) + qv_mid_col(:) = real(scalars(index_qv, :, i), kind_r8) tm_mid_col(:) = t_mid(:, i) * (1.0_kind_r8 + constant_rv / constant_rd * qv_mid_col(:)) ! Piecewise integrate hypsometric equation to derive `p_mid_col(1)`. ! The formulation used here is exact. p_mid_col(1) = p_by_hypsometric_equation( & p_sfc(i), & - zgrid(1, i), & + real(zgrid(1, i), kind_r8), & tm_mid_col(1) / (1.0_kind_r8 + qv_mid_col(1)), & - 0.5_kind_r8 * (zgrid(2, i) + zgrid(1, i))) + 0.5_kind_r8 * real(zgrid(2, i) + zgrid(1, i), kind_r8)) ! Piecewise integrate hypsometric equation to derive subsequent `p_mid_col(k)`. ! The formulation used here is exact. @@ -703,16 +703,16 @@ subroutine set_mpas_state_rho_theta() p_mid_col(k) = p_by_hypsometric_equation( & p_by_hypsometric_equation( & p_mid_col(k - 1), & - 0.5_kind_r8 * (zgrid(k, i) + zgrid(k - 1, i)), & + 0.5_kind_r8 * real(zgrid(k, i) + zgrid(k - 1, i), kind_r8), & tm_mid_col(k - 1) / (1.0_kind_r8 + qv_mid_col(k - 1)), & - zgrid(k, i)), & - zgrid(k, i), & + real(zgrid(k, i), kind_r8)), & + real(zgrid(k, i), kind_r8), & tm_mid_col(k) / (1.0_kind_r8 + qv_mid_col(k)), & - 0.5_kind_r8 * (zgrid(k + 1, i) + zgrid(k, i))) + 0.5_kind_r8 * real(zgrid(k + 1, i) + zgrid(k, i), kind_r8)) end do - rho(:, i) = p_mid_col(:) / (constant_rd * tm_mid_col(:)) - theta(:, i) = theta_by_poisson_equation(p_mid_col, t_mid(:, i), constant_p0) + rho(:, i) = real(p_mid_col(:) / (constant_rd * tm_mid_col(:)), kind_dyn_mpas) + theta(:, i) = real(theta_by_poisson_equation(p_mid_col, t_mid(:, i), constant_p0), kind_dyn_mpas) end do deallocate(p_mid_col) @@ -745,9 +745,9 @@ subroutine set_mpas_state_rho_base_theta_base() real(kind_r8), parameter :: t_base = 250.0_kind_r8 ! Base state temperature (K) of dry isothermal atmosphere. ! The value used here is identical to MPAS. real(kind_r8), allocatable :: p_base(:) ! Base state pressure (Pa) at layer midpoints of each column. - real(kind_r8), pointer :: rho_base(:, :) - real(kind_r8), pointer :: theta_base(:, :) - real(kind_r8), pointer :: zz(:, :) + real(kind_dyn_mpas), pointer :: rho_base(:, :) + real(kind_dyn_mpas), pointer :: theta_base(:, :) + real(kind_dyn_mpas), pointer :: zz(:, :) call dyn_debug_print('Setting MPAS state "rho_base" and "theta_base"') @@ -771,11 +771,11 @@ subroutine set_mpas_state_rho_base_theta_base() constant_p0, & 0.0_kind_r8, & t_base, & - 0.5_kind_r8 * (zgrid(k + 1, i) + zgrid(k, i))) + 0.5_kind_r8 * real(zgrid(k + 1, i) + zgrid(k, i), kind_r8)) end do - rho_base(:, i) = p_base(:) / (constant_rd * t_base * zz(:, i)) - theta_base(:, i) = theta_by_poisson_equation(p_base, t_base, constant_p0) + rho_base(:, i) = real(p_base(:) / (constant_rd * t_base * real(zz(:, i), kind_r8)), kind_dyn_mpas) + theta_base(:, i) = real(theta_by_poisson_equation(p_base, t_base, constant_p0), kind_dyn_mpas) end do deallocate(p_base) @@ -857,7 +857,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) logical, allocatable :: is_water_species(:) real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory. real(kind_r8), allocatable :: sigma_all_q(:) ! Summation of all water species mixing ratios. - real(kind_r8), pointer :: scalars(:, :, :) ! This points to MPAS memory. + real(kind_dyn_mpas), pointer :: scalars(:, :, :) ! This points to MPAS memory. select case (trim(adjustl(direction))) case ('e', 'export') @@ -932,13 +932,13 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) if (exchange) then ! Vertical index order is reversed between CAM-SIMA and MPAS. scalars(j, :, i) = & - reverse(constituents(i, :, mpas_dynamical_core % map_constituent_index(j))) + real(reverse(constituents(i, :, mpas_dynamical_core % map_constituent_index(j))), kind_dyn_mpas) end if if (conversion .and. is_conversion_needed(mpas_dynamical_core % map_constituent_index(j))) then ! Equation 8 in doi:10.1029/2017MS001257. scalars(j, :, i) = & - scalars(j, :, i) * sigma_all_q(:) + real(real(scalars(j, :, i), kind_r8) * sigma_all_q(:), kind_dyn_mpas) end if end do end do @@ -949,7 +949,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) do i = 1, ncells_solve if (conversion .and. any(is_conversion_needed)) then ! The summation term of equation 8 in doi:10.1029/2017MS001257. - sigma_all_q(:) = reverse(1.0_kind_r8 + sum(scalars(is_water_species_index, :, i), 1)) + sigma_all_q(:) = reverse(1.0_kind_r8 + sum(real(scalars(is_water_species_index, :, i), kind_r8), 1)) end if ! `j` is indexing into `constituents`, so it is regarded as constituent index. @@ -957,7 +957,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) if (exchange) then ! Vertical index order is reversed between CAM-SIMA and MPAS. constituents(i, :, j) = & - reverse(scalars(mpas_dynamical_core % map_mpas_scalar_index(j), :, i)) + reverse(real(scalars(mpas_dynamical_core % map_mpas_scalar_index(j), :, i), kind_r8)) end if if (conversion .and. is_conversion_needed(j)) then @@ -1084,7 +1084,7 @@ subroutine dyn_variable_dump() character(*), parameter :: subname = 'dyn_comp::dyn_variable_dump' integer :: ierr integer :: pio_ioformat, pio_iotype - real(kind_r8), pointer :: surface_pressure(:) + real(kind_dyn_mpas), pointer :: surface_pressure(:) type(file_desc_t), pointer :: pio_file type(iosystem_desc_t), pointer :: pio_iosystem @@ -1094,7 +1094,7 @@ subroutine dyn_variable_dump() call mpas_dynamical_core % get_variable_pointer(surface_pressure, 'diag', 'surface_pressure') - surface_pressure(1:ncells_solve) = phys_state % ps(:) + surface_pressure(1:ncells_solve) = real(phys_state % ps(:), kind_dyn_mpas) nullify(surface_pressure) diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 1b02d7a3..322e6a01 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -2,6 +2,8 @@ module dyn_coupling ! Module(s) from CESM Share. use shr_kind_mod, only: kind_r8 => shr_kind_r8, & len_cx => shr_kind_cx + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind implicit none @@ -43,13 +45,13 @@ subroutine dynamics_to_physics_coupling() real(kind_r8), allocatable :: u_mid_col(:), & ! Eastward wind velocity (m s-1). v_mid_col(:), & ! Northward wind velocity (m s-1). omega_mid_col(:) ! Vertical wind velocity (Pa s-1). - real(kind_r8), pointer :: exner(:, :) - real(kind_r8), pointer :: rho_zz(:, :) - real(kind_r8), pointer :: scalars(:, :, :) - real(kind_r8), pointer :: theta_m(:, :) - real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :), w(:, :) - real(kind_r8), pointer :: zgrid(:, :) - real(kind_r8), pointer :: zz(:, :) + real(kind_dyn_mpas), pointer :: exner(:, :) + real(kind_dyn_mpas), pointer :: rho_zz(:, :) + real(kind_dyn_mpas), pointer :: scalars(:, :, :) + real(kind_dyn_mpas), pointer :: theta_m(:, :) + real(kind_dyn_mpas), pointer :: ucellzonal(:, :), ucellmeridional(:, :), w(:, :) + real(kind_dyn_mpas), pointer :: zgrid(:, :) + real(kind_dyn_mpas), pointer :: zz(:, :) call init_shared_variables() @@ -198,10 +200,10 @@ subroutine update_shared_variables(i) ! Compute thermodynamic variables. ! By definition. - z_int_col(:) = zgrid(:, i) + z_int_col(:) = real(zgrid(:, i), kind_r8) dz_col(:) = z_int_col(2:pverp) - z_int_col(1:pver) - qv_mid_col(:) = scalars(index_qv, :, i) - rhod_mid_col(:) = rho_zz(:, i) * zz(:, i) + qv_mid_col(:) = real(scalars(index_qv, :, i), kind_r8) + rhod_mid_col(:) = real(rho_zz(:, i) * zz(:, i), kind_r8) ! Equation 5 in doi:10.1029/2017MS001257. rho_mid_col(:) = rhod_mid_col(:) * sigma_all_q_mid_col(:) @@ -211,7 +213,7 @@ subroutine update_shared_variables(i) dp_col(:) = -rho_mid_col(:) * constant_g * dz_col(:) ! By definition of Exner function. Also see below. - tm_mid_col(:) = theta_m(:, i) * exner(:, i) + tm_mid_col(:) = real(theta_m(:, i) * exner(:, i), kind_r8) ! The paragraph below equation 2.7 in doi:10.5065/1DFH-6P97. ! The paragraph below equation 2 in doi:10.1175/MWR-D-11-00215.1. @@ -239,9 +241,9 @@ subroutine update_shared_variables(i) ! Compute momentum variables. ! By definition. - u_mid_col(:) = ucellzonal(:, i) - v_mid_col(:) = ucellmeridional(:, i) - omega_mid_col(:) = -rhod_mid_col(:) * constant_g * 0.5_kind_r8 * (w(1:pver, i) + w(2:pverp, i)) + u_mid_col(:) = real(ucellzonal(:, i), kind_r8) + v_mid_col(:) = real(ucellmeridional(:, i), kind_r8) + omega_mid_col(:) = -rhod_mid_col(:) * constant_g * 0.5_kind_r8 * real(w(1:pver, i) + w(2:pverp, i), kind_r8) end subroutine update_shared_variables !> Set variables for the specific column, indicated by `i`, in the `physics_state` derived type. @@ -414,9 +416,9 @@ subroutine physics_to_dynamics_coupling() integer, pointer :: index_qv real(kind_r8), allocatable :: qv_prev(:, :) ! Water vapor mixing ratio (kg kg-1) ! before being updated by physics. - real(kind_r8), pointer :: rho_zz(:, :) - real(kind_r8), pointer :: scalars(:, :, :) - real(kind_r8), pointer :: zz(:, :) + real(kind_dyn_mpas), pointer :: rho_zz(:, :) + real(kind_dyn_mpas), pointer :: scalars(:, :, :) + real(kind_dyn_mpas), pointer :: zz(:, :) call init_shared_variables() @@ -458,7 +460,7 @@ subroutine init_shared_variables() ! Save water vapor mixing ratio before being updated by physics because `set_mpas_physics_tendency_rtheta` ! needs it. This must be done before calling `dyn_exchange_constituent_state`. - qv_prev(:, :) = scalars(index_qv, :, 1:ncells_solve) + qv_prev(:, :) = real(scalars(index_qv, :, 1:ncells_solve), kind_r8) end subroutine init_shared_variables !> Finalize variables that are shared and repeatedly used by the `set_mpas_physics_tendency_*` internal subroutines. @@ -484,7 +486,7 @@ subroutine set_mpas_physics_tendency_ru() character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_ru' integer :: i - real(kind_r8), pointer :: u_tendency(:, :), v_tendency(:, :) + real(kind_dyn_mpas), pointer :: u_tendency(:, :), v_tendency(:, :) call dyn_debug_print('Setting MPAS physics tendency "tend_ru_physics"') @@ -496,8 +498,8 @@ subroutine set_mpas_physics_tendency_ru() ! Vertical index order is reversed between CAM-SIMA and MPAS. ! Always call `reverse` when assigning anything to/from the `physics_tend` derived type. do i = 1, ncells_solve - u_tendency(:, i) = reverse(phys_tend % dudt_total(i, :)) * rho_zz(:, i) - v_tendency(:, i) = reverse(phys_tend % dvdt_total(i, :)) * rho_zz(:, i) + u_tendency(:, i) = real(reverse(phys_tend % dudt_total(i, :)) * real(rho_zz(:, i), kind_r8), kind_dyn_mpas) + v_tendency(:, i) = real(reverse(phys_tend % dvdt_total(i, :)) * real(rho_zz(:, i), kind_r8), kind_dyn_mpas) end do nullify(u_tendency, v_tendency) @@ -513,7 +515,7 @@ subroutine set_mpas_physics_tendency_rho() use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, ncells_solve character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling::set_mpas_physics_tendency_rho' - real(kind_r8), pointer :: rho_tendency(:, :) + real(kind_dyn_mpas), pointer :: rho_tendency(:, :) call dyn_debug_print('Setting MPAS physics tendency "tend_rho_physics"') @@ -522,7 +524,7 @@ subroutine set_mpas_physics_tendency_rho() call mpas_dynamical_core % get_variable_pointer(rho_tendency, 'tend_physics', 'tend_rho_physics') ! The material derivative of `rho` (i.e., dry air density) is zero for incompressible fluid. - rho_tendency(:, 1:ncells_solve) = 0.0_kind_r8 + rho_tendency(:, 1:ncells_solve) = 0.0_kind_dyn_mpas nullify(rho_tendency) @@ -553,8 +555,8 @@ subroutine set_mpas_physics_tendency_rtheta() real(kind_r8), allocatable :: t_col_prev(:), t_col_curr(:) ! Temperature (K). real(kind_r8), allocatable :: theta_col_prev(:), theta_col_curr(:) ! Potential temperature (K). real(kind_r8), allocatable :: thetam_col_prev(:), thetam_col_curr(:) ! Modified "moist" potential temperature (K). - real(kind_r8), pointer :: theta_m(:, :) - real(kind_r8), pointer :: theta_m_tendency(:, :) + real(kind_dyn_mpas), pointer :: theta_m(:, :) + real(kind_dyn_mpas), pointer :: theta_m_tendency(:, :) call dyn_debug_print('Setting MPAS physics tendency "tend_rtheta_physics"') @@ -591,11 +593,11 @@ subroutine set_mpas_physics_tendency_rtheta() ! Set `theta_m_tendency` column by column. This way, peak memory usage can be reduced. do i = 1, ncells_solve - qv_col_curr(:) = scalars(index_qv, :, i) + qv_col_curr(:) = real(scalars(index_qv, :, i), kind_r8) qv_col_prev(:) = qv_prev(:, i) - rhod_col(:) = rho_zz(:, i) * zz(:, i) + rhod_col(:) = real(rho_zz(:, i) * zz(:, i), kind_r8) - thetam_col_prev(:) = theta_m(:, i) + thetam_col_prev(:) = real(theta_m(:, i), kind_r8) theta_col_prev(:) = thetam_col_prev(:) / (1.0_kind_r8 + constant_rv / constant_rd * qv_col_prev(:)) t_col_prev(:) = t_of_theta_rhod_qv(theta_col_prev, rhod_col, qv_col_prev) @@ -605,7 +607,8 @@ subroutine set_mpas_physics_tendency_rtheta() theta_col_curr(:) = theta_of_t_rhod_qv(t_col_curr, rhod_col, qv_col_curr) thetam_col_curr(:) = theta_col_curr(:) * (1.0_kind_r8 + constant_rv / constant_rd * qv_col_curr(:)) - theta_m_tendency(:, i) = (thetam_col_curr(:) - thetam_col_prev(:)) * rho_zz(:, i) / dtime_phys + theta_m_tendency(:, i) = & + real((thetam_col_curr(:) - thetam_col_prev(:)) * real(rho_zz(:, i), kind_r8) / dtime_phys, kind_dyn_mpas) end do deallocate(qv_col_prev, qv_col_curr) diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 2356fcd4..e6780933 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -5,6 +5,8 @@ module dyn_grid use physics_column_type, only: kind_pcol ! Module(s) from CESM Share. use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind implicit none @@ -144,7 +146,7 @@ subroutine init_reference_pressure() ! `dzw` denotes the delta/difference between `zw`. ! `rdzw` denotes the reciprocal of `dzw`. real(kind_r8), allocatable :: dzw(:) - real(kind_r8), pointer :: rdzw(:) + real(kind_dyn_mpas), pointer :: rdzw(:) real(kind_r8), pointer :: zu(:) ! CANNOT be safely deallocated because `add_vert_coord` ! just uses pointers to point at it internally. real(kind_r8), pointer :: zw(:) ! CANNOT be safely deallocated because `add_vert_coord` @@ -160,7 +162,7 @@ subroutine init_reference_pressure() allocate(dzw(pver), stat=ierr) call check_allocate(ierr, subname, 'dzw(pver)', 'dyn_grid', __LINE__) - dzw(:) = 1.0_kind_r8 / rdzw + dzw(:) = 1.0_kind_r8 / real(rdzw(:), kind_r8) nullify(rdzw) @@ -245,10 +247,10 @@ subroutine init_physics_grid() integer :: hdim1_d, hdim2_d ! First and second horizontal dimensions of physics grid. integer :: i integer :: ierr - integer, pointer :: indextocellid(:) ! Global indexes of cell centers. - real(kind_r8), pointer :: areacell(:) ! Cell areas (square meters). - real(kind_r8), pointer :: latcell(:) ! Cell center latitudes (radians). - real(kind_r8), pointer :: loncell(:) ! Cell center longitudes (radians). + integer, pointer :: indextocellid(:) ! Global indexes of cell centers. + real(kind_dyn_mpas), pointer :: areacell(:) ! Cell areas (square meters). + real(kind_dyn_mpas), pointer :: latcell(:) ! Cell center latitudes (radians). + real(kind_dyn_mpas), pointer :: loncell(:) ! Cell center longitudes (radians). type(physics_column_t), allocatable :: dyn_column(:) ! Grid and mapping information between global and local indexes. nullify(areacell) @@ -273,12 +275,13 @@ subroutine init_physics_grid() ! Column information. dyn_column(i) % lat_rad = real(latcell(i), kind_pcol) dyn_column(i) % lon_rad = real(loncell(i), kind_pcol) - dyn_column(i) % lat_deg = real(latcell(i) * rad_to_deg, kind_pcol) - dyn_column(i) % lon_deg = real(loncell(i) * rad_to_deg, kind_pcol) + dyn_column(i) % lat_deg = real(real(latcell(i), kind_r8) * rad_to_deg, kind_pcol) + dyn_column(i) % lon_deg = real(real(loncell(i), kind_r8) * rad_to_deg, kind_pcol) ! Cell areas normalized to unit sphere. dyn_column(i) % area = real(areacell(i) / (sphere_radius ** 2), kind_pcol) ! Cell weights normalized to unity. - dyn_column(i) % weight = real(areacell(i) / (4.0_kind_r8 * constant_pi * sphere_radius ** 2), kind_pcol) + dyn_column(i) % weight = & + real(real(areacell(i), kind_r8) / (4.0_kind_r8 * constant_pi * real(sphere_radius, kind_r8) ** 2), kind_pcol) ! File decomposition. ! For unstructured grid, `coord_indices` is not used by `phys_grid_init`. @@ -337,13 +340,13 @@ subroutine define_cam_grid() integer, pointer :: indextocellid(:) ! Global indexes of cell centers. integer, pointer :: indextoedgeid(:) ! Global indexes of edge nodes. integer, pointer :: indextovertexid(:) ! Global indexes of vertex nodes. - real(kind_r8), pointer :: areacell(:) ! Cell areas (square meters). - real(kind_r8), pointer :: latcell(:) ! Cell center latitudes (radians). - real(kind_r8), pointer :: latedge(:) ! Edge node latitudes (radians). - real(kind_r8), pointer :: latvertex(:) ! Vertex node latitudes (radians). - real(kind_r8), pointer :: loncell(:) ! Cell center longitudes (radians). - real(kind_r8), pointer :: lonedge(:) ! Edge node longitudes (radians). - real(kind_r8), pointer :: lonvertex(:) ! Vertex node longitudes (radians). + real(kind_dyn_mpas), pointer :: areacell(:) ! Cell areas (square meters). + real(kind_dyn_mpas), pointer :: latcell(:) ! Cell center latitudes (radians). + real(kind_dyn_mpas), pointer :: latedge(:) ! Edge node latitudes (radians). + real(kind_dyn_mpas), pointer :: latvertex(:) ! Vertex node latitudes (radians). + real(kind_dyn_mpas), pointer :: loncell(:) ! Cell center longitudes (radians). + real(kind_dyn_mpas), pointer :: lonedge(:) ! Edge node longitudes (radians). + real(kind_dyn_mpas), pointer :: lonvertex(:) ! Vertex node longitudes (radians). ! Global grid indexes. CAN be safely deallocated because its values are copied internally by ! `cam_grid_attribute_register` and `horiz_coord_create`. @@ -390,9 +393,9 @@ subroutine define_cam_grid() global_grid_index(:) = int(indextocellid(1:ncells_solve), kind_imap) lat_coord => horiz_coord_create('latCell', 'nCells', ncells_global, 'latitude', 'degrees_north', & - 1, ncells_solve, latcell * rad_to_deg, map=global_grid_index) + 1, ncells_solve, real(latcell, kind_r8) * rad_to_deg, map=global_grid_index) lon_coord => horiz_coord_create('lonCell', 'nCells', ncells_global, 'longitude', 'degrees_east', & - 1, ncells_solve, loncell * rad_to_deg, map=global_grid_index) + 1, ncells_solve, real(loncell, kind_r8) * rad_to_deg, map=global_grid_index) allocate(cell_area(ncells_solve), stat=ierr) call check_allocate(ierr, subname, 'cell_area(ncells_solve)', 'dyn_grid', __LINE__) @@ -402,8 +405,8 @@ subroutine define_cam_grid() call check_allocate(ierr, subname, 'global_grid_map(3, ncells_solve)', 'dyn_grid', __LINE__) do i = 1, ncells_solve - cell_area(i) = areacell(i) - cell_weight(i) = areacell(i) / (4.0_kind_r8 * constant_pi * sphere_radius ** 2) + cell_area(i) = real(areacell(i), kind_r8) + cell_weight(i) = real(areacell(i), kind_r8) / (4.0_kind_r8 * constant_pi * real(sphere_radius, kind_r8) ** 2) global_grid_map(1, i) = int(i, kind_imap) global_grid_map(2, i) = int(1, kind_imap) @@ -426,9 +429,9 @@ subroutine define_cam_grid() ! Standard CAM-SIMA coordinate and dimension names are used. lat_coord => horiz_coord_create('lat', 'ncol', ncells_global, 'latitude', 'degrees_north', & - 1, ncells_solve, latcell * rad_to_deg, map=global_grid_index) + 1, ncells_solve, real(latcell, kind_r8) * rad_to_deg, map=global_grid_index) lon_coord => horiz_coord_create('lon', 'ncol', ncells_global, 'longitude', 'degrees_east', & - 1, ncells_solve, loncell * rad_to_deg, map=global_grid_index) + 1, ncells_solve, real(loncell, kind_r8) * rad_to_deg, map=global_grid_index) call dyn_debug_print('Registering grid "cam_cell" with id ' // stringify([dyn_grid_id('cam_cell')])) @@ -456,9 +459,9 @@ subroutine define_cam_grid() global_grid_index(:) = int(indextoedgeid(1:nedges_solve), kind_imap) lat_coord => horiz_coord_create('latEdge', 'nEdges', nedges_global, 'latitude', 'degrees_north', & - 1, nedges_solve, latedge * rad_to_deg, map=global_grid_index) + 1, nedges_solve, real(latedge, kind_r8) * rad_to_deg, map=global_grid_index) lon_coord => horiz_coord_create('lonEdge', 'nEdges', nedges_global, 'longitude', 'degrees_east', & - 1, nedges_solve, lonedge * rad_to_deg, map=global_grid_index) + 1, nedges_solve, real(lonedge, kind_r8) * rad_to_deg, map=global_grid_index) allocate(global_grid_map(3, nedges_solve), stat=ierr) call check_allocate(ierr, subname, 'global_grid_map(3, nedges_solve)', 'dyn_grid', __LINE__) @@ -494,9 +497,9 @@ subroutine define_cam_grid() global_grid_index(:) = int(indextovertexid(1:nvertices_solve), kind_imap) lat_coord => horiz_coord_create('latVertex', 'nVertices', nvertices_global, 'latitude', 'degrees_north', & - 1, nvertices_solve, latvertex * rad_to_deg, map=global_grid_index) + 1, nvertices_solve, real(latvertex, kind_r8) * rad_to_deg, map=global_grid_index) lon_coord => horiz_coord_create('lonVertex', 'nVertices', nvertices_global, 'longitude', 'degrees_east', & - 1, nvertices_solve, lonvertex * rad_to_deg, map=global_grid_index) + 1, nvertices_solve, real(lonvertex, kind_r8) * rad_to_deg, map=global_grid_index) allocate(global_grid_map(3, nvertices_solve), stat=ierr) call check_allocate(ierr, subname, 'global_grid_map(3, nvertices_solve)', 'dyn_grid', __LINE__) From 4491172914c8a3fb566db3ba1b5b12d1fa65f632 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Sat, 9 Nov 2024 04:05:39 -0700 Subject: [PATCH 10/24] Nullify pointers that are no longer needed --- src/dynamics/mpas/dyn_grid.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index e6780933..8f785f20 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -74,6 +74,8 @@ subroutine model_grid_init() ! Read time-invariant (e.g., grid/mesh) variables. call mpas_dynamical_core % read_write_stream(pio_file, 'r', 'invariant') + nullify(pio_file) + ! Compute local east, north and edge-normal unit vectors whenever time-invariant (e.g., grid/mesh) variables are read. call mpas_dynamical_core % compute_unit_vector() @@ -422,6 +424,7 @@ subroutine define_cam_grid() call cam_grid_attribute_register('mpas_cell', 'cell_weight', 'MPAS cell weight', 'nCells', cell_weight, & map=global_grid_index) + nullify(areacell) nullify(cell_area, cell_weight) nullify(lat_coord, lon_coord) @@ -438,7 +441,6 @@ subroutine define_cam_grid() call cam_grid_register('cam_cell', dyn_grid_id('cam_cell'), lat_coord, lon_coord, global_grid_map, & unstruct=.true., block_indexed=.false.) - nullify(areacell) nullify(indextocellid) nullify(latcell, loncell) From b32fe781a388e82a497cb235e9a9e333d57c0634 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 22 Nov 2024 13:53:30 -0700 Subject: [PATCH 11/24] Add `protected` attribute to grid/mesh variables Protect these variables from being modified externally. --- src/dynamics/mpas/dyn_comp.F90 | 33 ++++++++++++++++++++++++++++++--- src/dynamics/mpas/dyn_grid.F90 | 24 +++--------------------- 2 files changed, 33 insertions(+), 24 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 4e28bc8a..b8759a6c 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -18,6 +18,7 @@ module dyn_comp public :: dyn_debug_print public :: dyn_exchange_constituent_state + public :: dyn_inquire_mesh_dimensions public :: reverse public :: mpas_dynamical_core public :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels @@ -42,9 +43,9 @@ module dyn_comp type(mpas_dynamical_core_type) :: mpas_dynamical_core ! Local and global mesh dimensions of MPAS dynamical core. - integer :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels - integer :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max - real(kind_dyn_mpas) :: sphere_radius + integer, protected :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels + integer, protected :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max + real(kind_dyn_mpas), protected :: sphere_radius contains !> Print a debug message with optionally the value(s) of a variable. !> If `printer` is not supplied, the MPI root rank will print. Otherwise, the designated MPI rank will print instead. @@ -983,6 +984,32 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) end if end subroutine dyn_exchange_constituent_state + !> Inquire local and global mesh dimensions. Save them as protected module variables. + !> (KCW, 2024-11-21) + subroutine dyn_inquire_mesh_dimensions() + ! Module(s) from CAM-SIMA. + use string_utils, only: stringify + + character(*), parameter :: subname = 'dyn_comp::dyn_inquire_mesh_dimensions' + + call dyn_debug_print('Inquiring local and global mesh dimensions') + + call mpas_dynamical_core % get_local_mesh_dimension( & + ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels) + + call mpas_dynamical_core % get_global_mesh_dimension( & + ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max, & + sphere_radius) + + call dyn_debug_print('ncells_global = ' // stringify([ncells_global])) + call dyn_debug_print('nedges_global = ' // stringify([nedges_global])) + call dyn_debug_print('nvertices_global = ' // stringify([nvertices_global])) + call dyn_debug_print('nvertlevels = ' // stringify([nvertlevels])) + call dyn_debug_print('ncells_max = ' // stringify([ncells_max])) + call dyn_debug_print('nedges_max = ' // stringify([nedges_max])) + call dyn_debug_print('sphere_radius = ' // stringify([sphere_radius])) + end subroutine dyn_inquire_mesh_dimensions + !> Mark everything in the `physics_types` module along with constituents as initialized !> to prevent physics from attempting to read them from a file. !> (KCW, 2024-05-23) diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 8f785f20..2b4e3468 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -37,10 +37,7 @@ subroutine model_grid_init() use cam_abortutils, only: endrun use cam_constituents, only: num_advected use cam_initfiles, only: initial_file_get_id - use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, & - ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels, & - ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max, & - sphere_radius + use dyn_comp, only: dyn_debug_print, dyn_inquire_mesh_dimensions, mpas_dynamical_core, nvertlevels use dynconst, only: dynconst_init use string_utils, only: stringify use vert_coord, only: pver, vert_coord_init @@ -79,23 +76,8 @@ subroutine model_grid_init() ! Compute local east, north and edge-normal unit vectors whenever time-invariant (e.g., grid/mesh) variables are read. call mpas_dynamical_core % compute_unit_vector() - ! Inquire local and global mesh dimensions and save them as module variables. - call dyn_debug_print('Inquiring local and global mesh dimensions') - - call mpas_dynamical_core % get_local_mesh_dimension( & - ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels) - - call mpas_dynamical_core % get_global_mesh_dimension( & - ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max, & - sphere_radius) - - call dyn_debug_print('ncells_global = ' // stringify([ncells_global])) - call dyn_debug_print('nedges_global = ' // stringify([nedges_global])) - call dyn_debug_print('nvertices_global = ' // stringify([nvertices_global])) - call dyn_debug_print('nvertlevels = ' // stringify([nvertlevels])) - call dyn_debug_print('ncells_max = ' // stringify([ncells_max])) - call dyn_debug_print('nedges_max = ' // stringify([nedges_max])) - call dyn_debug_print('sphere_radius = ' // stringify([sphere_radius])) + ! Inquire local and global mesh dimensions. + call dyn_inquire_mesh_dimensions() ! Check for consistency in numbers of vertical layers. if (nvertlevels /= pver) then From bc45c6723a9a09d3d045e83477fddad3689aa6fe Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 2 Dec 2024 12:01:59 -0700 Subject: [PATCH 12/24] Change some wording to plural --- src/dynamics/mpas/dyn_comp.F90 | 22 +++++++++++----------- src/dynamics/mpas/dyn_coupling.F90 | 12 ++++++------ 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index b8759a6c..db651f1b 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -17,7 +17,7 @@ module dyn_comp public :: dyn_final public :: dyn_debug_print - public :: dyn_exchange_constituent_state + public :: dyn_exchange_constituent_states public :: dyn_inquire_mesh_dimensions public :: reverse public :: mpas_dynamical_core @@ -266,9 +266,9 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) ! Perform default initialization for all constituents. ! Subsequently, they can be overridden depending on the namelist option (below) and ! the actual availability (checked and handled by MPAS). - call dyn_debug_print('Calling dyn_exchange_constituent_state') + call dyn_debug_print('Calling dyn_exchange_constituent_states') - call dyn_exchange_constituent_state(direction='e', exchange=.true., conversion=.false.) + call dyn_exchange_constituent_states(direction='e', exchange=.true., conversion=.false.) ! Namelist option that controls if constituents are to be read from the file. if (readtrace) then @@ -287,7 +287,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) end if call clean_iodesc_list() - call mark_variable_as_initialized() + call mark_variables_as_initialized() nullify(pio_init_file) nullify(pio_topo_file) @@ -830,13 +830,13 @@ end subroutine set_analytic_initial_condition !> If `exchange` is `.true.` and `direction` is "i" or "import", set physics state `constituents` from MPAS state `scalars`. !> Think of it as "exporting/importing constituent states in CAM-SIMA to/from MPAS". !> Otherwise, if `exchange` is `.false.`, no exchange is performed at all. - !> If `conversion` is `.true.`, appropriate conversion is performed for constituent mixing ratio that has different + !> If `conversion` is `.true.`, appropriate conversion is performed for constituent mixing ratios that have different !> definitions between CAM-SIMA and MPAS (i.e., dry/moist). !> Otherwise, if `conversion` is `.false.`, no conversion is performed at all. !> This subroutine is intentionally designed to have these elaborate controls due to complications in CAM-SIMA. !> Some procedures in CAM-SIMA expect constituent states to be dry, while the others expect them to be moist. !> (KCW, 2024-09-26) - subroutine dyn_exchange_constituent_state(direction, exchange, conversion) + subroutine dyn_exchange_constituent_states(direction, exchange, conversion) ! Module(s) from CAM-SIMA. use cam_abortutils, only: check_allocate, endrun use cam_constituents, only: const_is_dry, const_is_water_species, num_advected @@ -850,7 +850,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) logical, intent(in) :: exchange logical, intent(in) :: conversion - character(*), parameter :: subname = 'dyn_comp::dyn_exchange_constituent_state' + character(*), parameter :: subname = 'dyn_comp::dyn_exchange_constituent_states' integer :: i, j integer :: ierr integer, allocatable :: is_water_species_index(:) @@ -982,7 +982,7 @@ subroutine dyn_exchange_constituent_state(direction, exchange, conversion) ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. call mpas_dynamical_core % exchange_halo('scalars') end if - end subroutine dyn_exchange_constituent_state + end subroutine dyn_exchange_constituent_states !> Inquire local and global mesh dimensions. Save them as protected module variables. !> (KCW, 2024-11-21) @@ -1013,13 +1013,13 @@ end subroutine dyn_inquire_mesh_dimensions !> Mark everything in the `physics_types` module along with constituents as initialized !> to prevent physics from attempting to read them from a file. !> (KCW, 2024-05-23) - subroutine mark_variable_as_initialized() + subroutine mark_variables_as_initialized() ! Module(s) from CAM-SIMA. use cam_constituents, only: const_name, num_advected ! Module(s) from CCPP. use phys_vars_init_check, only: mark_as_initialized - character(*), parameter :: subname = 'dyn_comp::mark_variable_as_initialized' + character(*), parameter :: subname = 'dyn_comp::mark_variables_as_initialized' integer :: i ! The variables below are managed by dynamics interface. @@ -1074,7 +1074,7 @@ subroutine mark_variable_as_initialized() call mark_as_initialized('vertically_integrated_total_energy_using_physics_energy_formula_at_start_of_physics_timestep') call mark_as_initialized('vertically_integrated_total_water') call mark_as_initialized('vertically_integrated_total_water_at_start_of_physics_timestep') - end subroutine mark_variable_as_initialized + end subroutine mark_variables_as_initialized !> Run MPAS dynamical core to integrate the dynamical states with time. !> (KCW, 2024-07-11) diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 322e6a01..40e37fb3 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -17,7 +17,7 @@ module dyn_coupling !> (KCW, 2024-07-31) subroutine dynamics_to_physics_coupling() ! Module(s) from CAM-SIMA. - use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_state, ncells_solve + use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_states, ncells_solve character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling' integer :: column_index @@ -55,7 +55,7 @@ subroutine dynamics_to_physics_coupling() call init_shared_variables() - call dyn_exchange_constituent_state(direction='i', exchange=.true., conversion=.false.) + call dyn_exchange_constituent_states(direction='i', exchange=.true., conversion=.false.) call dyn_debug_print('Setting physics state variables column by column') @@ -362,7 +362,7 @@ subroutine set_physics_state_external() end do ! Note that constituents become moist after this. - call dyn_exchange_constituent_state(direction='i', exchange=.false., conversion=.true.) + call dyn_exchange_constituent_states(direction='i', exchange=.false., conversion=.true.) ! Impose minimum limits on constituents. call qneg_run(subname, ncells_solve, pver, minimum_constituents, constituents, ierr, cerr) @@ -410,7 +410,7 @@ end subroutine dynamics_to_physics_coupling !> (KCW, 2024-09-20) subroutine physics_to_dynamics_coupling() ! Module(s) from CAM-SIMA. - use dyn_comp, only: dyn_exchange_constituent_state + use dyn_comp, only: dyn_exchange_constituent_states character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling' integer, pointer :: index_qv @@ -422,7 +422,7 @@ subroutine physics_to_dynamics_coupling() call init_shared_variables() - call dyn_exchange_constituent_state(direction='e', exchange=.true., conversion=.true.) + call dyn_exchange_constituent_states(direction='e', exchange=.true., conversion=.true.) call set_mpas_physics_tendency_ru() call set_mpas_physics_tendency_rho() @@ -459,7 +459,7 @@ subroutine init_shared_variables() 'dyn_coupling', __LINE__) ! Save water vapor mixing ratio before being updated by physics because `set_mpas_physics_tendency_rtheta` - ! needs it. This must be done before calling `dyn_exchange_constituent_state`. + ! needs it. This must be done before calling `dyn_exchange_constituent_states`. qv_prev(:, :) = real(scalars(index_qv, :, 1:ncells_solve), kind_r8) end subroutine init_shared_variables From be1449f569f6e6d61fc60ca43cd6128693e6b8e7 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 22 Nov 2024 14:10:35 -0700 Subject: [PATCH 13/24] Include argument keyword for better clarity --- src/dynamics/mpas/dyn_comp.F90 | 2 +- src/dynamics/mpas/dyn_coupling.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index db651f1b..453c6bf3 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -517,7 +517,7 @@ subroutine set_mpas_state_u() nullify(ucellzonal, ucellmeridional) - call mpas_dynamical_core % compute_edge_wind(.false.) + call mpas_dynamical_core % compute_edge_wind(wind_tendency=.false.) end subroutine set_mpas_state_u !> Set MPAS state `w` (i.e., vertical velocity at cell interfaces). diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 40e37fb3..831bcd02 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -504,7 +504,7 @@ subroutine set_mpas_physics_tendency_ru() nullify(u_tendency, v_tendency) - call mpas_dynamical_core % compute_edge_wind(.true.) + call mpas_dynamical_core % compute_edge_wind(wind_tendency=.true.) end subroutine set_mpas_physics_tendency_ru !> Set MPAS physics tendency `tend_rho_physics` (i.e., "coupled" tendency of dry air density due to physics). From efd27e8c75c776bea67269fb22a9798f9e3f29fd Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Tue, 5 Mar 2024 16:01:12 -0700 Subject: [PATCH 14/24] Reorganize directory structure for better clarity --- cime_config/buildlib | 6 +++--- src/dynamics/mpas/{ => assets}/Makefile | 0 src/dynamics/mpas/{ => assets}/Makefile.in.CESM | 0 3 files changed, 3 insertions(+), 3 deletions(-) rename src/dynamics/mpas/{ => assets}/Makefile (100%) rename src/dynamics/mpas/{ => assets}/Makefile.in.CESM (100%) diff --git a/cime_config/buildlib b/cime_config/buildlib index 214e6182..33f1d3ec 100755 --- a/cime_config/buildlib +++ b/cime_config/buildlib @@ -189,9 +189,9 @@ def _setup_mpas(case: Case) -> None: shutil.copytree(mpas_dycore_src_root, mpas_dycore_bld_root, copy_function=_copy2_as_needed, dirs_exist_ok=True) shutil.move(os.path.join(mpas_dycore_bld_root, "Makefile"), os.path.join(mpas_dycore_bld_root, "Makefile.CESM")) - _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "driver", "dyn_mpas_subdriver.F90")), os.path.join(mpas_dycore_bld_root, "driver")) - _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "Makefile")), mpas_dycore_bld_root) - _copy2_as_needed(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "Makefile.in.CESM")), mpas_dycore_bld_root) + + shutil.copytree(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "assets")), mpas_dycore_bld_root, copy_function=_copy2_as_needed, dirs_exist_ok=True) + shutil.copytree(os.path.normpath(os.path.join(mpas_dycore_src_root, os.pardir, os.pardir, "driver")), os.path.join(mpas_dycore_bld_root, "driver"), copy_function=_copy2_as_needed, dirs_exist_ok=True) def _copy2_as_needed(src: str, dst: str) -> None: """ diff --git a/src/dynamics/mpas/Makefile b/src/dynamics/mpas/assets/Makefile similarity index 100% rename from src/dynamics/mpas/Makefile rename to src/dynamics/mpas/assets/Makefile diff --git a/src/dynamics/mpas/Makefile.in.CESM b/src/dynamics/mpas/assets/Makefile.in.CESM similarity index 100% rename from src/dynamics/mpas/Makefile.in.CESM rename to src/dynamics/mpas/assets/Makefile.in.CESM From 37c86517b358c241e6e4ad0923442faa4c244183 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 4 Mar 2024 10:52:43 -0700 Subject: [PATCH 15/24] Remove hacky `sed` commands from `Makefile` Manipulating source code by running `sed` commands during building is neither elegant nor obvious. Downstream modifications like these are better achieved and communicated through patches. This is also the standard practice in open source software. --- src/dynamics/mpas/assets/Makefile.in.CESM | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/src/dynamics/mpas/assets/Makefile.in.CESM b/src/dynamics/mpas/assets/Makefile.in.CESM index b116c79e..ea33ba87 100644 --- a/src/dynamics/mpas/assets/Makefile.in.CESM +++ b/src/dynamics/mpas/assets/Makefile.in.CESM @@ -75,7 +75,7 @@ all: @echo ' `make libmpas-clean ESM="CESM" LIBROOT="..." SRCROOT="..."`' .PHONY: libmpas-prepare -libmpas-prepare: libmpas-archiver-script.txt libmpas-no-physics libmpas-prefix-namelist-groups libmpas-preview +libmpas-prepare: libmpas-archiver-script.txt libmpas-preview # Combine multiple static libraries into `libmpas.a` via archiver/MRI script. This requires GNU or GNU-like archiver (`ar`) program. libmpas-archiver-script.txt: @@ -87,17 +87,6 @@ libmpas-archiver-script.txt: @echo "save" >> $(@) @echo "end" >> $(@) -# Do not use built-in MPAS/WRF physics. -.PHONY: libmpas-no-physics -libmpas-no-physics: - @sed -E -i -e "s/^ *PHYSICS=.+$$/PHYSICS=/g" core_atmosphere/Makefile - -# Prefix `mpas_` to MPAS namelist groups to avoid naming collisions. -.PHONY: libmpas-prefix-namelist-groups -libmpas-prefix-namelist-groups: - @sed -E -i -e "s/(^ *< *nml_record.+name=)\"(mpas_)?(\w+)\"/\1\"mpas_\3\"/g" core_atmosphere/Registry.xml - @sed -E -i -e "s/(^ *< *nml_record.+name=)\"(mpas_)?(\w+)\"/\1\"mpas_\3\"/g" core_atmosphere/diagnostics/Registry_*.xml - .PHONY: libmpas-preview libmpas-preview: @echo "Previewing build options for $(LIBROOT)/libmpas.a:" From 40b1ef90a2fbaa54923f5b1e3bd776a6b749c9e8 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 12 Jul 2024 16:33:00 -0600 Subject: [PATCH 16/24] Convert CAM-SIMA-specific modifications to MPAS into patches These patches are applied before building. --- ...MPAS-namelist-group-and-option-names.patch | 190 ++++++++++++++++++ ...e-physics-for-MPAS-dycore-only-build.patch | 42 ++++ ...ants-at-the-native-precision-of-MPAS.patch | 94 +++++++++ 3 files changed, 326 insertions(+) create mode 100644 src/dynamics/mpas/assets/0001-Prefix-all-MPAS-namelist-group-and-option-names.patch create mode 100644 src/dynamics/mpas/assets/0002-Disable-physics-for-MPAS-dycore-only-build.patch create mode 100644 src/dynamics/mpas/assets/0003-Declare-constants-at-the-native-precision-of-MPAS.patch diff --git a/src/dynamics/mpas/assets/0001-Prefix-all-MPAS-namelist-group-and-option-names.patch b/src/dynamics/mpas/assets/0001-Prefix-all-MPAS-namelist-group-and-option-names.patch new file mode 100644 index 00000000..c0d4e1e5 --- /dev/null +++ b/src/dynamics/mpas/assets/0001-Prefix-all-MPAS-namelist-group-and-option-names.patch @@ -0,0 +1,190 @@ +From 0000000000000000000000000000000000000000 Mon Sep 17 00:00:00 2001 +From: Kuan-Chih Wang +Date: Mon, 1 Apr 2024 16:46:50 -0600 +Subject: [PATCH] Prefix all MPAS namelist group and option names + +The following transformations are performed for each MPAS namelist group and +option name: + +1. Leading `config_` is removed recursively from the name. Case-insensitive. +2. Leading `mpas_` is removed recursively from the name. Case-insensitive. +3. Prepend `mpas_` to the name. + +As a result, it is easier to distinguish MPAS namelist groups and options from +CAM-SIMA ones. Note that only namelist I/O is affected. Internally, MPAS still +refers to its namelist options by their original names due to compatibility reasons. + +This downstream patch is maintained by CAM-SIMA for its particular use case of MPAS. +--- + src/tools/registry/gen_inc.c | 80 +++++++++++++++++++++++++++++++++++- + src/tools/registry/gen_inc.h | 4 ++ + 2 files changed, 82 insertions(+), 2 deletions(-) + +diff --git a/src/tools/registry/gen_inc.c b/src/tools/registry/gen_inc.c +index 94f5f714..95b2502c 100644 +--- a/src/tools/registry/gen_inc.c ++++ b/src/tools/registry/gen_inc.c +@@ -15,6 +15,10 @@ + #include "fortprintf.h" + #include "utility.h" + ++#ifdef MPAS_CAM_DYCORE ++#include ++#endif ++ + void process_core_macro(const char *macro, const char *val, va_list ap); + void process_domain_macro(const char *macro, const char *val, va_list ap); + +@@ -696,8 +700,18 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + ezxml_t nmlrecs_xml, nmlopt_xml; + + const char *const_core; +- const char *nmlrecname, *nmlrecindef, *nmlrecinsub; +- const char *nmloptname, *nmlopttype, *nmloptval, *nmloptunits, *nmloptdesc, *nmloptposvals, *nmloptindef; ++#ifdef MPAS_CAM_DYCORE ++ const char *orinmlrecname; ++ const char *orinmloptname; ++ // Fortran variable names have a length limit of 63 characters. +1 for the terminating null character. ++ char nmlrecname[64]; ++ char nmloptname[64]; ++#else ++ const char *nmlrecname; ++ const char *nmloptname; ++#endif ++ const char *nmlrecindef, *nmlrecinsub; ++ const char *nmlopttype, *nmloptval, *nmloptunits, *nmloptdesc, *nmloptposvals, *nmloptindef; + + char pool_name[1024]; + char core_string[1024]; +@@ -743,7 +757,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + + // Parse Namelist Records + for (nmlrecs_xml = ezxml_child(registry, "nml_record"); nmlrecs_xml; nmlrecs_xml = nmlrecs_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmlrecname = ezxml_attr(nmlrecs_xml, "name"); ++ transform_name(nmlrecname, sizeof(nmlrecname), orinmlrecname); ++#else + nmlrecname = ezxml_attr(nmlrecs_xml, "name"); ++#endif + nmlrecindef = ezxml_attr(nmlrecs_xml, "in_defaults"); + nmlrecinsub = ezxml_attr(nmlrecs_xml, "in_subpool"); + +@@ -777,7 +796,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + + // Define variable definitions prior to reading the namelist in. + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); ++#endif + nmlopttype = ezxml_attr(nmlopt_xml, "type"); + nmloptval = ezxml_attr(nmlopt_xml, "default_value"); + nmloptunits = ezxml_attr(nmlopt_xml, "units"); +@@ -809,7 +833,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + // Define the namelist block, to read the namelist record in. + fortprintf(fd, " namelist /%s/ &\n", nmlrecname); + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); ++#endif + if(nmlopt_xml->next){ + fortprintf(fd, " %s, &\n", nmloptname); + } else { +@@ -840,7 +869,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + // Define broadcast calls for namelist values. + fortprintf(fd, " if (ierr <= 0) then\n"); + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); ++#endif + nmlopttype = ezxml_attr(nmlopt_xml, "type"); + + if(strncmp(nmlopttype, "real", 1024) == 0){ +@@ -858,7 +892,12 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + fortprintf(fd, " call mpas_log_write(' The following values will be used for variables in this record:')\n"); + fortprintf(fd, " call mpas_log_write(' ')\n"); + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); ++#endif + nmlopttype = ezxml_attr(nmlopt_xml, "type"); + + if (strncmp(nmlopttype, "character", 1024) == 0) { +@@ -885,10 +924,18 @@ int parse_namelist_records_from_registry(ezxml_t registry)/*{{{*/ + fortprintf(fd, "\n"); + + for (nmlopt_xml = ezxml_child(nmlrecs_xml, "nml_option"); nmlopt_xml; nmlopt_xml = nmlopt_xml->next){ ++#ifdef MPAS_CAM_DYCORE ++ orinmloptname = ezxml_attr(nmlopt_xml, "name"); ++ transform_name(nmloptname, sizeof(nmloptname), orinmloptname); ++ ++ fortprintf(fd, " call mpas_pool_add_config(%s, '%s', %s)\n", pool_name, orinmloptname, nmloptname); ++ fortprintf(fcg, " call mpas_pool_get_config(configPool, '%s', %s)\n", orinmloptname, nmloptname); ++#else + nmloptname = ezxml_attr(nmlopt_xml, "name"); + + fortprintf(fd, " call mpas_pool_add_config(%s, '%s', %s)\n", pool_name, nmloptname, nmloptname); + fortprintf(fcg, " call mpas_pool_get_config(configPool, '%s', %s)\n", nmloptname, nmloptname); ++#endif + } + fortprintf(fd, "\n"); + fortprintf(fcg, "\n"); +@@ -2532,3 +2579,32 @@ int parse_structs_from_registry(ezxml_t registry)/*{{{*/ + + return 0; + }/*}}}*/ ++ ++ ++#ifdef MPAS_CAM_DYCORE ++// Perform transformations for namelist group and option names. ++void transform_name(char *new_name, const size_t new_name_size, const char *old_name) { ++ const char *const new_prefix = "mpas_"; ++ const char *const old_prefix = "config_"; ++ size_t size = 0; ++ ++ if (!new_name || !old_name || new_name_size == 0) return; ++ ++ // Remove all leading whitespaces by moving pointer forward. ++ while (*old_name != '\0' && isspace((unsigned char) *old_name)) old_name++; ++ ++ // Remove all leading `config_` by moving pointer forward. ++ while (strncasecmp(old_name, old_prefix, strlen(old_prefix)) == 0) old_name += strlen(old_prefix); ++ ++ // Remove all leading `mpas_` by moving pointer forward. ++ while (strncasecmp(old_name, new_prefix, strlen(new_prefix)) == 0) old_name += strlen(new_prefix); ++ ++ *new_name = '\0'; ++ size = snprintf(NULL, 0, "%s%s", new_prefix, old_name) + 1; ++ snprintf(new_name, size > new_name_size ? new_name_size : size, "%s%s", new_prefix, old_name); ++ ++ // Remove all trailing whitespaces by zeroing (nulling) out. ++ new_name += strlen(new_name) - 1; ++ while (*new_name != '\0' && isspace((unsigned char) *new_name)) *new_name-- = '\0'; ++} ++#endif +diff --git a/src/tools/registry/gen_inc.h b/src/tools/registry/gen_inc.h +index fc94e78b..69a76ace 100644 +--- a/src/tools/registry/gen_inc.h ++++ b/src/tools/registry/gen_inc.h +@@ -38,3 +38,7 @@ int push_attributes(ezxml_t currentPosition); + int merge_structs_and_var_arrays(ezxml_t currentPosition); + int merge_streams(ezxml_t registry); + int parse_structs_from_registry(ezxml_t registry); ++ ++#ifdef MPAS_CAM_DYCORE ++void transform_name(char *new_name, const size_t new_name_size, const char *old_name); ++#endif +-- +2.43.0 + diff --git a/src/dynamics/mpas/assets/0002-Disable-physics-for-MPAS-dycore-only-build.patch b/src/dynamics/mpas/assets/0002-Disable-physics-for-MPAS-dycore-only-build.patch new file mode 100644 index 00000000..4f664d88 --- /dev/null +++ b/src/dynamics/mpas/assets/0002-Disable-physics-for-MPAS-dycore-only-build.patch @@ -0,0 +1,42 @@ +From 0000000000000000000000000000000000000000 Mon Sep 17 00:00:00 2001 +From: Kuan-Chih Wang +Date: Thu, 1 Aug 2024 17:09:58 -0600 +Subject: [PATCH] Disable physics for MPAS dycore-only build + +When building MPAS as a dycore, all physics-related components are disabled. +This build configuration is usually used by CAM/CAM-SIMA, and can be achieved by +defining the `MPAS_CAM_DYCORE` macro in `CPPFLAGS`. + +The `PHYSICS` variable controls whether physics are enabled in MPAS, but its logic +is decoupled from the `MPAS_CAM_DYCORE` macro. Disabling physics in MPAS currently +requires manual interventions. + +Therefore, automatically disable physics when the `MPAS_CAM_DYCORE` macro is found +in `CPPFLAGS`. + +This downstream patch is maintained by CAM-SIMA for its particular use case of MPAS. +--- + src/core_atmosphere/Makefile | 7 +++++-- + 1 file changed, 5 insertions(+), 2 deletions(-) + +diff --git a/src/core_atmosphere/Makefile b/src/core_atmosphere/Makefile +index 8d9f4f1a..cac8255e 100644 +--- a/src/core_atmosphere/Makefile ++++ b/src/core_atmosphere/Makefile +@@ -4,8 +4,11 @@ + # To build a dycore-only MPAS-Atmosphere model, comment-out or delete + # the definition of PHYSICS, below + # +-PHYSICS=-DDO_PHYSICS +- ++# If MPAS_CAM_DYCORE is found in CPPFLAGS, PHYSICS will become undefined automatically ++# ++ifeq ($(findstring MPAS_CAM_DYCORE,$(CPPFLAGS)),) ++ PHYSICS = -DDO_PHYSICS ++endif + + ifdef PHYSICS + PHYSCORE = physcore +-- +2.43.0 + diff --git a/src/dynamics/mpas/assets/0003-Declare-constants-at-the-native-precision-of-MPAS.patch b/src/dynamics/mpas/assets/0003-Declare-constants-at-the-native-precision-of-MPAS.patch new file mode 100644 index 00000000..d88d02ac --- /dev/null +++ b/src/dynamics/mpas/assets/0003-Declare-constants-at-the-native-precision-of-MPAS.patch @@ -0,0 +1,94 @@ +From 0000000000000000000000000000000000000000 Mon Sep 17 00:00:00 2001 +From: Kuan-Chih Wang +Date: Wed, 13 Nov 2024 17:19:10 -0700 +Subject: [PATCH] Declare constants at the native precision of MPAS + +When building MPAS as a dycore, certain constants in the `mpas_constants` module are +imported from the `physconst` module, which is a part of CAM/CAM-SIMA. However, +multiple issues arise if the precision of those constants differs from MPAS. + +For example, building MPAS in single precision mode with CAM-SIMA fails due to +multiple occurrences of type mismatch between actual and dummy arguments. + +mpas_geometry_utils.F:885:157: + + 885 | call mpas_log_write('$r', MPAS_LOG_ERR, realArgs=(/mpas_triangle_signed_area_sphere(a,b,c,sphereRadius) - pii/2.0_RKIND*sphereRadius*sphereRadius/)) + | 1 +Error: Type mismatch in argument 'realargs' at (1); passed REAL(8) to REAL(4) + +Here, `pii` is declared by CAM-SIMA to be double precision, and it causes unintended +floating-point promotion in the expression. + +The solution is to ensure that constants in the `mpas_constants` module are declared +at the native precision of MPAS. + +This downstream patch is maintained by CAM-SIMA for its particular use case of MPAS. +--- + src/framework/mpas_constants.F | 44 ++++++++++++++++++++++++++-------- + 1 file changed, 34 insertions(+), 10 deletions(-) + +diff --git a/src/framework/mpas_constants.F b/src/framework/mpas_constants.F +index c98cb810..68164bab 100644 +--- a/src/framework/mpas_constants.F ++++ b/src/framework/mpas_constants.F +@@ -23,16 +23,30 @@ module mpas_constants + use mpas_kind_types + + #ifdef MPAS_CAM_DYCORE +- use physconst, only : pii => pi +- use physconst, only : gravity => gravit +- use physconst, only : omega +- use physconst, only : a => rearth +- use physconst, only : cp => cpair +- use physconst, only : rgas => rair +- use physconst, only : rv => rh2o +- real (kind=RKIND) :: rvord = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived +- real (kind=RKIND) :: cv = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived +- real (kind=RKIND) :: cvpm = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ use physconst, only: external_pii => pi ++ use physconst, only: external_a => rearth ++ use physconst, only: external_omega => omega ++ use physconst, only: external_gravity => gravit ++ use physconst, only: external_rgas => rair ++ use physconst, only: external_rv => rh2o ++ use physconst, only: external_cp => cpair ++ private :: external_pii ++ private :: external_a ++ private :: external_omega ++ private :: external_gravity ++ private :: external_rgas ++ private :: external_rv ++ private :: external_cp ++ real (kind=RKIND), protected :: pii = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: a = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: omega = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: gravity = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: rgas = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: rv = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: cp = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: rvord = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: cv = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived ++ real (kind=RKIND), protected :: cvpm = huge(1.0_RKIND) ! Derived in mpas_constants_compute_derived + #else + real (kind=RKIND), parameter :: pii = 3.141592653589793_RKIND !< Constant: Pi + real (kind=RKIND), parameter :: a = 6371229.0_RKIND !< Constant: Spherical Earth radius [m] +@@ -77,6 +91,16 @@ module mpas_constants + implicit none + + #ifdef MPAS_CAM_DYCORE ++ ! Convert external constants to the native precision of MPAS (i.e., `RKIND`). ++ ++ pii = real(external_pii, RKIND) ++ a = real(external_a, RKIND) ++ omega = real(external_omega, RKIND) ++ gravity = real(external_gravity, RKIND) ++ rgas = real(external_rgas, RKIND) ++ rv = real(external_rv, RKIND) ++ cp = real(external_cp, RKIND) ++ + ! + ! In the case of CAM-MPAS, rgas may depend on a CAM namelist option, + ! so physical constants that depend on rgas must be computed here after +-- +2.43.0 + From cdc0b8e0661b846983eba14a49b35c1b5552de96 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 12 Jul 2024 16:33:49 -0600 Subject: [PATCH 17/24] Apply patches before building --- src/dynamics/mpas/assets/Makefile.in.CESM | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/dynamics/mpas/assets/Makefile.in.CESM b/src/dynamics/mpas/assets/Makefile.in.CESM index ea33ba87..4b1c5708 100644 --- a/src/dynamics/mpas/assets/Makefile.in.CESM +++ b/src/dynamics/mpas/assets/Makefile.in.CESM @@ -75,7 +75,17 @@ all: @echo ' `make libmpas-clean ESM="CESM" LIBROOT="..." SRCROOT="..."`' .PHONY: libmpas-prepare -libmpas-prepare: libmpas-archiver-script.txt libmpas-preview +libmpas-prepare: libmpas-apply-patch libmpas-archiver-script.txt libmpas-preview + +# Apply patches. +.PHONY: libmpas-apply-patch +libmpas-apply-patch: + @for file in *.patch; do \ + if git apply --check -p2 "$${file}"; then \ + echo "Applying $${file}"; \ + git apply -p2 "$${file}"; \ + fi; \ + done # Combine multiple static libraries into `libmpas.a` via archiver/MRI script. This requires GNU or GNU-like archiver (`ar`) program. libmpas-archiver-script.txt: @@ -88,7 +98,7 @@ libmpas-archiver-script.txt: @echo "end" >> $(@) .PHONY: libmpas-preview -libmpas-preview: +libmpas-preview: libmpas-apply-patch libmpas-archiver-script.txt @echo "Previewing build options for $(LIBROOT)/libmpas.a:" @echo "AR = $(AR)" @echo "ARFLAGS = $(ARFLAGS)" From a05b30f6f3c05b7469bf3659a09f9b71cbef8c8f Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 5 Apr 2024 11:38:30 -0600 Subject: [PATCH 18/24] Add Python script for generating namelist definition This Python script generates the XML namelist definition file of MPAS dynamical core from MPAS registry. Optionally, if an XML schema is provided, the generated file can also be validated for correctness. --- .../assets/generate_namelist_definition.py | 301 ++++++++++++++++++ 1 file changed, 301 insertions(+) create mode 100755 src/dynamics/mpas/assets/generate_namelist_definition.py diff --git a/src/dynamics/mpas/assets/generate_namelist_definition.py b/src/dynamics/mpas/assets/generate_namelist_definition.py new file mode 100755 index 00000000..f72aa69d --- /dev/null +++ b/src/dynamics/mpas/assets/generate_namelist_definition.py @@ -0,0 +1,301 @@ +#!/usr/bin/env python3 + +''' +Generate XML namelist definition file for MPAS dynamical core in CAM-SIMA. +''' + +import argparse +import textwrap +import xml.etree.ElementTree as ET + +# These namelist groups are irrelevant to CAM-SIMA for its particular use case of MPAS. +# Hide them to prevent users from inadvertently modifying them. +EXCLUDED_NAMELIST_GROUP = [ + 'assimilation', + 'iau', + 'limited_area', + 'physics', + 'restart' +] +# These namelist options are forcefully controlled by CAM-SIMA at run-time. +# Hide them to prevent users from inadvertently modifying them. +EXCLUDED_NAMELIST_OPTION = [ + 'config_calendar_type', + 'config_do_restart', + 'config_run_duration', + 'config_start_time', + 'config_stop_time' +] +# List all overridden namelist options below. +# `OVERRIDDEN_NAMELIST_OPTION` is a dictionary with `str` as keys and `list[str]` as values. +OVERRIDDEN_NAMELIST_OPTION = { + 'mpas_block_decomp_file_prefix': [ + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa30.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa15.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa12.graph.info.part.', + '${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa15-3.graph.info.part.' + ], + 'mpas_coef_3rd_order': [ + '1.0' + ], + 'mpas_dt': [ + '1800.0', + '900.0', + '450.0', + '225.0' + ], + 'mpas_len_disp': [ + '480000.0', + '120000.0', + '60000.0', + '30000.0' + ], + 'mpas_number_cam_damping_levels': [ + '0' + ], + 'mpas_number_rayleigh_damp_u_levels': [ + '5' + ], + 'mpas_print_detailed_minmax_vel': [ + '.true.' + ], + 'mpas_rayleigh_damp_u': [ + '.true.' + ] +} + +INDENT_PER_LEVEL = ' ' * 4 +NEW_PREFIX = 'mpas_' +OLD_PREFIX = 'config_' + +def parse_argument() -> argparse.Namespace: + ''' + Parse command line arguments. + ''' + + parser = argparse.ArgumentParser( + description='Generate XML namelist definition file for MPAS dynamical core in CAM-SIMA.' + ) + + parser.add_argument( + '-r', '--registry', + default='Registry.xml', + type=str, + required=False, + help='XML MPAS registry file.', + dest='reg_xml' + ) + parser.add_argument( + '-n', '--namelist', + default='Namelist.xml', + type=str, + required=False, + help='XML namelist definition file.', + dest='nml_xml' + ) + parser.add_argument( + '-s', '--schema', + default=None, + type=str, + required=False, + help='XML schema for namelist definition file.', + dest='nml_xsd' + ) + + argument = parser.parse_args() + + return argument + +def parse_xml(xml_file: str) -> ET.ElementTree: + ''' + Parse XML file into element tree. + ''' + + xml_et = ET.parse(xml_file) + + return xml_et + +def validate_xml(xml_file: str, xsd_file: str) -> bool: + ''' + Validate XML file against XSD file. + ''' + + # Only import `xmlschema` if XML validation is requested. Shush pylint about it. + import xmlschema # pylint: disable=import-outside-toplevel + + xml_schema = xmlschema.XMLSchema(xsd_file) + + return xml_schema.is_valid(xml_file) + +def transform_name(name: str) -> str: + ''' + Change prefix of namelist option/group name. + ''' + + while name.startswith(OLD_PREFIX): + name = name[len(OLD_PREFIX):] + + while name.startswith(NEW_PREFIX): + name = name[len(NEW_PREFIX):] + + name = NEW_PREFIX + name + + return name + +def translate_element_tree(reg_xml_et: ET.ElementTree) -> ET.ElementTree: + ''' + Translate MPAS registry into namelist definition. + ''' + + # `entry_id_pg` is the root element in namelist definition. + entry_id_pg_element = ET.Element('entry_id_pg', {'version': '0.1'}) + + comment_element = ET.Comment( + '\n' + + INDENT_PER_LEVEL * 2 + 'MPAS dycore' + '\n' + + '\n' + + INDENT_PER_LEVEL * 2 + 'Note to developers/maintainers:' + '\n' + + INDENT_PER_LEVEL * 2 + 'This file is auto-generated from MPAS registry. Do not edit directly.' + '\n' + + INDENT_PER_LEVEL * 2 + 'Instead, use the Python script at `src/dynamics/mpas/assets/generate_namelist_definition.py`.' + '\n' + + INDENT_PER_LEVEL + ) + entry_id_pg_element.append(comment_element) + + for namelist_group in reg_xml_et.findall('nml_record'): + if namelist_group.attrib['name'].strip().lower() in EXCLUDED_NAMELIST_GROUP: + continue + + for namelist_option in namelist_group.findall('nml_option'): + if namelist_option.attrib['name'].strip().lower() in EXCLUDED_NAMELIST_OPTION: + continue + + # The `entry_id_pg` root element contains many `entry` elements. + # Each `entry` element describes a namelist option, indicated by its `id` attribute. + entry_element = ET.SubElement(entry_id_pg_element, 'entry', {'id': transform_name(namelist_option.attrib['name'].strip().lower())}) + + # The `category` element. + category_element = ET.SubElement(entry_element, 'category') + category_element.text = 'mpas' + + # The `desc` element. + desc_text = ' '.join(namelist_option.attrib['description'].strip(' .').split()) + desc_text = desc_text[0].upper() + desc_text[1:] + desc_text = '\n' + textwrap.fill(desc_text, 80, initial_indent=INDENT_PER_LEVEL * 3, subsequent_indent=INDENT_PER_LEVEL * 3) + '\n' + INDENT_PER_LEVEL * 2 + + desc_element = ET.SubElement(entry_element, 'desc') + desc_element.text = desc_text + + # The `group` element. + group_element = ET.SubElement(entry_element, 'group') + group_element.text = transform_name(namelist_group.attrib['name'].strip().lower()) + + # The `type` element. + # The `values` element and its containing `value` element. + type_text = namelist_option.attrib['type'].strip().lower() + value_text = namelist_option.attrib['default_value'] + + # Do some sanitization. + if type_text.startswith('ch'): + type_text = 'char*256' + value_text = value_text.strip() + + if not value_text.isascii(): + raise ValueError('"' + value_text + '" is not ASCII') + elif type_text.startswith('in'): + type_text = 'integer' + value_text = canonicalize_int(value_text) + elif type_text.startswith('lo'): + type_text = 'logical' + value_text = value_text.strip().lower() + + if value_text.startswith(('t', '.t')): + value_text = '.true.' + elif value_text.startswith(('f', '.f')): + value_text = '.false.' + else: + raise ValueError('"' + value_text + '" does not represent a logical') + elif type_text.startswith('re'): + type_text = 'real' + value_text = canonicalize_real(value_text) + else: + raise ValueError('"' + type_text + '" is not a supported type') + + type_element = ET.SubElement(entry_element, 'type') + type_element.text = type_text + + values_element = ET.SubElement(entry_element, 'values') + + if entry_element.attrib['id'] in OVERRIDDEN_NAMELIST_OPTION: + for value_string in OVERRIDDEN_NAMELIST_OPTION[entry_element.attrib['id']]: + value_element = ET.fromstring(value_string) + values_element.append(value_element) + else: + value_element = ET.SubElement(values_element, 'value') + value_element.text = value_text + + # Sort the `entry` elements for result stability except for the comment element at index 0. + entry_id_pg_element[1:] = sorted(entry_id_pg_element.iterfind('entry'), key=lambda entry: entry.get('id')) + nml_xml_et = ET.ElementTree(entry_id_pg_element) + + return nml_xml_et + +def canonicalize_int(s: str) -> str: + ''' + Canonicalize a string that represents an integer. + ''' + + try: + i = int(s) + except Exception as e: + raise ValueError('"' + s + '" does not represent an integer') from e + + return str(i) + +def canonicalize_real(s: str) -> str: + ''' + Canonicalize a string that represents a real. + ''' + + try: + f = float(s) + except Exception as e: + raise ValueError('"' + s + '" does not represent a real') from e + + return str(f) + +def write_element_tree(xml_file: str, xml_et: ET.ElementTree) -> None: + ''' + Write element tree into XML file. + ''' + + ET.indent(xml_et, space=INDENT_PER_LEVEL) + + xml_et.write( + xml_file, + encoding='UTF-8', + xml_declaration=True, + default_namespace=None, + method='xml', + short_empty_elements=False + ) + + # The file written by `ElementTree.write()` contains no newline at end of file. Add it manually. + with open(xml_file, 'a', encoding='utf-8') as nml_xml_file: + nml_xml_file.write('\n') + +if __name__ == '__main__': + arg = parse_argument() + reg_xml_element_tree = parse_xml(arg.reg_xml) + nml_xml_element_tree = translate_element_tree(reg_xml_element_tree) + + write_element_tree(arg.nml_xml, nml_xml_element_tree) + print('Generated ' + arg.nml_xml) + + if arg.nml_xsd is not None: + if validate_xml(arg.nml_xml, arg.nml_xsd): + print('Successfully validated ' + arg.nml_xml + ' against ' + arg.nml_xsd) + else: + print('Failed to validate ' + arg.nml_xml + ' against ' + arg.nml_xsd) From 6c43d37437b339375b990f7f26e060cd54cffe53 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 5 Apr 2024 11:40:32 -0600 Subject: [PATCH 19/24] Update namelist definition for MPAS dynamical core From now on, it can be generated automatically from MPAS registry. --- .../mpas/namelist_definition_mpas_dycore.xml | 1334 ++++++++--------- 1 file changed, 626 insertions(+), 708 deletions(-) diff --git a/src/dynamics/mpas/namelist_definition_mpas_dycore.xml b/src/dynamics/mpas/namelist_definition_mpas_dycore.xml index 874b6d86..650ab0e8 100644 --- a/src/dynamics/mpas/namelist_definition_mpas_dycore.xml +++ b/src/dynamics/mpas/namelist_definition_mpas_dycore.xml @@ -1,710 +1,628 @@ - + - - - - char*256 - mpas - mpas_nhyd_model - - Time integration scheme in MPAS. - - Possible values: `SRK3' - Default: SRK3 - - - SRK3 - - - - integer - mpas - mpas_nhyd_model - - Order for RK time integration in MPAS. - - Possible values: 2 or 3 - Default: 2 - - - 2 - - - - real - mpas - mpas_nhyd_model - - Model time step, seconds in MPAS. - - Possible values: Positive real values - Default: 720.0 - - - 720.0 - - - - logical - mpas - mpas_nhyd_model - - Whether to super-cycle scalar transport in MPAS. - - Possible values: Logical values - Default: true - - - true - - - - integer - mpas - mpas_nhyd_model - - Number of acoustic steps per full RK step in MPAS. - - Possible values: Positive, even integer values, typically 2 or 6 depending - on transport splitting - Default: 2 - - - 2 - - - - integer - mpas - mpas_nhyd_model - - When config_split_dynamics_transport = T, the number of RK steps per - transport step in MPAS. - - Possible values: Positive integer values - Default: 3 - - - 3 - - - - real - mpas - mpas_nhyd_model - - $\nabla^2$ eddy viscosity for horizontal diffusion of momentum in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of momentum in - MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^2$ eddy viscosity for vertical diffusion of momentum in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^2$ eddy viscosity for horizontal diffusion of theta in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of theta in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - $\nabla^2$ eddy viscosity for vertical diffusion of theta in MPAS. - - Possible values: Positive real values - Default: 0.0 - - - 0.0 - - - - char*256 - mpas - mpas_nhyd_model - - Formulation of horizontal mixing in MPAS. - - Possible values: `2d_fixed' or `2d_smagorinsky' - Default: 2d_smagorinsky - - - 2d_smagorinsky - - - - real - mpas - mpas_nhyd_model - - Horizontal length scale, used by the Smagorinsky formulation of horizontal - diffusion and by 3-d divergence damping in MPAS. - - Possible values: Positive real values. A zero value implies that the - length scale is prescribed by the nominalMinDc value in the input file. - Default: 0.0 - - - 0.0 - - - - real - mpas - mpas_nhyd_model - - Scaling coefficient of $\delta x^3$ to obtain $\nabla^4$ diffusion - coefficient in MPAS. - - Possible values: Non-negative real values - Default: 0.05 - - - 0.05 - - - - real - mpas - mpas_nhyd_model - - Scaling factor for the divergent component of $\nabla^4 u$ calculation in - MPAS. - - Possible values: Positive real values - Default: 10.0 - - - 10.0 - - - - integer - mpas - mpas_nhyd_model - - Horizontal advection order for w in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Horizontal advection order for theta in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Horizontal advection order for scalars in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Vertical advection order for normal velocities (u) in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Vertical advection order for w in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Vertical advection order for theta in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - integer - mpas - mpas_nhyd_model - - Vertical advection order for scalars in MPAS. - - Possible values: 2, 3, or 4 - Default: 3 - - - 3 - - - - logical - mpas - mpas_nhyd_model - - Whether to advect scalar fields in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - logical - mpas - mpas_nhyd_model - - Whether to enable positive-definite advection of scalars in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - logical - mpas - mpas_nhyd_model - - Whether to enable monotonic limiter in scalar advection in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - real - mpas - mpas_nhyd_model - - Upwinding coefficient in the 3rd order advection scheme in MPAS. - - Possible values: 0 $\leq$ config_coef_3rd_order $\leq$ 1 - Default: 0.25 - - - 0.25 - - - - real - mpas - mpas_nhyd_model - - Dimensionless empirical parameter relating the strain tensor to the eddy - viscosity in the Smagorinsky turbulence model in MPAS. - - Possible values: Real values typically in the range 0.1 to 0.4 - Default: 0.125 - - - 0.125 - - - - logical - mpas - mpas_nhyd_model - - Mix full $\theta$ and $u$ fields, or mix perturbation from intitial state - in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - real - mpas - mpas_nhyd_model - - Off-centering parameter for the vertically implicit acoustic integration - in MPAS. - - Possible values: Positive real values - Default: 0.1 - - - 0.1 - - - - real - mpas - mpas_nhyd_model - - 3-d divergence damping coefficient in MPAS. - - Possible values: Positive real values - Default: 0.1 - - - 0.1 - - - - real - mpas - mpas_nhyd_model - - Amount of upwinding in APVM in MPAS. - - Possible values: 0 $\leq$ config_apvm_upwinding $\leq$ 1 - Default: 0.5 - - - 0.5 - - - - logical - mpas - mpas_nhyd_model - - Scale eddy viscosities with mesh-density function for horizontal diffusion - in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - real - mpas - mpas_nhyd_model - - Coefficient for the divergent component of the Laplacian filter of - momentum in the relaxation zone in MPAS. - - Possible values: Positive real values - Default: 6.0 - - - 6.0 - - - - real - mpas - mpas_damping - - Height MSL to begin w-damping profile in MPAS. - - Possible values: Positive real values - Default: 22000.0 - - - 22000.0 - - - - real - mpas - mpas_damping - - Maximum w-damping coefficient at model top in MPAS. - - Possible values: 0 $\leq$ config_xnutr $\leq$ 1 - Default: 0.2 - - - 0.2 - - - - real - mpas - mpas_damping - - Coefficient for scaling the 2nd-order horizontal mixing in the mpas_cam - absorbing layer in MPAS. - - Possible values: 0 $\leq$ config_mpas_cam_coef $\leq$ 1, standard value is - 0.2 - Default: 0.0 - - - 0.0 - - - - integer - mpas - mpas_damping - - Number of layers in which to apply cam 2nd-order horizontal filter top of - model; viscosity linearly ramps to zero by layer number from the top in - MPAS. - - Possible values: Positive integer values - Default: 4 - - - 4 - - - - logical - mpas - mpas_damping - - Whether to apply Rayleigh damping on horizontal velocity in the top-most - model levels. The number of levels is specified by the - config_number_rayleigh_damp_u_levels option, and the damping timescale is - specified by the config_rayleigh_damp_u_timescale_days option. in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - real - mpas - mpas_damping - - Timescale, in days (86400 s), for the Rayleigh damping on horizontal - velocity in the top-most model levels. in MPAS. - - Possible values: Positive real values - Default: 5.0 - - - 5.0 - - - - integer - mpas - mpas_damping - - Number of layers in which to apply Rayleigh damping on horizontal velocity - at top of model; damping linearly ramps to zero by layer number from the - top in MPAS. - - Possible values: Positive integer values - Default: 6 - - - 6 - - - - logical - mpas - mpas_limited_area - - Whether to apply lateral boundary conditions in MPAS. - - Possible values: true or false; this option must be set to true for - limited-area simulations and false for global simulations - Default: false - - - false - - - - char*256 - mpas - mpas_decomposition - - Prefix of graph decomposition file, to be suffixed with the MPI task count - in MPAS. - - Possible values: Any valid filename - Default: x1.40962.graph.info.part. - - - x1.40962.graph.info.part. - - - - logical - mpas - mpas_restart - - Whether this run of the model is to restart from a previous restart file - or not in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - logical - mpas - mpas_printout - - Whether to print the global min/max of horizontal normal velocity and - vertical velocity each timestep in MPAS. - - Possible values: .true. or .false. - Default: true - - - true - - - - logical - mpas - mpas_printout - - Whether to print the global min/max of horizontal normal velocity and - vertical velocity each timestep, along with the location in the domain - where those extrema occurred in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - logical - mpas - mpas_printout - - Whether to print the global min/max of scalar fields each timestep in - MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - - - logical - mpas - mpas_assimilation - - Whether this run is within the JEDI data assimilation framework; used to - add temperature and specific humidity as diagnostics in MPAS. - - Possible values: .true. or .false. - Default: false - - - false - - + + + mpas + + Amount of upwinding in APVM + + mpas_nhyd_model + real + + 0.5 + + + + mpas + + Prefix of graph decomposition file, to be suffixed with the MPI task + count + + mpas_decomposition + char*256 + + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa30.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa15.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa12.graph.info.part. + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa15-3.graph.info.part. + + + + mpas + + Coefficient for scaling the 2nd-order horizontal mixing in the + mpas_cam absorbing layer + + mpas_damping + real + + 0.0 + + + + mpas + + Upwinding coefficient in the 3rd order advection scheme + + mpas_nhyd_model + real + + 1.0 + + + + mpas + + Scaling factor for the divergent component of $\nabla^4 u$ + calculation + + mpas_nhyd_model + real + + 10.0 + + + + mpas + + Model time step, seconds + + mpas_nhyd_model + real + + 1800.0 + 900.0 + 450.0 + 225.0 + + + + mpas + + When config_split_dynamics_transport = T, the number of RK steps per + transport step + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Off-centering parameter for the vertically implicit acoustic + integration + + mpas_nhyd_model + real + + 0.1 + + + + mpas + + Whether to use an explicit mapping of blocks to MPI tasks + + mpas_decomposition + logical + + .false. + + + + mpas + + $\nabla^2$ eddy viscosity for horizontal diffusion of momentum + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of momentum + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + Scale eddy viscosities with mesh-density function for horizontal + diffusion + + mpas_nhyd_model + logical + + .true. + + + + mpas + + $\nabla^2$ eddy viscosity for horizontal diffusion of theta + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + $\nabla^4$ eddy hyper-viscosity for horizontal diffusion of theta + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + Method to use for exchanging halos + + mpas_development + char*256 + + mpas_halo + + + + mpas + + Formulation of horizontal mixing + + mpas_nhyd_model + char*256 + + 2d_smagorinsky + + + + mpas + + Horizontal length scale, used by the Smagorinsky formulation of + horizontal diffusion and by 3-d divergence damping + + mpas_nhyd_model + real + + 480000.0 + 120000.0 + 60000.0 + 30000.0 + + + + mpas + + Mix full $\theta$ and $u$ fields, or mix perturbation from intitial + state + + mpas_nhyd_model + logical + + .true. + + + + mpas + + Whether to enable monotonic limiter in scalar advection + + mpas_nhyd_model + logical + + .true. + + + + mpas + + Number of halo layers for fields + + mpas_nhyd_model + integer + + 2 + + + + mpas + + Number of layers in which to apply cam 2nd-order horizontal filter + top of model; viscosity linearly ramps to zero by layer number from + the top + + mpas_damping + integer + + 0 + + + + mpas + + Number of blocks to assign to each MPI task + + mpas_decomposition + integer + + 0 + + + + mpas + + Number of acoustic steps per full RK step + + mpas_nhyd_model + integer + + 2 + + + + mpas + + Number of layers in which to apply Rayleigh damping on horizontal + velocity at top of model; damping linearly ramps to zero by layer + number from the top + + mpas_damping + integer + + 5 + + + + mpas + + Number of tasks to perform file I/O + + mpas_io + integer + + 0 + + + + mpas + + Stride between file I/O tasks + + mpas_io + integer + + 1 + + + + mpas + + Whether to enable positive-definite advection of scalars + + mpas_nhyd_model + logical + + .false. + + + + mpas + + Whether to print the global min/max of horizontal normal velocity + and vertical velocity each timestep, along with the location in the + domain where those extrema occurred + + mpas_printout + logical + + .true. + + + + mpas + + Whether to print the global min/max of scalar fields each timestep + + mpas_printout + logical + + .false. + + + + mpas + + Whether to print the global min/max of horizontal normal velocity + and vertical velocity each timestep + + mpas_printout + logical + + .true. + + + + mpas + + Prefix of block mapping file + + mpas_decomposition + char*256 + + graph.info.part. + + + + mpas + + Whether to apply Rayleigh damping on horizontal velocity in the top- + most model levels. The number of levels is specified by the + config_number_rayleigh_damp_u_levels option, and the damping + timescale is specified by the config_rayleigh_damp_u_timescale_days + option + + mpas_damping + logical + + .true. + + + + mpas + + Timescale, in days (86400 s), for the Rayleigh damping on horizontal + velocity in the top-most model levels + + mpas_damping + real + + 5.0 + + + + mpas + + Coefficient for the divergent component of the Laplacian filter of + momentum in the relaxation zone + + mpas_nhyd_model + real + + 6.0 + + + + mpas + + Filename used to store most recent restart time stamp + + mpas_io + char*256 + + restart_timestamp + + + + mpas + + Horizontal advection order for scalars + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Whether to advect scalar fields + + mpas_nhyd_model + logical + + .true. + + + + mpas + + Vertical advection order for scalars + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Dimensionless empirical parameter relating the strain tensor to the + eddy viscosity in the Smagorinsky turbulence model + + mpas_nhyd_model + real + + 0.125 + + + + mpas + + 3-d divergence damping coefficient + + mpas_nhyd_model + real + + 0.1 + + + + mpas + + Whether to super-cycle scalar transport + + mpas_nhyd_model + logical + + .true. + + + + mpas + + Horizontal advection order for theta + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Vertical advection order for theta + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Time integration scheme + + mpas_nhyd_model + char*256 + + SRK3 + + + + mpas + + Order for RK time integration + + mpas_nhyd_model + integer + + 2 + + + + mpas + + Vertical advection order for normal velocities (u) + + mpas_nhyd_model + integer + + 3 + + + + mpas + + $\nabla^2$ eddy viscosity for vertical diffusion of momentum + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + $\nabla^2$ eddy viscosity for vertical diffusion of theta + + mpas_nhyd_model + real + + 0.0 + + + + mpas + + Scaling coefficient of $\delta x^3$ to obtain $\nabla^4$ diffusion + coefficient + + mpas_nhyd_model + real + + 0.05 + + + + mpas + + Horizontal advection order for w + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Vertical advection order for w + + mpas_nhyd_model + integer + + 3 + + + + mpas + + Maximum w-damping coefficient at model top + + mpas_damping + real + + 0.2 + + + + mpas + + Height MSL to begin w-damping profile + + mpas_damping + real + + 22000.0 + + From e36fcdbf562c4849d73975dfa72ffd6185067530 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Sat, 16 Nov 2024 06:08:19 -0700 Subject: [PATCH 20/24] Update namelist definition for CAM-SIMA Add MPAS initial files for running analytic cases. --- cime_config/namelist_definition_cam.xml | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/cime_config/namelist_definition_cam.xml b/cime_config/namelist_definition_cam.xml index fce70728..3054ebb6 100644 --- a/cime_config/namelist_definition_cam.xml +++ b/cime_config/namelist_definition_cam.xml @@ -149,8 +149,16 @@ ${DIN_LOC_ROOT}/atm/waccm/ic/waccm5_1850_ne30np4_L70_0001-01-11-00000_c151217.nc ${DIN_LOC_ROOT}/atm/waccm/ic/fw2000_ne30np4_L70_c181221.nc ${DIN_LOC_ROOT}/atm/cam/inic/gaus/cami_0000-09-01_64x128_L30_c031210.nc - ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120_L32_notopo_coords_c201216.nc - ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480_L32_notopo_coords_c201125.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480_L32_notopo_coords_c240507.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120_L32_notopo_coords_c240507.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60_L32_notopo_coords_c240507.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa30_L32_notopo_coords_c240507.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480_L58_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120_L58_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60_L58_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa480_L93_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa120_L93_notopo_coords_c240814.nc + ${DIN_LOC_ROOT}/atm/cam/inic/mpas/mpasa60_L93_notopo_coords_c240814.nc From 4212d77ed5525ca5f6fbc3266566dbaa93207e53 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Thu, 2 Jan 2025 11:31:45 -0700 Subject: [PATCH 21/24] Adjust some wording --- src/dynamics/mpas/assets/generate_namelist_definition.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dynamics/mpas/assets/generate_namelist_definition.py b/src/dynamics/mpas/assets/generate_namelist_definition.py index f72aa69d..12dd0c54 100755 --- a/src/dynamics/mpas/assets/generate_namelist_definition.py +++ b/src/dynamics/mpas/assets/generate_namelist_definition.py @@ -93,7 +93,7 @@ def parse_argument() -> argparse.Namespace: default='Namelist.xml', type=str, required=False, - help='XML namelist definition file.', + help='XML CAM-SIMA namelist definition file.', dest='nml_xml' ) parser.add_argument( @@ -101,7 +101,7 @@ def parse_argument() -> argparse.Namespace: default=None, type=str, required=False, - help='XML schema for namelist definition file.', + help='XML schema for CAM-SIMA namelist definition file.', dest='nml_xsd' ) @@ -158,7 +158,7 @@ def translate_element_tree(reg_xml_et: ET.ElementTree) -> ET.ElementTree: INDENT_PER_LEVEL * 2 + 'MPAS dycore' + '\n' + '\n' + INDENT_PER_LEVEL * 2 + 'Note to developers/maintainers:' + '\n' + - INDENT_PER_LEVEL * 2 + 'This file is auto-generated from MPAS registry. Do not edit directly.' + '\n' + + INDENT_PER_LEVEL * 2 + 'This file is auto-generated from the MPAS registry. Do not edit directly.' + '\n' + INDENT_PER_LEVEL * 2 + 'Instead, use the Python script at `src/dynamics/mpas/assets/generate_namelist_definition.py`.' + '\n' + INDENT_PER_LEVEL ) From a9c253f9ae5e8b24140b25924724fc8e136ba29c Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Thu, 2 Jan 2025 11:34:42 -0700 Subject: [PATCH 22/24] Regenerate CAM-SIMA namelist definition file --- src/dynamics/mpas/namelist_definition_mpas_dycore.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamics/mpas/namelist_definition_mpas_dycore.xml b/src/dynamics/mpas/namelist_definition_mpas_dycore.xml index 650ab0e8..c59570bf 100644 --- a/src/dynamics/mpas/namelist_definition_mpas_dycore.xml +++ b/src/dynamics/mpas/namelist_definition_mpas_dycore.xml @@ -4,7 +4,7 @@ MPAS dycore Note to developers/maintainers: - This file is auto-generated from MPAS registry. Do not edit directly. + This file is auto-generated from the MPAS registry. Do not edit directly. Instead, use the Python script at `src/dynamics/mpas/assets/generate_namelist_definition.py`. --> From 155ef57c88fd33a5cbf18b5f1421d9b64ee82e0e Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Thu, 2 Jan 2025 13:28:05 -0700 Subject: [PATCH 23/24] Move `use` statements to the smallest scope possible ... Also for kind and length parameters. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 3 ++- src/dynamics/mpas/dyn_comp.F90 | 14 ++++++++++--- src/dynamics/mpas/dyn_coupling.F90 | 16 ++++++++------ src/dynamics/mpas/dyn_grid.F90 | 21 ++++++++++++------- src/dynamics/mpas/stepon.F90 | 7 ++++--- 5 files changed, 41 insertions(+), 20 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index ea793dcf..b4ec0921 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -17,7 +17,7 @@ module dyn_mpas_subdriver #endif ! Module(s) from MPAS. use mpas_derived_types, only: core_type, domain_type - use mpas_kind_types, only: rkind, r4kind, r8kind, strkind + use mpas_kind_types, only: rkind, strkind implicit none @@ -1897,6 +1897,7 @@ subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compa use mpas_derived_types, only: field0dchar, field1dchar, & field0dinteger, field1dinteger, field2dinteger, field3dinteger, & field0dreal, field1dreal, field2dreal, field3dreal, field4dreal, field5dreal + use mpas_kind_types, only: r4kind, r8kind use mpas_pool_routines, only: mpas_pool_get_field class(mpas_dynamical_core_type), intent(in) :: self diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 453c6bf3..f6e3a10a 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -1,7 +1,4 @@ module dyn_comp - ! Module(s) from CESM Share. - use shr_kind_mod, only: kind_r8 => shr_kind_r8, & - len_cs => shr_kind_cs ! Module(s) from MPAS. use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind, mpas_dynamical_core_type @@ -99,6 +96,7 @@ subroutine dyn_readnl(namelist_path) use time_manager, only: get_start_date, get_stop_date, get_run_duration, timemgr_get_calendar_cf ! Module(s) from CESM Share. use shr_file_mod, only: shr_file_getunit + use shr_kind_mod, only: len_cs => shr_kind_cs use shr_pio_mod, only: shr_pio_getiosys ! Module(s) from external libraries. use pio, only: iosystem_desc_t @@ -312,6 +310,8 @@ subroutine check_topography_data(pio_file) use cam_field_read, only: cam_read_field use cam_logfile, only: debug_output, debugout_none use dynconst, only: constant_g => gravit + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 ! Module(s) from external libraries. use pio, only: file_desc_t, pio_file_is_open @@ -376,6 +376,9 @@ end subroutine check_topography_data !> Set analytic initial condition for MPAS. !> (KCW, 2024-05-22) subroutine set_analytic_initial_condition() + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition' integer, allocatable :: global_grid_index(:) real(kind_r8), allocatable :: buffer_2d_real(:, :), buffer_3d_real(:, :, :) @@ -845,6 +848,8 @@ subroutine dyn_exchange_constituent_states(direction, exchange, conversion) ! Module(s) from CCPP. use cam_ccpp_cap, only: cam_constituents_array use ccpp_kinds, only: kind_phys + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 character(*), intent(in) :: direction logical, intent(in) :: exchange @@ -1156,6 +1161,9 @@ end subroutine dyn_variable_dump !> Helper function for reversing the order of elements in `array`. !> (KCW, 2024-07-17) pure function reverse(array) + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + real(kind_r8), intent(in) :: array(:) real(kind_r8) :: reverse(size(array)) diff --git a/src/dynamics/mpas/dyn_coupling.F90 b/src/dynamics/mpas/dyn_coupling.F90 index 831bcd02..aa176a73 100644 --- a/src/dynamics/mpas/dyn_coupling.F90 +++ b/src/dynamics/mpas/dyn_coupling.F90 @@ -1,10 +1,4 @@ module dyn_coupling - ! Module(s) from CESM Share. - use shr_kind_mod, only: kind_r8 => shr_kind_r8, & - len_cx => shr_kind_cx - ! Module(s) from MPAS. - use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind - implicit none private @@ -18,6 +12,10 @@ module dyn_coupling subroutine dynamics_to_physics_coupling() ! Module(s) from CAM-SIMA. use dyn_comp, only: dyn_debug_print, dyn_exchange_constituent_states, ncells_solve + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling' integer :: column_index @@ -308,6 +306,8 @@ subroutine set_physics_state_external() use geopotential_temp, only: geopotential_temp_run use qneg, only: qneg_run use static_energy, only: update_dry_static_energy_run + ! Module(s) from CESM Share. + use shr_kind_mod, only: len_cx => shr_kind_cx character(*), parameter :: subname = 'dyn_coupling::dynamics_to_physics_coupling::set_physics_state_external' character(len_cx) :: cerr @@ -411,6 +411,10 @@ end subroutine dynamics_to_physics_coupling subroutine physics_to_dynamics_coupling() ! Module(s) from CAM-SIMA. use dyn_comp, only: dyn_exchange_constituent_states + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind character(*), parameter :: subname = 'dyn_coupling::physics_to_dynamics_coupling' integer, pointer :: index_qv diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 2b4e3468..874b751d 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -1,12 +1,6 @@ module dyn_grid ! Module(s) from CAM-SIMA. use cam_grid_support, only: max_hcoordname_len - use cam_map_utils, only: kind_imap => imap - use physics_column_type, only: kind_pcol - ! Module(s) from CESM Share. - use shr_kind_mod, only: kind_r8 => shr_kind_r8 - ! Module(s) from MPAS. - use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind implicit none @@ -115,6 +109,10 @@ subroutine init_reference_pressure() use std_atm_profile, only: std_atm_pres use string_utils, only: stringify use vert_coord, only: pver, pverp + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind character(*), parameter :: subname = 'dyn_grid::init_reference_pressure' ! Number of pure pressure levels at model top. @@ -221,10 +219,14 @@ subroutine init_physics_grid() use cam_abortutils, only: check_allocate use dyn_comp, only: mpas_dynamical_core, ncells_global, ncells_solve, sphere_radius use dynconst, only: constant_pi => pi, rad_to_deg - use physics_column_type, only: physics_column_t + use physics_column_type, only: kind_pcol, physics_column_t use physics_grid, only: phys_grid_init use spmd_utils, only: iam use string_utils, only: stringify + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind character(*), parameter :: subname = 'dyn_grid::init_physics_grid' character(max_hcoordname_len), allocatable :: dyn_attribute_name(:) @@ -311,12 +313,17 @@ subroutine define_cam_grid() use cam_abortutils, only: check_allocate use cam_grid_support, only: cam_grid_attribute_register, cam_grid_register, & horiz_coord_create, horiz_coord_t + use cam_map_utils, only: kind_imap => imap use dyn_comp, only: dyn_debug_print, mpas_dynamical_core, & ncells_global, nedges_global, nvertices_global, & ncells_solve, nedges_solve, nvertices_solve, & sphere_radius use dynconst, only: constant_pi => pi, rad_to_deg use string_utils, only: stringify + ! Module(s) from CESM Share. + use shr_kind_mod, only: kind_r8 => shr_kind_r8 + ! Module(s) from MPAS. + use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind character(*), parameter :: subname = 'dyn_grid::define_cam_grid' integer :: i diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90 index dfaf2c18..a89c7c75 100644 --- a/src/dynamics/mpas/stepon.F90 +++ b/src/dynamics/mpas/stepon.F90 @@ -1,7 +1,4 @@ module stepon - ! Module(s) from CCPP. - use ccpp_kinds, only: kind_phys - implicit none private @@ -31,6 +28,8 @@ subroutine stepon_timestep_init(dtime_phys, cam_runtime_opts, phys_state, phys_t use physics_types, only: physics_state, physics_tend use runtime_obj, only: runtime_options use time_manager, only: get_step_size + ! Module(s) from CCPP. + use ccpp_kinds, only: kind_phys real(kind_phys), intent(out) :: dtime_phys type(runtime_options), intent(in) :: cam_runtime_opts @@ -69,6 +68,8 @@ subroutine stepon_run3(dtime_phys, cam_runtime_opts, cam_out, phys_state, dyn_in use dyn_comp, only: dyn_export_t, dyn_import_t, dyn_run use physics_types, only: physics_state use runtime_obj, only: runtime_options + ! Module(s) from CCPP. + use ccpp_kinds, only: kind_phys real(kind_phys), intent(in) :: dtime_phys type(runtime_options), intent(in) :: cam_runtime_opts From 9b500857b7784775a347b19b6f118aeabf9bb4a5 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 13 Jan 2025 17:57:14 -0700 Subject: [PATCH 24/24] Existing test failures are all clear --- test/existing-test-failures.txt | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/existing-test-failures.txt b/test/existing-test-failures.txt index 975b8681..2328b59d 100644 --- a/test/existing-test-failures.txt +++ b/test/existing-test-failures.txt @@ -1,3 +1 @@ -SMS_Ln2.mpasa480_mpasa480.FKESSLER.derecho_intel.cam-outfrq_kessler_mpas_derecho (Overall: FAIL) -SMS_Ln2.mpasa480_mpasa480.FKESSLER.derecho_gnu.cam-outfrq_kessler_mpas_derecho (Overall: FAIL) - - will fail until MPAS is fully integrated +All clear!