From 01af7c752ee39f879cdec8622210f5c390a3149b Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 29 May 2024 11:12:43 -0600 Subject: [PATCH 01/31] Remove relative humidity from MPAS input stream Relative humidity (`relhum`) is a diagnostic variable. It should not be included in MPAS input stream. --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 9d2085f1..8b65d20a 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -231,7 +231,6 @@ end subroutine model_error_if type(var_info_type), parameter :: input_var_info_list(*) = [ & var_info_type('Time' , 'real' , 1), & var_info_type('initial_time' , 'character' , 0), & - var_info_type('relhum' , 'real' , 3), & var_info_type('rho' , 'real' , 3), & var_info_type('rho_base' , 'real' , 3), & var_info_type('scalars' , 'real' , 3), & From 6b69ea90bdf1f734439b6d3516d377283bc21b80 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 29 May 2024 11:22:39 -0600 Subject: [PATCH 02/31] Fix style --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 8b65d20a..339ca7d3 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -483,7 +483,7 @@ subroutine dyn_mpas_init_phase1(self, mpi_comm, model_error_impl, log_unit, mpas ! We need: ! 1) `domain_ptr` to be allocated; ! 2) `dmpar_init` to be completed for accessing `dminfo`; - ! 3) `*_setup_core` to assign the `setup_log` function pointer. + ! 3) `*_setup_core` to assign the `setup_log` procedure pointer. ierr = self % domain_ptr % core % setup_log(self % domain_ptr % loginfo, self % domain_ptr, unitnumbers=mpas_log_unit) if (ierr /= 0) then @@ -1193,7 +1193,7 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file call add_stream_attribute('is_periodic', self % domain_ptr % is_periodic) call add_stream_attribute('mesh_spec', self % domain_ptr % mesh_spec) call add_stream_attribute('on_a_sphere', self % domain_ptr % on_a_sphere) - call add_stream_attribute('parent_id', self % domain_ptr % parent_id) + call add_stream_attribute('parent_id', self % domain_ptr % parent_id) call add_stream_attribute('sphere_radius', self % domain_ptr % sphere_radius) call add_stream_attribute('x_period', self % domain_ptr % x_period) call add_stream_attribute('y_period', self % domain_ptr % y_period) From 61ccb16a7221ba7b13a2e306ead936ac288f4cb0 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 29 May 2024 11:37:49 -0600 Subject: [PATCH 03/31] Handle the allocation of constituents more carefully More elaborate error checking. Also use a more descriptive name for number of constituents. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 55 ++++++++++++++++--- 1 file changed, 46 insertions(+), 9 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 339ca7d3..0089c642 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -35,7 +35,7 @@ module dyn_mpas_subdriver 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_get_dimension, & + mpas_pool_add_dimension, mpas_pool_get_dimension, & mpas_pool_get_field, mpas_pool_get_field_info use mpas_stream_manager, only: postread_reindex, prewrite_reindex, postwrite_reindex @@ -61,6 +61,7 @@ end subroutine model_error_if logical, public :: debug_output = .false. + ! Initialized by `dyn_mpas_init_phase1`. integer :: log_unit = output_unit integer :: mpi_comm = mpi_comm_null integer :: mpi_rank = 0 @@ -71,6 +72,9 @@ end subroutine model_error_if type(core_type), pointer :: corelist => null() type(domain_type), pointer :: domain_ptr => null() + + ! Initialized by `dyn_mpas_init_phase3`. + integer :: number_of_constituents = 0 contains private @@ -702,37 +706,50 @@ end subroutine dyn_mpas_init_phase2 !------------------------------------------------------------------------------- ! subroutine dyn_mpas_init_phase3 ! - !> \brief Tracks `mpas_init` up to the point of calling `core_init` + !> \brief Tracks `mpas_init` up to the point of calling `atm_core_init` !> \author Michael Duda !> \date 19 April 2019 !> \details !> This subroutine follows the stand-alone MPAS subdriver after the check on !> the existence of the `streams.` file up to, but not including, - !> the point where `core_init` is called. It completes MPAS framework + !> the point where `atm_core_init` is called. It completes MPAS framework !> initialization, including the allocation of all blocks and fields managed - !> by MPAS. + !> by MPAS. However, scalars are allocated but not yet defined. + !> `dyn_mpas_define_scalar` must be called afterwards. Also note that MPAS uses + !> the term "scalar", but CAM-SIMA calls it "constituent". !> \addenda !> Ported and refactored for CAM-SIMA. (KCW, 2024-03-06) ! !------------------------------------------------------------------------------- - subroutine dyn_mpas_init_phase3(self, cam_pcnst, pio_file) - class(mpas_dynamical_core_type), intent(in) :: self - integer, intent(in) :: cam_pcnst + subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file) + class(mpas_dynamical_core_type), intent(inout) :: self + integer, intent(in) :: number_of_constituents type(file_desc_t), pointer, intent(in) :: pio_file character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase3' character(strkind) :: mesh_filename integer :: mesh_format + integer, pointer :: num_scalars => null() + type(mpas_pool_type), pointer :: mpas_pool => null() call self % debug_print(subname // ' entered') - call self % debug_print('Number of constituents is ', [cam_pcnst]) + ! In MPAS, there must be at least one constituent, `qv`, which denotes water vapor mixing ratio. + ! Because MPAS has some hard-coded array accesses through the `index_qv` index, it will crash + ! (i.e., segmentation fault due to invalid memory access) if `qv` is not allocated. + if (number_of_constituents < 1) then + self % number_of_constituents = 1 + else + self % number_of_constituents = number_of_constituents + end if + + call self % debug_print('Number of constituents is ', [self % number_of_constituents]) ! Adding a config named `cam_pcnst` with the number of constituents will indicate to MPAS that ! it is operating as a dynamical core, and therefore it needs to allocate scalars separately ! from other Registry-defined fields. The special logic is located in `atm_setup_block`. ! This must be done before calling `mpas_bootstrap_framework_phase1`. - call mpas_pool_add_config(self % domain_ptr % configs, 'cam_pcnst', cam_pcnst) + call mpas_pool_add_config(self % domain_ptr % configs, 'cam_pcnst', self % number_of_constituents) ! Not actually used because a PIO file descriptor is directly supplied. mesh_filename = 'external mesh' @@ -758,6 +775,26 @@ subroutine dyn_mpas_init_phase3(self, cam_pcnst, pio_file) ! Finish setting up fields. call mpas_bootstrap_framework_phase2(self % domain_ptr, pio_file_desc=pio_file) + ! `num_scalars` is a dimension variable, but it only exists in MPAS `state` pool. + ! Fix this inconsistency by also adding it to MPAS `dimension` pool. + call self % get_pool_pointer(mpas_pool, 'state') + + call mpas_pool_get_dimension(mpas_pool, 'num_scalars', num_scalars) + + if (.not. associated(num_scalars)) then + call self % model_error('Failed to find variable "num_scalars"', subname, __LINE__) + end if + + ! While we are at it, check if its value is consistent. + if (num_scalars /= self % number_of_constituents) then + call self % model_error('Failed to allocate constituents', subname, __LINE__) + end if + + call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'num_scalars', num_scalars) + + nullify(mpas_pool) + nullify(num_scalars) + call self % debug_print(subname // ' completed') end subroutine dyn_mpas_init_phase3 From bf2a21d2763f9e129bb2f50b8f986a81989693d4 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 29 May 2024 11:40:49 -0600 Subject: [PATCH 04/31] Implement MPAS subdriver * Finish MPAS dynamical core initialization. * Support for deferring the definition of constituents until run-time. * Accessor functions/subroutines for querying: * Local mesh dimensions. * Constituent names and indexes. * Mapping between constituent indexes and MPAS scalar indexes. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 900 +++++++++++++++++- 1 file changed, 899 insertions(+), 1 deletion(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 0089c642..09a5ff57 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -15,12 +15,20 @@ module dyn_mpas_subdriver pio_inq_varid, pio_inq_varndims, pio_noerr ! Modules from MPAS. + use atm_core, only: 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 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_bootstrapping, only: mpas_bootstrap_framework_phase1, mpas_bootstrap_framework_phase2 + use mpas_constants, only: mpas_constants_compute_derived 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_start_time, & mpas_io_native_precision, mpas_io_pnetcdf, mpas_io_read, mpas_io_write, & field0dchar, field1dchar, & field0dinteger, field1dinteger, field2dinteger, field3dinteger, & @@ -36,8 +44,12 @@ module dyn_mpas_subdriver 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_get_field, mpas_pool_get_field_info, & + mpas_pool_initialize_time_levels use mpas_stream_manager, only: postread_reindex, prewrite_reindex, postwrite_reindex + use mpas_string_utils, only: mpas_string_replace + use mpas_timekeeping, only: mpas_get_clock_time, mpas_get_time + use mpas_vector_operations, only: mpas_initialize_vectors implicit none @@ -75,6 +87,15 @@ end subroutine model_error_if ! Initialized by `dyn_mpas_init_phase3`. integer :: number_of_constituents = 0 + + ! Initialized by `dyn_mpas_define_scalar`. + character(strkind), allocatable :: constituent_name(:) + integer, allocatable :: index_constituent_to_mpas_scalar(:) + integer, allocatable :: index_mpas_scalar_to_constituent(:) + logical, allocatable :: is_water_species(:) + + ! Initialized by `dyn_mpas_init_phase4`. + integer :: coupling_time_interval contains private @@ -83,12 +104,23 @@ end subroutine model_error_if procedure, pass, public :: read_namelist => dyn_mpas_read_namelist procedure, pass, public :: init_phase2 => dyn_mpas_init_phase2 procedure, pass, public :: init_phase3 => dyn_mpas_init_phase3 + procedure, pass, public :: define_scalar => dyn_mpas_define_scalar procedure, pass, public :: read_write_stream => dyn_mpas_read_write_stream procedure, pass :: init_stream_with_pool => dyn_mpas_init_stream_with_pool procedure, pass, public :: exchange_halo => dyn_mpas_exchange_halo + procedure, pass, public :: compute_unit_vector => dyn_mpas_compute_unit_vector + procedure, pass, public :: compute_edge_wind => dyn_mpas_compute_edge_wind + procedure, pass, public :: init_phase4 => dyn_mpas_init_phase4 ! Accessor subroutines for users to access internal states of MPAS dynamical core. + procedure, pass, public :: get_constituent_name => dyn_mpas_get_constituent_name + procedure, pass, public :: get_constituent_index => dyn_mpas_get_constituent_index + + procedure, pass, public :: map_mpas_scalar_index => dyn_mpas_map_mpas_scalar_index + procedure, pass, public :: map_constituent_index => dyn_mpas_map_constituent_index + + procedure, pass, public :: get_local_mesh_dimension => dyn_mpas_get_local_mesh_dimension procedure, pass, public :: get_global_mesh_dimension => dyn_mpas_get_global_mesh_dimension procedure, pass :: get_pool_pointer => dyn_mpas_get_pool_pointer @@ -798,6 +830,267 @@ subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file) call self % debug_print(subname // ' completed') end subroutine dyn_mpas_init_phase3 + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_define_scalar + ! + !> \brief Defines the names of constituents at run-time + !> \author Michael Duda + !> \date 21 May 2020 + !> \details + !> Given arrays of constituent names and their corresponding waterness, which + !> must have sizes equal to the number of constituents used to call + !> `dyn_mpas_init_phase3`, this subroutine defines the scalars inside MPAS. + !> Note that MPAS uses the term "scalar", but CAM-SIMA calls it "constituent". + !> Furthermore, because MPAS expects all water scalars to appear in a + !> contiguous index range, this subroutine may reorder the scalars to satisfy + !> this constrain. Index mapping between MPAS scalars and constituent names + !> can be looked up through `index_constituent_to_mpas_scalar` and + !> `index_mpas_scalar_to_constituent`. + !> \addenda + !> Ported and refactored for CAM-SIMA. (KCW, 2024-05-19) + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species) + class(mpas_dynamical_core_type), intent(inout) :: self + character(*), intent(in) :: constituent_name(:) + logical, intent(in) :: is_water_species(:) + + ! Possible CCPP standard names of `qv`, which denotes water vapor mixing ratio. + ! They are hard-coded here because MPAS needs to know where `qv` is. + ! Index 1 is exactly what MPAS wants. Others also work, but need to be converted. + character(*), parameter :: mpas_scalar_qv_standard_name(*) = [ character(strkind) :: & + 'water_vapor_mixing_ratio_wrt_dry_air', & + 'water_vapor_mixing_ratio_wrt_moist_air', & + 'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water' & + ] + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_define_scalar' + integer :: i, j, ierr + integer :: index_qv, index_water_start, index_water_end + integer :: time_level + type(field3dreal), pointer :: field_3d_real => null() + type(mpas_pool_type), pointer :: mpas_pool => null() + + call self % debug_print(subname // ' entered') + + if (self % number_of_constituents == 0) then + call self % model_error('Constituents must be allocated before being defined', subname, __LINE__) + end if + + ! Input sanitization. + + if (size(constituent_name) /= size(is_water_species)) then + call self % model_error('Mismatch between numbers of constituent names and their waterness', subname, __LINE__) + end if + + if (size(constituent_name) == 0 .and. self % number_of_constituents == 1) then + ! If constituent definitions are empty, `qv` is the only constituent per MPAS requirements. + ! See `dyn_mpas_init_phase3` for details. + allocate(self % constituent_name(1), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate constituent_name', subname, __LINE__) + end if + + allocate(self % is_water_species(1), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate is_water_species', subname, __LINE__) + end if + + self % constituent_name(1) = mpas_scalar_qv_standard_name(1) + self % is_water_species(1) = .true. + else + if (size(constituent_name) /= self % number_of_constituents) then + call self % model_error('Mismatch between numbers of constituents and their names', subname, __LINE__) + end if + + if (any(len_trim(adjustl(constituent_name)) > len(self % constituent_name))) then + call self % model_error('Constituent names are too long', subname, __LINE__) + end if + + allocate(self % constituent_name(self % number_of_constituents), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate constituent_name', subname, __LINE__) + end if + + self % constituent_name(:) = adjustl(constituent_name) + + allocate(self % is_water_species(self % number_of_constituents), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate is_water_species', subname, __LINE__) + end if + + self % is_water_species(:) = is_water_species + + if (size(self % constituent_name) /= size(index_unique(self % constituent_name))) then + call self % model_error('Constituent names must be unique', subname, __LINE__) + end if + + ! `qv` must be present in constituents per MPAS requirements. It is a water species by definition. + ! See `dyn_mpas_init_phase3` for details. + index_qv = 0 + + ! Lower index in `mpas_scalar_qv_standard_name` has higher precedence, with index 1 being exactly what MPAS wants. + set_index_qv: do i = 1, size(mpas_scalar_qv_standard_name) + do j = 1, self % number_of_constituents + if (self % constituent_name(j) == mpas_scalar_qv_standard_name(i) .and. self % is_water_species(j)) then + index_qv = j + + ! The best candidate of `qv` has been found. Exit prematurely. + exit set_index_qv + end if + end do + end do set_index_qv + + if (index_qv == 0) then + call self % model_error('Constituent names must contain one of: ' // & + stringify(mpas_scalar_qv_standard_name) // ' and it must be a water species', subname, __LINE__) + end if + end if + + ! Create index mapping between MPAS scalars and constituent names. For example, + ! MPAS scalar index `i` corresponds to constituent index `index_mpas_scalar_to_constituent(i)`. + + allocate(self % index_mpas_scalar_to_constituent(self % number_of_constituents), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate index_mpas_scalar_to_constituent', subname, __LINE__) + end if + + self % index_mpas_scalar_to_constituent(:) = 0 + j = 1 + + ! Place water species first per MPAS requirements. + do i = 1, self % number_of_constituents + if (self % is_water_species(i)) then + self % index_mpas_scalar_to_constituent(j) = i + j = j + 1 + end if + end do + + index_water_start = 1 + index_water_end = count(self % is_water_species) + + ! Place non-water species second per MPAS requirements. + do i = 1, self % number_of_constituents + if (.not. self % is_water_species(i)) then + self % index_mpas_scalar_to_constituent(j) = i + j = j + 1 + end if + end do + + ! Create inverse index mapping between MPAS scalars and constituent names. For example, + ! Constituent index `i` corresponds to MPAS scalar index `index_constituent_to_mpas_scalar(i)`. + + allocate(self % index_constituent_to_mpas_scalar(self % number_of_constituents), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate index_constituent_to_mpas_scalar', subname, __LINE__) + end if + + self % index_constituent_to_mpas_scalar(:) = 0 + + do i = 1, self % number_of_constituents + self % index_constituent_to_mpas_scalar(self % index_mpas_scalar_to_constituent(i)) = i + end do + + ! Set the index of `qv` in terms of MPAS scalars. + index_qv = self % index_constituent_to_mpas_scalar(index_qv) + + ! Print information about constituents. + do i = 1, self % number_of_constituents + call self % debug_print('Constituent index ' // stringify([i])) + call self % debug_print(' Constituent name: ' // & + stringify([self % constituent_name(i)])) + call self % debug_print(' Is water species: ' // & + stringify([self % is_water_species(i)])) + call self % debug_print(' Index mapping from constituent to MPAS scalar: ' // & + stringify([i]) // ' -> ' // stringify([self % index_constituent_to_mpas_scalar(i)])) + end do + + ! Define "scalars" for MPAS. + + call self % get_pool_pointer(mpas_pool, 'state') + + call mpas_pool_add_dimension(mpas_pool, 'index_qv', index_qv) + call mpas_pool_add_dimension(mpas_pool, 'moist_start', index_water_start) + call mpas_pool_add_dimension(mpas_pool, 'moist_end', index_water_end) + + ! MPAS `state` pool has two time levels. + time_level = 2 + + do i = 1, time_level + call mpas_pool_get_field(mpas_pool, 'scalars', field_3d_real, timelevel=i) + + if (.not. associated(field_3d_real)) then + call self % model_error('Failed to find variable "scalars"', subname, __LINE__) + end if + + do j = 1, self % number_of_constituents + field_3d_real % constituentnames(j) = & + trim(adjustl(self % constituent_name(self % index_mpas_scalar_to_constituent(j)))) + + ! Print information about MPAS scalars. Only do it once. + if (i == 1) then + call self % debug_print('MPAS scalar index ' // stringify([j])) + call self % debug_print(' MPAS scalar name: ' // & + stringify([field_3d_real % constituentnames(j)])) + call self % debug_print(' Is water species: ' // & + stringify([self % is_water_species(self % index_mpas_scalar_to_constituent(j))])) + call self % debug_print(' Index mapping from MPAS scalar to constituent: ' // & + stringify([j]) // ' -> ' // stringify([self % index_mpas_scalar_to_constituent(j)])) + end if + end do + + nullify(field_3d_real) + end do + + nullify(mpas_pool) + + ! Define "scalars_tend" for MPAS. + + call self % get_pool_pointer(mpas_pool, 'tend') + + call mpas_pool_add_dimension(mpas_pool, 'index_qv', index_qv) + call mpas_pool_add_dimension(mpas_pool, 'moist_start', index_water_start) + call mpas_pool_add_dimension(mpas_pool, 'moist_end', index_water_end) + + ! MPAS `tend` pool only has one time level. + time_level = 1 + + do i = 1, time_level + call mpas_pool_get_field(mpas_pool, 'scalars_tend', field_3d_real, timelevel=i) + + if (.not. associated(field_3d_real)) then + call self % model_error('Failed to find variable "scalars_tend"', subname, __LINE__) + end if + + do j = 1, self % number_of_constituents + field_3d_real % constituentnames(j) = & + 'tendency_of_' // trim(adjustl(self % constituent_name(self % index_mpas_scalar_to_constituent(j)))) + end do + + nullify(field_3d_real) + end do + + nullify(mpas_pool) + + ! For consistency, also add dimension variables to MPAS `dimension` pool. + + call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'index_qv', index_qv) + call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'moist_start', index_water_start) + call mpas_pool_add_dimension(self % domain_ptr % blocklist % dimensions, 'moist_end', index_water_end) + + call self % debug_print('index_qv = ' // stringify([index_qv])) + call self % debug_print('moist_start = ' // stringify([index_water_start])) + call self % debug_print('moist_end = ' // stringify([index_water_end])) + + call self % debug_print(subname // ' completed') + end subroutine dyn_mpas_define_scalar + !------------------------------------------------------------------------------- ! subroutine dyn_mpas_read_write_stream ! @@ -1609,6 +1902,611 @@ subroutine dyn_mpas_exchange_halo(self, field_name) call self % debug_print(subname // ' completed') end subroutine dyn_mpas_exchange_halo + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_compute_unit_vector + ! + !> \brief Computes local east, north and edge-normal unit vectors + !> \author Michael Duda + !> \date 15 January 2020 + !> \details + !> This subroutine computes the local east and north unit vectors at all cells, + !> storing the results in MPAS `mesh` pool as `east` and `north`, respectively. + !> It also computes the edge-normal unit vectors at all edges by calling + !> `mpas_initialize_vectors`. Before calling this subroutine, MPAS `mesh` pool + !> must contain `latCell` and `lonCell` that are valid for all cells (not just + !> solve cells), plus any additional variables that are required by + !> `mpas_initialize_vectors`. + !> For stand-alone MPAS, the whole deal is handled by `init_dirs_forphys` + !> during physics initialization. However, MPAS as a dynamical core does + !> not have physics, hence this subroutine. + !> \addenda + !> Ported and refactored for CAM-SIMA. (KCW, 2024-04-23) + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_compute_unit_vector(self) + class(mpas_dynamical_core_type), intent(in) :: self + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_unit_vector' + integer :: i + integer, pointer :: ncells => null() + real(rkind), pointer :: latcell(:) => null() + real(rkind), pointer :: loncell(:) => null() + real(rkind), pointer :: east(:, :) => null() + real(rkind), pointer :: north(:, :) => null() + type(mpas_pool_type), pointer :: mpas_pool => null() + + call self % debug_print(subname // ' entered') + + ! Input. + call self % get_variable_pointer(ncells, 'dim', 'nCells') + call self % get_variable_pointer(latcell, 'mesh', 'latCell') + call self % get_variable_pointer(loncell, 'mesh', 'lonCell') + + ! Output. + call self % get_variable_pointer(east, 'mesh', 'east') + call self % get_variable_pointer(north, 'mesh', 'north') + + do i = 1, ncells + east(1, i) = -sin(loncell(i)) + east(2, i) = cos(loncell(i)) + east(3, i) = 0.0_rkind + ! `r3_normalize` has been inlined below. + east(1:3, i) = east(1:3, i) / sqrt(sum(east(1:3, i) * east(1:3, i))) + + north(1, i) = -cos(loncell(i)) * sin(latcell(i)) + north(2, i) = -sin(loncell(i)) * sin(latcell(i)) + north(3, i) = cos(latcell(i)) + ! `r3_normalize` has been inlined below. + north(1:3, i) = north(1:3, i) / sqrt(sum(north(1:3, i) * north(1:3, i))) + end do + + nullify(ncells) + nullify(latcell) + nullify(loncell) + + nullify(east) + nullify(north) + + call self % get_pool_pointer(mpas_pool, 'mesh') + call mpas_initialize_vectors(mpas_pool) + + nullify(mpas_pool) + + call self % debug_print(subname // ' completed') + end subroutine dyn_mpas_compute_unit_vector + + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_compute_edge_wind + ! + !> \brief Computes the edge-normal wind vectors at edge points + !> \author Michael Duda + !> \date 16 January 2020 + !> \details + !> This subroutine computes the edge-normal wind vectors at edge points + !> (i.e., `u` in MPAS `state` pool) from wind components at cell points + !> (i.e., `uReconstruct{Zonal,Meridional}` in MPAS `diag` pool). In MPAS, the + !> former are PROGNOSTIC variables, while the latter are DIAGNOSTIC variables + !> that are "reconstructed" from the former. This subroutine is essentially the + !> inverse function of that reconstruction. The purpose is to provide an + !> alternative way for MPAS to initialize from zonal and meridional wind + !> components at cell points. + !> \addenda + !> Ported and refactored for CAM-SIMA. (KCW, 2024-05-08) + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_compute_edge_wind(self) + class(mpas_dynamical_core_type), intent(in) :: self + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_edge_wind' + integer :: cell1, cell2, i + integer, pointer :: cellsonedge(:, :) => null() + integer, pointer :: nedges => null() + real(rkind), pointer :: east(:, :) => null() + real(rkind), pointer :: edgenormalvectors(:, :) => null() + real(rkind), pointer :: north(:, :) => null() + real(rkind), pointer :: ucellmeridional(:, :) => null() + real(rkind), pointer :: ucellzonal(:, :) => null() + real(rkind), pointer :: uedge(:, :) => null() + + call self % debug_print(subname // ' entered') + + ! Make sure halo layers are up-to-date before computation. + call self % exchange_halo('uReconstructZonal') + call self % exchange_halo('uReconstructMeridional') + + ! Input. + call self % get_variable_pointer(nedges, 'dim', 'nEdges') + + call self % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal') + call self % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional') + + call self % get_variable_pointer(cellsonedge, 'mesh', 'cellsOnEdge') + call self % get_variable_pointer(east, 'mesh', 'east') + call self % get_variable_pointer(edgenormalvectors, 'mesh', 'edgeNormalVectors') + call self % get_variable_pointer(north, 'mesh', 'north') + + ! Output. + call self % get_variable_pointer(uedge, 'state', 'u', time_level=1) + + do i = 1, nedges + cell1 = cellsonedge(1, i) + cell2 = cellsonedge(2, i) + + uedge(:, i) = ucellzonal(:, cell1) * 0.5_rkind * (edgenormalvectors(1, i) * east(1, cell1) + & + edgenormalvectors(2, i) * east(2, cell1) + & + edgenormalvectors(3, i) * east(3, cell1)) + & + ucellmeridional(:, cell1) * 0.5_rkind * (edgenormalvectors(1, i) * north(1, cell1) + & + edgenormalvectors(2, i) * north(2, cell1) + & + edgenormalvectors(3, i) * north(3, cell1)) + & + ucellzonal(:, cell2) * 0.5_rkind * (edgenormalvectors(1, i) * east(1, cell2) + & + edgenormalvectors(2, i) * east(2, cell2) + & + edgenormalvectors(3, i) * east(3, cell2)) + & + ucellmeridional(:, cell2) * 0.5_rkind * (edgenormalvectors(1, i) * north(1, cell2) + & + edgenormalvectors(2, i) * north(2, cell2) + & + edgenormalvectors(3, i) * north(3, cell2)) + end do + + nullify(nedges) + + nullify(ucellzonal) + nullify(ucellmeridional) + + nullify(cellsonedge) + nullify(east) + nullify(edgenormalvectors) + nullify(north) + + nullify(uedge) + + ! Make sure halo layers are up-to-date after computation. + call self % exchange_halo('u') + + call self % debug_print(subname // ' completed') + end subroutine dyn_mpas_compute_edge_wind + + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_init_phase4 + ! + !> \brief Tracks `atm_core_init` to finish MPAS dynamical core initialization + !> \author Michael Duda + !> \date 29 February 2020 + !> \details + !> This subroutine completes MPAS dynamical core initialization. + !> Essentially, it closely follows what is done in `atm_core_init`, but without + !> any calls to MPAS diagnostics manager or MPAS stream manager. + !> \addenda + !> Ported and refactored for CAM-SIMA. (KCW, 2024-05-25) + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_init_phase4(self, coupling_time_interval) + 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. + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase4' + character(strkind) :: date_time + character(strkind), pointer :: initial_time_1_pointer => null(), initial_time_2_pointer => null() + character(strkind), pointer :: xtime_pointer => null() + integer :: ierr + integer, pointer :: nvertlevels => null(), maxedges => null(), maxedges2 => null(), num_scalars => null() + logical, pointer :: config_do_restart => null() + real(rkind), pointer :: config_dt => null() + type(field0dreal), pointer :: field_0d_real => null() + type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_time_type) :: mpas_time + + call self % debug_print(subname // ' entered') + + if (real(coupling_time_interval, rkind) < 1.0_rkind) then + call self % model_error('Invalid coupling time interval ' // stringify([real(coupling_time_interval, rkind)]), & + subname, __LINE__) + end if + + call self % get_variable_pointer(config_dt, 'cfg', 'config_dt') + + if (config_dt < 1.0_rkind) then + call self % model_error('Invalid time step ' // stringify([config_dt]), & + subname, __LINE__) + end if + + ! `config_dt` in MPAS is a floating-point number. Testing floating-point numbers for divisibility is not trivial and + ! should be done carefully. + if (.not. almost_divisible(real(coupling_time_interval, rkind), config_dt)) then + call self % model_error('Coupling time interval ' // stringify([real(coupling_time_interval, rkind)]) // & + ' must be divisible by time step ' // stringify([config_dt]), subname, __LINE__) + end if + + self % coupling_time_interval = coupling_time_interval + + call self % debug_print('Coupling time interval is ' // stringify([real(self % coupling_time_interval, rkind)]) // & + ' seconds') + call self % debug_print('Time step is ' // stringify([config_dt]) // ' seconds') + + nullify(config_dt) + + ! Compute derived constants. + call mpas_constants_compute_derived() + + ! Set up OpenMP threading. + call self % debug_print('Setting up OpenMP threading') + + call mpas_atm_threading_init(self % domain_ptr % blocklist, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('OpenMP threading setup failed for core ' // trim(self % domain_ptr % core % corename), & + subname, __LINE__) + end if + + ! Set up inner dimensions used by arrays in optimized dynamics subroutines. + call self % debug_print('Setting up dimensions') + + call self % get_variable_pointer(nvertlevels, 'dim', 'nVertLevels') + call self % get_variable_pointer(maxedges, 'dim', 'maxEdges') + call self % get_variable_pointer(maxedges2, 'dim', 'maxEdges2') + call self % get_variable_pointer(num_scalars, 'dim', 'num_scalars') + + call mpas_atm_set_dims(nvertlevels, maxedges, maxedges2, num_scalars) + + nullify(nvertlevels) + nullify(maxedges) + nullify(maxedges2) + nullify(num_scalars) + + ! Build halo exchange groups and set the `exchange_halo_group` procedure pointer, which is used to + ! exchange the halo layers of all fields in the named group. + nullify(exchange_halo_group) + + call atm_build_halo_groups(self % domain_ptr, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to build halo exchange groups', subname, __LINE__) + end if + + if (.not. associated(exchange_halo_group)) then + call self % model_error('Failed to build halo exchange groups', subname, __LINE__) + end if + + ! Variables in MPAS `state` pool have more than one time level. Copy the values from the first time level of + ! such variables into all subsequent time levels to initialize them. + call self % get_variable_pointer(config_do_restart, 'cfg', 'config_do_restart') + + if (.not. config_do_restart) then + ! Run type is initial run. + call self % debug_print('Initializing time levels') + + call self % get_pool_pointer(mpas_pool, 'state') + + call mpas_pool_initialize_time_levels(mpas_pool) + + nullify(mpas_pool) + end if + + nullify(config_do_restart) + + call exchange_halo_group(self % domain_ptr, 'initialization:u', ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to exchange halo layers for group "initialization:u"', subname, __LINE__) + end if + + ! Initialize atmospheric variables (e.g., momentum, thermodynamic... variables in governing equations) + ! as well as various aspects of time in MPAS. + + call self % debug_print('Initializing atmospheric variables') + + ! Controlled by `config_start_time` in namelist. + mpas_time = mpas_get_clock_time(self % domain_ptr % clock, mpas_start_time, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to get time for "mpas_start_time"', subname, __LINE__) + end if + + call mpas_get_time(mpas_time, datetimestring=date_time, ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to get time for "mpas_start_time"', subname, __LINE__) + end if + + ! Controlled by `config_dt` in namelist. + call self % get_pool_pointer(mpas_pool, 'mesh') + call self % get_variable_pointer(config_dt, 'cfg', 'config_dt') + + call atm_mpas_init_block(self % domain_ptr % dminfo, self % domain_ptr % streammanager, self % domain_ptr % blocklist, & + mpas_pool, config_dt) + + nullify(mpas_pool) + nullify(config_dt) + + call self % get_variable_pointer(xtime_pointer, 'state', 'xtime', time_level=1) + + xtime_pointer = date_time + + nullify(xtime_pointer) + + ! Initialize `initial_time` in the second time level. We need to do this manually because initial states + ! are read into time level 1, and if we write anything from time level 2, `initial_time` will be invalid. + call self % get_variable_pointer(initial_time_1_pointer, 'state', 'initial_time', time_level=1) + call self % get_variable_pointer(initial_time_2_pointer, 'state', 'initial_time', time_level=2) + + initial_time_2_pointer = initial_time_1_pointer + + ! Set time units to CF-compliant "seconds since ". + call self % get_pool_pointer(mpas_pool, 'state') + + call mpas_pool_get_field(mpas_pool, 'Time', field_0d_real, timelevel=1) + + if (.not. associated(field_0d_real)) then + call self % model_error('Failed to find variable "Time"', subname, __LINE__) + end if + + call mpas_modify_att(field_0d_real % attlists(1) % attlist, 'units', & + 'seconds since ' // mpas_string_replace(initial_time_1_pointer, '_', ' '), ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to set time units', subname, __LINE__) + end if + + nullify(initial_time_1_pointer) + nullify(initial_time_2_pointer) + nullify(mpas_pool) + nullify(field_0d_real) + + call exchange_halo_group(self % domain_ptr, 'initialization:pv_edge,ru,rw', ierr=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to exchange halo layers for group "initialization:pv_edge,ru,rw"', subname, __LINE__) + end if + + call self % debug_print('Initializing dynamics') + + ! Prepare dynamics for time integration. + call mpas_atm_dynamics_init(self % domain_ptr) + + call self % debug_print(subname // ' completed') + + call self % debug_print('Successful initialization of MPAS dynamical core') + contains + !> Test if `a` is divisible by `b`, where `a` and `b` are both reals. + !> (KCW, 2024-05-25) + pure function almost_divisible(a, b) + real(rkind), intent(in) :: a, b + logical :: almost_divisible + + real(rkind) :: error_tolerance + + error_tolerance = epsilon(1.0_rkind) * max(abs(a), abs(b)) + + if (almost_equal(mod(abs(a), abs(b)), 0.0_rkind, absolute_tolerance=error_tolerance) .or. & + almost_equal(mod(abs(a), abs(b)), abs(b), absolute_tolerance=error_tolerance)) then + almost_divisible = .true. + + return + end if + + almost_divisible = .false. + end function almost_divisible + + !> Test `a` and `b` for approximate equality, where `a` and `b` are both reals. + !> (KCW, 2024-05-25) + pure function almost_equal(a, b, absolute_tolerance, relative_tolerance) + real(rkind), intent(in) :: a, b + real(rkind), optional, intent(in) :: absolute_tolerance, relative_tolerance + logical :: almost_equal + + real(rkind) :: error_tolerance + + if (present(relative_tolerance)) then + error_tolerance = relative_tolerance * max(abs(a), abs(b)) + else + error_tolerance = epsilon(1.0_rkind) * max(abs(a), abs(b)) + end if + + if (present(absolute_tolerance)) then + error_tolerance = max(absolute_tolerance, error_tolerance) + end if + + if (abs(a - b) <= error_tolerance) then + almost_equal = .true. + + return + end if + + almost_equal = .false. + end function almost_equal + end subroutine dyn_mpas_init_phase4 + + !------------------------------------------------------------------------------- + ! function dyn_mpas_get_constituent_name + ! + !> \brief Query constituent name by its index + !> \author Kuan-Chih Wang + !> \date 2024-05-16 + !> \details + !> This function returns the constituent name that corresponds to the given + !> constituent index. In case of errors, an empty character string is produced. + ! + !------------------------------------------------------------------------------- + pure function dyn_mpas_get_constituent_name(self, constituent_index) result(constituent_name) + class(mpas_dynamical_core_type), intent(in) :: self + integer, intent(in) :: constituent_index + + character(:), allocatable :: constituent_name + + ! Catch segmentation fault. + if (.not. allocated(self % constituent_name)) then + constituent_name = '' + + return + end if + + if (constituent_index < lbound(self % constituent_name, 1) .or. & + constituent_index > ubound(self % constituent_name, 1)) then + constituent_name = '' + + return + end if + + constituent_name = trim(adjustl(self % constituent_name(constituent_index))) + end function dyn_mpas_get_constituent_name + + !------------------------------------------------------------------------------- + ! function dyn_mpas_get_constituent_index + ! + !> \brief Query constituent index by its name + !> \author Kuan-Chih Wang + !> \date 2024-05-16 + !> \details + !> This function returns the constituent index that corresponds to the given + !> constituent name. In case of errors, zero is produced. + ! + !------------------------------------------------------------------------------- + pure function dyn_mpas_get_constituent_index(self, constituent_name) result(constituent_index) + class(mpas_dynamical_core_type), intent(in) :: self + character(*), intent(in) :: constituent_name + + integer :: i + integer :: constituent_index + + ! Catch segmentation fault. + if (.not. allocated(self % constituent_name)) then + constituent_index = 0 + + return + end if + + do i = 1, self % number_of_constituents + if (trim(adjustl(constituent_name)) == trim(adjustl(self % constituent_name(i)))) then + constituent_index = i + + return + end if + end do + + constituent_index = 0 + end function dyn_mpas_get_constituent_index + + !------------------------------------------------------------------------------- + ! function dyn_mpas_map_mpas_scalar_index + ! + !> \brief Map MPAS scalar index from constituent index + !> \author Kuan-Chih Wang + !> \date 2024-05-16 + !> \details + !> This function returns the MPAS scalar index that corresponds to the given + !> constituent index. In case of errors, zero is produced. + ! + !------------------------------------------------------------------------------- + pure function dyn_mpas_map_mpas_scalar_index(self, constituent_index) result(mpas_scalar_index) + class(mpas_dynamical_core_type), intent(in) :: self + integer, intent(in) :: constituent_index + + integer :: mpas_scalar_index + + ! Catch segmentation fault. + if (.not. allocated(self % index_constituent_to_mpas_scalar)) then + mpas_scalar_index = 0 + + return + end if + + if (constituent_index < lbound(self % index_constituent_to_mpas_scalar, 1) .or. & + constituent_index > ubound(self % index_constituent_to_mpas_scalar, 1)) then + mpas_scalar_index = 0 + + return + end if + + mpas_scalar_index = self % index_constituent_to_mpas_scalar(constituent_index) + end function dyn_mpas_map_mpas_scalar_index + + !------------------------------------------------------------------------------- + ! function dyn_mpas_map_constituent_index + ! + !> \brief Map constituent index from MPAS scalar index + !> \author Kuan-Chih Wang + !> \date 2024-05-16 + !> \details + !> This function returns the constituent index that corresponds to the given + !> MPAS scalar index. In case of errors, zero is produced. + ! + !------------------------------------------------------------------------------- + pure function dyn_mpas_map_constituent_index(self, mpas_scalar_index) result(constituent_index) + class(mpas_dynamical_core_type), intent(in) :: self + integer, intent(in) :: mpas_scalar_index + + integer :: constituent_index + + ! Catch segmentation fault. + if (.not. allocated(self % index_mpas_scalar_to_constituent)) then + constituent_index = 0 + + return + end if + + if (mpas_scalar_index < lbound(self % index_mpas_scalar_to_constituent, 1) .or. & + mpas_scalar_index > ubound(self % index_mpas_scalar_to_constituent, 1)) then + constituent_index = 0 + + return + end if + + constituent_index = self % index_mpas_scalar_to_constituent(mpas_scalar_index) + end function dyn_mpas_map_constituent_index + + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_get_local_mesh_dimension + ! + !> \brief Returns local mesh dimensions + !> \author Kuan-Chih Wang + !> \date 2024-05-09 + !> \details + !> This subroutine returns local mesh dimensions, including: + !> * Numbers of local mesh cells, edges, vertices and vertical levels + !> on each individual task, both with/without halo points. + ! + !------------------------------------------------------------------------------- + subroutine dyn_mpas_get_local_mesh_dimension(self, & + ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels) + class(mpas_dynamical_core_type), intent(in) :: self + integer, intent(out) :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_local_mesh_dimension' + integer, pointer :: ncells_pointer => null() + integer, pointer :: ncellssolve_pointer => null() + integer, pointer :: nedges_pointer => null() + integer, pointer :: nedgessolve_pointer => null() + integer, pointer :: nvertices_pointer => null() + integer, pointer :: nverticessolve_pointer => null() + integer, pointer :: nvertlevels_pointer => null() + + call self % get_variable_pointer(ncells_pointer, 'dim', 'nCells') + call self % get_variable_pointer(ncellssolve_pointer, 'dim', 'nCellsSolve') + call self % get_variable_pointer(nedges_pointer, 'dim', 'nEdges') + call self % get_variable_pointer(nedgessolve_pointer, 'dim', 'nEdgesSolve') + call self % get_variable_pointer(nvertices_pointer, 'dim', 'nVertices') + call self % get_variable_pointer(nverticessolve_pointer, 'dim', 'nVerticesSolve') + call self % get_variable_pointer(nvertlevels_pointer, 'dim', 'nVertLevels') + + ncells = ncells_pointer ! Number of cells, including halo cells. + ncells_solve = ncellssolve_pointer ! Number of cells, excluding halo cells. + nedges = nedges_pointer ! Number of edges, including halo edges. + nedges_solve = nedgessolve_pointer ! Number of edges, excluding halo edges. + nvertices = nvertices_pointer ! Number of vertices, including halo vertices. + nvertices_solve = nverticessolve_pointer ! Number of vertices, excluding halo vertices. + + ! Vertical levels are not decomposed. + ! All tasks have the same number of vertical levels. + nvertlevels = nvertlevels_pointer + + nullify(ncells_pointer) + nullify(ncellssolve_pointer) + nullify(nedges_pointer) + nullify(nedgessolve_pointer) + nullify(nvertices_pointer) + nullify(nverticessolve_pointer) + nullify(nvertlevels_pointer) + end subroutine dyn_mpas_get_local_mesh_dimension + !------------------------------------------------------------------------------- ! subroutine dyn_mpas_get_global_mesh_dimension ! From 94355357a66b0f0b5912587950eab2ac62c84aad Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Sun, 2 Jun 2024 21:41:05 -0600 Subject: [PATCH 05/31] Reimplement `parse_stream_name` for greater flexibility --- .../mpas/driver/dyn_mpas_subdriver.F90 | 65 ------------------- 1 file changed, 65 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 09a5ff57..584ce0f4 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -1593,71 +1593,6 @@ subroutine add_stream_attribute_1d(attribute_name, attribute_value) end subroutine add_stream_attribute_1d end subroutine dyn_mpas_init_stream_with_pool - !> Parse one or more stream names and return variable information contained in those streams as a list of `var_info_type`. - !> Multiple stream names should be separated by `+` (i.e., a plus). - !> Duplicate variable names in the resulting list are discarded. - !> (KCW, 2024-03-15) - pure recursive function parse_stream_name(stream_name) result(var_info_list) - character(*), intent(in) :: stream_name - type(var_info_type), allocatable :: var_info_list(:) - - character(64), allocatable :: var_name_list(:) - integer :: i, n, offset - type(var_info_type), allocatable :: var_info_list_append(:) - - allocate(var_info_list(0)) - - n = len(stream_name) - offset = 0 - - if (offset + 1 > n) then - return - end if - - i = index(stream_name(offset + 1:), '+') - - do while (i > 0) - if (i > 1) then - var_info_list_append = parse_stream_name(stream_name(offset + 1:offset + i - 1)) - var_info_list = [var_info_list, var_info_list_append] - - deallocate(var_info_list_append) - end if - - offset = offset + i - - if (offset + 1 > n) then - exit - end if - - i = index(stream_name(offset + 1:), '+') - end do - - if (offset + 1 > n) then - return - end if - - select case (trim(adjustl(stream_name(offset + 1:)))) - case ('invariant') - allocate(var_info_list_append, source=invariant_var_info_list) - case ('input') - allocate(var_info_list_append, source=input_var_info_list) - case ('restart') - allocate(var_info_list_append, source=restart_var_info_list) - case default - allocate(var_info_list_append(0)) - end select - - var_info_list = [var_info_list, var_info_list_append] - - ! Discard duplicate variable information by names. - var_name_list = var_info_list(:) % name - var_info_list = var_info_list(index_unique(var_name_list)) - - deallocate(var_info_list_append) - deallocate(var_name_list) - end function parse_stream_name - !> Return the index of unique elements in `array`, which can be any intrinsic data types, as an integer array. !> If `array` contains zero element or is of unsupported data types, an empty integer array is produced. !> (KCW, 2024-03-22) From be95a75e46bb33e26d411ad9118cf820f0214f45 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Sun, 2 Jun 2024 21:43:01 -0600 Subject: [PATCH 06/31] Reimplement `parse_stream_name` for greater flexibility It is now possible to construct an arbitrary list of variables to input/output data to/from MPAS dynamical core. This functionality may be overengineered, but it supports what CAM-SIMA needs. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 136 ++++++++++++++++++ 1 file changed, 136 insertions(+) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 584ce0f4..85fc31c9 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -1593,6 +1593,142 @@ subroutine add_stream_attribute_1d(attribute_name, attribute_value) end subroutine add_stream_attribute_1d end subroutine dyn_mpas_init_stream_with_pool + !> Parse a stream name, which consists of one or more stream name fragments, and return the corresponding variable information + !> as a list of `var_info_type`. Multiple stream name fragments should be separated by `+` (i.e., a plus, meaning "addition" + !> operation) or `-` (i.e., a minus, meaning "subtraction" operation). + !> A stream name fragment can be a predefined stream name (e.g., "invariant", "input", "restart") or a single variable name. + !> Duplicate variable names in the resulting list are discarded. + !> (KCW, 2024-06-01) + pure function parse_stream_name(stream_name) result(var_info_list) + character(*), intent(in) :: stream_name + type(var_info_type), allocatable :: var_info_list(:) + + character(*), parameter :: supported_stream_name_operator = '+-' + character(1) :: stream_name_operator + character(:), allocatable :: stream_name_fragment + character(len(invariant_var_info_list % name)), allocatable :: var_name_list(:) + integer :: i, j, n, offset + type(var_info_type), allocatable :: var_info_list_buffer(:) + + n = len_trim(stream_name) + + if (n == 0) then + ! Empty character string means empty list. + var_info_list = parse_stream_name_fragment('') + + return + end if + + i = scan(stream_name, supported_stream_name_operator) + + if (i == 0) then + ! No operators are present in the stream name. It is just a single stream name fragment. + stream_name_fragment = stream_name + var_info_list = parse_stream_name_fragment(stream_name_fragment) + + return + end if + + offset = 0 + var_info_list = parse_stream_name_fragment('') + + do while (.true.) + ! Extract operator from the stream name. + if (offset > 0) then + stream_name_operator = stream_name(offset:offset) + else + stream_name_operator = '+' + end if + + ! Extract stream name fragment from the stream name. + if (i > 1) then + stream_name_fragment = stream_name(offset + 1:offset + i - 1) + else + stream_name_fragment = '' + end if + + ! Process the stream name fragment according to the operator. + if (len_trim(stream_name_fragment) > 0) then + var_info_list_buffer = parse_stream_name_fragment(stream_name_fragment) + + select case (stream_name_operator) + case ('+') + var_info_list = [var_info_list, var_info_list_buffer] + case ('-') + do j = 1, size(var_info_list_buffer) + var_name_list = var_info_list % name + var_info_list = pack(var_info_list, var_name_list /= var_info_list_buffer(j) % name) + end do + case default + ! Do nothing for unknown operators. Should not happen at all. + end select + end if + + offset = offset + i + + ! Terminate loop when everything in the stream name has been processed. + if (offset + 1 > n) then + exit + end if + + i = scan(stream_name(offset + 1:), supported_stream_name_operator) + + ! Run the loop one last time for the remaining stream name fragment. + if (i == 0) then + i = n - offset + 1 + end if + end do + + ! Discard duplicate variable information by names. + var_name_list = var_info_list % name + var_info_list = var_info_list(index_unique(var_name_list)) + end function parse_stream_name + + !> Parse a stream name fragment and return the corresponding variable information as a list of `var_info_type`. + !> A stream name fragment can be a predefined stream name (e.g., "invariant", "input", "restart") or a single variable name. + !> (KCW, 2024-06-01) + pure function parse_stream_name_fragment(stream_name_fragment) result(var_info_list) + character(*), intent(in) :: stream_name_fragment + type(var_info_type), allocatable :: var_info_list(:) + + character(len(invariant_var_info_list % name)), allocatable :: var_name_list(:) + type(var_info_type), allocatable :: var_info_list_buffer(:) + + select case (trim(adjustl(stream_name_fragment))) + case ('') + allocate(var_info_list(0)) + case ('invariant') + allocate(var_info_list, source=invariant_var_info_list) + case ('input') + allocate(var_info_list, source=input_var_info_list) + case ('restart') + allocate(var_info_list, source=restart_var_info_list) + case default + allocate(var_info_list(0)) + + var_name_list = invariant_var_info_list % name + + if (any(var_name_list == trim(adjustl(stream_name_fragment)))) then + var_info_list_buffer = pack(invariant_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 = input_var_info_list % name + + if (any(var_name_list == trim(adjustl(stream_name_fragment)))) then + var_info_list_buffer = pack(input_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 = restart_var_info_list % name + + if (any(var_name_list == trim(adjustl(stream_name_fragment)))) then + 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 + end select + end function parse_stream_name_fragment + !> Return the index of unique elements in `array`, which can be any intrinsic data types, as an integer array. !> If `array` contains zero element or is of unsupported data types, an empty integer array is produced. !> (KCW, 2024-03-22) From 6c848ce222a4dc68f72ae235e69c96917f23a900 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 5 Jun 2024 14:56:24 -0600 Subject: [PATCH 07/31] Fix and enhance handling of variables during input There are two kinds of variables in MPAS: ordinary variables and variable arrays. The old implementation would always skip the latter entirely when encountering them during input. The new implementation correctly distinguishes between the two, and checks whether each variable is eligible to be input. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 434 +++++++++++++++++- 1 file changed, 416 insertions(+), 18 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 85fc31c9..046ba6fe 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -11,8 +11,9 @@ module dyn_mpas_subdriver ! Modules from external libraries. use mpi, only: mpi_comm_null, mpi_comm_rank, mpi_success + use netcdf, only: nf90_char, nf90_int, nf90_float, nf90_double use pio, only: file_desc_t, iosystem_desc_t, pio_file_is_open, pio_iosystem_is_active, & - pio_inq_varid, pio_inq_varndims, pio_noerr + pio_inq_varid, pio_inq_varndims, pio_inq_vartype, pio_noerr ! Modules from MPAS. use atm_core, only: atm_mpas_init_block @@ -39,7 +40,7 @@ module dyn_mpas_subdriver use mpas_framework, only: mpas_framework_init_phase1, mpas_framework_init_phase2 use mpas_io_streams, only: mpas_createstream, mpas_closestream, mpas_streamaddfield, & mpas_readstream, mpas_writestream, mpas_writestreamatt - use mpas_kind_types, only: rkind, strkind + use mpas_kind_types, only: rkind, r4kind, r8kind, strkind 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, & @@ -107,6 +108,7 @@ end subroutine model_error_if procedure, pass, public :: define_scalar => dyn_mpas_define_scalar procedure, pass, public :: read_write_stream => dyn_mpas_read_write_stream procedure, pass :: init_stream_with_pool => dyn_mpas_init_stream_with_pool + procedure, pass :: check_variable_status => dyn_mpas_check_variable_status procedure, pass, public :: exchange_halo => dyn_mpas_exchange_halo procedure, pass, public :: compute_unit_vector => dyn_mpas_compute_unit_vector procedure, pass, public :: compute_edge_wind => dyn_mpas_compute_edge_wind @@ -1223,7 +1225,9 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_stream_with_pool' character(strkind) :: stream_filename - integer :: i, ierr, ndims, stream_format, varid + integer :: i, ierr, stream_format + logical, allocatable :: var_is_present(:) + logical, allocatable :: var_is_tkr_compatible(:) type(field0dchar), pointer :: field_0d_char => null() type(field1dchar), pointer :: field_1d_char => null() type(field0dinteger), pointer :: field_0d_integer => null() @@ -1299,28 +1303,23 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file stringify([var_info_list(i) % rank])) if (trim(adjustl(stream_mode)) == 'r' .or. trim(adjustl(stream_mode)) == 'read') then - ! Check if "" is present. - ierr = pio_inq_varid(pio_file, trim(adjustl(var_info_list(i) % name)), varid) + call self % check_variable_status(var_is_present, var_is_tkr_compatible, pio_file, var_info_list(i)) ! Do not hard crash the model if a variable is missing and cannot be read. ! This can happen if users attempt to initialize/restart the model with data generated by ! older versions of MPAS. Print a debug message to let users decide if this is acceptable. - if (ierr /= pio_noerr) then - call self % debug_print('Skipping variable "' // trim(adjustl(var_info_list(i) % name)) // '"') + if (.not. any(var_is_present)) then + call self % debug_print('Skipping variable "' // trim(adjustl(var_info_list(i) % name)) // & + '" due to not present') cycle end if - ierr = pio_inq_varndims(pio_file, varid, ndims) + if (any(var_is_present .and. .not. var_is_tkr_compatible)) then + call self % debug_print('Skipping variable "' // trim(adjustl(var_info_list(i) % name)) // & + '" due to not TKR compatible') - if (ierr /= pio_noerr) then - call self % model_error('Failed to inquire variable rank for "' // trim(adjustl(var_info_list(i) % name)) // & - '"', subname, __LINE__) - end if - - if (ndims /= var_info_list(i) % rank) then - call self % model_error('Variable rank mismatch for "' // trim(adjustl(var_info_list(i) % name)) // & - '"', subname, __LINE__) + cycle end if end if @@ -1506,8 +1505,6 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file end if end do - deallocate(var_info_list) - if (trim(adjustl(stream_mode)) == 'w' .or. trim(adjustl(stream_mode)) == 'write') then ! Add MPAS-specific attributes to stream. call self % debug_print('Adding attributes to stream') @@ -1807,6 +1804,407 @@ pure function index_unique(array) index_unique = pack([(i, i = 1, n)], mask_unique) end function index_unique + !------------------------------------------------------------------------------- + ! subroutine dyn_mpas_check_variable_status + ! + !> \brief Check and return variable status on the given file. + !> \author Kuan-Chih Wang + !> \date 2024-06-04 + !> \details + !> On the given file (i.e., `pio_file`), this subroutine checks whether the + !> given variable (i.e., `var_info`) is present, and whether it is "TKR" + !> compatible with what MPAS expects. "TKR" means type, kind and rank. + !> This subroutine can handle both ordinary variables and variable arrays. + !> They are indicated by the `var` and `var_array` elements, respectively, + !> in MPAS registry. For an ordinary variable, the checks are performed on + !> itself. Otherwise, for a variable array, the checks are performed on its + !> constituent parts instead. + !------------------------------------------------------------------------------- + subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compatible, pio_file, var_info) + class(mpas_dynamical_core_type), intent(in) :: self + logical, allocatable, intent(out) :: var_is_present(:) + logical, allocatable, intent(out) :: var_is_tkr_compatible(:) + type(file_desc_t), pointer, intent(in) :: pio_file + type(var_info_type), intent(in) :: var_info + + character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_check_variable_status' + character(strkind), allocatable :: var_name_list(:) + integer :: i, ierr, varid, varndims, vartype + type(field0dchar), pointer :: field_0d_char => null() + type(field1dchar), pointer :: field_1d_char => null() + type(field0dinteger), pointer :: field_0d_integer => null() + type(field1dinteger), pointer :: field_1d_integer => null() + type(field2dinteger), pointer :: field_2d_integer => null() + type(field3dinteger), pointer :: field_3d_integer => null() + type(field0dreal), pointer :: field_0d_real => null() + type(field1dreal), pointer :: field_1d_real => null() + type(field2dreal), pointer :: field_2d_real => null() + type(field3dreal), pointer :: field_3d_real => null() + type(field4dreal), pointer :: field_4d_real => null() + type(field5dreal), pointer :: field_5d_real => null() + + call self % debug_print(subname // ' entered') + + ! Extract a list of variable names to check on the file. + ! For an ordinary variable, this list just contains its name. + ! For a variable array, this list contains the names of its constituent parts. + select case (trim(adjustl(var_info % type))) + case ('character') + select case (var_info % rank) + case (0) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_0d_char, timelevel=1) + + if (.not. associated(field_0d_char)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_0d_char % isvararray .and. associated(field_0d_char % constituentnames)) then + allocate(var_name_list(size(field_0d_char % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_0d_char % constituentnames(:) + end if + + nullify(field_0d_char) + case (1) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_1d_char, timelevel=1) + + if (.not. associated(field_1d_char)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_1d_char % isvararray .and. associated(field_1d_char % constituentnames)) then + allocate(var_name_list(size(field_1d_char % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_1d_char % constituentnames(:) + end if + + nullify(field_1d_char) + case default + call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]), & + subname, __LINE__) + end select + case ('integer') + select case (var_info % rank) + case (0) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_0d_integer, timelevel=1) + + if (.not. associated(field_0d_integer)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_0d_integer % isvararray .and. associated(field_0d_integer % constituentnames)) then + allocate(var_name_list(size(field_0d_integer % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_0d_integer % constituentnames(:) + end if + + nullify(field_0d_integer) + case (1) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_1d_integer, timelevel=1) + + if (.not. associated(field_1d_integer)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_1d_integer % isvararray .and. associated(field_1d_integer % constituentnames)) then + allocate(var_name_list(size(field_1d_integer % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_1d_integer % constituentnames(:) + end if + + nullify(field_1d_integer) + case (2) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_2d_integer, timelevel=1) + + if (.not. associated(field_2d_integer)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_2d_integer % isvararray .and. associated(field_2d_integer % constituentnames)) then + allocate(var_name_list(size(field_2d_integer % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_2d_integer % constituentnames(:) + end if + + nullify(field_2d_integer) + case (3) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_3d_integer, timelevel=1) + + if (.not. associated(field_3d_integer)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_3d_integer % isvararray .and. associated(field_3d_integer % constituentnames)) then + allocate(var_name_list(size(field_3d_integer % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_3d_integer % constituentnames(:) + end if + + nullify(field_3d_integer) + case default + call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]), & + subname, __LINE__) + end select + case ('real') + select case (var_info % rank) + case (0) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_0d_real, timelevel=1) + + if (.not. associated(field_0d_real)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_0d_real % isvararray .and. associated(field_0d_real % constituentnames)) then + allocate(var_name_list(size(field_0d_real % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_0d_real % constituentnames(:) + end if + + nullify(field_0d_real) + case (1) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_1d_real, timelevel=1) + + if (.not. associated(field_1d_real)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_1d_real % isvararray .and. associated(field_1d_real % constituentnames)) then + allocate(var_name_list(size(field_1d_real % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_1d_real % constituentnames(:) + end if + + nullify(field_1d_real) + case (2) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_2d_real, timelevel=1) + + if (.not. associated(field_2d_real)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_2d_real % isvararray .and. associated(field_2d_real % constituentnames)) then + allocate(var_name_list(size(field_2d_real % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_2d_real % constituentnames(:) + end if + + nullify(field_2d_real) + case (3) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_3d_real, timelevel=1) + + if (.not. associated(field_3d_real)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_3d_real % isvararray .and. associated(field_3d_real % constituentnames)) then + allocate(var_name_list(size(field_3d_real % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_3d_real % constituentnames(:) + end if + + nullify(field_3d_real) + case (4) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_4d_real, timelevel=1) + + if (.not. associated(field_4d_real)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_4d_real % isvararray .and. associated(field_4d_real % constituentnames)) then + allocate(var_name_list(size(field_4d_real % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_4d_real % constituentnames(:) + end if + + nullify(field_4d_real) + case (5) + call mpas_pool_get_field(self % domain_ptr % blocklist % allfields, & + trim(adjustl(var_info % name)), field_5d_real, timelevel=1) + + if (.not. associated(field_5d_real)) then + call self % model_error('Failed to find variable "' // trim(adjustl(var_info % name)) // & + '"', subname, __LINE__) + end if + + if (field_5d_real % isvararray .and. associated(field_5d_real % constituentnames)) then + allocate(var_name_list(size(field_5d_real % constituentnames)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(:) = field_5d_real % constituentnames(:) + end if + + nullify(field_5d_real) + case default + call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]), & + subname, __LINE__) + end select + case default + call self % model_error('Unsupported variable type "' // trim(adjustl(var_info % type)) // & + '"', subname, __LINE__) + end select + + if (.not. allocated(var_name_list)) then + allocate(var_name_list(1), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_name_list', subname, __LINE__) + end if + + var_name_list(1) = var_info % name + end if + + allocate(var_is_present(size(var_name_list)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_is_present', subname, __LINE__) + end if + + var_is_present(:) = .false. + + allocate(var_is_tkr_compatible(size(var_name_list)), stat=ierr) + + if (ierr /= 0) then + call self % model_error('Failed to allocate var_is_tkr_compatible', subname, __LINE__) + end if + + var_is_tkr_compatible(:) = .false. + + if (.not. associated(pio_file)) then + return + end if + + if (.not. pio_file_is_open(pio_file)) then + return + end if + + do i = 1, size(var_name_list) + ! Check if the variable is present on the file. + ierr = pio_inq_varid(pio_file, trim(adjustl(var_name_list(i))), varid) + + if (ierr /= pio_noerr) then + cycle + end if + + var_is_present(i) = .true. + + ! Check if the variable is "TK"R compatible between MPAS and the file. + ierr = pio_inq_vartype(pio_file, varid, vartype) + + if (ierr /= pio_noerr) then + cycle + end if + + select case (trim(adjustl(var_info % type))) + case ('character') + if (vartype /= nf90_char) then + cycle + end if + case ('integer') + if (vartype /= nf90_int) then + cycle + end if + case ('real') + if (rkind == r4kind .and. vartype /= nf90_float) then + cycle + end if + + if (rkind == r8kind .and. vartype /= nf90_double) then + cycle + end if + case default + cycle + end select + + ! Check if the variable is TK"R" compatible between MPAS and the file. + ierr = pio_inq_varndims(pio_file, varid, varndims) + + if (ierr /= pio_noerr) then + cycle + end if + + if (varndims /= var_info % rank) then + cycle + end if + + var_is_tkr_compatible(i) = .true. + end do + + call self % debug_print('var_name_list = ' // stringify(var_name_list)) + call self % debug_print('var_is_present = ' // stringify(var_is_present)) + call self % debug_print('var_is_tkr_compatible = ' // stringify(var_is_tkr_compatible)) + + call self % debug_print(subname // ' completed') + end subroutine dyn_mpas_check_variable_status + !------------------------------------------------------------------------------- ! subroutine dyn_mpas_exchange_halo ! From c38a9b3aeb2dcb27848b7f8364ba7694d9a6b123 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 29 Apr 2024 16:52:25 -0600 Subject: [PATCH 08/31] Add an additional include path for building MPAS subdriver --- src/dynamics/mpas/Makefile.in.CESM | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamics/mpas/Makefile.in.CESM b/src/dynamics/mpas/Makefile.in.CESM index a0bff7ed..f4016f0e 100644 --- a/src/dynamics/mpas/Makefile.in.CESM +++ b/src/dynamics/mpas/Makefile.in.CESM @@ -130,5 +130,5 @@ externals: $(AUTOCLEAN_DEPS) subdrv: driver/dyn_mpas_subdriver.o %.o: %.F90 dycore frame ops - ( cd $( Date: Wed, 29 May 2024 11:57:37 -0600 Subject: [PATCH 09/31] Move all shared variables to `dyn_comp` due to circular dependency issue The instance of MPAS dynamical core resides in `dyn_comp`. `dyn_grid` depends on `dyn_comp`, but not vice versa. --- src/dynamics/mpas/dyn_comp.F90 | 12 ++++++++++++ src/dynamics/mpas/dyn_grid.F90 | 30 +++++++----------------------- 2 files changed, 19 insertions(+), 23 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 5b4af9cd..cb0dcc1a 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -32,6 +32,10 @@ module dyn_comp public :: dyn_debug_print public :: mpas_dynamical_core + public :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels + public :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max + public :: sphere_radius + public :: deg_to_rad, rad_to_deg type :: dyn_import_t end type dyn_import_t @@ -41,6 +45,14 @@ module dyn_comp !> The "instance/object" of MPAS dynamical core. 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_r8) :: sphere_radius + + real(kind_r8), parameter :: deg_to_rad = constant_pi / 180.0_kind_r8 ! Convert degrees to radians. + real(kind_r8), parameter :: rad_to_deg = 180.0_kind_r8 / constant_pi ! Convert radians to degrees. 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. diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 1effa78c..0074df22 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -7,7 +7,11 @@ module dyn_grid max_hcoordname_len use cam_initfiles, only: initial_file_get_id use cam_map_utils, only: kind_imap => imap - use dyn_comp, only: dyn_debug_print, mpas_dynamical_core + 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, & + deg_to_rad, rad_to_deg use dynconst, only: dynconst_init, pi use physics_column_type, only: kind_pcol, physics_column_t use physics_grid, only: phys_decomp, phys_grid_init @@ -40,14 +44,6 @@ module dyn_grid 'mpas_edge', & 'mpas_vertex' & ] - - real(kind_r8), parameter :: deg_to_rad = pi / 180.0_kind_r8 ! Convert degrees to radians. - real(kind_r8), parameter :: rad_to_deg = 180.0_kind_r8 / pi ! Convert radians to degrees. - - ! Local and global mesh dimensions. - integer :: ncells_solve, nedges_solve, nvertices_solve - integer :: ncells_global, nedges_global, nvertices_global, nvertlevels, ncells_max, nedges_max - real(kind_r8) :: sphere_radius contains !> Initialize various model grids (e.g., dynamics, physics) in terms of dynamics decomposition. !> Additionally, MPAS framework initialization and reading time-invariant (e.g., grid/mesh) variables @@ -57,9 +53,6 @@ module dyn_grid ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine model_grid_init() character(*), parameter :: subname = 'dyn_grid::model_grid_init' - integer, pointer :: ncellssolve => null() - integer, pointer :: nedgessolve => null() - integer, pointer :: nverticessolve => null() type(file_desc_t), pointer :: pio_file => null() ! Initialize mathematical and physical constants for dynamics. @@ -87,17 +80,8 @@ subroutine model_grid_init() ! 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_variable_pointer(ncellssolve, 'dim', 'nCellsSolve') - call mpas_dynamical_core % get_variable_pointer(nedgessolve, 'dim', 'nEdgesSolve') - call mpas_dynamical_core % get_variable_pointer(nverticessolve, 'dim', 'nVerticesSolve') - - ncells_solve = ncellssolve - nedges_solve = nedgessolve - nvertices_solve = nverticessolve - - nullify(ncellssolve) - nullify(nedgessolve) - nullify(nverticessolve) + 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, & From a0ac03b87863ab1d68b0fcaf77fadd0117096905 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 29 May 2024 12:01:48 -0600 Subject: [PATCH 10/31] Implement `dyn_init` --- src/dynamics/mpas/dyn_comp.F90 | 608 ++++++++++++++++++++++++++++++++- src/dynamics/mpas/dyn_grid.F90 | 9 +- 2 files changed, 603 insertions(+), 14 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index cb0dcc1a..00ff9825 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -2,22 +2,41 @@ module dyn_comp use dyn_mpas_subdriver, only: mpas_dynamical_core_type ! Modules from CAM. - use cam_abortutils, only: endrun + 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_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_g => gravit, constant_pi => pi + use inic_analytic, only: analytic_ic_active, dyn_set_inic_col + use physconst, only: constant_cpd => cpair, constant_p0 => pref, constant_rd => rair, constant_rv => rh2o 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_run_duration, timemgr_get_calendar_cf + 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 CESM Share. use shr_file_mod, only: shr_file_getunit - use shr_kind_mod, only: shr_kind_cs + 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 CCPP. + use phys_vars_init_check, only: mark_as_initialized, std_name_len + ! Modules from external libraries. - use pio, only: iosystem_desc_t + use pio, only: file_desc_t, iosystem_desc_t, pio_file_is_open implicit none @@ -37,9 +56,17 @@ module dyn_comp public :: sphere_radius public :: deg_to_rad, rad_to_deg + !> NOTE: + !> This derived type is not used by MPAS dynamical core. It exists only as a placeholder because CAM Control requires it. + !> Developers/Maintainers/Users who wish to interact with MPAS dynamical core may do so by using the "instance/object" + !> below. type :: dyn_import_t end type dyn_import_t + !> NOTE: + !> This derived type is not used by MPAS dynamical core. It exists only as a placeholder because CAM Control requires it. + !> Developers/Maintainers/Users who wish to interact with MPAS dynamical core may do so by using the "instance/object" + !> below. type :: dyn_export_t end type dyn_export_t @@ -95,7 +122,7 @@ subroutine dyn_readnl(namelist_path) character(*), intent(in) :: namelist_path character(*), parameter :: subname = 'dyn_comp::dyn_readnl' - character(shr_kind_cs) :: cam_calendar + character(kind_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. @@ -111,7 +138,6 @@ subroutine dyn_readnl(namelist_path) log_unit(2) = shr_file_getunit() ! Initialize MPAS framework with supplied MPI communicator group and log units. - ! See comment blocks in `src/dynamics/mpas/driver/dyn_mpas_subdriver.F90` for details. call mpas_dynamical_core % init_phase1(mpicom, endrun, iulog, log_unit) cam_calendar = timemgr_get_calendar_cf() @@ -126,14 +152,12 @@ subroutine dyn_readnl(namelist_path) run_duration(2:4) = sec_to_hour_min_sec(sec_since_midnight) ! Read MPAS-related namelist variables from `namelist_path`, e.g., `atm_in`. - ! See comment blocks in `src/dynamics/mpas/driver/dyn_mpas_subdriver.F90` for details. call mpas_dynamical_core % read_namelist(namelist_path, & cam_calendar, start_date_time, stop_date_time, run_duration, initial_run) pio_iosystem => shr_pio_getiosys(atm_id) ! Initialize MPAS framework with supplied PIO system descriptor. - ! See comment blocks in `src/dynamics/mpas/driver/dyn_mpas_subdriver.F90` for details. call mpas_dynamical_core % init_phase2(pio_iosystem) nullify(pio_iosystem) @@ -151,13 +175,575 @@ pure function sec_to_hour_min_sec(sec) result(hour_min_sec) end function sec_to_hour_min_sec end subroutine dyn_readnl + !> Initialize MPAS dynamical core by one of the following: + !> 1. Setting analytic initial condition; + !> 2. Reading initial condition from a file; + !> 3. Restarting from a file. + !> (KCW, 2024-05-28) + ! ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) - type(runtime_options), intent(in) :: cam_runtime_opts - type(dyn_import_t), intent(out) :: dyn_in - type(dyn_export_t), intent(out) :: 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 + + character(*), parameter :: subname = 'dyn_comp::dyn_init' + character(std_name_len), allocatable :: constituent_name(:) + integer :: coupling_time_interval + integer :: i + integer :: ierr + logical, allocatable :: is_water_species(:) + type(file_desc_t), pointer :: pio_init_file => null() + type(file_desc_t), pointer :: pio_topo_file => null() + + allocate(constituent_name(num_advected), stat=ierr) + call check_allocate(ierr, 'dyn_init', 'constituent_name', 'dyn_comp', __LINE__) + + allocate(is_water_species(num_advected), stat=ierr) + call check_allocate(ierr, 'dyn_init', 'is_water_species', 'dyn_comp', __LINE__) + + do i = 1, num_advected + constituent_name(i) = const_name(i) + is_water_species(i) = const_is_water_species(i) + end do + + ! Inform MPAS about constituent names and their corresponding waterness. + call mpas_dynamical_core % define_scalar(constituent_name, is_water_species) + + ! Provide mapping information between MPAS scalars and constituent names to CAM-SIMA. + do i = 1, thermodynamic_active_species_num + thermodynamic_active_species_idx_dycore(i) = & + mpas_dynamical_core % map_mpas_scalar_index(thermodynamic_active_species_idx(i)) + end do + + do i = 1, thermodynamic_active_species_liq_num + thermodynamic_active_species_liq_idx_dycore(i) = & + mpas_dynamical_core % map_mpas_scalar_index(thermodynamic_active_species_liq_idx(i)) + end do + + do i = 1, thermodynamic_active_species_ice_num + thermodynamic_active_species_ice_idx_dycore(i) = & + mpas_dynamical_core % map_mpas_scalar_index(thermodynamic_active_species_ice_idx(i)) + end do + + pio_init_file => initial_file_get_id() + pio_topo_file => topo_file_get_id() + + if (initial_run) then + ! Run type is initial run. + + call dyn_debug_print('Calling check_topography_data') + + call check_topography_data(pio_topo_file) + + if (analytic_ic_active()) then + call dyn_debug_print('Calling set_analytic_initial_condition') + + call set_analytic_initial_condition() + else + ! Namelist option that controls if constituents are to be read from the file. + if (readtrace) then + ! Read variables that belong to the "input" stream in MPAS. + call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input') + + ! TODO: + ! `cnst_init_default` or its replacement is not yet implemented in CAM-SIMA. + ! Default constituent initialization should be performed here + ! for constituents that fail to be read from the file. + else + ! Read variables that belong to the "input" stream in MPAS, excluding constituents. + call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input-scalars') + + ! TODO: + ! `cnst_init_default` or its replacement is not yet implemented in CAM-SIMA. + ! Default constituent initialization should be performed here + ! because constituents are not read from the file. + end if + end if + else + ! Run type is branch or restart run. + + ! Read variables that belong to the "input" and "restart" streams in MPAS. + call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input+restart') + end if + + call clean_iodesc_list() + call mark_variable_as_initialized(constituent_name) + + deallocate(constituent_name) + deallocate(is_water_species) + + nullify(pio_init_file) + nullify(pio_topo_file) + + ! This is the time interval for dynamics-physics coupling in CAM-SIMA. + ! Each time MPAS dynamical core is called to run, it will integrate with time for this specific interval, + ! then yield control back to the caller. + coupling_time_interval = get_step_size() + + ! Finish MPAS dynamical core initialization. After this point, MPAS dynamical core is ready for time integration. + call mpas_dynamical_core % init_phase4(coupling_time_interval) end subroutine dyn_init + !> Check for consistency in topography data. The presence of topography file is inferred from the `pio_file` pointer. + !> If topography file is used, check that the `PHIS` variable, which denotes surface geopotential, + !> is consistent with the surface geometric height in MPAS. + !> 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) + type(file_desc_t), pointer, intent(in) :: pio_file + + character(*), parameter :: subname = 'dyn_comp::check_topography_data' + integer :: ierr + logical :: success + 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(:, :) => null() ! From MPAS. Geometric height (meters) at layer interfaces. + + call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid') + + if (associated(pio_file)) then + call dyn_debug_print('Topography file is used') + + if (.not. pio_file_is_open(pio_file)) then + call endrun('Invalid PIO file descriptor', subname, __LINE__) + end if + + allocate(surface_geopotential(ncells_solve), stat=ierr) + call check_allocate(ierr, 'check_topography_data', 'surface_geopotential', 'dyn_comp', __LINE__) + + allocate(surface_geometric_height(ncells_solve), stat=ierr) + call check_allocate(ierr, 'check_topography_data', 'surface_geometric_height', 'dyn_comp', __LINE__) + + surface_geopotential(:) = 0.0_kind_r8 + surface_geometric_height(:) = 0.0_kind_r8 + + call cam_read_field('PHIS', pio_file, surface_geopotential, success, & + gridname='cam_cell', timelevel=1, log_output=(debug_output > debugout_none)) + + if (success) then + surface_geometric_height(:) = surface_geopotential(:) / constant_g + else + call endrun('Failed to find variable "PHIS"', subname, __LINE__) + end if + + ! 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 + call endrun('Surface geometric height in MPAS is not consistent with topography data', subname, __LINE__) + end if + + deallocate(surface_geopotential) + deallocate(surface_geometric_height) + else + 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 + call endrun('Surface geometric height in MPAS is not zero', subname, __LINE__) + end if + end if + + nullify(zgrid) + end subroutine check_topography_data + + !> Set analytic initial condition for MPAS. + !> (KCW, 2024-05-22) + subroutine set_analytic_initial_condition() + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition' + integer, allocatable :: global_grid_index(:) + real(kind_r8), allocatable :: buffer_1d_real(:), 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(:, :) => null() ! Geometric height (meters) at layer interfaces. + ! Dimension and vertical index orders follow MPAS convention. + + call init_shared_variable() + + call set_mpas_state_u() + call set_mpas_state_w() + call set_mpas_state_scalars() + call set_mpas_state_rho_theta() + call set_mpas_state_rho_base_theta_base() + + deallocate(global_grid_index) + deallocate(lat_rad) + deallocate(lon_rad) + deallocate(z_int) + + nullify(zgrid) + contains + !> Initialize variables that are shared and repeatedly used by the `set_mpas_state_*` internal subroutines. + !> (KCW, 2024-05-13) + subroutine init_shared_variable() + integer :: ierr + integer :: k + integer, pointer :: indextocellid(:) => null() + real(kind_r8), pointer :: lat_deg(:) => null(), lon_deg(:) => null() + + call dyn_debug_print('Preparing to set analytic initial condition') + + allocate(global_grid_index(ncells_solve), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'global_grid_index', 'dyn_comp', __LINE__) + + call mpas_dynamical_core % get_variable_pointer(indextocellid, 'mesh', 'indexToCellID') + + global_grid_index(:) = indextocellid(1:ncells_solve) + + nullify(indextocellid) + + allocate(lat_rad(ncells_solve), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'lat_rad', 'dyn_comp', __LINE__) + + allocate(lon_rad(ncells_solve), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'lon_rad', 'dyn_comp', __LINE__) + + ! "mpas_cell" is a registered grid name that is defined in `dyn_grid`. + lat_deg => cam_grid_get_latvals(cam_grid_id('mpas_cell')) + lon_deg => cam_grid_get_lonvals(cam_grid_id('mpas_cell')) + + if (.not. associated(lat_deg)) then + call endrun('Failed to find variable "lat_deg"', subname, __LINE__) + end if + + if (.not. associated(lon_deg)) then + call endrun('Failed to find variable "lon_deg"', subname, __LINE__) + end if + + lat_rad(:) = lat_deg(:) * deg_to_rad + lon_rad(:) = lon_deg(:) * deg_to_rad + + nullify(lat_deg) + nullify(lon_deg) + + allocate(z_int(ncells_solve, pverp), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'z_int', 'dyn_comp', __LINE__) + + call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid') + + ! Vertical index order is reversed between CAM-SIMA and MPAS. + do k = 1, pverp + z_int(:, k) = zgrid(pverp - k + 1, 1:ncells_solve) + end do + end subroutine init_shared_variable + + !> Set MPAS state `u` (i.e., horizontal velocity at edge interfaces). + !> (KCW, 2024-05-13) + subroutine set_mpas_state_u() + integer :: ierr + integer :: k + real(kind_r8), pointer :: ucellmeridional(:, :) => null() + real(kind_r8), pointer :: ucellzonal(:, :) => null() + + call dyn_debug_print('Setting MPAS state "u"') + + allocate(buffer_2d_real(ncells_solve, pver), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_2d_real', 'dyn_comp', __LINE__) + + call mpas_dynamical_core % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal') + call mpas_dynamical_core % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional') + + buffer_2d_real(:, :) = 0.0_kind_r8 + + call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, u=buffer_2d_real) + + ! Vertical index order is reversed between CAM-SIMA and MPAS. + do k = 1, pver + ucellzonal(k, 1:ncells_solve) = buffer_2d_real(:, pver - k + 1) + end do + + buffer_2d_real(:, :) = 0.0_kind_r8 + + call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, v=buffer_2d_real) + + ! Vertical index order is reversed between CAM-SIMA and MPAS. + do k = 1, pver + ucellmeridional(k, 1:ncells_solve) = buffer_2d_real(:, pver - k + 1) + end do + + deallocate(buffer_2d_real) + + nullify(ucellzonal) + nullify(ucellmeridional) + + call mpas_dynamical_core % compute_edge_wind() + end subroutine set_mpas_state_u + + !> Set MPAS state `w` (i.e., vertical velocity at cell interfaces). + !> (KCW, 2024-05-13) + subroutine set_mpas_state_w() + real(kind_r8), pointer :: w(:, :) => null() + + call dyn_debug_print('Setting MPAS state "w"') + + call mpas_dynamical_core % get_variable_pointer(w, 'state', 'w', time_level=1) + + w(:, 1:ncells_solve) = 0.0_kind_r8 + + nullify(w) + + ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. + call mpas_dynamical_core % exchange_halo('w') + end subroutine set_mpas_state_w + + !> Set MPAS state `scalars` (i.e., constituents). + !> (KCW, 2024-05-17) + subroutine set_mpas_state_scalars() + ! CCPP standard name of `qv`, which denotes water vapor mixing ratio. + character(*), parameter :: constituent_qv_standard_name = & + 'water_vapor_mixing_ratio_wrt_dry_air' + + integer :: i, k + integer :: ierr + integer, allocatable :: constituent_index(:) + integer, pointer :: index_qv => null() + real(kind_r8), pointer :: scalars(:, :, :) => null() + + call dyn_debug_print('Setting MPAS state "scalars"') + + allocate(buffer_3d_real(ncells_solve, pver, num_advected), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_3d_real', 'dyn_comp', __LINE__) + + allocate(constituent_index(num_advected), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'constituent_index', 'dyn_comp', __LINE__) + + call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv') + call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) + + buffer_3d_real(:, :, :) = 0.0_kind_r8 + constituent_index(:) = [(i, i = 1, num_advected)] + + call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, q=buffer_3d_real, & + m_cnst=constituent_index) + + ! `i` is indexing into `scalars`, so it is regarded as MPAS scalar index. + do i = 1, num_advected + ! Vertical index order is reversed between CAM-SIMA and MPAS. + do k = 1, pver + scalars(i, k, 1:ncells_solve) = & + buffer_3d_real(:, pver - k + 1, mpas_dynamical_core % map_constituent_index(i)) + end do + end do + + if (mpas_dynamical_core % get_constituent_name(mpas_dynamical_core % map_constituent_index(index_qv)) == & + constituent_qv_standard_name) then + ! The definition of `qv` matches exactly what MPAS wants. No conversion is needed. + call dyn_debug_print('No conversion is needed for water vapor mixing ratio') + else + ! The definition of `qv` actually represents specific humidity. Conversion is needed. + call dyn_debug_print('Conversion is needed and applied for water vapor mixing ratio') + + ! 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)) + end if + + deallocate(buffer_3d_real) + deallocate(constituent_index) + + nullify(index_qv) + nullify(scalars) + + ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. + call mpas_dynamical_core % exchange_halo('scalars') + 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() + integer :: i, k + integer :: ierr + integer, pointer :: index_qv => null() + real(kind_r8), allocatable :: p_mid_col(:) ! Pressure (Pa) at layer midpoints of each column. This is full pressure, + ! which also accounts for water vapor. + real(kind_r8), allocatable :: p_sfc(:) ! Pressure (Pa) at surface. This is full pressure, + ! which also accounts for water vapor. + real(kind_r8), allocatable :: qv_mid_col(:) ! Water vapor mixing ratio (kg/kg) at layer midpoints of each column. + real(kind_r8), allocatable :: t_mid(:, :) ! Temperature (K) at layer midpoints. + real(kind_r8), allocatable :: tm_mid_col(:) ! Modified "moist" temperature (K) at layer midpoints of each column. + real(kind_r8), pointer :: rho(:, :) => null() + real(kind_r8), pointer :: theta(:, :) => null() + real(kind_r8), pointer :: scalars(:, :, :) => null() + + call dyn_debug_print('Setting MPAS state "rho" and "theta"') + + allocate(buffer_1d_real(ncells_solve), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_1d_real', 'dyn_comp', __LINE__) + + allocate(p_sfc(ncells_solve), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'p_sfc', 'dyn_comp', __LINE__) + + buffer_1d_real(:) = 0.0_kind_r8 + + call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, ps=buffer_1d_real) + + p_sfc(:) = buffer_1d_real(:) + + deallocate(buffer_1d_real) + + allocate(buffer_2d_real(ncells_solve, pver), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_2d_real', 'dyn_comp', __LINE__) + + allocate(t_mid(pver, ncells_solve), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 't_mid', 'dyn_comp', __LINE__) + + buffer_2d_real(:, :) = 0.0_kind_r8 + + call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, t=buffer_2d_real) + + ! Vertical index order is reversed between CAM-SIMA and MPAS. + do k = 1, pver + t_mid(k, :) = buffer_2d_real(:, pver - k + 1) + end do + + deallocate(buffer_2d_real) + + allocate(p_mid_col(pver), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'p_mid_col', 'dyn_comp', __LINE__) + + allocate(qv_mid_col(pver), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'qv_mid_col', 'dyn_comp', __LINE__) + + allocate(tm_mid_col(pver), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'tm_mid_col', 'dyn_comp', __LINE__) + + call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv') + call mpas_dynamical_core % get_variable_pointer(rho, 'diag', 'rho') + call mpas_dynamical_core % get_variable_pointer(theta, 'diag', 'theta') + call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) + + ! 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) + 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_sfc(i) * & + exp(-0.5_kind_r8 * (zgrid(2, i) - zgrid(1, i)) * & + constant_g / (constant_rd * tm_mid_col(1)) * (1.0_kind_r8 + qv_mid_col(1))) + + ! Piecewise integrate hypsometric equation to derive subsequent `p_mid_col(k)`. + ! The formulation used here is exact. + do k = 2, pver + p_mid_col(k) = p_mid_col(k - 1) * & + exp(-0.5_kind_r8 * (zgrid(k , i) - zgrid(k - 1, i)) * & + constant_g / (constant_rd * tm_mid_col(k - 1)) * (1.0_kind_r8 + qv_mid_col(k - 1))) * & + exp(-0.5_kind_r8 * (zgrid(k + 1, i) - zgrid(k , i)) * & + constant_g / (constant_rd * tm_mid_col(k )) * (1.0_kind_r8 + qv_mid_col(k ))) + end do + + rho(:, i) = p_mid_col(:) / (constant_rd * tm_mid_col(:)) + theta(:, i) = t_mid(:, i) * ((constant_p0 / p_mid_col(:)) ** (constant_rd / constant_cpd)) + end do + + deallocate(p_mid_col) + deallocate(p_sfc) + deallocate(qv_mid_col) + deallocate(t_mid) + deallocate(tm_mid_col) + + nullify(index_qv) + nullify(rho) + nullify(theta) + nullify(scalars) + + ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. + call mpas_dynamical_core % exchange_halo('rho') + call mpas_dynamical_core % exchange_halo('theta') + 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() + integer :: i, k + integer :: ierr + 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(:, :) => null() + real(kind_r8), pointer :: theta_base(:, :) => null() + real(kind_r8), pointer :: zz(:, :) => null() + + call dyn_debug_print('Setting MPAS state "rho_base" and "theta_base"') + + allocate(p_base(pver), stat=ierr) + call check_allocate(ierr, 'set_analytic_initial_condition', 'p_base', 'dyn_comp', __LINE__) + + call mpas_dynamical_core % get_variable_pointer(rho_base, 'diag', 'rho_base') + call mpas_dynamical_core % get_variable_pointer(theta_base, 'diag', 'theta_base') + call mpas_dynamical_core % get_variable_pointer(zz, 'mesh', 'zz') + + ! Set `rho_base` and `theta_base` column by column. This way, peak memory usage can be reduced. + do i = 1, ncells_solve + do k = 1, pver + ! Derive `p_base` by hypsometric equation. + ! The formulation used here is exact and identical to MPAS. + p_base(k) = constant_p0 * exp(-0.5_kind_r8 * (zgrid(k, i) + zgrid(k + 1, i)) / & + (constant_rd * t_base / constant_g)) + end do + + rho_base(:, i) = p_base(:) / (constant_rd * t_base * zz(:, i)) + theta_base(:, i) = t_base * ((constant_p0 / p_base(:)) ** (constant_rd / constant_cpd)) + end do + + deallocate(p_base) + + nullify(rho_base) + nullify(theta_base) + nullify(zz) + + ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. + call mpas_dynamical_core % exchange_halo('rho_base') + call mpas_dynamical_core % exchange_halo('theta_base') + end subroutine set_mpas_state_rho_base_theta_base + end subroutine set_analytic_initial_condition + + !> 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. + !> (KCW, 2024-05-23) + subroutine mark_variable_as_initialized(constituent_name) + character(*), intent(in) :: constituent_name(:) + + character(*), parameter :: subname = 'dyn_comp::mark_variable_as_initialized' + integer :: i + + ! CCPP standard names of physical quantities in the `physics_{state,tend}` derived types. + call mark_as_initialized('air_pressure') + call mark_as_initialized('air_pressure_at_interface') + call mark_as_initialized('air_pressure_of_dry_air') + call mark_as_initialized('air_pressure_of_dry_air_at_interface') + call mark_as_initialized('air_pressure_thickness') + call mark_as_initialized('air_pressure_thickness_of_dry_air') + call mark_as_initialized('air_temperature') + call mark_as_initialized('dry_static_energy') + call mark_as_initialized('eastward_wind') + call mark_as_initialized('geopotential_height_wrt_surface') + call mark_as_initialized('geopotential_height_wrt_surface_at_interface') + call mark_as_initialized('inverse_exner_function_wrt_surface_pressure') + call mark_as_initialized('lagrangian_tendency_of_air_pressure') + call mark_as_initialized('ln_air_pressure') + call mark_as_initialized('ln_air_pressure_at_interface') + call mark_as_initialized('ln_air_pressure_of_dry_air') + call mark_as_initialized('ln_air_pressure_of_dry_air_at_interface') + call mark_as_initialized('northward_wind') + 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('surface_air_pressure') + call mark_as_initialized('surface_geopotential') + call mark_as_initialized('surface_pressure_of_dry_air') + call mark_as_initialized('tendency_of_air_temperature_due_to_model_physics') + call mark_as_initialized('tendency_of_eastward_wind_due_to_model_physics') + call mark_as_initialized('tendency_of_northward_wind_due_to_model_physics') + + ! CCPP standard names of constituents. + do i = 1, num_advected + call mark_as_initialized(trim(adjustl(constituent_name(i)))) + end do + end subroutine mark_variable_as_initialized + ! Not used for now. Intended to be called by `stepon_run*` in `src/dynamics/mpas/stepon.F90`. ! subroutine dyn_run() ! end subroutine dyn_run diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 0074df22..f65d657f 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -12,7 +12,7 @@ module dyn_grid ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max, & sphere_radius, & deg_to_rad, rad_to_deg - use dynconst, only: dynconst_init, pi + use dynconst, only: constant_pi => pi, 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 @@ -77,6 +77,9 @@ subroutine model_grid_init() ! Read time-invariant (e.g., grid/mesh) variables. call mpas_dynamical_core % read_write_stream(pio_file, 'r', 'invariant') + ! 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') @@ -232,7 +235,7 @@ subroutine init_physics_grid() ! 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 * pi * sphere_radius ** 2), kind_pcol) + dyn_column(i) % weight = real(areacell(i) / (4.0_kind_r8 * constant_pi * sphere_radius ** 2), kind_pcol) ! File decomposition. ! For unstructured grid, `coord_indices` is not used by `phys_grid_init`. @@ -331,7 +334,7 @@ subroutine define_cam_grid() do i = 1, ncells_solve cell_area(i) = areacell(i) - cell_weight(i) = areacell(i) / (4.0_kind_r8 * pi * sphere_radius ** 2) + cell_weight(i) = areacell(i) / (4.0_kind_r8 * constant_pi * sphere_radius ** 2) global_grid_map(1, i) = int(i, kind_imap) global_grid_map(2, i) = int(1, kind_imap) From 0c564c0eb34d95cef3a90fb9b8841b4f14c9e851 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 29 May 2024 17:11:10 -0600 Subject: [PATCH 11/31] Replace mentions of "CAM" with "CAM-SIMA" --- src/dynamics/mpas/dyn_comp.F90 | 8 ++++---- src/dynamics/mpas/dyn_grid.F90 | 4 ++-- src/dynamics/mpas/stepon.F90 | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 00ff9825..c064a8fc 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -1,7 +1,7 @@ module dyn_comp use dyn_mpas_subdriver, only: mpas_dynamical_core_type - ! Modules from CAM. + ! Modules from CAM-SIMA. use air_composition, only: thermodynamic_active_species_num, & thermodynamic_active_species_liq_num, & thermodynamic_active_species_ice_num, & @@ -41,7 +41,7 @@ module dyn_comp implicit none private - ! Provide APIs required by CAM Control. + ! Provide APIs required by CAM-SIMA. public :: dyn_import_t public :: dyn_export_t public :: dyn_readnl @@ -57,14 +57,14 @@ module dyn_comp public :: deg_to_rad, rad_to_deg !> NOTE: - !> This derived type is not used by MPAS dynamical core. It exists only as a placeholder because CAM Control requires it. + !> This derived type is not used by MPAS dynamical core. It exists only as a placeholder because CAM-SIMA requires it. !> Developers/Maintainers/Users who wish to interact with MPAS dynamical core may do so by using the "instance/object" !> below. type :: dyn_import_t end type dyn_import_t !> NOTE: - !> This derived type is not used by MPAS dynamical core. It exists only as a placeholder because CAM Control requires it. + !> This derived type is not used by MPAS dynamical core. It exists only as a placeholder because CAM-SIMA requires it. !> Developers/Maintainers/Users who wish to interact with MPAS dynamical core may do so by using the "instance/object" !> below. type :: dyn_export_t diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index f65d657f..fd8d49b2 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -1,5 +1,5 @@ module dyn_grid - ! Modules from CAM. + ! 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, & @@ -30,7 +30,7 @@ module dyn_grid implicit none private - ! Provide APIs required by CAM Control. + ! Provide APIs required by CAM-SIMA. public :: model_grid_init public :: dyn_grid_id diff --git a/src/dynamics/mpas/stepon.F90 b/src/dynamics/mpas/stepon.F90 index 49784c5e..de99739e 100644 --- a/src/dynamics/mpas/stepon.F90 +++ b/src/dynamics/mpas/stepon.F90 @@ -8,7 +8,7 @@ module stepon implicit none private - ! Provide APIs required by CAM Control. + ! Provide APIs required by CAM-SIMA. public :: stepon_init public :: stepon_run1 public :: stepon_run2 From a87368ac1125ca20e2ce73bbf6ed7c6a030afefc Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 24 Jun 2024 16:11:03 -0600 Subject: [PATCH 12/31] Add more code comments --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 046ba6fe..cd8c9c4d 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -1226,7 +1226,9 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_stream_with_pool' character(strkind) :: stream_filename integer :: i, ierr, stream_format + ! Whether a variable is present on the file (i.e., `pio_file`). logical, allocatable :: var_is_present(:) + ! Whether a variable is type, kind and rank compatible with what MPAS expects on the file (i.e., `pio_file`). logical, allocatable :: var_is_tkr_compatible(:) type(field0dchar), pointer :: field_0d_char => null() type(field1dchar), pointer :: field_1d_char => null() @@ -1807,7 +1809,7 @@ end function index_unique !------------------------------------------------------------------------------- ! subroutine dyn_mpas_check_variable_status ! - !> \brief Check and return variable status on the given file. + !> \brief Check and return variable status on the given file !> \author Kuan-Chih Wang !> \date 2024-06-04 !> \details From b7dccf3e5e708c46edd7062d58396bcdc5ef59fe Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 24 Jun 2024 16:24:28 -0600 Subject: [PATCH 13/31] Include more constants from `physconst` in `dynconst` --- src/dynamics/utils/dynconst.F90 | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/dynamics/utils/dynconst.F90 b/src/dynamics/utils/dynconst.F90 index c774cf82..84574749 100644 --- a/src/dynamics/utils/dynconst.F90 +++ b/src/dynamics/utils/dynconst.F90 @@ -30,6 +30,10 @@ module dynconst real(kind_dyn), protected, public :: cpair ! Dry air gas constant [J/K/kg] real(kind_dyn), protected, public :: rair + ! water vapor gas constant [J/K/kg] + real(kind_dyn), protected, public :: rh2o + ! reference surface pressure [Pa] + real(kind_dyn), protected, public :: pref ! reference temperature [K] real(kind_dyn), protected, public :: tref ! reference lapse rate [K/m] @@ -57,10 +61,12 @@ subroutine dynconst_init use physconst, only: phys_omega=>omega use physconst, only: phys_cpair=>cpair use physconst, only: phys_gravit=>gravit + use physconst, only: phys_pref=>pref use physconst, only: phys_tref=>tref use physconst, only: phys_lapse_rate=>lapse_rate use physconst, only: phys_cappa=>cappa use physconst, only: phys_rair=>rair + use physconst, only: phys_rh2o=>rh2o !Set constants used by the dynamics: @@ -69,7 +75,9 @@ subroutine dynconst_init omega = real(phys_omega, kind_dyn) cpair = real(phys_cpair, kind_dyn) rair = real(phys_rair, kind_dyn) + rh2o = real(phys_rh2o, kind_dyn) gravit = real(phys_gravit, kind_dyn) + pref = real(phys_pref, kind_dyn) tref = real(phys_tref, kind_dyn) lapse_rate = real(phys_lapse_rate, kind_dyn) cappa = real(phys_cappa, kind_dyn) From 03cf4e001348d372c7362d2db6162aafbd6d00b4 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 24 Jun 2024 16:25:17 -0600 Subject: [PATCH 14/31] Use constants from `dynconst` only for consistent kind parameters --- src/dynamics/mpas/dyn_comp.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index c064a8fc..ead7bf0a 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -18,9 +18,9 @@ module dyn_comp 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_g => gravit, constant_pi => pi + use dynconst, only: constant_cpd => cpair, constant_g => gravit, constant_p0 => pref, constant_pi => pi, & + constant_rd => rair, constant_rv => rh2o use inic_analytic, only: analytic_ic_active, dyn_set_inic_col - use physconst, only: constant_cpd => cpair, constant_p0 => pref, constant_rd => rair, constant_rv => rh2o use runtime_obj, only: runtime_options use spmd_utils, only: iam, masterproc, mpicom use string_utils, only: stringify From ef3f0761d6714e5328c99f39b8af48c59325f169 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 24 Jun 2024 17:27:42 -0600 Subject: [PATCH 15/31] Include array shapes in allocation error messages --- src/dynamics/mpas/dyn_comp.F90 | 45 ++++++++++++++++++++-------------- src/dynamics/mpas/dyn_grid.F90 | 33 +++++++++++++------------ 2 files changed, 43 insertions(+), 35 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index ead7bf0a..8a6cc3ff 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -197,10 +197,10 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) type(file_desc_t), pointer :: pio_topo_file => null() allocate(constituent_name(num_advected), stat=ierr) - call check_allocate(ierr, 'dyn_init', 'constituent_name', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'dyn_init', 'constituent_name(num_advected)', 'dyn_comp', __LINE__) allocate(is_water_species(num_advected), stat=ierr) - call check_allocate(ierr, 'dyn_init', 'is_water_species', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'dyn_init', 'is_water_species(num_advected)', 'dyn_comp', __LINE__) do i = 1, num_advected constituent_name(i) = const_name(i) @@ -311,10 +311,10 @@ subroutine check_topography_data(pio_file) end if allocate(surface_geopotential(ncells_solve), stat=ierr) - call check_allocate(ierr, 'check_topography_data', 'surface_geopotential', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'check_topography_data', 'surface_geopotential(ncells_solve)', 'dyn_comp', __LINE__) allocate(surface_geometric_height(ncells_solve), stat=ierr) - call check_allocate(ierr, 'check_topography_data', 'surface_geometric_height', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'check_topography_data', 'surface_geometric_height(ncells_solve)', 'dyn_comp', __LINE__) surface_geopotential(:) = 0.0_kind_r8 surface_geometric_height(:) = 0.0_kind_r8 @@ -385,7 +385,7 @@ subroutine init_shared_variable() call dyn_debug_print('Preparing to set analytic initial condition') allocate(global_grid_index(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'global_grid_index', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'global_grid_index(ncells_solve)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(indextocellid, 'mesh', 'indexToCellID') @@ -394,10 +394,10 @@ subroutine init_shared_variable() nullify(indextocellid) allocate(lat_rad(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'lat_rad', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'lat_rad(ncells_solve)', 'dyn_comp', __LINE__) allocate(lon_rad(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'lon_rad', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'lon_rad(ncells_solve)', 'dyn_comp', __LINE__) ! "mpas_cell" is a registered grid name that is defined in `dyn_grid`. lat_deg => cam_grid_get_latvals(cam_grid_id('mpas_cell')) @@ -418,7 +418,7 @@ subroutine init_shared_variable() nullify(lon_deg) allocate(z_int(ncells_solve, pverp), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'z_int', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'z_int(ncells_solve, pverp)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid') @@ -439,7 +439,8 @@ subroutine set_mpas_state_u() call dyn_debug_print('Setting MPAS state "u"') allocate(buffer_2d_real(ncells_solve, pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_2d_real', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_2d_real(ncells_solve, pver)', & + 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal') call mpas_dynamical_core % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional') @@ -503,10 +504,12 @@ subroutine set_mpas_state_scalars() call dyn_debug_print('Setting MPAS state "scalars"') allocate(buffer_3d_real(ncells_solve, pver, num_advected), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_3d_real', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_3d_real(ncells_solve, pver, num_advected)', & + 'dyn_comp', __LINE__) allocate(constituent_index(num_advected), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'constituent_index', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'constituent_index(num_advected)', & + 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv') call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) @@ -569,10 +572,12 @@ subroutine set_mpas_state_rho_theta() call dyn_debug_print('Setting MPAS state "rho" and "theta"') allocate(buffer_1d_real(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_1d_real', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_1d_real(ncells_solve)', & + 'dyn_comp', __LINE__) allocate(p_sfc(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'p_sfc', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'p_sfc(ncells_solve)', & + 'dyn_comp', __LINE__) buffer_1d_real(:) = 0.0_kind_r8 @@ -583,10 +588,12 @@ subroutine set_mpas_state_rho_theta() deallocate(buffer_1d_real) allocate(buffer_2d_real(ncells_solve, pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_2d_real', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_2d_real(ncells_solve, pver)', & + 'dyn_comp', __LINE__) allocate(t_mid(pver, ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 't_mid', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 't_mid(pver, ncells_solve)', & + 'dyn_comp', __LINE__) buffer_2d_real(:, :) = 0.0_kind_r8 @@ -600,13 +607,13 @@ subroutine set_mpas_state_rho_theta() deallocate(buffer_2d_real) allocate(p_mid_col(pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'p_mid_col', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'p_mid_col(pver)', 'dyn_comp', __LINE__) allocate(qv_mid_col(pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'qv_mid_col', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'qv_mid_col(pver)', 'dyn_comp', __LINE__) allocate(tm_mid_col(pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'tm_mid_col', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'tm_mid_col(pver)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv') call mpas_dynamical_core % get_variable_pointer(rho, 'diag', 'rho') @@ -669,7 +676,7 @@ subroutine set_mpas_state_rho_base_theta_base() call dyn_debug_print('Setting MPAS state "rho_base" and "theta_base"') allocate(p_base(pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'p_base', 'dyn_comp', __LINE__) + call check_allocate(ierr, 'set_analytic_initial_condition', 'p_base(pver)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(rho_base, 'diag', 'rho_base') call mpas_dynamical_core % get_variable_pointer(theta_base, 'diag', 'theta_base') diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index fd8d49b2..e10d6eb5 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -144,16 +144,16 @@ subroutine init_reference_pressure() call mpas_dynamical_core % get_variable_pointer(rdzw, 'mesh', 'rdzw') allocate(dzw(pver), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'dzw', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'init_reference_pressure', 'dzw(pver)', 'dyn_grid', __LINE__) dzw(:) = 1.0_kind_r8 / rdzw nullify(rdzw) allocate(zw(pverp), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'zw', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'init_reference_pressure', 'zw(pverp)', 'dyn_grid', __LINE__) allocate(zu(pver), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'zu', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'init_reference_pressure', 'zu(pver)', 'dyn_grid', __LINE__) ! In MPAS, zeta coordinates are stored in increasing order (i.e., bottom to top of atmosphere). ! In CAM-SIMA, however, index order is reversed (i.e., top to bottom of atmosphere). @@ -167,12 +167,12 @@ subroutine init_reference_pressure() ! Compute reference pressure from reference height. allocate(p_ref_int(pverp), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'p_ref_int', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'init_reference_pressure', 'p_ref_int(pverp)', 'dyn_grid', __LINE__) call std_atm_pres(zw, p_ref_int) allocate(p_ref_mid(pver), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'p_ref_mid', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'init_reference_pressure', 'p_ref_mid(pver)', 'dyn_grid', __LINE__) p_ref_mid(:) = 0.5_kind_r8 * (p_ref_int(1:pver) + p_ref_int(2:pverp)) @@ -224,7 +224,7 @@ subroutine init_physics_grid() call mpas_dynamical_core % get_variable_pointer(loncell, 'mesh', 'lonCell') allocate(dyn_column(ncells_solve), stat=ierr) - call check_allocate(ierr, 'init_physics_grid', 'dyn_column', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'init_physics_grid', 'dyn_column(ncells_solve)', 'dyn_grid', __LINE__) do i = 1, ncells_solve ! Column information. @@ -249,7 +249,8 @@ subroutine init_physics_grid() dyn_column(i) % local_dyn_block = i ! `dyn_block_index` is not used due to no dynamics block offset, but it still needs to be allocated. allocate(dyn_column(i) % dyn_block_index(0), stat=ierr) - call check_allocate(ierr, 'init_physics_grid', 'dyn_column % dyn_block_index', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'init_physics_grid', 'dyn_column(' // stringify([i]) // ') % dyn_block_index(0)', & + 'dyn_grid', __LINE__) end do nullify(areacell) @@ -260,7 +261,7 @@ subroutine init_physics_grid() ! `phys_grid_init` expects to receive the `area` attribute from dynamics. ! However, do not let it because dynamics grid is different from physics grid. allocate(dyn_attribute_name(0), stat=ierr) - call check_allocate(ierr, 'init_physics_grid', 'dyn_attribute_name', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'init_physics_grid', 'dyn_attribute_name(0)', 'dyn_grid', __LINE__) call phys_grid_init(hdim1_d, hdim2_d, 'mpas', dyn_column, 'mpas_cell', dyn_attribute_name) end subroutine init_physics_grid @@ -316,7 +317,7 @@ subroutine define_cam_grid() call mpas_dynamical_core % get_variable_pointer(loncell, 'mesh', 'lonCell') allocate(global_grid_index(ncells_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_index', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'define_cam_grid', 'global_grid_index(ncells_solve)', 'dyn_grid', __LINE__) global_grid_index(:) = int(indextocellid(1:ncells_solve), kind_imap) @@ -326,11 +327,11 @@ subroutine define_cam_grid() 1, ncells_solve, loncell * rad_to_deg, map=global_grid_index) allocate(cell_area(ncells_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'cell_area', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'define_cam_grid', 'cell_area(ncells_solve)', 'dyn_grid', __LINE__) allocate(cell_weight(ncells_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'cell_weight', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'define_cam_grid', 'cell_weight(ncells_solve)', 'dyn_grid', __LINE__) allocate(global_grid_map(3, ncells_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_map', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'define_cam_grid', 'global_grid_map(3, ncells_solve)', 'dyn_grid', __LINE__) do i = 1, ncells_solve cell_area(i) = areacell(i) @@ -387,7 +388,7 @@ subroutine define_cam_grid() call mpas_dynamical_core % get_variable_pointer(lonedge, 'mesh', 'lonEdge') allocate(global_grid_index(nedges_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_index', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'define_cam_grid', 'global_grid_index(nedges_solve)', 'dyn_grid', __LINE__) global_grid_index(:) = int(indextoedgeid(1:nedges_solve), kind_imap) @@ -397,7 +398,7 @@ subroutine define_cam_grid() 1, nedges_solve, lonedge * rad_to_deg, map=global_grid_index) allocate(global_grid_map(3, nedges_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_map', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'define_cam_grid', 'global_grid_map(3, nedges_solve)', 'dyn_grid', __LINE__) do i = 1, nedges_solve global_grid_map(1, i) = int(i, kind_imap) @@ -428,7 +429,7 @@ subroutine define_cam_grid() call mpas_dynamical_core % get_variable_pointer(lonvertex, 'mesh', 'lonVertex') allocate(global_grid_index(nvertices_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_index', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'define_cam_grid', 'global_grid_index(nvertices_solve)', 'dyn_grid', __LINE__) global_grid_index(:) = int(indextovertexid(1:nvertices_solve), kind_imap) @@ -438,7 +439,7 @@ subroutine define_cam_grid() 1, nvertices_solve, lonvertex * rad_to_deg, map=global_grid_index) allocate(global_grid_map(3, nvertices_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_map', 'dyn_grid', __LINE__) + call check_allocate(ierr, 'define_cam_grid', 'global_grid_map(3, nvertices_solve)', 'dyn_grid', __LINE__) do i = 1, nvertices_solve global_grid_map(1, i) = int(i, kind_imap) From 7a6451c936d57dbcaea05cd504a26614b7787bad Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 24 Jun 2024 17:51:40 -0600 Subject: [PATCH 16/31] Avoid unnecessary array allocation and copy --- src/dynamics/mpas/dyn_comp.F90 | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 8a6cc3ff..6432480f 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -352,7 +352,7 @@ end subroutine check_topography_data subroutine set_analytic_initial_condition() character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition' integer, allocatable :: global_grid_index(:) - real(kind_r8), allocatable :: buffer_1d_real(:), buffer_2d_real(:, :), buffer_3d_real(:, :, :) + 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. @@ -571,21 +571,13 @@ subroutine set_mpas_state_rho_theta() call dyn_debug_print('Setting MPAS state "rho" and "theta"') - allocate(buffer_1d_real(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_1d_real(ncells_solve)', & - 'dyn_comp', __LINE__) - allocate(p_sfc(ncells_solve), stat=ierr) call check_allocate(ierr, 'set_analytic_initial_condition', 'p_sfc(ncells_solve)', & 'dyn_comp', __LINE__) - buffer_1d_real(:) = 0.0_kind_r8 - - call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, ps=buffer_1d_real) - - p_sfc(:) = buffer_1d_real(:) + p_sfc(:) = 0.0_kind_r8 - deallocate(buffer_1d_real) + call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, ps=p_sfc) allocate(buffer_2d_real(ncells_solve, pver), stat=ierr) call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_2d_real(ncells_solve, pver)', & From b3f1937bb868856590c9380d97c5897ee6600ae0 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 28 Jun 2024 18:02:16 -0600 Subject: [PATCH 17/31] Include variable names in error messages --- .../mpas/driver/dyn_mpas_subdriver.F90 | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index cd8c9c4d..d1456c22 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -1361,8 +1361,8 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file nullify(field_1d_char) case default - call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]), & - subname, __LINE__) + call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]) // & + ' for "' // trim(adjustl(var_info_list(i) % name)) // '"', subname, __LINE__) end select case ('integer') select case (var_info_list(i) % rank) @@ -1415,8 +1415,8 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file nullify(field_3d_integer) case default - call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]), & - subname, __LINE__) + call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]) // & + ' for "' // trim(adjustl(var_info_list(i) % name)) // '"', subname, __LINE__) end select case ('real') select case (var_info_list(i) % rank) @@ -1493,12 +1493,12 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file nullify(field_5d_real) case default - call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]), & - subname, __LINE__) + call self % model_error('Unsupported variable rank ' // stringify([var_info_list(i) % rank]) // & + ' for "' // trim(adjustl(var_info_list(i) % name)) // '"', subname, __LINE__) end select case default call self % model_error('Unsupported variable type "' // trim(adjustl(var_info_list(i) % type)) // & - '"', subname, __LINE__) + '" for "' // trim(adjustl(var_info_list(i) % name)) // '"', subname, __LINE__) end select if (ierr /= mpas_stream_noerr) then @@ -1894,8 +1894,8 @@ subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compa nullify(field_1d_char) case default - call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]), & - subname, __LINE__) + call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]) // & + ' for "' // trim(adjustl(var_info % name)) // '"', subname, __LINE__) end select case ('integer') select case (var_info % rank) @@ -1980,8 +1980,8 @@ subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compa nullify(field_3d_integer) case default - call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]), & - subname, __LINE__) + call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]) // & + ' for "' // trim(adjustl(var_info % name)) // '"', subname, __LINE__) end select case ('real') select case (var_info % rank) @@ -2106,12 +2106,12 @@ subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compa nullify(field_5d_real) case default - call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]), & - subname, __LINE__) + call self % model_error('Unsupported variable rank ' // stringify([var_info % rank]) // & + ' for "' // trim(adjustl(var_info % name)) // '"', subname, __LINE__) end select case default call self % model_error('Unsupported variable type "' // trim(adjustl(var_info % type)) // & - '"', subname, __LINE__) + '" for "' // trim(adjustl(var_info % name)) // '"', subname, __LINE__) end select if (.not. allocated(var_name_list)) then From 77cf091c24a89ee229b2ba291f928b85cb9997d5 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 5 Jul 2024 17:06:30 -0600 Subject: [PATCH 18/31] Extract and separate equations into their own --- src/dynamics/mpas/dyn_comp.F90 | 73 ++++++++++++++++++++++++++++------ 1 file changed, 61 insertions(+), 12 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 6432480f..c0a52a4f 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -619,22 +619,37 @@ subroutine set_mpas_state_rho_theta() ! Piecewise integrate hypsometric equation to derive `p_mid_col(1)`. ! The formulation used here is exact. - p_mid_col(1) = p_sfc(i) * & - exp(-0.5_kind_r8 * (zgrid(2, i) - zgrid(1, i)) * & - constant_g / (constant_rd * tm_mid_col(1)) * (1.0_kind_r8 + qv_mid_col(1))) + ! p_mid_col(1) = p_sfc(i) * & + ! exp(-0.5_kind_r8 * (zgrid(2, i) - zgrid(1, i)) * & + ! constant_g / (constant_rd * tm_mid_col(1)) * (1.0_kind_r8 + qv_mid_col(1))) + p_mid_col(1) = p_by_hypsometric_equation( & + p_sfc(i), & + zgrid(1, i), & + tm_mid_col(1) / (1.0_kind_r8 + qv_mid_col(1)), & + 0.5_kind_r8 * (zgrid(2, i) + zgrid(1, i))) ! Piecewise integrate hypsometric equation to derive subsequent `p_mid_col(k)`. ! The formulation used here is exact. do k = 2, pver - p_mid_col(k) = p_mid_col(k - 1) * & - exp(-0.5_kind_r8 * (zgrid(k , i) - zgrid(k - 1, i)) * & - constant_g / (constant_rd * tm_mid_col(k - 1)) * (1.0_kind_r8 + qv_mid_col(k - 1))) * & - exp(-0.5_kind_r8 * (zgrid(k + 1, i) - zgrid(k , i)) * & - constant_g / (constant_rd * tm_mid_col(k )) * (1.0_kind_r8 + qv_mid_col(k ))) + ! p_mid_col(k) = p_mid_col(k - 1) * & + ! exp(-0.5_kind_r8 * (zgrid(k , i) - zgrid(k - 1, i)) * & + ! constant_g / (constant_rd * tm_mid_col(k - 1)) * (1.0_kind_r8 + qv_mid_col(k - 1))) * & + ! exp(-0.5_kind_r8 * (zgrid(k + 1, i) - zgrid(k , i)) * & + ! constant_g / (constant_rd * tm_mid_col(k )) * (1.0_kind_r8 + qv_mid_col(k ))) + 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)), & + tm_mid_col(k - 1) / (1.0_kind_r8 + qv_mid_col(k - 1)), & + zgrid(k, i)), & + zgrid(k, i), & + tm_mid_col(k) / (1.0_kind_r8 + qv_mid_col(k)), & + 0.5_kind_r8 * (zgrid(k + 1, i) + zgrid(k, i))) end do rho(:, i) = p_mid_col(:) / (constant_rd * tm_mid_col(:)) - theta(:, i) = t_mid(:, i) * ((constant_p0 / p_mid_col(:)) ** (constant_rd / constant_cpd)) + ! theta(:, i) = t_mid(:, i) * ((constant_p0 / p_mid_col(:)) ** (constant_rd / constant_cpd)) + theta(:, i) = theta_by_poisson_equation(p_mid_col, t_mid(:, i), constant_p0) end do deallocate(p_mid_col) @@ -679,12 +694,18 @@ subroutine set_mpas_state_rho_base_theta_base() do k = 1, pver ! Derive `p_base` by hypsometric equation. ! The formulation used here is exact and identical to MPAS. - p_base(k) = constant_p0 * exp(-0.5_kind_r8 * (zgrid(k, i) + zgrid(k + 1, i)) / & - (constant_rd * t_base / constant_g)) + ! p_base(k) = constant_p0 * exp(-0.5_kind_r8 * (zgrid(k, i) + zgrid(k + 1, i)) / & + ! (constant_rd * t_base / constant_g)) + p_base(k) = p_by_hypsometric_equation( & + constant_p0, & + 0.0_kind_r8, & + t_base, & + 0.5_kind_r8 * (zgrid(k + 1, i) + zgrid(k, i))) end do rho_base(:, i) = p_base(:) / (constant_rd * t_base * zz(:, i)) - theta_base(:, i) = t_base * ((constant_p0 / p_base(:)) ** (constant_rd / constant_cpd)) + ! theta_base(:, i) = t_base * ((constant_p0 / p_base(:)) ** (constant_rd / constant_cpd)) + theta_base(:, i) = theta_by_poisson_equation(p_base, t_base, constant_p0) end do deallocate(p_base) @@ -697,6 +718,34 @@ subroutine set_mpas_state_rho_base_theta_base() call mpas_dynamical_core % exchange_halo('rho_base') call mpas_dynamical_core % exchange_halo('theta_base') end subroutine set_mpas_state_rho_base_theta_base + + ! ----- p_2, z_2 ----- (Layer 2) + ! t_v + ! ----- p_1, z_1 ----- (Layer 1) + ! + !> Compute the pressure `p_2` at height `z_2` from the pressure `p_1` at height `z_1` by hypsometric equation. + !> `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) + real(kind_r8), intent(in) :: p_1, z_1, t_v, z_2 + real(kind_r8) :: p_2 + + p_2 = p_1 * exp(-(z_2 - z_1) * constant_g / (constant_rd * t_v)) + end function p_by_hypsometric_equation + + ! ----- p_1, t_1 ----- (Arbitrary layer) + ! + ! ----- p_0, t_0 ----- (Reference layer) + ! + !> Compute the potential temperature `t_0` at reference pressure `p_0` from the temperature `t_1` at pressure `p_1` by + !> Poisson equation. + !> (KCW, 2024-07-02) + pure elemental function theta_by_poisson_equation(p_1, t_1, p_0) result(t_0) + real(kind_r8), intent(in) :: p_1, t_1, p_0 + real(kind_r8) :: t_0 + + t_0 = t_1 * ((p_0 / p_1) ** (constant_rd / constant_cpd)) + end function theta_by_poisson_equation end subroutine set_analytic_initial_condition !> Mark everything in the `physics_{state,tend}` derived types along with constituents as initialized From d3567bf2793e631053314022dca16613ba83dfd7 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Tue, 9 Jul 2024 17:19:07 -0600 Subject: [PATCH 19/31] Implement default constituent initialization --- src/dynamics/mpas/dyn_comp.F90 | 53 +++++++++++++++++++++++++++------- 1 file changed, 43 insertions(+), 10 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index c0a52a4f..fb3f8602 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -33,6 +33,8 @@ module dyn_comp use shr_pio_mod, only: shr_pio_getiosys ! 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 ! Modules from external libraries. @@ -241,23 +243,20 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) call set_analytic_initial_condition() else + ! 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 set_default_constituent') + + call set_default_constituent() + ! Namelist option that controls if constituents are to be read from the file. if (readtrace) then ! Read variables that belong to the "input" stream in MPAS. call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input') - - ! TODO: - ! `cnst_init_default` or its replacement is not yet implemented in CAM-SIMA. - ! Default constituent initialization should be performed here - ! for constituents that fail to be read from the file. else ! Read variables that belong to the "input" stream in MPAS, excluding constituents. call mpas_dynamical_core % read_write_stream(pio_init_file, 'r', 'input-scalars') - - ! TODO: - ! `cnst_init_default` or its replacement is not yet implemented in CAM-SIMA. - ! Default constituent initialization should be performed here - ! because constituents are not read from the file. end if end if else @@ -748,6 +747,40 @@ pure elemental function theta_by_poisson_equation(p_1, t_1, p_0) result(t_0) end function theta_by_poisson_equation end subroutine set_analytic_initial_condition + !> Set default MPAS state `scalars` (i.e., constituents) in accordance with CCPP, which is a component of CAM-SIMA. + !> (KCW, 2024-07-09) + subroutine set_default_constituent() + character(*), parameter :: subname = 'dyn_comp::set_default_constituent' + integer :: i, k + real(kind_phys), pointer :: constituents(:, :, :) => null() ! This points to CCPP memory. + real(kind_r8), pointer :: scalars(:, :, :) => null() ! This points to MPAS memory. + + call dyn_debug_print('Setting default MPAS state "scalars"') + + constituents => cam_constituents_array() + + if (.not. associated(constituents)) then + call endrun('Failed to find variable "constituents"', subname, __LINE__) + end if + + call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) + + ! `i` is indexing into `scalars`, so it is regarded as MPAS scalar index. + do i = 1, num_advected + ! Vertical index order is reversed between CAM-SIMA and MPAS. + do k = 1, pver + scalars(i, k, 1:ncells_solve) = & + constituents(:, pver - k + 1, mpas_dynamical_core % map_constituent_index(i)) + end do + end do + + nullify(constituents) + nullify(scalars) + + ! Because we are injecting data directly into MPAS memory, halo layers need to be updated manually. + call mpas_dynamical_core % exchange_halo('scalars') + end subroutine set_default_constituent + !> 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. From 724d5c859556d4ed9bc13de039b574f0810b7252 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Tue, 9 Jul 2024 17:51:25 -0600 Subject: [PATCH 20/31] Use accessor functions from `cam_constituents` directly --- src/dynamics/mpas/dyn_comp.F90 | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index fb3f8602..10864129 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -212,6 +212,9 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) ! Inform MPAS about constituent names and their corresponding waterness. call mpas_dynamical_core % define_scalar(constituent_name, is_water_species) + deallocate(constituent_name) + deallocate(is_water_species) + ! Provide mapping information between MPAS scalars and constituent names to CAM-SIMA. do i = 1, thermodynamic_active_species_num thermodynamic_active_species_idx_dycore(i) = & @@ -267,10 +270,7 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) end if call clean_iodesc_list() - call mark_variable_as_initialized(constituent_name) - - deallocate(constituent_name) - deallocate(is_water_species) + call mark_variable_as_initialized() nullify(pio_init_file) nullify(pio_topo_file) @@ -785,9 +785,7 @@ end subroutine set_default_constituent !> to prevent physics from attempting to read them from a file. These variables are to be exchanged later !> during dynamics-physics coupling. !> (KCW, 2024-05-23) - subroutine mark_variable_as_initialized(constituent_name) - character(*), intent(in) :: constituent_name(:) - + subroutine mark_variable_as_initialized() character(*), parameter :: subname = 'dyn_comp::mark_variable_as_initialized' integer :: i @@ -821,7 +819,7 @@ subroutine mark_variable_as_initialized(constituent_name) ! CCPP standard names of constituents. do i = 1, num_advected - call mark_as_initialized(trim(adjustl(constituent_name(i)))) + call mark_as_initialized(trim(adjustl(const_name(i)))) end do end subroutine mark_variable_as_initialized From d57d286b7208a14129bd42bf8b450a637f452845 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Thu, 11 Jul 2024 10:42:53 -0600 Subject: [PATCH 21/31] Add a blank line in comments Maintain consistent style with others. --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index d1456c22..ed7d11f1 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -1821,6 +1821,7 @@ end function index_unique !> in MPAS registry. For an ordinary variable, the checks are performed on !> itself. Otherwise, for a variable array, the checks are performed on its !> constituent parts instead. + ! !------------------------------------------------------------------------------- subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compatible, pio_file, var_info) class(mpas_dynamical_core_type), intent(in) :: self From 3ca360bb44c2935ae58afce307c82a2829cc51b6 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 22 Jul 2024 13:31:15 -0600 Subject: [PATCH 22/31] Remove defunct code --- src/dynamics/mpas/dyn_comp.F90 | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 10864129..5d862b63 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -618,9 +618,6 @@ subroutine set_mpas_state_rho_theta() ! Piecewise integrate hypsometric equation to derive `p_mid_col(1)`. ! The formulation used here is exact. - ! p_mid_col(1) = p_sfc(i) * & - ! exp(-0.5_kind_r8 * (zgrid(2, i) - zgrid(1, i)) * & - ! constant_g / (constant_rd * tm_mid_col(1)) * (1.0_kind_r8 + qv_mid_col(1))) p_mid_col(1) = p_by_hypsometric_equation( & p_sfc(i), & zgrid(1, i), & @@ -630,11 +627,6 @@ subroutine set_mpas_state_rho_theta() ! Piecewise integrate hypsometric equation to derive subsequent `p_mid_col(k)`. ! The formulation used here is exact. do k = 2, pver - ! p_mid_col(k) = p_mid_col(k - 1) * & - ! exp(-0.5_kind_r8 * (zgrid(k , i) - zgrid(k - 1, i)) * & - ! constant_g / (constant_rd * tm_mid_col(k - 1)) * (1.0_kind_r8 + qv_mid_col(k - 1))) * & - ! exp(-0.5_kind_r8 * (zgrid(k + 1, i) - zgrid(k , i)) * & - ! constant_g / (constant_rd * tm_mid_col(k )) * (1.0_kind_r8 + qv_mid_col(k ))) p_mid_col(k) = p_by_hypsometric_equation( & p_by_hypsometric_equation( & p_mid_col(k - 1), & @@ -647,7 +639,6 @@ subroutine set_mpas_state_rho_theta() end do rho(:, i) = p_mid_col(:) / (constant_rd * tm_mid_col(:)) - ! theta(:, i) = t_mid(:, i) * ((constant_p0 / p_mid_col(:)) ** (constant_rd / constant_cpd)) theta(:, i) = theta_by_poisson_equation(p_mid_col, t_mid(:, i), constant_p0) end do @@ -693,8 +684,6 @@ subroutine set_mpas_state_rho_base_theta_base() do k = 1, pver ! Derive `p_base` by hypsometric equation. ! The formulation used here is exact and identical to MPAS. - ! p_base(k) = constant_p0 * exp(-0.5_kind_r8 * (zgrid(k, i) + zgrid(k + 1, i)) / & - ! (constant_rd * t_base / constant_g)) p_base(k) = p_by_hypsometric_equation( & constant_p0, & 0.0_kind_r8, & @@ -703,7 +692,6 @@ subroutine set_mpas_state_rho_base_theta_base() end do rho_base(:, i) = p_base(:) / (constant_rd * t_base * zz(:, i)) - ! theta_base(:, i) = t_base * ((constant_p0 / p_base(:)) ** (constant_rd / constant_cpd)) theta_base(:, i) = theta_by_poisson_equation(p_base, t_base, constant_p0) end do From 9cc38a476cc6b27bd50d4123e225bb954bfa5112 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Mon, 22 Jul 2024 14:08:03 -0600 Subject: [PATCH 23/31] Move degree-radian conversion constants to `dynconst` --- src/dynamics/mpas/dyn_comp.F90 | 7 ++----- src/dynamics/mpas/dyn_grid.F90 | 5 ++--- src/dynamics/utils/dynconst.F90 | 5 +++++ 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 5d862b63..e783f519 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -19,7 +19,8 @@ module dyn_comp 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 + 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 @@ -56,7 +57,6 @@ module dyn_comp public :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels public :: ncells_global, nedges_global, nvertices_global, ncells_max, nedges_max public :: sphere_radius - public :: deg_to_rad, rad_to_deg !> NOTE: !> This derived type is not used by MPAS dynamical core. It exists only as a placeholder because CAM-SIMA requires it. @@ -79,9 +79,6 @@ module dyn_comp 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_r8), parameter :: deg_to_rad = constant_pi / 180.0_kind_r8 ! Convert degrees to radians. - real(kind_r8), parameter :: rad_to_deg = 180.0_kind_r8 / constant_pi ! Convert radians to degrees. 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. diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index e10d6eb5..d32e1dbb 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -10,9 +10,8 @@ module dyn_grid 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, & - deg_to_rad, rad_to_deg - use dynconst, only: constant_pi => pi, dynconst_init + sphere_radius + use dynconst, only: 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 diff --git a/src/dynamics/utils/dynconst.F90 b/src/dynamics/utils/dynconst.F90 index 84574749..78c5b04c 100644 --- a/src/dynamics/utils/dynconst.F90 +++ b/src/dynamics/utils/dynconst.F90 @@ -18,6 +18,11 @@ module dynconst !circle's circumference/diameter [unitless] real(kind_dyn), parameter, public :: pi = real(phys_pi, kind_dyn) + ! Convert degrees to radians + real(kind_dyn), parameter, public :: deg_to_rad = pi / 180.0_kind_dyn + ! Convert radians to degrees + real(kind_dyn), parameter, public :: rad_to_deg = 180.0_kind_dyn / pi + ! radius of earth [m] real(kind_dyn), protected, public :: rearth ! reciprocal of earth's radius [1/m] From 155b51c4e3f01e2b6114979d2c3051ff7cb27f55 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 24 Jul 2024 11:53:07 -0600 Subject: [PATCH 24/31] Use `max` intrinsic when determining the number of constituents Same logic, but the statement is now more compact. --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index ed7d11f1..dbe65730 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -771,11 +771,7 @@ subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file) ! In MPAS, there must be at least one constituent, `qv`, which denotes water vapor mixing ratio. ! Because MPAS has some hard-coded array accesses through the `index_qv` index, it will crash ! (i.e., segmentation fault due to invalid memory access) if `qv` is not allocated. - if (number_of_constituents < 1) then - self % number_of_constituents = 1 - else - self % number_of_constituents = number_of_constituents - end if + self % number_of_constituents = max(1, number_of_constituents) call self % debug_print('Number of constituents is ', [self % number_of_constituents]) From b2ef93e070812e3b2f9a42ec779865762985e97b Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 24 Jul 2024 12:26:25 -0600 Subject: [PATCH 25/31] Use whole array section when assigning arrays The statement is now clearer on its intent. Other similar instances in the code already follow this style. --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index dbe65730..8daa07da 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -921,7 +921,7 @@ subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species) call self % model_error('Failed to allocate is_water_species', subname, __LINE__) end if - self % is_water_species(:) = is_water_species + self % is_water_species(:) = is_water_species(:) if (size(self % constituent_name) /= size(index_unique(self % constituent_name))) then call self % model_error('Constituent names must be unique', subname, __LINE__) From 03768738e0dbcc08112a8cdca9682ed59e7d18fe Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 24 Jul 2024 12:35:58 -0600 Subject: [PATCH 26/31] Make status check on reading surface geopotential more intuitive --- src/dynamics/mpas/dyn_comp.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index e783f519..e6776216 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -318,12 +318,12 @@ subroutine check_topography_data(pio_file) call cam_read_field('PHIS', pio_file, surface_geopotential, success, & gridname='cam_cell', timelevel=1, log_output=(debug_output > debugout_none)) - if (success) then - surface_geometric_height(:) = surface_geopotential(:) / constant_g - else - call endrun('Failed to find variable "PHIS"', subname, __LINE__) + if (.not. success) then + call endrun('Failed to find variable "PHIS"', subname, __LINE__) end if + 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 call endrun('Surface geometric height in MPAS is not consistent with topography data', subname, __LINE__) From 3308a26ac86d5162b0369cecca487bef6454dfa1 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 24 Jul 2024 12:45:51 -0600 Subject: [PATCH 27/31] Switch to use PIO type constants to check variable types --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 8daa07da..0bede481 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -11,8 +11,8 @@ module dyn_mpas_subdriver ! Modules from external libraries. use mpi, only: mpi_comm_null, mpi_comm_rank, mpi_success - use netcdf, only: nf90_char, nf90_int, nf90_float, nf90_double - use pio, only: file_desc_t, iosystem_desc_t, pio_file_is_open, pio_iosystem_is_active, & + 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. @@ -2164,19 +2164,19 @@ subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compa select case (trim(adjustl(var_info % type))) case ('character') - if (vartype /= nf90_char) then + if (vartype /= pio_char) then cycle end if case ('integer') - if (vartype /= nf90_int) then + if (vartype /= pio_int) then cycle end if case ('real') - if (rkind == r4kind .and. vartype /= nf90_float) then + if (rkind == r4kind .and. vartype /= pio_real) then cycle end if - if (rkind == r8kind .and. vartype /= nf90_double) then + if (rkind == r8kind .and. vartype /= pio_double) then cycle end if case default From 149cbca853d77f7ba5dbc64ac4c14dc5c2b9d8e0 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Wed, 24 Jul 2024 14:06:39 -0600 Subject: [PATCH 28/31] Directly pass subroutine name constant to `check_allocate` --- src/dynamics/mpas/dyn_comp.F90 | 42 +++++++++++++++------------------- src/dynamics/mpas/dyn_grid.F90 | 35 +++++++++++++++------------- 2 files changed, 37 insertions(+), 40 deletions(-) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index e6776216..5eaf6fed 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -196,10 +196,10 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) type(file_desc_t), pointer :: pio_topo_file => null() allocate(constituent_name(num_advected), stat=ierr) - call check_allocate(ierr, 'dyn_init', 'constituent_name(num_advected)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__) allocate(is_water_species(num_advected), stat=ierr) - call check_allocate(ierr, 'dyn_init', 'is_water_species(num_advected)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'is_water_species(num_advected)', 'dyn_comp', __LINE__) do i = 1, num_advected constituent_name(i) = const_name(i) @@ -307,10 +307,10 @@ subroutine check_topography_data(pio_file) end if allocate(surface_geopotential(ncells_solve), stat=ierr) - call check_allocate(ierr, 'check_topography_data', 'surface_geopotential(ncells_solve)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'surface_geopotential(ncells_solve)', 'dyn_comp', __LINE__) allocate(surface_geometric_height(ncells_solve), stat=ierr) - call check_allocate(ierr, 'check_topography_data', 'surface_geometric_height(ncells_solve)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'surface_geometric_height(ncells_solve)', 'dyn_comp', __LINE__) surface_geopotential(:) = 0.0_kind_r8 surface_geometric_height(:) = 0.0_kind_r8 @@ -381,7 +381,7 @@ subroutine init_shared_variable() call dyn_debug_print('Preparing to set analytic initial condition') allocate(global_grid_index(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'global_grid_index(ncells_solve)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'global_grid_index(ncells_solve)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(indextocellid, 'mesh', 'indexToCellID') @@ -390,10 +390,10 @@ subroutine init_shared_variable() nullify(indextocellid) allocate(lat_rad(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'lat_rad(ncells_solve)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'lat_rad(ncells_solve)', 'dyn_comp', __LINE__) allocate(lon_rad(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'lon_rad(ncells_solve)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'lon_rad(ncells_solve)', 'dyn_comp', __LINE__) ! "mpas_cell" is a registered grid name that is defined in `dyn_grid`. lat_deg => cam_grid_get_latvals(cam_grid_id('mpas_cell')) @@ -414,7 +414,7 @@ subroutine init_shared_variable() nullify(lon_deg) allocate(z_int(ncells_solve, pverp), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'z_int(ncells_solve, pverp)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'z_int(ncells_solve, pverp)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid') @@ -435,8 +435,7 @@ subroutine set_mpas_state_u() call dyn_debug_print('Setting MPAS state "u"') allocate(buffer_2d_real(ncells_solve, pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_2d_real(ncells_solve, pver)', & - 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'buffer_2d_real(ncells_solve, pver)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(ucellzonal, 'diag', 'uReconstructZonal') call mpas_dynamical_core % get_variable_pointer(ucellmeridional, 'diag', 'uReconstructMeridional') @@ -500,12 +499,10 @@ subroutine set_mpas_state_scalars() call dyn_debug_print('Setting MPAS state "scalars"') allocate(buffer_3d_real(ncells_solve, pver, num_advected), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_3d_real(ncells_solve, pver, num_advected)', & - 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'buffer_3d_real(ncells_solve, pver, num_advected)', 'dyn_comp', __LINE__) allocate(constituent_index(num_advected), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'constituent_index(num_advected)', & - 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'constituent_index(num_advected)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv') call mpas_dynamical_core % get_variable_pointer(scalars, 'state', 'scalars', time_level=1) @@ -568,20 +565,17 @@ subroutine set_mpas_state_rho_theta() call dyn_debug_print('Setting MPAS state "rho" and "theta"') allocate(p_sfc(ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'p_sfc(ncells_solve)', & - 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'p_sfc(ncells_solve)', 'dyn_comp', __LINE__) p_sfc(:) = 0.0_kind_r8 call dyn_set_inic_col(vc_height, lat_rad, lon_rad, global_grid_index, zint=z_int, ps=p_sfc) allocate(buffer_2d_real(ncells_solve, pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'buffer_2d_real(ncells_solve, pver)', & - 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'buffer_2d_real(ncells_solve, pver)', 'dyn_comp', __LINE__) allocate(t_mid(pver, ncells_solve), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 't_mid(pver, ncells_solve)', & - 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 't_mid(pver, ncells_solve)', 'dyn_comp', __LINE__) buffer_2d_real(:, :) = 0.0_kind_r8 @@ -595,13 +589,13 @@ subroutine set_mpas_state_rho_theta() deallocate(buffer_2d_real) allocate(p_mid_col(pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'p_mid_col(pver)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'p_mid_col(pver)', 'dyn_comp', __LINE__) allocate(qv_mid_col(pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'qv_mid_col(pver)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'qv_mid_col(pver)', 'dyn_comp', __LINE__) allocate(tm_mid_col(pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'tm_mid_col(pver)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'tm_mid_col(pver)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(index_qv, 'dim', 'index_qv') call mpas_dynamical_core % get_variable_pointer(rho, 'diag', 'rho') @@ -670,7 +664,7 @@ subroutine set_mpas_state_rho_base_theta_base() call dyn_debug_print('Setting MPAS state "rho_base" and "theta_base"') allocate(p_base(pver), stat=ierr) - call check_allocate(ierr, 'set_analytic_initial_condition', 'p_base(pver)', 'dyn_comp', __LINE__) + call check_allocate(ierr, subname, 'p_base(pver)', 'dyn_comp', __LINE__) call mpas_dynamical_core % get_variable_pointer(rho_base, 'diag', 'rho_base') call mpas_dynamical_core % get_variable_pointer(theta_base, 'diag', 'theta_base') diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index d32e1dbb..16598eea 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -124,6 +124,7 @@ end subroutine model_grid_init !> Initialize reference pressure by computing necessary variables and calling `ref_pres_init`. !> (KCW, 2024-03-25) subroutine init_reference_pressure() + character(*), parameter :: subname = 'dyn_grid::init_reference_pressure' ! Number of pure pressure levels at model top. integer, parameter :: num_pure_p_lev = 0 integer :: ierr @@ -143,16 +144,16 @@ subroutine init_reference_pressure() call mpas_dynamical_core % get_variable_pointer(rdzw, 'mesh', 'rdzw') allocate(dzw(pver), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'dzw(pver)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'dzw(pver)', 'dyn_grid', __LINE__) dzw(:) = 1.0_kind_r8 / rdzw nullify(rdzw) allocate(zw(pverp), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'zw(pverp)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'zw(pverp)', 'dyn_grid', __LINE__) allocate(zu(pver), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'zu(pver)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'zu(pver)', 'dyn_grid', __LINE__) ! In MPAS, zeta coordinates are stored in increasing order (i.e., bottom to top of atmosphere). ! In CAM-SIMA, however, index order is reversed (i.e., top to bottom of atmosphere). @@ -166,12 +167,12 @@ subroutine init_reference_pressure() ! Compute reference pressure from reference height. allocate(p_ref_int(pverp), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'p_ref_int(pverp)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'p_ref_int(pverp)', 'dyn_grid', __LINE__) call std_atm_pres(zw, p_ref_int) allocate(p_ref_mid(pver), stat=ierr) - call check_allocate(ierr, 'init_reference_pressure', 'p_ref_mid(pver)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'p_ref_mid(pver)', 'dyn_grid', __LINE__) p_ref_mid(:) = 0.5_kind_r8 * (p_ref_int(1:pver) + p_ref_int(2:pverp)) @@ -202,6 +203,7 @@ 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() + character(*), parameter :: subname = 'dyn_grid::init_physics_grid' character(max_hcoordname_len), allocatable :: dyn_attribute_name(:) integer :: hdim1_d, hdim2_d integer :: i @@ -223,7 +225,7 @@ subroutine init_physics_grid() call mpas_dynamical_core % get_variable_pointer(loncell, 'mesh', 'lonCell') allocate(dyn_column(ncells_solve), stat=ierr) - call check_allocate(ierr, 'init_physics_grid', 'dyn_column(ncells_solve)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'dyn_column(ncells_solve)', 'dyn_grid', __LINE__) do i = 1, ncells_solve ! Column information. @@ -248,7 +250,7 @@ subroutine init_physics_grid() dyn_column(i) % local_dyn_block = i ! `dyn_block_index` is not used due to no dynamics block offset, but it still needs to be allocated. allocate(dyn_column(i) % dyn_block_index(0), stat=ierr) - call check_allocate(ierr, 'init_physics_grid', 'dyn_column(' // stringify([i]) // ') % dyn_block_index(0)', & + call check_allocate(ierr, subname, 'dyn_column(' // stringify([i]) // ') % dyn_block_index(0)', & 'dyn_grid', __LINE__) end do @@ -260,7 +262,7 @@ subroutine init_physics_grid() ! `phys_grid_init` expects to receive the `area` attribute from dynamics. ! However, do not let it because dynamics grid is different from physics grid. allocate(dyn_attribute_name(0), stat=ierr) - call check_allocate(ierr, 'init_physics_grid', 'dyn_attribute_name(0)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'dyn_attribute_name(0)', 'dyn_grid', __LINE__) call phys_grid_init(hdim1_d, hdim2_d, 'mpas', dyn_column, 'mpas_cell', dyn_attribute_name) end subroutine init_physics_grid @@ -273,6 +275,7 @@ end subroutine init_physics_grid !> * "mpas_vertex": Grid that is centered at MPAS "vertex" points. !> (KCW, 2024-03-28) subroutine define_cam_grid() + character(*), parameter :: subname = 'dyn_grid::define_cam_grid' integer :: i integer :: ierr integer, pointer :: indextocellid(:) => null() ! Global indexes of cell centers. @@ -316,7 +319,7 @@ subroutine define_cam_grid() call mpas_dynamical_core % get_variable_pointer(loncell, 'mesh', 'lonCell') allocate(global_grid_index(ncells_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_index(ncells_solve)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'global_grid_index(ncells_solve)', 'dyn_grid', __LINE__) global_grid_index(:) = int(indextocellid(1:ncells_solve), kind_imap) @@ -326,11 +329,11 @@ subroutine define_cam_grid() 1, ncells_solve, loncell * rad_to_deg, map=global_grid_index) allocate(cell_area(ncells_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'cell_area(ncells_solve)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'cell_area(ncells_solve)', 'dyn_grid', __LINE__) allocate(cell_weight(ncells_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'cell_weight(ncells_solve)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'cell_weight(ncells_solve)', 'dyn_grid', __LINE__) allocate(global_grid_map(3, ncells_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_map(3, ncells_solve)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'global_grid_map(3, ncells_solve)', 'dyn_grid', __LINE__) do i = 1, ncells_solve cell_area(i) = areacell(i) @@ -387,7 +390,7 @@ subroutine define_cam_grid() call mpas_dynamical_core % get_variable_pointer(lonedge, 'mesh', 'lonEdge') allocate(global_grid_index(nedges_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_index(nedges_solve)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'global_grid_index(nedges_solve)', 'dyn_grid', __LINE__) global_grid_index(:) = int(indextoedgeid(1:nedges_solve), kind_imap) @@ -397,7 +400,7 @@ subroutine define_cam_grid() 1, nedges_solve, lonedge * rad_to_deg, map=global_grid_index) allocate(global_grid_map(3, nedges_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_map(3, nedges_solve)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'global_grid_map(3, nedges_solve)', 'dyn_grid', __LINE__) do i = 1, nedges_solve global_grid_map(1, i) = int(i, kind_imap) @@ -428,7 +431,7 @@ subroutine define_cam_grid() call mpas_dynamical_core % get_variable_pointer(lonvertex, 'mesh', 'lonVertex') allocate(global_grid_index(nvertices_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_index(nvertices_solve)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'global_grid_index(nvertices_solve)', 'dyn_grid', __LINE__) global_grid_index(:) = int(indextovertexid(1:nvertices_solve), kind_imap) @@ -438,7 +441,7 @@ subroutine define_cam_grid() 1, nvertices_solve, lonvertex * rad_to_deg, map=global_grid_index) allocate(global_grid_map(3, nvertices_solve), stat=ierr) - call check_allocate(ierr, 'define_cam_grid', 'global_grid_map(3, nvertices_solve)', 'dyn_grid', __LINE__) + call check_allocate(ierr, subname, 'global_grid_map(3, nvertices_solve)', 'dyn_grid', __LINE__) do i = 1, nvertices_solve global_grid_map(1, i) = int(i, kind_imap) From a731e9769498b5a3b374c2329aaf97eaec2ff853 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 26 Jul 2024 12:33:58 -0600 Subject: [PATCH 29/31] Avoid explicit initialization of pointers According to "8.4 Initialization", it states "Explicit initialization of a variable that is not in a common block implies the SAVE attribute". This also applies to pointer variables like here. Avoid the implied SAVE attribute by removing explicit initialization. Uphold the principle of least astonishment. --- .../mpas/driver/dyn_mpas_subdriver.F90 | 342 +++++++++++------- src/dynamics/mpas/dyn_comp.F90 | 80 ++-- src/dynamics/mpas/dyn_grid.F90 | 96 ++--- 3 files changed, 326 insertions(+), 192 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 0bede481..2e309671 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -555,12 +555,15 @@ subroutine dyn_mpas_read_namelist(self, namelist_path, & character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_read_namelist' character(strkind) :: mpas_calendar - character(strkind), pointer :: config_pointer_c => null() + character(strkind), pointer :: config_pointer_c integer :: ierr - logical, pointer :: config_pointer_l => null() + logical, pointer :: config_pointer_l call self % debug_print(subname // ' entered') + nullify(config_pointer_c) + nullify(config_pointer_l) + call self % debug_print('Reading namelist at ', [namelist_path]) ! Override namelist filename so that we can rely on upstream MPAS functionality for reading its own namelist. @@ -763,11 +766,14 @@ subroutine dyn_mpas_init_phase3(self, number_of_constituents, pio_file) character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase3' character(strkind) :: mesh_filename integer :: mesh_format - integer, pointer :: num_scalars => null() - type(mpas_pool_type), pointer :: mpas_pool => null() + integer, pointer :: num_scalars + type(mpas_pool_type), pointer :: mpas_pool call self % debug_print(subname // ' entered') + nullify(mpas_pool) + nullify(num_scalars) + ! In MPAS, there must be at least one constituent, `qv`, which denotes water vapor mixing ratio. ! Because MPAS has some hard-coded array accesses through the `index_qv` index, it will crash ! (i.e., segmentation fault due to invalid memory access) if `qv` is not allocated. @@ -866,11 +872,14 @@ subroutine dyn_mpas_define_scalar(self, constituent_name, is_water_species) integer :: i, j, ierr integer :: index_qv, index_water_start, index_water_end integer :: time_level - type(field3dreal), pointer :: field_3d_real => null() - type(mpas_pool_type), pointer :: mpas_pool => null() + type(field3dreal), pointer :: field_3d_real + type(mpas_pool_type), pointer :: mpas_pool call self % debug_print(subname // ' entered') + nullify(field_3d_real) + nullify(mpas_pool) + if (self % number_of_constituents == 0) then call self % model_error('Constituents must be allocated before being defined', subname, __LINE__) end if @@ -1111,12 +1120,15 @@ subroutine dyn_mpas_read_write_stream(self, pio_file, stream_mode, stream_name) character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_read_write_stream' integer :: i, ierr - type(mpas_pool_type), pointer :: mpas_pool => null() - type(mpas_stream_type), pointer :: mpas_stream => null() + type(mpas_pool_type), pointer :: mpas_pool + type(mpas_stream_type), pointer :: mpas_stream type(var_info_type), allocatable :: var_info_list(:) call self % debug_print(subname // ' entered') + nullify(mpas_pool) + nullify(mpas_stream) + call self % debug_print('Initializing stream "' // trim(adjustl(stream_name)) // '"') call self % init_stream_with_pool(mpas_pool, mpas_stream, pio_file, stream_mode, stream_name) @@ -1226,22 +1238,35 @@ subroutine dyn_mpas_init_stream_with_pool(self, mpas_pool, mpas_stream, pio_file logical, allocatable :: var_is_present(:) ! Whether a variable is type, kind and rank compatible with what MPAS expects on the file (i.e., `pio_file`). logical, allocatable :: var_is_tkr_compatible(:) - type(field0dchar), pointer :: field_0d_char => null() - type(field1dchar), pointer :: field_1d_char => null() - type(field0dinteger), pointer :: field_0d_integer => null() - type(field1dinteger), pointer :: field_1d_integer => null() - type(field2dinteger), pointer :: field_2d_integer => null() - type(field3dinteger), pointer :: field_3d_integer => null() - type(field0dreal), pointer :: field_0d_real => null() - type(field1dreal), pointer :: field_1d_real => null() - type(field2dreal), pointer :: field_2d_real => null() - type(field3dreal), pointer :: field_3d_real => null() - type(field4dreal), pointer :: field_4d_real => null() - type(field5dreal), pointer :: field_5d_real => null() + type(field0dchar), pointer :: field_0d_char + type(field1dchar), pointer :: field_1d_char + type(field0dinteger), pointer :: field_0d_integer + type(field1dinteger), pointer :: field_1d_integer + type(field2dinteger), pointer :: field_2d_integer + type(field3dinteger), pointer :: field_3d_integer + type(field0dreal), pointer :: field_0d_real + type(field1dreal), pointer :: field_1d_real + type(field2dreal), pointer :: field_2d_real + type(field3dreal), pointer :: field_3d_real + type(field4dreal), pointer :: field_4d_real + type(field5dreal), pointer :: field_5d_real type(var_info_type), allocatable :: var_info_list(:) call self % debug_print(subname // ' entered') + nullify(field_0d_char) + nullify(field_1d_char) + nullify(field_0d_integer) + nullify(field_1d_integer) + nullify(field_2d_integer) + nullify(field_3d_integer) + nullify(field_0d_real) + nullify(field_1d_real) + nullify(field_2d_real) + nullify(field_3d_real) + nullify(field_4d_real) + nullify(field_5d_real) + call mpas_pool_create_pool(mpas_pool) allocate(mpas_stream, stat=ierr) @@ -1829,21 +1854,34 @@ subroutine dyn_mpas_check_variable_status(self, var_is_present, var_is_tkr_compa character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_check_variable_status' character(strkind), allocatable :: var_name_list(:) integer :: i, ierr, varid, varndims, vartype - type(field0dchar), pointer :: field_0d_char => null() - type(field1dchar), pointer :: field_1d_char => null() - type(field0dinteger), pointer :: field_0d_integer => null() - type(field1dinteger), pointer :: field_1d_integer => null() - type(field2dinteger), pointer :: field_2d_integer => null() - type(field3dinteger), pointer :: field_3d_integer => null() - type(field0dreal), pointer :: field_0d_real => null() - type(field1dreal), pointer :: field_1d_real => null() - type(field2dreal), pointer :: field_2d_real => null() - type(field3dreal), pointer :: field_3d_real => null() - type(field4dreal), pointer :: field_4d_real => null() - type(field5dreal), pointer :: field_5d_real => null() + type(field0dchar), pointer :: field_0d_char + type(field1dchar), pointer :: field_1d_char + type(field0dinteger), pointer :: field_0d_integer + type(field1dinteger), pointer :: field_1d_integer + type(field2dinteger), pointer :: field_2d_integer + type(field3dinteger), pointer :: field_3d_integer + type(field0dreal), pointer :: field_0d_real + type(field1dreal), pointer :: field_1d_real + type(field2dreal), pointer :: field_2d_real + type(field3dreal), pointer :: field_3d_real + type(field4dreal), pointer :: field_4d_real + type(field5dreal), pointer :: field_5d_real call self % debug_print(subname // ' entered') + nullify(field_0d_char) + nullify(field_1d_char) + nullify(field_0d_integer) + nullify(field_1d_integer) + nullify(field_2d_integer) + nullify(field_3d_integer) + nullify(field_0d_real) + nullify(field_1d_real) + nullify(field_2d_real) + nullify(field_3d_real) + nullify(field_4d_real) + nullify(field_5d_real) + ! Extract a list of variable names to check on the file. ! For an ordinary variable, this list just contains its name. ! For a variable array, this list contains the names of its constituent parts. @@ -2222,18 +2260,27 @@ subroutine dyn_mpas_exchange_halo(self, field_name) character(*), intent(in) :: field_name character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_exchange_halo' - type(field1dinteger), pointer :: field_1d_integer => null() - type(field2dinteger), pointer :: field_2d_integer => null() - type(field3dinteger), pointer :: field_3d_integer => null() - type(field1dreal), pointer :: field_1d_real => null() - type(field2dreal), pointer :: field_2d_real => null() - type(field3dreal), pointer :: field_3d_real => null() - type(field4dreal), pointer :: field_4d_real => null() - type(field5dreal), pointer :: field_5d_real => null() + type(field1dinteger), pointer :: field_1d_integer + type(field2dinteger), pointer :: field_2d_integer + type(field3dinteger), pointer :: field_3d_integer + type(field1dreal), pointer :: field_1d_real + type(field2dreal), pointer :: field_2d_real + type(field3dreal), pointer :: field_3d_real + type(field4dreal), pointer :: field_4d_real + type(field5dreal), pointer :: field_5d_real type(mpas_pool_field_info_type) :: mpas_pool_field_info call self % debug_print(subname // ' entered') + nullify(field_1d_integer) + nullify(field_2d_integer) + nullify(field_3d_integer) + nullify(field_1d_real) + nullify(field_2d_real) + nullify(field_3d_real) + nullify(field_4d_real) + nullify(field_5d_real) + call self % debug_print('Inquiring field information for "' // trim(adjustl(field_name)) // '"') call mpas_pool_get_field_info(self % domain_ptr % blocklist % allfields, & @@ -2396,15 +2443,20 @@ subroutine dyn_mpas_compute_unit_vector(self) character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_unit_vector' integer :: i - integer, pointer :: ncells => null() - real(rkind), pointer :: latcell(:) => null() - real(rkind), pointer :: loncell(:) => null() - real(rkind), pointer :: east(:, :) => null() - real(rkind), pointer :: north(:, :) => null() - type(mpas_pool_type), pointer :: mpas_pool => null() + integer, pointer :: ncells + real(rkind), pointer :: latcell(:), loncell(:) + real(rkind), pointer :: east(:, :), north(:, :) + type(mpas_pool_type), pointer :: mpas_pool call self % debug_print(subname // ' entered') + nullify(ncells) + nullify(latcell, loncell) + + nullify(east, north) + + nullify(mpas_pool) + ! Input. call self % get_variable_pointer(ncells, 'dim', 'nCells') call self % get_variable_pointer(latcell, 'mesh', 'latCell') @@ -2429,11 +2481,9 @@ subroutine dyn_mpas_compute_unit_vector(self) end do nullify(ncells) - nullify(latcell) - nullify(loncell) + nullify(latcell, loncell) - nullify(east) - nullify(north) + nullify(east, north) call self % get_pool_pointer(mpas_pool, 'mesh') call mpas_initialize_vectors(mpas_pool) @@ -2467,17 +2517,25 @@ subroutine dyn_mpas_compute_edge_wind(self) character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_compute_edge_wind' integer :: cell1, cell2, i - integer, pointer :: cellsonedge(:, :) => null() - integer, pointer :: nedges => null() - real(rkind), pointer :: east(:, :) => null() - real(rkind), pointer :: edgenormalvectors(:, :) => null() - real(rkind), pointer :: north(:, :) => null() - real(rkind), pointer :: ucellmeridional(:, :) => null() - real(rkind), pointer :: ucellzonal(:, :) => null() - real(rkind), pointer :: uedge(:, :) => null() + integer, pointer :: cellsonedge(:, :) + integer, pointer :: nedges + real(rkind), pointer :: east(:, :), north(:, :) + real(rkind), pointer :: edgenormalvectors(:, :) + real(rkind), pointer :: ucellzonal(:, :), ucellmeridional(:, :) + real(rkind), pointer :: uedge(:, :) call self % debug_print(subname // ' entered') + nullify(nedges) + + nullify(ucellzonal, ucellmeridional) + + nullify(cellsonedge) + nullify(east, north) + nullify(edgenormalvectors) + + nullify(uedge) + ! Make sure halo layers are up-to-date before computation. call self % exchange_halo('uReconstructZonal') call self % exchange_halo('uReconstructMeridional') @@ -2490,8 +2548,8 @@ subroutine dyn_mpas_compute_edge_wind(self) call self % get_variable_pointer(cellsonedge, 'mesh', 'cellsOnEdge') call self % get_variable_pointer(east, 'mesh', 'east') - call self % get_variable_pointer(edgenormalvectors, 'mesh', 'edgeNormalVectors') call self % get_variable_pointer(north, 'mesh', 'north') + call self % get_variable_pointer(edgenormalvectors, 'mesh', 'edgeNormalVectors') ! Output. call self % get_variable_pointer(uedge, 'state', 'u', time_level=1) @@ -2516,13 +2574,11 @@ subroutine dyn_mpas_compute_edge_wind(self) nullify(nedges) - nullify(ucellzonal) - nullify(ucellmeridional) + nullify(ucellzonal, ucellmeridional) nullify(cellsonedge) - nullify(east) + nullify(east, north) nullify(edgenormalvectors) - nullify(north) nullify(uedge) @@ -2553,18 +2609,26 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval) character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_init_phase4' character(strkind) :: date_time - character(strkind), pointer :: initial_time_1_pointer => null(), initial_time_2_pointer => null() - character(strkind), pointer :: xtime_pointer => null() + character(strkind), pointer :: initial_time_1, initial_time_2 + character(strkind), pointer :: xtime integer :: ierr - integer, pointer :: nvertlevels => null(), maxedges => null(), maxedges2 => null(), num_scalars => null() - logical, pointer :: config_do_restart => null() - real(rkind), pointer :: config_dt => null() - type(field0dreal), pointer :: field_0d_real => null() - type(mpas_pool_type), pointer :: mpas_pool => null() + integer, pointer :: nvertlevels, maxedges, maxedges2, num_scalars + logical, pointer :: config_do_restart + real(rkind), pointer :: config_dt + type(field0dreal), pointer :: field_0d_real + type(mpas_pool_type), pointer :: mpas_pool type(mpas_time_type) :: mpas_time call self % debug_print(subname // ' entered') + nullify(initial_time_1, initial_time_2) + nullify(xtime) + nullify(nvertlevels, maxedges, maxedges2, num_scalars) + nullify(config_do_restart) + nullify(config_dt) + nullify(field_0d_real) + nullify(mpas_pool) + if (real(coupling_time_interval, rkind) < 1.0_rkind) then call self % model_error('Invalid coupling time interval ' // stringify([real(coupling_time_interval, rkind)]), & subname, __LINE__) @@ -2615,10 +2679,7 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval) call mpas_atm_set_dims(nvertlevels, maxedges, maxedges2, num_scalars) - nullify(nvertlevels) - nullify(maxedges) - nullify(maxedges2) - nullify(num_scalars) + nullify(nvertlevels, maxedges, maxedges2, num_scalars) ! Build halo exchange groups and set the `exchange_halo_group` procedure pointer, which is used to ! exchange the halo layers of all fields in the named group. @@ -2685,18 +2746,18 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval) nullify(mpas_pool) nullify(config_dt) - call self % get_variable_pointer(xtime_pointer, 'state', 'xtime', time_level=1) + call self % get_variable_pointer(xtime, 'state', 'xtime', time_level=1) - xtime_pointer = date_time + xtime = date_time - nullify(xtime_pointer) + nullify(xtime) ! Initialize `initial_time` in the second time level. We need to do this manually because initial states ! are read into time level 1, and if we write anything from time level 2, `initial_time` will be invalid. - call self % get_variable_pointer(initial_time_1_pointer, 'state', 'initial_time', time_level=1) - call self % get_variable_pointer(initial_time_2_pointer, 'state', 'initial_time', time_level=2) + call self % get_variable_pointer(initial_time_1, 'state', 'initial_time', time_level=1) + call self % get_variable_pointer(initial_time_2, 'state', 'initial_time', time_level=2) - initial_time_2_pointer = initial_time_1_pointer + initial_time_2 = initial_time_1 ! Set time units to CF-compliant "seconds since ". call self % get_pool_pointer(mpas_pool, 'state') @@ -2708,14 +2769,13 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval) end if call mpas_modify_att(field_0d_real % attlists(1) % attlist, 'units', & - 'seconds since ' // mpas_string_replace(initial_time_1_pointer, '_', ' '), ierr=ierr) + 'seconds since ' // mpas_string_replace(initial_time_1, '_', ' '), ierr=ierr) if (ierr /= 0) then call self % model_error('Failed to set time units', subname, __LINE__) end if - nullify(initial_time_1_pointer) - nullify(initial_time_2_pointer) + nullify(initial_time_1, initial_time_2) nullify(mpas_pool) nullify(field_0d_real) @@ -2939,13 +2999,21 @@ subroutine dyn_mpas_get_local_mesh_dimension(self, & integer, intent(out) :: ncells, ncells_solve, nedges, nedges_solve, nvertices, nvertices_solve, nvertlevels character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_local_mesh_dimension' - integer, pointer :: ncells_pointer => null() - integer, pointer :: ncellssolve_pointer => null() - integer, pointer :: nedges_pointer => null() - integer, pointer :: nedgessolve_pointer => null() - integer, pointer :: nvertices_pointer => null() - integer, pointer :: nverticessolve_pointer => null() - integer, pointer :: nvertlevels_pointer => null() + integer, pointer :: ncells_pointer + integer, pointer :: ncellssolve_pointer + integer, pointer :: nedges_pointer + integer, pointer :: nedgessolve_pointer + integer, pointer :: nvertices_pointer + integer, pointer :: nverticessolve_pointer + integer, pointer :: nvertlevels_pointer + + nullify(ncells_pointer) + nullify(ncellssolve_pointer) + nullify(nedges_pointer) + nullify(nedgessolve_pointer) + nullify(nvertices_pointer) + nullify(nverticessolve_pointer) + nullify(nvertlevels_pointer) call self % get_variable_pointer(ncells_pointer, 'dim', 'nCells') call self % get_variable_pointer(ncellssolve_pointer, 'dim', 'nCellsSolve') @@ -2999,11 +3067,17 @@ subroutine dyn_mpas_get_global_mesh_dimension(self, & real(rkind), intent(out) :: sphere_radius character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_global_mesh_dimension' - integer, pointer :: maxedges_pointer => null() - integer, pointer :: ncellssolve_pointer => null() - integer, pointer :: nedgessolve_pointer => null() - integer, pointer :: nverticessolve_pointer => null() - integer, pointer :: nvertlevels_pointer => null() + integer, pointer :: maxedges_pointer + integer, pointer :: ncellssolve_pointer + integer, pointer :: nedgessolve_pointer + integer, pointer :: nverticessolve_pointer + integer, pointer :: nvertlevels_pointer + + nullify(maxedges_pointer) + nullify(ncellssolve_pointer) + nullify(nedgessolve_pointer) + nullify(nverticessolve_pointer) + nullify(nvertlevels_pointer) call self % get_variable_pointer(maxedges_pointer, 'dim', 'maxEdges') call self % get_variable_pointer(ncellssolve_pointer, 'dim', 'nCellsSolve') @@ -3087,8 +3161,9 @@ subroutine dyn_mpas_get_variable_pointer_c0(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_c0' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) @@ -3114,8 +3189,9 @@ subroutine dyn_mpas_get_variable_pointer_c1(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_c1' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level) @@ -3135,8 +3211,9 @@ subroutine dyn_mpas_get_variable_pointer_i0(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_i0' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) @@ -3165,8 +3242,9 @@ subroutine dyn_mpas_get_variable_pointer_i1(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_i1' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) @@ -3192,8 +3270,9 @@ subroutine dyn_mpas_get_variable_pointer_i2(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_i2' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level) @@ -3213,8 +3292,9 @@ subroutine dyn_mpas_get_variable_pointer_i3(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_i3' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level) @@ -3234,8 +3314,9 @@ subroutine dyn_mpas_get_variable_pointer_l0(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_l0' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) @@ -3259,8 +3340,9 @@ subroutine dyn_mpas_get_variable_pointer_r0(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r0' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) @@ -3286,8 +3368,9 @@ subroutine dyn_mpas_get_variable_pointer_r1(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r1' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level) @@ -3307,8 +3390,9 @@ subroutine dyn_mpas_get_variable_pointer_r2(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r2' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level) @@ -3328,8 +3412,9 @@ subroutine dyn_mpas_get_variable_pointer_r3(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r3' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level) @@ -3349,8 +3434,9 @@ subroutine dyn_mpas_get_variable_pointer_r4(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r4' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level) @@ -3370,8 +3456,9 @@ subroutine dyn_mpas_get_variable_pointer_r5(self, variable_pointer, pool_name, v integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_pointer_r5' - type(mpas_pool_type), pointer :: mpas_pool => null() + type(mpas_pool_type), pointer :: mpas_pool + nullify(mpas_pool) call self % get_pool_pointer(mpas_pool, pool_name) nullify(variable_pointer) call mpas_pool_get_array(mpas_pool, trim(adjustl(variable_name)), variable_pointer, timelevel=time_level) @@ -3404,9 +3491,10 @@ subroutine dyn_mpas_get_variable_value_c0(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_c0' - character(strkind), pointer :: variable_pointer => null() + character(strkind), pointer :: variable_pointer integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3425,9 +3513,10 @@ subroutine dyn_mpas_get_variable_value_c1(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_c1' - character(strkind), pointer :: variable_pointer(:) => null() + character(strkind), pointer :: variable_pointer(:) integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3446,9 +3535,10 @@ subroutine dyn_mpas_get_variable_value_i0(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_i0' - integer, pointer :: variable_pointer => null() + integer, pointer :: variable_pointer integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3467,9 +3557,10 @@ subroutine dyn_mpas_get_variable_value_i1(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_i1' - integer, pointer :: variable_pointer(:) => null() + integer, pointer :: variable_pointer(:) integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3488,9 +3579,10 @@ subroutine dyn_mpas_get_variable_value_i2(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_i2' - integer, pointer :: variable_pointer(:, :) => null() + integer, pointer :: variable_pointer(:, :) integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3509,9 +3601,10 @@ subroutine dyn_mpas_get_variable_value_i3(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_i3' - integer, pointer :: variable_pointer(:, :, :) => null() + integer, pointer :: variable_pointer(:, :, :) integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3530,9 +3623,10 @@ subroutine dyn_mpas_get_variable_value_l0(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_l0' - logical, pointer :: variable_pointer => null() + logical, pointer :: variable_pointer integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3551,9 +3645,10 @@ subroutine dyn_mpas_get_variable_value_r0(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r0' - real(rkind), pointer :: variable_pointer => null() + real(rkind), pointer :: variable_pointer integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3572,9 +3667,10 @@ subroutine dyn_mpas_get_variable_value_r1(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r1' - real(rkind), pointer :: variable_pointer(:) => null() + real(rkind), pointer :: variable_pointer(:) integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3593,9 +3689,10 @@ subroutine dyn_mpas_get_variable_value_r2(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r2' - real(rkind), pointer :: variable_pointer(:, :) => null() + real(rkind), pointer :: variable_pointer(:, :) integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3614,9 +3711,10 @@ subroutine dyn_mpas_get_variable_value_r3(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r3' - real(rkind), pointer :: variable_pointer(:, :, :) => null() + real(rkind), pointer :: variable_pointer(:, :, :) integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3635,9 +3733,10 @@ subroutine dyn_mpas_get_variable_value_r4(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r4' - real(rkind), pointer :: variable_pointer(:, :, :, :) => null() + real(rkind), pointer :: variable_pointer(:, :, :, :) integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) @@ -3656,9 +3755,10 @@ subroutine dyn_mpas_get_variable_value_r5(self, variable_value, pool_name, varia integer, optional, intent(in) :: time_level character(*), parameter :: subname = 'dyn_mpas_subdriver::dyn_mpas_get_variable_value_r5' - real(rkind), pointer :: variable_pointer(:, :, :, :, :) => null() + real(rkind), pointer :: variable_pointer(:, :, :, :, :) integer :: ierr + nullify(variable_pointer) call self % get_variable_pointer(variable_pointer, pool_name, variable_name, time_level=time_level) allocate(variable_value, source=variable_pointer, stat=ierr) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index 5eaf6fed..aec4a265 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -127,7 +127,9 @@ subroutine dyn_readnl(namelist_path) stop_date_time(6), & ! YYYY, MM, DD, hh, mm, ss. run_duration(4), & ! DD, hh, mm, ss. sec_since_midnight ! Second(s) since midnight. - type(iosystem_desc_t), pointer :: pio_iosystem => null() + type(iosystem_desc_t), pointer :: pio_iosystem + + nullify(pio_iosystem) ! Enable/disable the debug output of MPAS dynamical core according to the debug verbosity level of CAM-SIMA. mpas_dynamical_core % debug_output = (debug_output > debugout_none) @@ -192,8 +194,11 @@ subroutine dyn_init(cam_runtime_opts, dyn_in, dyn_out) integer :: i integer :: ierr logical, allocatable :: is_water_species(:) - type(file_desc_t), pointer :: pio_init_file => null() - type(file_desc_t), pointer :: pio_topo_file => null() + type(file_desc_t), pointer :: pio_init_file + type(file_desc_t), pointer :: pio_topo_file + + nullify(pio_init_file) + nullify(pio_topo_file) allocate(constituent_name(num_advected), stat=ierr) call check_allocate(ierr, subname, 'constituent_name(num_advected)', 'dyn_comp', __LINE__) @@ -295,7 +300,9 @@ 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(:, :) => null() ! From MPAS. Geometric height (meters) at layer interfaces. + real(kind_r8), pointer :: zgrid(:, :) ! From MPAS. Geometric height (meters) at layer interfaces. + + nullify(zgrid) call mpas_dynamical_core % get_variable_pointer(zgrid, 'mesh', 'zgrid') @@ -352,7 +359,7 @@ subroutine set_analytic_initial_condition() 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(:, :) => null() ! Geometric height (meters) at layer interfaces. + real(kind_r8), pointer :: zgrid(:, :) ! Geometric height (meters) at layer interfaces. ! Dimension and vertical index orders follow MPAS convention. call init_shared_variable() @@ -364,8 +371,7 @@ subroutine set_analytic_initial_condition() call set_mpas_state_rho_base_theta_base() deallocate(global_grid_index) - deallocate(lat_rad) - deallocate(lon_rad) + deallocate(lat_rad, lon_rad) deallocate(z_int) nullify(zgrid) @@ -375,11 +381,15 @@ subroutine set_analytic_initial_condition() subroutine init_shared_variable() integer :: ierr integer :: k - integer, pointer :: indextocellid(:) => null() - real(kind_r8), pointer :: lat_deg(:) => null(), lon_deg(:) => null() + integer, pointer :: indextocellid(:) + real(kind_r8), pointer :: lat_deg(:), lon_deg(:) call dyn_debug_print('Preparing to set analytic initial condition') + nullify(zgrid) + nullify(indextocellid) + nullify(lat_deg, lon_deg) + allocate(global_grid_index(ncells_solve), stat=ierr) call check_allocate(ierr, subname, 'global_grid_index(ncells_solve)', 'dyn_comp', __LINE__) @@ -410,8 +420,7 @@ subroutine init_shared_variable() lat_rad(:) = lat_deg(:) * deg_to_rad lon_rad(:) = lon_deg(:) * deg_to_rad - nullify(lat_deg) - nullify(lon_deg) + nullify(lat_deg, lon_deg) allocate(z_int(ncells_solve, pverp), stat=ierr) call check_allocate(ierr, subname, 'z_int(ncells_solve, pverp)', 'dyn_comp', __LINE__) @@ -429,11 +438,12 @@ end subroutine init_shared_variable subroutine set_mpas_state_u() integer :: ierr integer :: k - real(kind_r8), pointer :: ucellmeridional(:, :) => null() - real(kind_r8), pointer :: ucellzonal(:, :) => null() + real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :) call dyn_debug_print('Setting MPAS state "u"') + nullify(ucellzonal, ucellmeridional) + allocate(buffer_2d_real(ncells_solve, pver), stat=ierr) call check_allocate(ierr, subname, 'buffer_2d_real(ncells_solve, pver)', 'dyn_comp', __LINE__) @@ -460,8 +470,7 @@ subroutine set_mpas_state_u() deallocate(buffer_2d_real) - nullify(ucellzonal) - nullify(ucellmeridional) + nullify(ucellzonal, ucellmeridional) call mpas_dynamical_core % compute_edge_wind() end subroutine set_mpas_state_u @@ -469,10 +478,12 @@ end subroutine set_mpas_state_u !> Set MPAS state `w` (i.e., vertical velocity at cell interfaces). !> (KCW, 2024-05-13) subroutine set_mpas_state_w() - real(kind_r8), pointer :: w(:, :) => null() + real(kind_r8), pointer :: w(:, :) call dyn_debug_print('Setting MPAS state "w"') + nullify(w) + call mpas_dynamical_core % get_variable_pointer(w, 'state', 'w', time_level=1) w(:, 1:ncells_solve) = 0.0_kind_r8 @@ -493,11 +504,14 @@ subroutine set_mpas_state_scalars() integer :: i, k integer :: ierr integer, allocatable :: constituent_index(:) - integer, pointer :: index_qv => null() - real(kind_r8), pointer :: scalars(:, :, :) => null() + integer, pointer :: index_qv + real(kind_r8), pointer :: scalars(:, :, :) call dyn_debug_print('Setting MPAS state "scalars"') + nullify(index_qv) + nullify(scalars) + allocate(buffer_3d_real(ncells_solve, pver, num_advected), stat=ierr) call check_allocate(ierr, subname, 'buffer_3d_real(ncells_solve, pver, num_advected)', 'dyn_comp', __LINE__) @@ -550,7 +564,7 @@ end subroutine set_mpas_state_scalars subroutine set_mpas_state_rho_theta() integer :: i, k integer :: ierr - integer, pointer :: index_qv => null() + integer, pointer :: index_qv real(kind_r8), allocatable :: p_mid_col(:) ! Pressure (Pa) at layer midpoints of each column. This is full pressure, ! which also accounts for water vapor. real(kind_r8), allocatable :: p_sfc(:) ! Pressure (Pa) at surface. This is full pressure, @@ -558,12 +572,17 @@ subroutine set_mpas_state_rho_theta() real(kind_r8), allocatable :: qv_mid_col(:) ! Water vapor mixing ratio (kg/kg) at layer midpoints of each column. real(kind_r8), allocatable :: t_mid(:, :) ! Temperature (K) at layer midpoints. real(kind_r8), allocatable :: tm_mid_col(:) ! Modified "moist" temperature (K) at layer midpoints of each column. - real(kind_r8), pointer :: rho(:, :) => null() - real(kind_r8), pointer :: theta(:, :) => null() - real(kind_r8), pointer :: scalars(:, :, :) => null() + real(kind_r8), pointer :: rho(:, :) + real(kind_r8), pointer :: theta(:, :) + real(kind_r8), pointer :: scalars(:, :, :) call dyn_debug_print('Setting MPAS state "rho" and "theta"') + nullify(index_qv) + nullify(rho) + nullify(theta) + nullify(scalars) + allocate(p_sfc(ncells_solve), stat=ierr) call check_allocate(ierr, subname, 'p_sfc(ncells_solve)', 'dyn_comp', __LINE__) @@ -657,12 +676,16 @@ 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(:, :) => null() - real(kind_r8), pointer :: theta_base(:, :) => null() - real(kind_r8), pointer :: zz(:, :) => null() + real(kind_r8), pointer :: rho_base(:, :) + real(kind_r8), pointer :: theta_base(:, :) + real(kind_r8), pointer :: zz(:, :) call dyn_debug_print('Setting MPAS state "rho_base" and "theta_base"') + nullify(rho_base) + nullify(theta_base) + nullify(zz) + allocate(p_base(pver), stat=ierr) call check_allocate(ierr, subname, 'p_base(pver)', 'dyn_comp', __LINE__) @@ -731,11 +754,14 @@ end subroutine set_analytic_initial_condition subroutine set_default_constituent() character(*), parameter :: subname = 'dyn_comp::set_default_constituent' integer :: i, k - real(kind_phys), pointer :: constituents(:, :, :) => null() ! This points to CCPP memory. - real(kind_r8), pointer :: scalars(:, :, :) => null() ! This points to MPAS memory. + real(kind_phys), pointer :: constituents(:, :, :) ! This points to CCPP memory. + real(kind_r8), pointer :: scalars(:, :, :) ! This points to MPAS memory. call dyn_debug_print('Setting default MPAS state "scalars"') + nullify(constituents) + nullify(scalars) + constituents => cam_constituents_array() if (.not. associated(constituents)) then diff --git a/src/dynamics/mpas/dyn_grid.F90 b/src/dynamics/mpas/dyn_grid.F90 index 16598eea..26b0e345 100644 --- a/src/dynamics/mpas/dyn_grid.F90 +++ b/src/dynamics/mpas/dyn_grid.F90 @@ -52,7 +52,9 @@ module dyn_grid ! Called by `cam_init` in `src/control/cam_comp.F90`. subroutine model_grid_init() character(*), parameter :: subname = 'dyn_grid::model_grid_init' - type(file_desc_t), pointer :: pio_file => null() + type(file_desc_t), pointer :: pio_file + + nullify(pio_file) ! Initialize mathematical and physical constants for dynamics. call dyn_debug_print('Calling dynconst_init') @@ -138,7 +140,9 @@ subroutine init_reference_pressure() ! `dzw` denotes the delta/difference between `zw`. ! `rdzw` denotes the reciprocal of `dzw`. real(kind_r8), allocatable :: zu(:), zw(:), dzw(:) - real(kind_r8), pointer :: rdzw(:) => null() + real(kind_r8), pointer :: rdzw(:) + + nullify(rdzw) ! Compute reference height. call mpas_dynamical_core % get_variable_pointer(rdzw, 'mesh', 'rdzw') @@ -208,12 +212,17 @@ subroutine init_physics_grid() integer :: hdim1_d, hdim2_d integer :: i integer :: ierr - integer, pointer :: indextocellid(:) => null() ! Global indexes of cell centers. - real(kind_r8), pointer :: areacell(:) => null() ! Cell areas (square meters). - real(kind_r8), pointer :: latcell(:) => null() ! Cell center latitudes (radians). - real(kind_r8), pointer :: loncell(:) => null() ! Cell center longitudes (radians). + 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). type(physics_column_t), allocatable :: dyn_column(:) ! Grid and mapping information between global and local indexes. + nullify(areacell) + nullify(indextocellid) + nullify(latcell) + nullify(loncell) + hdim1_d = ncells_global ! Setting `hdim2_d` to `1` indicates unstructured grid. @@ -278,37 +287,47 @@ subroutine define_cam_grid() character(*), parameter :: subname = 'dyn_grid::define_cam_grid' integer :: i integer :: ierr - integer, pointer :: indextocellid(:) => null() ! Global indexes of cell centers. - integer, pointer :: indextoedgeid(:) => null() ! Global indexes of edge nodes. - integer, pointer :: indextovertexid(:) => null() ! Global indexes of vertex nodes. - real(kind_r8), pointer :: areacell(:) => null() ! Cell areas (square meters). - real(kind_r8), pointer :: latcell(:) => null() ! Cell center latitudes (radians). - real(kind_r8), pointer :: latedge(:) => null() ! Edge node latitudes (radians). - real(kind_r8), pointer :: latvertex(:) => null() ! Vertex node latitudes (radians). - real(kind_r8), pointer :: loncell(:) => null() ! Cell center longitudes (radians). - real(kind_r8), pointer :: lonedge(:) => null() ! Edge node longitudes (radians). - real(kind_r8), pointer :: lonvertex(:) => null() ! Vertex node longitudes (radians). + 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). ! Global grid indexes. CAN be safely deallocated because its values are copied into ! `cam_grid_attribute_*_t` and `horiz_coord_t`. ! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`. - integer(kind_imap), pointer :: global_grid_index(:) => null() + integer(kind_imap), pointer :: global_grid_index(:) ! Global grid maps. CANNOT be safely deallocated because `cam_filemap_t` ! just uses pointers to point at it. ! `kind_imap` is an integer kind of `PIO_OFFSET_KIND`. - integer(kind_imap), pointer :: global_grid_map(:, :) => null() + integer(kind_imap), pointer :: global_grid_map(:, :) ! Cell areas (square meters). CANNOT be safely deallocated because `cam_grid_attribute_*_t` ! just uses pointers to point at it. - real(kind_r8), pointer :: cell_area(:) => null() + real(kind_r8), pointer :: cell_area(:) ! Cell weights normalized to unity. CANNOT be safely deallocated because `cam_grid_attribute_*_t` ! just uses pointers to point at it. - real(kind_r8), pointer :: cell_weight(:) => null() + real(kind_r8), pointer :: cell_weight(:) ! Latitude coordinates. CANNOT be safely deallocated because `cam_grid_t` ! just uses pointers to point at it. - type(horiz_coord_t), pointer :: lat_coord => null() + type(horiz_coord_t), pointer :: lat_coord ! Longitude coordinates. CANNOT be safely deallocated because `cam_grid_t` ! just uses pointers to point at it. - type(horiz_coord_t), pointer :: lon_coord => null() + type(horiz_coord_t), pointer :: lon_coord + + nullify(indextocellid, indextoedgeid, indextovertexid) + nullify(areacell) + nullify(latcell, loncell) + nullify(latedge, lonedge) + nullify(latvertex, lonvertex) + + nullify(global_grid_index, global_grid_map) + nullify(cell_area, cell_weight) + nullify(lat_coord, lon_coord) ! Construct coordinate and grid objects for cell center grid (i.e., "mpas_cell"). ! Standard MPAS coordinate and dimension names are used. @@ -353,10 +372,8 @@ subroutine define_cam_grid() call cam_grid_attribute_register('mpas_cell', 'cell_weight', 'MPAS cell weight', 'nCells', cell_weight, & map=global_grid_index) - nullify(cell_area) - nullify(cell_weight) - nullify(lat_coord) - nullify(lon_coord) + nullify(cell_area, cell_weight) + nullify(lat_coord, lon_coord) ! Construct coordinate and grid objects for cell center grid (i.e., "cam_cell"). ! Standard CAM-SIMA coordinate and dimension names are used. @@ -373,14 +390,11 @@ subroutine define_cam_grid() nullify(areacell) nullify(indextocellid) - nullify(latcell) - nullify(loncell) + nullify(latcell, loncell) deallocate(global_grid_index) - nullify(global_grid_index) - nullify(global_grid_map) - nullify(lat_coord) - nullify(lon_coord) + nullify(global_grid_index, global_grid_map) + nullify(lat_coord, lon_coord) ! Construct coordinate and grid objects for edge node grid (i.e., "mpas_edge"). ! Standard MPAS coordinate and dimension names are used. @@ -414,14 +428,11 @@ subroutine define_cam_grid() unstruct=.true., block_indexed=.false.) nullify(indextoedgeid) - nullify(latedge) - nullify(lonedge) + nullify(latedge, lonedge) deallocate(global_grid_index) - nullify(global_grid_index) - nullify(global_grid_map) - nullify(lat_coord) - nullify(lon_coord) + nullify(global_grid_index, global_grid_map) + nullify(lat_coord, lon_coord) ! Construct coordinate and grid objects for vertex node grid (i.e., "mpas_vertex"). ! Standard MPAS coordinate and dimension names are used. @@ -455,14 +466,11 @@ subroutine define_cam_grid() unstruct=.true., block_indexed=.false.) nullify(indextovertexid) - nullify(latvertex) - nullify(lonvertex) + nullify(latvertex, lonvertex) deallocate(global_grid_index) - nullify(global_grid_index) - nullify(global_grid_map) - nullify(lat_coord) - nullify(lon_coord) + nullify(global_grid_index, global_grid_map) + nullify(lat_coord, lon_coord) end subroutine define_cam_grid !> Helper function for returning grid id given its name. From d18264ba2467694d053246789215e23c2048f176 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 26 Jul 2024 13:27:13 -0600 Subject: [PATCH 30/31] Also add `subname` to internal subroutines Easier to pinpoint where error occurs. --- src/dynamics/mpas/dyn_comp.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/dynamics/mpas/dyn_comp.F90 b/src/dynamics/mpas/dyn_comp.F90 index aec4a265..5c2d8f0b 100644 --- a/src/dynamics/mpas/dyn_comp.F90 +++ b/src/dynamics/mpas/dyn_comp.F90 @@ -379,6 +379,7 @@ 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_variable() + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::init_shared_variable' integer :: ierr integer :: k integer, pointer :: indextocellid(:) @@ -436,6 +437,7 @@ end subroutine init_shared_variable !> Set MPAS state `u` (i.e., horizontal velocity at edge interfaces). !> (KCW, 2024-05-13) subroutine set_mpas_state_u() + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_u' integer :: ierr integer :: k real(kind_r8), pointer :: ucellzonal(:, :), ucellmeridional(:, :) @@ -478,6 +480,7 @@ end subroutine set_mpas_state_u !> Set MPAS state `w` (i.e., vertical velocity at cell interfaces). !> (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(:, :) call dyn_debug_print('Setting MPAS state "w"') @@ -501,6 +504,7 @@ subroutine set_mpas_state_scalars() character(*), parameter :: constituent_qv_standard_name = & 'water_vapor_mixing_ratio_wrt_dry_air' + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_scalars' integer :: i, k integer :: ierr integer, allocatable :: constituent_index(:) @@ -562,6 +566,7 @@ 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() + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_rho_theta' integer :: i, k integer :: ierr integer, pointer :: index_qv @@ -671,6 +676,7 @@ 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() + character(*), parameter :: subname = 'dyn_comp::set_analytic_initial_condition::set_mpas_state_rho_base_theta_base' integer :: i, k integer :: ierr real(kind_r8), parameter :: t_base = 250.0_kind_r8 ! Base state temperature (K) of dry isothermal atmosphere. From 6428560b42e2eee6e8db89d1ffe04f1016212cd9 Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Fri, 26 Jul 2024 16:41:31 -0600 Subject: [PATCH 31/31] Adjust validity check for coupling time interval and time step Zero or negative values are not allowed. Others are just fine. --- src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 index 2e309671..b21cdcfe 100644 --- a/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 +++ b/src/dynamics/mpas/driver/dyn_mpas_subdriver.F90 @@ -2629,14 +2629,14 @@ subroutine dyn_mpas_init_phase4(self, coupling_time_interval) nullify(field_0d_real) nullify(mpas_pool) - if (real(coupling_time_interval, rkind) < 1.0_rkind) then + if (coupling_time_interval <= 0) then call self % model_error('Invalid coupling time interval ' // stringify([real(coupling_time_interval, rkind)]), & subname, __LINE__) end if call self % get_variable_pointer(config_dt, 'cfg', 'config_dt') - if (config_dt < 1.0_rkind) then + if (config_dt <= 0.0_rkind) then call self % model_error('Invalid time step ' // stringify([config_dt]), & subname, __LINE__) end if