diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 49cfedc8b8..ce0343f714 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -23,8 +23,8 @@ module MOM use MOM_diag_mediator, only : diag_grid_storage, diag_grid_storage_init use MOM_diag_mediator, only : diag_save_grids, diag_restore_grids use MOM_diag_mediator, only : diag_copy_storage_to_diag, diag_copy_diag_to_storage -use MOM_domains, only : MOM_domains_init, clone_MOM_domain -use MOM_domains, only : sum_across_PEs, pass_var, pass_vector +use MOM_domain_init, only : MOM_domains_init +use MOM_domains, only : sum_across_PEs, pass_var, pass_vector, clone_MOM_domain use MOM_domains, only : To_North, To_East, To_South, To_West use MOM_domains, only : To_All, Omit_corners, CGRID_NE, SCALAR_PAIR use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 2be0a978c1..a89514cde4 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -16,7 +16,6 @@ module MOM_dynamics_split_RK2 use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids -use MOM_domains, only : MOM_domains_init use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR use MOM_domains, only : To_North, To_East, Omit_Corners use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index a8de99df47..10b1f2e857 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -61,9 +61,8 @@ module MOM_dynamics_unsplit use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids -use MOM_domains, only : MOM_domains_init, pass_var, pass_vector -use MOM_domains, only : pass_var_start, pass_var_complete -use MOM_domains, only : pass_vector_start, pass_vector_complete +use MOM_domains, only : pass_var, pass_var_start, pass_var_complete +use MOM_domains, only : pass_vector, pass_vector_start, pass_vector_complete use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index c9da85fda9..8ca671d463 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -59,9 +59,8 @@ module MOM_dynamics_unsplit_RK2 use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr use MOM_diag_mediator, only : register_diag_field, register_static_field use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl -use MOM_domains, only : MOM_domains_init, pass_var, pass_vector -use MOM_domains, only : pass_var_start, pass_var_complete -use MOM_domains, only : pass_vector_start, pass_vector_complete +use MOM_domains, only : pass_var, pass_var_start, pass_var_complete +use MOM_domains, only : pass_vector, pass_vector_start, pass_vector_complete use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : MOM_set_verbosity diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 8844c65f40..9ca98adf71 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -5,7 +5,7 @@ module MOM_grid use MOM_hor_index, only : hor_index_type, hor_index_init use MOM_domains, only : MOM_domain_type, get_domain_extent, compute_block_extent -use MOM_domains, only : get_global_shape, get_domain_extent_dsamp2 +use MOM_domains, only : get_global_shape, get_domain_extent_dsamp2, deallocate_MOM_domain use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_unit_scaling, only : unit_scale_type @@ -630,8 +630,9 @@ subroutine MOM_grid_end(G) deallocate(G%gridLonT) ; deallocate(G%gridLatT) deallocate(G%gridLonB) ; deallocate(G%gridLatB) - deallocate(G%Domain%mpp_domain) - deallocate(G%Domain) + ! The cursory flag avoids doing any deallocation of memory in the underlying + ! infrastructure to avoid problems due to shared pointers. + call deallocate_MOM_domain(G%Domain, cursory=.true.) end subroutine MOM_grid_end diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 550f4d65e8..0746a120f2 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -4,7 +4,7 @@ module MOM_sum_output ! This file is part of MOM6. See LICENSE.md for the license. use iso_fortran_env, only : int64 -use MOM_coms, only : sum_across_PEs, PE_here, root_PE, num_PEs, max_across_PEs +use MOM_coms, only : sum_across_PEs, PE_here, root_PE, num_PEs, max_across_PEs, field_chksum use MOM_coms, only : reproducing_sum, reproducing_sum_EFP, EFP_to_real, real_to_EFP use MOM_coms, only : EFP_type, operator(+), operator(-), assignment(=), EFP_sum_across_PEs use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe, MOM_mesg @@ -25,7 +25,6 @@ module MOM_sum_output use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use mpp_mod, only : mpp_chksum use netcdf @@ -1511,13 +1510,13 @@ subroutine get_depth_list_checksums(G, depth_chksum, area_chksum) do j=G%jsc,G%jec ; do i=G%isc,G%iec field(i,j) = G%bathyT(i,j) enddo ; enddo - write(depth_chksum, '(Z16)') mpp_chksum(field(:,:)) + write(depth_chksum, '(Z16)') field_chksum(field(:,:)) ! Area checksum do j=G%jsc,G%jec ; do i=G%isc,G%iec field(i,j) = G%mask2dT(i,j) * G%US%L_to_m**2*G%areaT(i,j) enddo ; enddo - write(area_chksum, '(Z16)') mpp_chksum(field(:,:)) + write(area_chksum, '(Z16)') field_chksum(field(:,:)) deallocate(field) end subroutine get_depth_list_checksums diff --git a/src/framework/MOM_checksums.F90 b/src/framework/MOM_checksums.F90 index c3174dbe7b..5c503836f0 100644 --- a/src/framework/MOM_checksums.F90 +++ b/src/framework/MOM_checksums.F90 @@ -3,10 +3,11 @@ module MOM_checksums ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_array_transform, only: rotate_array, rotate_array_pair, rotate_vector +use MOM_array_transform, only : rotate_array, rotate_array_pair, rotate_vector +use MOM_array_transform, only : allocate_rotated_array use MOM_coms, only : PE_here, root_PE, num_PEs, sum_across_PEs use MOM_coms, only : min_across_PEs, max_across_PEs -use MOM_coms, only : reproducing_sum +use MOM_coms, only : reproducing_sum, field_chksum use MOM_error_handler, only : MOM_error, FATAL, is_root_pe use MOM_file_parser, only : log_version, param_file_type use MOM_hor_index, only : hor_index_type, rotate_hor_index @@ -15,7 +16,7 @@ module MOM_checksums implicit none ; private -public :: chksum0, zchksum +public :: chksum0, zchksum, rotated_field_chksum public :: hchksum, Bchksum, uchksum, vchksum, qchksum, is_NaN, chksum public :: hchksum_pair, uvchksum, Bchksum_pair public :: MOM_checksums_init @@ -75,6 +76,15 @@ module MOM_checksums module procedure is_NaN_0d, is_NaN_1d, is_NaN_2d, is_NaN_3d end interface +!> Rotate and compute the checksum of a field +interface rotated_field_chksum + module procedure rotated_field_chksum_real_0d + module procedure rotated_field_chksum_real_1d + module procedure rotated_field_chksum_real_2d + module procedure rotated_field_chksum_real_3d + module procedure rotated_field_chksum_real_4d +end interface rotated_field_chksum + integer, parameter :: bc_modulus = 1000000000 !< Modulus of checksum bitcount integer, parameter :: default_shift=0 !< The default array shift logical :: calculateStatistics=.true. !< If true, report min, max and mean. @@ -2021,16 +2031,16 @@ function is_NaN_1d(x, skip_mpp) logical :: is_NaN_1d integer :: i, n - logical :: call_mpp + logical :: global_check n = 0 do i = LBOUND(x,1), UBOUND(x,1) if (is_NaN_0d(x(i))) n = n + 1 enddo - call_mpp = .true. - if (present(skip_mpp)) call_mpp = .not.skip_mpp + global_check = .true. + if (present(skip_mpp)) global_check = .not.skip_mpp - if (call_mpp) call sum_across_PEs(n) + if (global_check) call sum_across_PEs(n) is_NaN_1d = .false. if (n>0) is_NaN_1d = .true. @@ -2072,6 +2082,121 @@ function is_NaN_3d(x) end function is_NaN_3d +! The following set of routines do a checksum across the computational domain of +! a field, with the potential for rotation of this field and masking. + +!> Compute the field checksum of a scalar. +function rotated_field_chksum_real_0d(field, pelist, mask_val, turns) & + result(chksum) + real, intent(in) :: field !< Input scalar + integer, optional, intent(in) :: pelist(:) !< PE list of ranks to checksum + real, optional, intent(in) :: mask_val !< FMS mask value + integer, optional, intent(in) :: turns !< Number of quarter turns + integer :: chksum !< checksum of scalar + + if (present(turns)) call MOM_error(FATAL, "Rotation not supported for 0d fields.") + + chksum = field_chksum(field, pelist=pelist, mask_val=mask_val) +end function rotated_field_chksum_real_0d + + +!> Compute the field checksum of a 1d field. +function rotated_field_chksum_real_1d(field, pelist, mask_val, turns) & + result(chksum) + real, dimension(:), intent(in) :: field !< Input array + integer, optional, intent(in) :: pelist(:) !< PE list of ranks to checksum + real, optional, intent(in) :: mask_val !< FMS mask value + integer, optional, intent(in) :: turns !< Number of quarter turns + integer :: chksum !< checksum of array + + if (present(turns)) call MOM_error(FATAL, "Rotation not supported for 1d fields.") + + chksum = field_chksum(field, pelist=pelist, mask_val=mask_val) +end function rotated_field_chksum_real_1d + + +!> Compute the field checksum of a rotated 2d field. +function rotated_field_chksum_real_2d(field, pelist, mask_val, turns) & + result(chksum) + real, dimension(:,:), intent(in) :: field !< Unrotated input field + integer, optional, intent(in) :: pelist(:) !< PE list of ranks to checksum + real, optional, intent(in) :: mask_val !< FMS mask value + integer, optional, intent(in) :: turns !< Number of quarter turns + integer :: chksum !< checksum of array + + ! Local variables + real, allocatable :: field_rot(:,:) ! A rotated version of field, with the same units + integer :: qturns ! The number of quarter turns through which to rotate field + + qturns = 0 + if (present(turns)) & + qturns = modulo(turns, 4) + + if (qturns == 0) then + chksum = field_chksum(field, pelist=pelist, mask_val=mask_val) + else + call allocate_rotated_array(field, [1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + chksum = field_chksum(field_rot, pelist=pelist, mask_val=mask_val) + deallocate(field_rot) + endif +end function rotated_field_chksum_real_2d + +!> Compute the field checksum of a rotated 3d field. +function rotated_field_chksum_real_3d(field, pelist, mask_val, turns) & + result(chksum) + real, dimension(:,:,:), intent(in) :: field !< Unrotated input field + integer, optional, intent(in) :: pelist(:) !< PE list of ranks to checksum + real, optional, intent(in) :: mask_val !< FMS mask value + integer, optional, intent(in) :: turns !< Number of quarter turns + integer :: chksum !< checksum of array + + ! Local variables + real, allocatable :: field_rot(:,:,:) ! A rotated version of field, with the same units + integer :: qturns ! The number of quarter turns through which to rotate field + + qturns = 0 + if (present(turns)) & + qturns = modulo(turns, 4) + + if (qturns == 0) then + chksum = field_chksum(field, pelist=pelist, mask_val=mask_val) + else + call allocate_rotated_array(field, [1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + chksum = field_chksum(field_rot, pelist=pelist, mask_val=mask_val) + deallocate(field_rot) + endif +end function rotated_field_chksum_real_3d + +!> Compute the field checksum of a rotated 4d field. +function rotated_field_chksum_real_4d(field, pelist, mask_val, turns) & + result(chksum) + real, dimension(:,:,:,:), intent(in) :: field !< Unrotated input field + integer, optional, intent(in) :: pelist(:) !< PE list of ranks to checksum + real, optional, intent(in) :: mask_val !< FMS mask value + integer, optional, intent(in) :: turns !< Number of quarter turns + integer :: chksum !< checksum of array + + ! Local variables + real, allocatable :: field_rot(:,:,:,:) ! A rotated version of field, with the same units + integer :: qturns ! The number of quarter turns through which to rotate field + + qturns = 0 + if (present(turns)) & + qturns = modulo(turns, 4) + + if (qturns == 0) then + chksum = field_chksum(field, pelist=pelist, mask_val=mask_val) + else + call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + chksum = field_chksum(field_rot, pelist=pelist, mask_val=mask_val) + deallocate(field_rot) + endif +end function rotated_field_chksum_real_4d + + !> Write a message including the checksum of the non-shifted array subroutine chk_sum_msg1(fmsg, bc0, mesg, iounit) character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble diff --git a/src/framework/MOM_coms.F90 b/src/framework/MOM_coms.F90 index 0c6b948980..04ed46ad22 100644 --- a/src/framework/MOM_coms.F90 +++ b/src/framework/MOM_coms.F90 @@ -9,13 +9,13 @@ module MOM_coms use memutils_mod, only : print_memuse_stats use mpp_mod, only : PE_here => mpp_pe, root_PE => mpp_root_pe, num_PEs => mpp_npes use mpp_mod, only : Set_PElist => mpp_set_current_pelist, Get_PElist => mpp_get_current_pelist -use mpp_mod, only : broadcast => mpp_broadcast +use mpp_mod, only : broadcast => mpp_broadcast, field_chksum => mpp_chksum use mpp_mod, only : sum_across_PEs => mpp_sum, max_across_PEs => mpp_max, min_across_PEs => mpp_min implicit none ; private public :: PE_here, root_PE, num_PEs, MOM_infra_init, MOM_infra_end -public :: broadcast, sum_across_PEs, min_across_PEs, max_across_PEs +public :: broadcast, sum_across_PEs, min_across_PEs, max_across_PEs, field_chksum public :: reproducing_sum, reproducing_sum_EFP, EFP_sum_across_PEs, EFP_list_sum_across_PEs public :: EFP_plus, EFP_minus, EFP_to_real, real_to_EFP, EFP_real_diff public :: operator(+), operator(-), assignment(=) diff --git a/src/framework/MOM_domain_init.F90 b/src/framework/MOM_domain_init.F90 new file mode 100644 index 0000000000..25064cf24e --- /dev/null +++ b/src/framework/MOM_domain_init.F90 @@ -0,0 +1,330 @@ +!> Describes the decomposed MOM domain and has routines for communications across PEs +module MOM_domain_init + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_coms, only : num_PEs +use MOM_domains, only : MOM_domain_type, create_MOM_domain, MOM_define_layout +use MOM_domains, only : MOM_thread_affinity_set, set_MOM_thread_affinity +use MOM_error_handler, only : MOM_error, MOM_mesg, NOTE, WARNING, FATAL +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_io, only : file_exists +use MOM_string_functions, only : slasher + +implicit none ; private + +public :: MOM_domains_init, MOM_domain_type + +contains + +!> MOM_domains_init initializes a MOM_domain_type variable, based on the information +!! read in from a param_file_type, and optionally returns data describing various' +!! properties of the domain type. +subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & + NIHALO, NJHALO, NIGLOBAL, NJGLOBAL, NIPROC, NJPROC, & + min_halo, domain_name, include_name, param_suffix) + type(MOM_domain_type), pointer :: MOM_dom !< A pointer to the MOM_domain_type + !! being defined here. + type(param_file_type), intent(in) :: param_file !< A structure to parse for + !! run-time parameters + logical, optional, intent(in) :: symmetric !< If present, this specifies + !! whether this domain is symmetric, regardless of + !! whether the macro SYMMETRIC_MEMORY_ is defined. + logical, optional, intent(in) :: static_memory !< If present and true, this + !! domain type is set up for static memory and + !! error checking of various input values is + !! performed against those in the input file. + integer, optional, intent(in) :: NIHALO !< Default halo sizes, required + !! with static memory. + integer, optional, intent(in) :: NJHALO !< Default halo sizes, required + !! with static memory. + integer, optional, intent(in) :: NIGLOBAL !< Total domain sizes, required + !! with static memory. + integer, optional, intent(in) :: NJGLOBAL !< Total domain sizes, required + !! with static memory. + integer, optional, intent(in) :: NIPROC !< Processor counts, required with + !! static memory. + integer, optional, intent(in) :: NJPROC !< Processor counts, required with + !! static memory. + integer, dimension(2), optional, intent(inout) :: min_halo !< If present, this sets the + !! minimum halo size for this domain in the i- and j- + !! directions, and returns the actual halo size used. + character(len=*), optional, intent(in) :: domain_name !< A name for this domain, "MOM" + !! if missing. + character(len=*), optional, intent(in) :: include_name !< A name for model's include file, + !! "MOM_memory.h" if missing. + character(len=*), optional, intent(in) :: param_suffix !< A suffix to apply to + !! layout-specific parameters. + + ! Local variables + integer, dimension(2) :: layout ! The number of logical processors in the i- and j- directions + integer, dimension(2) :: io_layout ! The layout of logical processors for input and output + !$ integer :: ocean_nthreads ! Number of openMP threads + !$ logical :: ocean_omp_hyper_thread ! If true use openMP hyper-threads + integer, dimension(2) :: n_global ! The number of i- and j- points in the global computational domain + integer, dimension(2) :: n_halo ! The number of i- and j- points in the halos + integer :: nihalo_dflt, njhalo_dflt ! The default halo sizes + integer :: PEs_used ! The number of processors used + logical, dimension(2) :: reentrant ! True if the x- and y- directions are periodic. + logical, dimension(2,2) :: tripolar ! A set of flag indicating whether there is tripolar + ! connectivity for any of the four logical edges of the grid. + ! Currently only tripolar_N is implemented. + logical :: is_static ! If true, static memory is being used for this domain. + logical :: is_symmetric ! True if the domain being set up will use symmetric memory. + logical :: nonblocking ! If true, nonblocking halo updates will be used. + logical :: thin_halos ! If true, If true, optional arguments may be used to specify the + ! width of the halos that are updated with each call. + logical :: mask_table_exists ! True if there is a mask table file + character(len=128) :: inputdir ! The directory in which to find the diag table + character(len=200) :: mask_table ! The file name and later the full path to the diag table + character(len=64) :: inc_nm ! The name of the memory include file + character(len=200) :: mesg ! A string to use for error messages + + integer :: nip_parsed, njp_parsed + character(len=8) :: char_xsiz, char_ysiz, char_niglobal, char_njglobal + character(len=40) :: nihalo_nm, njhalo_nm, layout_nm, io_layout_nm, masktable_nm + character(len=40) :: niproc_nm, njproc_nm + ! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=40) :: mdl ! This module's name. + + PEs_used = num_PEs() + + mdl = "MOM_domains" !### Change this to "MOM_domain_init" + + is_symmetric = .true. ; if (present(symmetric)) is_symmetric = symmetric + if (present(min_halo)) mdl = trim(mdl)//" min_halo" + + inc_nm = "MOM_memory.h" ; if (present(include_name)) inc_nm = trim(include_name) + + nihalo_nm = "NIHALO" ; njhalo_nm = "NJHALO" + layout_nm = "LAYOUT" ; io_layout_nm = "IO_LAYOUT" ; masktable_nm = "MASKTABLE" + niproc_nm = "NIPROC" ; njproc_nm = "NJPROC" + if (present(param_suffix)) then ; if (len(trim(adjustl(param_suffix))) > 0) then + nihalo_nm = "NIHALO"//(trim(adjustl(param_suffix))) + njhalo_nm = "NJHALO"//(trim(adjustl(param_suffix))) + layout_nm = "LAYOUT"//(trim(adjustl(param_suffix))) + io_layout_nm = "IO_LAYOUT"//(trim(adjustl(param_suffix))) + masktable_nm = "MASKTABLE"//(trim(adjustl(param_suffix))) + niproc_nm = "NIPROC"//(trim(adjustl(param_suffix))) + njproc_nm = "NJPROC"//(trim(adjustl(param_suffix))) + endif ; endif + + is_static = .false. ; if (present(static_memory)) is_static = static_memory + if (is_static) then + if (.not.present(NIHALO)) call MOM_error(FATAL, "NIHALO must be "// & + "present in the call to MOM_domains_init with static memory.") + if (.not.present(NJHALO)) call MOM_error(FATAL, "NJHALO must be "// & + "present in the call to MOM_domains_init with static memory.") + if (.not.present(NIGLOBAL)) call MOM_error(FATAL, "NIGLOBAL must be "// & + "present in the call to MOM_domains_init with static memory.") + if (.not.present(NJGLOBAL)) call MOM_error(FATAL, "NJGLOBAL must be "// & + "present in the call to MOM_domains_init with static memory.") + if (.not.present(NIPROC)) call MOM_error(FATAL, "NIPROC must be "// & + "present in the call to MOM_domains_init with static memory.") + if (.not.present(NJPROC)) call MOM_error(FATAL, "NJPROC must be "// & + "present in the call to MOM_domains_init with static memory.") + endif + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "", log_to_all=.true., layout=.true.) + call get_param(param_file, mdl, "REENTRANT_X", reentrant(1), & + "If true, the domain is zonally reentrant.", default=.true.) + call get_param(param_file, mdl, "REENTRANT_Y", reentrant(2), & + "If true, the domain is meridionally reentrant.", & + default=.false.) + tripolar(1:2,1:2) = .false. + call get_param(param_file, mdl, "TRIPOLAR_N", tripolar(2,2), & + "Use tripolar connectivity at the northern edge of the "//& + "domain. With TRIPOLAR_N, NIGLOBAL must be even.", & + default=.false.) + +# ifndef NOT_SET_AFFINITY + !$ if (.not.MOM_thread_affinity_set()) then + !$ call get_param(param_file, mdl, "OCEAN_OMP_THREADS", ocean_nthreads, & + !$ "The number of OpenMP threads that MOM6 will use.", & + !$ default = 1, layoutParam=.true.) + !$ call get_param(param_file, mdl, "OCEAN_OMP_HYPER_THREAD", ocean_omp_hyper_thread, & + !$ "If True, use hyper-threading.", default = .false., layoutParam=.true.) + !$ call set_MOM_thread_affinity(ocean_nthreads, ocean_omp_hyper_thread) + !$ endif +# endif + + call log_param(param_file, mdl, "!SYMMETRIC_MEMORY_", is_symmetric, & + "If defined, the velocity point data domain includes every face of the "//& + "thickness points. In other words, some arrays are larger than others, "//& + "depending on where they are on the staggered grid. Also, the starting "//& + "index of the velocity-point arrays is usually 0, not 1. "//& + "This can only be set at compile time.",& + layoutParam=.true.) + call get_param(param_file, mdl, "NONBLOCKING_UPDATES", nonblocking, & + "If true, non-blocking halo updates may be used.", & + default=.false., layoutParam=.true.) + !### Note the duplicated "the the" in the following description, which should be fixed as a part + ! of a larger commit that also changes other MOM_parameter_doc file messages, but for now + ! reproduces the existing output files. + call get_param(param_file, mdl, "THIN_HALO_UPDATES", thin_halos, & + "If true, optional arguments may be used to specify the the width of the "//& + "halos that are updated with each call.", & + default=.true., layoutParam=.true.) + + nihalo_dflt = 4 ; njhalo_dflt = 4 + if (present(NIHALO)) nihalo_dflt = NIHALO + if (present(NJHALO)) njhalo_dflt = NJHALO + + call log_param(param_file, mdl, "!STATIC_MEMORY_", is_static, & + "If STATIC_MEMORY_ is defined, the principle variables will have sizes that "//& + "are statically determined at compile time. Otherwise the sizes are not "//& + "determined until run time. The STATIC option is substantially faster, but "//& + "does not allow the PE count to be changed at run time. This can only be "//& + "set at compile time.", layoutParam=.true.) + + if (is_static) then + call get_param(param_file, mdl, "NIGLOBAL", n_global(1), & + "The total number of thickness grid points in the x-direction in the physical "//& + "domain. With STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.", & + static_value=NIGLOBAL) + call get_param(param_file, mdl, "NJGLOBAL", n_global(2), & + "The total number of thickness grid points in the y-direction in the physical "//& + "domain. With STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.", & + static_value=NJGLOBAL) + if (n_global(1) /= NIGLOBAL) call MOM_error(FATAL,"MOM_domains_init: " // & + "static mismatch for NIGLOBAL_ domain size. Header file does not match input namelist") + if (n_global(2) /= NJGLOBAL) call MOM_error(FATAL,"MOM_domains_init: " // & + "static mismatch for NJGLOBAL_ domain size. Header file does not match input namelist") + + ! Check the requirement of equal sized compute domains when STATIC_MEMORY_ is used. + if ((MOD(NIGLOBAL, NIPROC) /= 0) .OR. (MOD(NJGLOBAL, NJPROC) /= 0)) then + write( char_xsiz, '(i4)' ) NIPROC + write( char_ysiz, '(i4)' ) NJPROC + write( char_niglobal, '(i4)' ) NIGLOBAL + write( char_njglobal, '(i4)' ) NJGLOBAL + call MOM_error(WARNING, 'MOM_domains: Processor decomposition (NIPROC_,NJPROC_) = ('//& + trim(char_xsiz)//','//trim(char_ysiz)//') does not evenly divide size '//& + 'set by preprocessor macro ('//trim(char_niglobal)//','//trim(char_njglobal)//').') + call MOM_error(FATAL,'MOM_domains: #undef STATIC_MEMORY_ in '//trim(inc_nm)//' to use '//& + 'dynamic allocation, or change processor decomposition to evenly divide the domain.') + endif + else + call get_param(param_file, mdl, "NIGLOBAL", n_global(1), & + "The total number of thickness grid points in the x-direction in the physical "//& + "domain. With STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.", & + fail_if_missing=.true.) + call get_param(param_file, mdl, "NJGLOBAL", n_global(2), & + "The total number of thickness grid points in the y-direction in the physical "//& + "domain. With STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.", & + fail_if_missing=.true.) + endif + + call get_param(param_file, mdl, trim(nihalo_nm), n_halo(1), & + "The number of halo points on each side in the x-direction. How this is set "//& + "varies with the calling component and static or dynamic memory configuration.", & + default=nihalo_dflt, static_value=nihalo_dflt) + call get_param(param_file, mdl, trim(njhalo_nm), n_halo(2), & + "The number of halo points on each side in the y-direction. How this is set "//& + "varies with the calling component and static or dynamic memory configuration.", & + default=njhalo_dflt, static_value=njhalo_dflt) + if (present(min_halo)) then + n_halo(1) = max(n_halo(1), min_halo(1)) + min_halo(1) = n_halo(1) + n_halo(2) = max(n_halo(2), min_halo(2)) + min_halo(2) = n_halo(2) + ! These are generally used only with static memory, so they are considerd layout params. + call log_param(param_file, mdl, "!NIHALO min_halo", n_halo(1), layoutParam=.true.) + call log_param(param_file, mdl, "!NJHALO min_halo", n_halo(2), layoutParam=.true.) + endif + if (is_static .and. .not.present(min_halo)) then + if (n_halo(1) /= NIHALO) call MOM_error(FATAL,"MOM_domains_init: " // & + "static mismatch for "//trim(nihalo_nm)//" domain size") + if (n_halo(2) /= NJHALO) call MOM_error(FATAL,"MOM_domains_init: " // & + "static mismatch for "//trim(njhalo_nm)//" domain size") + endif + + call get_param(param_file, mdl, "INPUTDIR", inputdir, do_not_log=.true., default=".") + inputdir = slasher(inputdir) + + call get_param(param_file, mdl, trim(masktable_nm), mask_table, & + "A text file to specify n_mask, layout and mask_list. This feature masks out "//& + "processors that contain only land points. The first line of mask_table is the "//& + "number of regions to be masked out. The second line is the layout of the "//& + "model and must be consistent with the actual model layout. The following "//& + "(n_mask) lines give the logical positions of the processors that are masked "//& + "out. The mask_table can be created by tools like check_mask. The following "//& + "example of mask_table masks out 2 processors, (1,2) and (3,6), out of the 24 "//& + "in a 4x6 layout: \n 2\n 4,6\n 1,2\n 3,6\n", default="MOM_mask_table", & + layoutParam=.true.) + mask_table = trim(inputdir)//trim(mask_table) + mask_table_exists = file_exists(mask_table) + + if (is_static) then + layout(1) = NIPROC ; layout(2) = NJPROC + else + call get_param(param_file, mdl, trim(layout_nm), layout, & + "The processor layout to be used, or 0, 0 to automatically set the layout "//& + "based on the number of processors.", default=0, do_not_log=.true.) + call get_param(param_file, mdl, trim(niproc_nm), nip_parsed, & + "The number of processors in the x-direction.", default=-1, do_not_log=.true.) + call get_param(param_file, mdl, trim(njproc_nm), njp_parsed, & + "The number of processors in the y-direction.", default=-1, do_not_log=.true.) + if (nip_parsed > -1) then + if ((layout(1) > 0) .and. (layout(1) /= nip_parsed)) & + call MOM_error(FATAL, trim(layout_nm)//" and "//trim(niproc_nm)//" set inconsistently. "//& + "Only LAYOUT should be used.") + layout(1) = nip_parsed + call MOM_mesg(trim(niproc_nm)//" used to set "//trim(layout_nm)//" in dynamic mode. "//& + "Shift to using "//trim(layout_nm)//" instead.") + endif + if (njp_parsed > -1) then + if ((layout(2) > 0) .and. (layout(2) /= njp_parsed)) & + call MOM_error(FATAL, trim(layout_nm)//" and "//trim(njproc_nm)//" set inconsistently. "//& + "Only "//trim(layout_nm)//" should be used.") + layout(2) = njp_parsed + call MOM_mesg(trim(njproc_nm)//" used to set "//trim(layout_nm)//" in dynamic mode. "//& + "Shift to using "//trim(layout_nm)//" instead.") + endif + + if ( (layout(1) == 0) .and. (layout(2) == 0) ) & + call MOM_define_layout( (/ 1, n_global(1), 1, n_global(2) /), PEs_used, layout) + if ( (layout(1) /= 0) .and. (layout(2) == 0) ) layout(2) = PEs_used / layout(1) + if ( (layout(1) == 0) .and. (layout(2) /= 0) ) layout(1) = PEs_used / layout(2) + + if (layout(1)*layout(2) /= PEs_used .and. (.not. mask_table_exists) ) then + write(mesg,'("MOM_domains_init: The product of the two components of layout, ", & + & 2i4,", is not the number of PEs used, ",i5,".")') & + layout(1), layout(2), PEs_used + call MOM_error(FATAL, mesg) + endif + endif + call log_param(param_file, mdl, trim(niproc_nm), layout(1), & + "The number of processors in the x-direction. With STATIC_MEMORY_ this "//& + "is set in "//trim(inc_nm)//" at compile time.", layoutParam=.true.) + call log_param(param_file, mdl, trim(njproc_nm), layout(2), & + "The number of processors in the y-direction. With STATIC_MEMORY_ this "//& + "is set in "//trim(inc_nm)//" at compile time.", layoutParam=.true.) + call log_param(param_file, mdl, trim(layout_nm), layout, & + "The processor layout that was actually used.", layoutParam=.true.) + + ! Idiot check that fewer PEs than columns have been requested + if (layout(1)*layout(2) > n_global(1)*n_global(2)) then + write(mesg,'(a,2(i5,x,a))') 'You requested to use',layout(1)*layout(2), & + 'PEs but there are only', n_global(1)*n_global(2), 'columns in the model' + call MOM_error(FATAL, mesg) + endif + + if (mask_table_exists) & + call MOM_error(NOTE, 'MOM_domains_init: reading maskmap information from '//trim(mask_table)) + + ! Set up the I/O layout, it will be checked later that it uses an even multiple of the number of + ! PEs in each direction. + io_layout(:) = (/ 1, 1 /) + call get_param(param_file, mdl, trim(io_layout_nm), io_layout, & + "The processor layout to be used, or 0,0 to automatically set the io_layout "//& + "to be the same as the layout.", default=1, layoutParam=.true.) + + call create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar, layout, & + io_layout=io_layout, domain_name=domain_name, mask_table=mask_table, & + symmetric=symmetric, thin_halos=thin_halos, nonblocking=nonblocking) + +end subroutine MOM_domains_init + +end module MOM_domain_init diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 46cc9c526a..56ac0b3ccf 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -12,11 +12,12 @@ module MOM_domains use MOM_file_parser, only : param_file_type use MOM_string_functions, only : slasher -use mpp_domains_mod, only : mpp_define_layout, mpp_get_boundary +use mpp_domains_mod, only : MOM_define_layout => mpp_define_layout, mpp_get_boundary use mpp_domains_mod, only : MOM_define_io_domain => mpp_define_io_domain use mpp_domains_mod, only : MOM_define_domain => mpp_define_domains -use mpp_domains_mod, only : domain2D, domain1D, mpp_get_data_domain +use mpp_domains_mod, only : domain2D, domain1D, mpp_get_data_domain, mpp_get_domain_components use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_global_domain +use mpp_domains_mod, only : mpp_get_domain_extents, mpp_deallocate_domain use mpp_domains_mod, only : global_field_sum => mpp_global_sum use mpp_domains_mod, only : mpp_update_domains, CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE use mpp_domains_mod, only : mpp_start_update_domains, mpp_complete_update_domains @@ -36,7 +37,9 @@ module MOM_domains implicit none ; private public :: MOM_domains_init, MOM_infra_init, MOM_infra_end, get_domain_extent, get_domain_extent_dsamp2 -public :: MOM_define_domain, MOM_define_io_domain, clone_MOM_domain +public :: create_MOM_domain, clone_MOM_domain, get_domain_components +public :: deallocate_MOM_domain, deallocate_domain_contents +public :: MOM_define_domain, MOM_define_layout, MOM_define_io_domain public :: pass_var, pass_vector, PE_here, root_PE, num_PEs public :: pass_var_start, pass_var_complete, fill_symmetric_edges, broadcast public :: pass_vector_start, pass_vector_complete @@ -46,9 +49,10 @@ module MOM_domains public :: To_East, To_West, To_North, To_South, To_All, Omit_Corners public :: create_group_pass, do_group_pass, group_pass_type public :: start_group_pass, complete_group_pass -public :: compute_block_extent, get_global_shape +public :: compute_block_extent, get_global_shape, get_layout_extents +public :: MOM_thread_affinity_set, set_MOM_thread_affinity public :: get_simple_array_i_ind, get_simple_array_j_ind -public :: domain2D +public :: domain2D, domain1D !> Do a halo update on an array interface pass_var @@ -100,7 +104,12 @@ module MOM_domains module procedure clone_MD_to_MD, clone_MD_to_d2D end interface clone_MOM_domain -!> The MOM_domain_type contains information about the domain decompositoin. +!> Extract the 1-d domain components from a MOM_domain or domain2d +interface get_domain_components + module procedure get_domain_components_MD, get_domain_components_d2D +end interface get_domain_components + +!> The MOM_domain_type contains information about the domain decomposition. type, public :: MOM_domain_type type(domain2D), pointer :: mpp_domain => NULL() !< The FMS domain with halos !! on this processor, centered at h points. @@ -1169,7 +1178,7 @@ subroutine complete_group_pass(group, MOM_dom, clock) end subroutine complete_group_pass -!> MOM_domains_init initalizes a MOM_domain_type variable, based on the information +!> MOM_domains_init initializes a MOM_domain_type variable, based on the information !! read in from a param_file_type, and optionally returns data describing various' !! properties of the domain type. subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & @@ -1183,8 +1192,9 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & !! whether this domain is symmetric, regardless of !! whether the macro SYMMETRIC_MEMORY_ is defined. logical, optional, intent(in) :: static_memory !< If present and true, this - !! domain type is set up for static memory and error checking of - !! various input values is performed against those in the input file. + !! domain type is set up for static memory and error + !! checking of various input values is performed against + !! those in the input file. integer, optional, intent(in) :: NIHALO !< Default halo sizes, required !! with static memory. integer, optional, intent(in) :: NJHALO !< Default halo sizes, required @@ -1198,8 +1208,8 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & integer, optional, intent(in) :: NJPROC !< Processor counts, required with !! static memory. integer, dimension(2), optional, intent(inout) :: min_halo !< If present, this sets the - !! minimum halo size for this domain in the i- and j- - !! directions, and returns the actual halo size used. + !! minimum halo size for this domain in the i- and j- + !! directions, and returns the actual halo size used. character(len=*), optional, intent(in) :: domain_name !< A name for this domain, "MOM" !! if missing. character(len=*), optional, intent(in) :: include_name !< A name for model's include file, @@ -1211,46 +1221,43 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & integer, dimension(2) :: layout = (/ 1, 1 /) integer, dimension(2) :: io_layout = (/ 0, 0 /) integer, dimension(4) :: global_indices -!$ integer :: ocean_nthreads ! Number of Openmp threads -!$ integer :: get_cpu_affinity, omp_get_thread_num, omp_get_num_threads -!$ logical :: ocean_omp_hyper_thread + !$ integer :: ocean_nthreads ! Number of Openmp threads + !$ logical :: ocean_omp_hyper_thread + integer, dimension(2) :: n_global ! The number of i- and j- points in the global computational domain. + integer, dimension(2) :: n_halo ! The number of i- and j- points in the halos. integer :: nihalo_dflt, njhalo_dflt integer :: pe, proc_used - integer :: X_FLAGS, Y_FLAGS - logical :: reentrant_x, reentrant_y, tripolar_N, is_static + logical, dimension(2) :: reentrant ! True if the x- and y- directions are periodic. + logical, dimension(2,2) :: tripolar ! A set of flag indicating whether there is tripolar + ! connectivity for any of the four logical edges of the grid. + ! Currently only tripolar_N is implemented. + logical :: is_static ! If true, static memory is being used for this domain. + logical :: is_symmetric ! True if the domainn being set up will use symmetric memory. + logical :: nonblocking ! If true, nonblocking halo updates will be used. + logical :: thin_halos ! If true, If true, optional arguments may be used to specify the + ! width of the halos that are updated with each call. logical :: mask_table_exists character(len=128) :: mask_table, inputdir - character(len=64) :: dom_name, inc_nm + character(len=64) :: inc_nm character(len=200) :: mesg - integer :: xsiz, ysiz, nip_parsed, njp_parsed - integer :: isc,iec,jsc,jec ! The bounding indices of the computational domain. + integer :: nip_parsed, njp_parsed character(len=8) :: char_xsiz, char_ysiz, char_niglobal, char_njglobal character(len=40) :: nihalo_nm, njhalo_nm, layout_nm, io_layout_nm, masktable_nm character(len=40) :: niproc_nm, njproc_nm - integer :: xhalo_d2,yhalo_d2 -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl ! This module's name. - if (.not.associated(MOM_dom)) then - allocate(MOM_dom) - allocate(MOM_dom%mpp_domain) - allocate(MOM_dom%mpp_domain_d2) - endif - pe = PE_here() proc_used = num_PEs() mdl = "MOM_domains" - MOM_dom%symmetric = .true. - if (present(symmetric)) then ; MOM_dom%symmetric = symmetric ; endif + is_symmetric = .true. ; if (present(symmetric)) is_symmetric = symmetric if (present(min_halo)) mdl = trim(mdl)//" min_halo" - dom_name = "MOM" ; inc_nm = "MOM_memory.h" - if (present(domain_name)) dom_name = trim(domain_name) - if (present(include_name)) inc_nm = trim(include_name) + inc_nm = "MOM_memory.h" ; if (present(include_name)) inc_nm = trim(include_name) nihalo_nm = "NIHALO" ; njhalo_nm = "NJHALO" layout_nm = "LAYOUT" ; io_layout_nm = "IO_LAYOUT" ; masktable_nm = "MASKTABLE" @@ -1283,36 +1290,29 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version, "", log_to_all=.true., layout=.true.) - call get_param(param_file, mdl, "REENTRANT_X", reentrant_x, & + call get_param(param_file, mdl, "REENTRANT_X", reentrant(1), & "If true, the domain is zonally reentrant.", default=.true.) - call get_param(param_file, mdl, "REENTRANT_Y", reentrant_y, & + call get_param(param_file, mdl, "REENTRANT_Y", reentrant(2), & "If true, the domain is meridionally reentrant.", & default=.false.) - call get_param(param_file, mdl, "TRIPOLAR_N", tripolar_N, & + tripolar(1:2,1:2) = .false. + call get_param(param_file, mdl, "TRIPOLAR_N", tripolar(2,2), & "Use tripolar connectivity at the northern edge of the "//& "domain. With TRIPOLAR_N, NIGLOBAL must be even.", & default=.false.) #ifndef NOT_SET_AFFINITY -!$ call fms_affinity_init -!$OMP PARALLEL -!$OMP master -!$ ocean_nthreads = omp_get_num_threads() -!$OMP END MASTER -!$OMP END PARALLEL -!$ if(ocean_nthreads < 2 ) then -!$ call get_param(param_file, mdl, "OCEAN_OMP_THREADS", ocean_nthreads, & -!$ "The number of OpenMP threads that MOM6 will use.", & -!$ default = 1, layoutParam=.true.) -!$ call get_param(param_file, mdl, "OCEAN_OMP_HYPER_THREAD", ocean_omp_hyper_thread, & -!$ "If True, use hyper-threading.", default = .false., layoutParam=.true.) -!$ call fms_affinity_set('OCEAN', ocean_omp_hyper_thread, ocean_nthreads) -!$ call omp_set_num_threads(ocean_nthreads) -!$ write(6,*) "MOM_domains_mod OMPthreading ", fms_affinity_get(), omp_get_thread_num(), omp_get_num_threads() -!$ flush(6) -!$ endif + !$ if (.not.MOM_thread_affinity_set()) then + !$ call get_param(param_file, mdl, "OCEAN_OMP_THREADS", ocean_nthreads, & + !$ "The number of OpenMP threads that MOM6 will use.", & + !$ default = 1, layoutParam=.true.) + !$ call get_param(param_file, mdl, "OCEAN_OMP_HYPER_THREAD", ocean_omp_hyper_thread, & + !$ "If True, use hyper-threading.", default = .false., layoutParam=.true.) + !$ call set_MOM_thread_affinity(ocean_nthreads, ocean_omp_hyper_thread) + !$ endif #endif - call log_param(param_file, mdl, "!SYMMETRIC_MEMORY_", MOM_dom%symmetric, & + + call log_param(param_file, mdl, "!SYMMETRIC_MEMORY_", is_symmetric, & "If defined, the velocity point data domain includes "//& "every face of the thickness points. In other words, "//& "some arrays are larger than others, depending on where "//& @@ -1320,10 +1320,10 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & "index of the velocity-point arrays is usually 0, not 1. "//& "This can only be set at compile time.",& layoutParam=.true.) - call get_param(param_file, mdl, "NONBLOCKING_UPDATES", MOM_dom%nonblocking_updates, & + call get_param(param_file, mdl, "NONBLOCKING_UPDATES", nonblocking, & "If true, non-blocking halo updates may be used.", & default=.false., layoutParam=.true.) - call get_param(param_file, mdl, "THIN_HALO_UPDATES", MOM_dom%thin_halo_updates, & + call get_param(param_file, mdl, "THIN_HALO_UPDATES", thin_halos, & "If true, optional arguments may be used to specify the "//& "the width of the halos that are updated with each call.", & default=.true., layoutParam=.true.) @@ -1342,60 +1342,72 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & layoutParam=.true.) if (is_static) then - call get_param(param_file, mdl, "NIGLOBAL", MOM_dom%niglobal, & + call get_param(param_file, mdl, "NIGLOBAL", n_global(1), & "The total number of thickness grid points in the "//& "x-direction in the physical domain. With STATIC_MEMORY_ "//& "this is set in "//trim(inc_nm)//" at compile time.", & static_value=NIGLOBAL) - call get_param(param_file, mdl, "NJGLOBAL", MOM_dom%njglobal, & + call get_param(param_file, mdl, "NJGLOBAL", n_global(2), & "The total number of thickness grid points in the "//& "y-direction in the physical domain. With STATIC_MEMORY_ "//& "this is set in "//trim(inc_nm)//" at compile time.", & static_value=NJGLOBAL) - if (MOM_dom%niglobal /= NIGLOBAL) call MOM_error(FATAL,"MOM_domains_init: " // & + if (n_global(1) /= NIGLOBAL) call MOM_error(FATAL,"MOM_domains_init: " // & "static mismatch for NIGLOBAL_ domain size. Header file does not match input namelist") - if (MOM_dom%njglobal /= NJGLOBAL) call MOM_error(FATAL,"MOM_domains_init: " // & + if (n_global(2) /= NJGLOBAL) call MOM_error(FATAL,"MOM_domains_init: " // & "static mismatch for NJGLOBAL_ domain size. Header file does not match input namelist") + ! Check the requirement of equal sized compute domains when STATIC_MEMORY_ is used. + if ((MOD(NIGLOBAL, NIPROC) /= 0) .OR. (MOD(NJGLOBAL, NJPROC) /= 0)) then + write( char_xsiz, '(i4)' ) NIPROC + write( char_ysiz, '(i4)' ) NJPROC + write( char_niglobal, '(i4)' ) NIGLOBAL + write( char_njglobal, '(i4)' ) NJGLOBAL + call MOM_error(WARNING, 'MOM_domains: Processor decomposition (NIPROC_,NJPROC_) = (' & + //trim(char_xsiz)//','//trim(char_ysiz)//') does not evenly divide size '//& + 'set by preprocessor macro ('//trim(char_niglobal)//','//trim(char_njglobal)//').') + call MOM_error(FATAL,'MOM_domains: #undef STATIC_MEMORY_ in '//trim(inc_nm)//' to use '//& + 'dynamic allocation, or change processor decomposition to evenly divide the domain.') + endif else - call get_param(param_file, mdl, "NIGLOBAL", MOM_dom%niglobal, & + call get_param(param_file, mdl, "NIGLOBAL", n_global(1), & "The total number of thickness grid points in the "//& "x-direction in the physical domain. With STATIC_MEMORY_ "//& "this is set in "//trim(inc_nm)//" at compile time.", & fail_if_missing=.true.) - call get_param(param_file, mdl, "NJGLOBAL", MOM_dom%njglobal, & + call get_param(param_file, mdl, "NJGLOBAL", n_global(2), & "The total number of thickness grid points in the "//& "y-direction in the physical domain. With STATIC_MEMORY_ "//& "this is set in "//trim(inc_nm)//" at compile time.", & fail_if_missing=.true.) endif - call get_param(param_file, mdl, trim(nihalo_nm), MOM_dom%nihalo, & + call get_param(param_file, mdl, trim(nihalo_nm), n_halo(1), & "The number of halo points on each side in the x-direction. How this is set "//& "varies with the calling component and static or dynamic memory configuration.", & default=nihalo_dflt, static_value=nihalo_dflt) - call get_param(param_file, mdl, trim(njhalo_nm), MOM_dom%njhalo, & + call get_param(param_file, mdl, trim(njhalo_nm), n_halo(2), & "The number of halo points on each side in the y-direction. How this is set "//& "varies with the calling component and static or dynamic memory configuration.", & default=njhalo_dflt, static_value=njhalo_dflt) if (present(min_halo)) then - MOM_dom%nihalo = max(MOM_dom%nihalo, min_halo(1)) - min_halo(1) = MOM_dom%nihalo - MOM_dom%njhalo = max(MOM_dom%njhalo, min_halo(2)) - min_halo(2) = MOM_dom%njhalo + n_halo(1) = max(n_halo(1), min_halo(1)) + min_halo(1) = n_halo(1) + n_halo(2) = max(n_halo(2), min_halo(2)) + min_halo(2) = n_halo(2) ! These are generally used only with static memory, so they are considerd layout params. - call log_param(param_file, mdl, "!NIHALO min_halo", MOM_dom%nihalo, layoutParam=.true.) - call log_param(param_file, mdl, "!NJHALO min_halo", MOM_dom%nihalo, layoutParam=.true.) + call log_param(param_file, mdl, "!NIHALO min_halo", n_halo(1), layoutParam=.true.) + call log_param(param_file, mdl, "!NJHALO min_halo", n_halo(2), layoutParam=.true.) endif if (is_static .and. .not.present(min_halo)) then - if (MOM_dom%nihalo /= NIHALO) call MOM_error(FATAL,"MOM_domains_init: " // & + if (n_halo(1) /= NIHALO) call MOM_error(FATAL,"MOM_domains_init: " // & "static mismatch for "//trim(nihalo_nm)//" domain size") - if (MOM_dom%njhalo /= NJHALO) call MOM_error(FATAL,"MOM_domains_init: " // & + if (n_halo(2) /= NJHALO) call MOM_error(FATAL,"MOM_domains_init: " // & "static mismatch for "//trim(njhalo_nm)//" domain size") endif - global_indices(1) = 1 ; global_indices(2) = MOM_dom%niglobal - global_indices(3) = 1 ; global_indices(4) = MOM_dom%njglobal + global_indices(1) = 1 ; global_indices(2) = n_global(1) + global_indices(3) = 1 ; global_indices(4) = n_global(2) call get_param(param_file, mdl, "INPUTDIR", inputdir, do_not_log=.true., default=".") inputdir = slasher(inputdir) @@ -1447,7 +1459,7 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & endif if ( layout(1)==0 .and. layout(2)==0 ) & - call mpp_define_layout(global_indices, proc_used, layout) + call MOM_define_layout(global_indices, proc_used, layout) if ( layout(1)/=0 .and. layout(2)==0 ) layout(2) = proc_used/layout(1) if ( layout(1)==0 .and. layout(2)/=0 ) layout(1) = proc_used/layout(2) @@ -1471,63 +1483,125 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & layoutParam=.true.) ! Idiot check that fewer PEs than columns have been requested - if (layout(1)*layout(2)>MOM_dom%niglobal*MOM_dom%njglobal) then + if (layout(1)*layout(2)>n_global(1)*n_global(2)) then write(mesg,'(a,2(i5,x,a))') 'You requested to use',layout(1)*layout(2), & - 'PEs but there are only',MOM_dom%niglobal*MOM_dom%njglobal,'columns in the model' + 'PEs but there are only', n_global(1)*n_global(2), 'columns in the model' call MOM_error(FATAL, mesg) endif - if (mask_table_exists) then - call MOM_error(NOTE, 'MOM_domains_init: reading maskmap information from '//& - trim(mask_table)) - allocate(MOM_dom%maskmap(layout(1), layout(2))) - call parse_mask_table(mask_table, MOM_dom%maskmap, dom_name) - endif + if (mask_table_exists) & + call MOM_error(NOTE, 'MOM_domains_init: reading maskmap information from '//trim(mask_table)) - ! Set up the I/O layout, and check that it uses an even multiple of the - ! number of PEs in each direction. + ! Set up the I/O layout, it will be checked later that it uses an even multiple of the number of + ! PEs in each direction. io_layout(:) = (/ 1, 1 /) call get_param(param_file, mdl, trim(io_layout_nm), io_layout, & "The processor layout to be used, or 0,0 to automatically "//& "set the io_layout to be the same as the layout.", default=1, & layoutParam=.true.) - if (io_layout(1) < 0) then - write(mesg,'("MOM_domains_init: IO_LAYOUT(1) = ",i4,". Negative values "//& - &"are not allowed in ")') io_layout(1) - call MOM_error(FATAL, mesg//trim(IO_layout_nm)) - elseif (io_layout(1) > 0) then ; if (modulo(layout(1), io_layout(1)) /= 0) then - write(mesg,'("MOM_domains_init: The i-direction I/O-layout, IO_LAYOUT(1)=",i4, & - &", does not evenly divide the i-direction layout, NIPROC=,",i4,".")') & - io_layout(1),layout(1) - call MOM_error(FATAL, mesg) - endif ; endif + call create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar, layout, io_layout=io_layout, & + domain_name=domain_name, mask_table=mask_table, symmetric=symmetric, & + thin_halos=thin_halos, nonblocking=nonblocking) - if (io_layout(2) < 0) then - write(mesg,'("MOM_domains_init: IO_LAYOUT(2) = ",i4,". Negative values "//& - &"are not allowed in ")') io_layout(2) - call MOM_error(FATAL, mesg//trim(IO_layout_nm)) - elseif (io_layout(2) /= 0) then ; if (modulo(layout(2), io_layout(2)) /= 0) then - write(mesg,'("MOM_domains_init: The j-direction I/O-layout, IO_LAYOUT(2)=",i4, & - &", does not evenly divide the j-direction layout, NJPROC=,",i4,".")') & - io_layout(2),layout(2) - call MOM_error(FATAL, mesg) - endif ; endif +end subroutine MOM_domains_init + +!> create_MOM_domain creates and initializes a MOM_domain_type variables, based on the information +!! provided in arguments. +subroutine create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar, layout, io_layout, & + domain_name, mask_table, symmetric, thin_halos, nonblocking) + type(MOM_domain_type), pointer :: MOM_dom !< A pointer to the MOM_domain_type being defined here. + integer, dimension(2), intent(in) :: n_global !< The number of points on the global grid in + !! the i- and j-directions + integer, dimension(2), intent(in) :: n_halo !< The number of halo points on each processor + logical, dimension(2), intent(in) :: reentrant !< If true the grid is periodic in the i- and j- directions + logical, dimension(2,2), intent(in) :: tripolar !< If true the grid uses tripolar connectivity on the two + !! ends (first index) of the i- and j-grids (second index) + integer, dimension(2), intent(in) :: layout !< The layout of logical PEs in the i- and j-directions. + integer, dimension(2), optional, intent(in) :: io_layout !< The layout for parallel input and output. + character(len=*), optional, intent(in) :: domain_name !< A name for this domain, "MOM" if missing. + character(len=*), optional, intent(in) :: mask_table !< The full relative or absolute path to the mask table. + logical, optional, intent(in) :: symmetric !< If present, this specifies whether this domain + !! uses symmetric memory, or true if missing. + logical, optional, intent(in) :: thin_halos !< If present, this specifies whether to permit the use of + !! thin halo updates, or true if missing. + logical, optional, intent(in) :: nonblocking !< If present, this specifies whether to permit the use of + !! nonblocking halo updates, or false if missing. + + ! local variables + integer, dimension(4) :: global_indices ! The lower and upper global i- and j-index bounds + integer :: X_FLAGS ! A combination of integers encoding the x-direction grid connectivity. + integer :: Y_FLAGS ! A combination of integers encoding the y-direction grid connectivity. + integer :: xhalo_d2, yhalo_d2 + character(len=200) :: mesg ! A string for use in error messages + character(len=64) :: dom_name ! The domain name + logical :: mask_table_exists ! Mask_table is present and the file it points to exists + + if (.not.associated(MOM_dom)) then + allocate(MOM_dom) + allocate(MOM_dom%mpp_domain) + allocate(MOM_dom%mpp_domain_d2) + endif - if (io_layout(2) == 0) io_layout(2) = layout(2) - if (io_layout(1) == 0) io_layout(1) = layout(1) + dom_name = "MOM" ; if (present(domain_name)) dom_name = trim(domain_name) X_FLAGS = 0 ; Y_FLAGS = 0 - if (reentrant_x) X_FLAGS = CYCLIC_GLOBAL_DOMAIN - if (reentrant_y) Y_FLAGS = CYCLIC_GLOBAL_DOMAIN - if (tripolar_N) then + if (reentrant(1)) X_FLAGS = CYCLIC_GLOBAL_DOMAIN + if (reentrant(2)) Y_FLAGS = CYCLIC_GLOBAL_DOMAIN + if (tripolar(2,2)) then Y_FLAGS = FOLD_NORTH_EDGE - if (reentrant_y) call MOM_error(FATAL,"MOM_domains: "// & - "TRIPOLAR_N and REENTRANT_Y may not be defined together.") + if (reentrant(2)) call MOM_error(FATAL,"MOM_domains: "// & + "TRIPOLAR_N and REENTRANT_Y may not be used together.") endif - global_indices(1) = 1 ; global_indices(2) = MOM_dom%niglobal - global_indices(3) = 1 ; global_indices(4) = MOM_dom%njglobal + MOM_dom%nonblocking_updates = nonblocking + MOM_dom%thin_halo_updates = thin_halos + MOM_dom%symmetric = .true. ; if (present(symmetric)) MOM_dom%symmetric = symmetric + MOM_dom%niglobal = n_global(1) ; MOM_dom%njglobal = n_global(2) + MOM_dom%nihalo = n_halo(1) ; MOM_dom%njhalo = n_halo(2) + + ! Save the extra data for creating other domains of different resolution that overlay this domain. + MOM_dom%X_FLAGS = X_FLAGS + MOM_dom%Y_FLAGS = Y_FLAGS + MOM_dom%layout(:) = layout(:) + + ! Set up the io_layout, with error handling. + MOM_dom%io_layout(:) = (/ 1, 1 /) + if (present(io_layout)) then + if (io_layout(1) == 0) then + MOM_dom%io_layout(1) = layout(1) + elseif (io_layout(1) > 1) then + MOM_dom%io_layout(1) = io_layout(1) + if (modulo(layout(1), io_layout(1)) /= 0) then + write(mesg,'("MOM_domains_init: The i-direction I/O-layout, IO_LAYOUT(1)=",i4, & + &", does not evenly divide the i-direction layout, NIPROC=,",i4,".")') io_layout(1), layout(1) + call MOM_error(FATAL, mesg) + endif + endif + + if (io_layout(2) == 0) then + MOM_dom%io_layout(2) = layout(2) + elseif (io_layout(2) > 1) then + MOM_dom%io_layout(2) = io_layout(2) + if (modulo(layout(2), io_layout(2)) /= 0) then + write(mesg,'("MOM_domains_init: The j-direction I/O-layout, IO_LAYOUT(2)=",i4, & + &", does not evenly divide the j-direction layout, NJPROC=,",i4,".")') io_layout(2), layout(2) + call MOM_error(FATAL, mesg) + endif + endif + endif + + global_indices(1:4) = (/ 1, MOM_dom%niglobal, 1, MOM_dom%njglobal /) + + if (present(mask_table)) then + mask_table_exists = file_exist(mask_table) + if (mask_table_exists) then + allocate(MOM_dom%maskmap(layout(1), layout(2))) + call parse_mask_table(mask_table, MOM_dom%maskmap, dom_name) + endif + else + mask_table_exists = .false. + endif if (mask_table_exists) then call MOM_define_domain( global_indices, layout, MOM_dom%mpp_domain, & @@ -1542,44 +1616,16 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & symmetry = MOM_dom%symmetric, name=dom_name) endif - if ((io_layout(1) > 0) .and. (io_layout(2) > 0) .and. & - (layout(1)*layout(2) > 1)) then - call MOM_define_io_domain(MOM_dom%mpp_domain, io_layout) - endif - -! Save the extra data for creating other domains of different resolution that overlay this domain - MOM_dom%X_FLAGS = X_FLAGS - MOM_dom%Y_FLAGS = Y_FLAGS - MOM_dom%layout = layout - MOM_dom%io_layout = io_layout - - if (is_static) then - ! A requirement of equal sized compute domains is necessary when STATIC_MEMORY_ - ! is used. - call mpp_get_compute_domain(MOM_dom%mpp_domain,isc,iec,jsc,jec) - xsiz = iec - isc + 1 - ysiz = jec - jsc + 1 - if (xsiz*NIPROC /= MOM_dom%niglobal .OR. ysiz*NJPROC /= MOM_dom%njglobal) then - write( char_xsiz,'(i4)' ) NIPROC - write( char_ysiz,'(i4)' ) NJPROC - write( char_niglobal,'(i4)' ) MOM_dom%niglobal - write( char_njglobal,'(i4)' ) MOM_dom%njglobal - call MOM_error(WARNING,'MOM_domains: Processor decomposition (NIPROC_,NJPROC_) = (' & - //trim(char_xsiz)//','//trim(char_ysiz)// & - ') does not evenly divide size set by preprocessor macro ('& - //trim(char_niglobal)//','//trim(char_njglobal)// '). ') - call MOM_error(FATAL,'MOM_domains: #undef STATIC_MEMORY_ in "//trim(inc_nm)//" to use & - &dynamic allocation, or change processor decomposition to evenly divide the domain.') - endif + if ((MOM_dom%io_layout(1) > 0) .and. (MOM_dom%io_layout(2) > 0) .and. (layout(1)*layout(2) > 1)) then + call MOM_define_io_domain(MOM_dom%mpp_domain, MOM_dom%io_layout) endif - global_indices(1) = 1 ; global_indices(2) = int(MOM_dom%niglobal/2) - global_indices(3) = 1 ; global_indices(4) = int(MOM_dom%njglobal/2) !For downsampled domain, recommend a halo of 1 (or 0?) since we're not doing wide-stencil computations. !But that does not work because the downsampled field would not have the correct size to pass the checks, e.g., we get !error: downsample_diag_indices_get: peculiar size 28 in i-direction\ndoes not match one of 24 25 26 27 xhalo_d2 = int(MOM_dom%nihalo/2) yhalo_d2 = int(MOM_dom%njhalo/2) + global_indices(1:4) = (/ 1, int(MOM_dom%niglobal/2), 1, int(MOM_dom%njglobal/2) /) if (mask_table_exists) then call MOM_define_domain( global_indices, layout, MOM_dom%mpp_domain_d2, & xflags=X_FLAGS, yflags=Y_FLAGS, & @@ -1593,12 +1639,98 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & symmetry = MOM_dom%symmetric, name=trim("MOMc")) endif - if ((io_layout(1) > 0) .and. (io_layout(2) > 0) .and. & + if ((MOM_dom%io_layout(1) > 0) .and. (MOM_dom%io_layout(2) > 0) .and. & (layout(1)*layout(2) > 1)) then - call MOM_define_io_domain(MOM_dom%mpp_domain_d2, io_layout) + call MOM_define_io_domain(MOM_dom%mpp_domain_d2, MOM_dom%io_layout) endif -end subroutine MOM_domains_init +end subroutine create_MOM_domain + +!> dealloc_MOM_domain deallocates memory associated with a pointer to a MOM_domain_type +!! and all of its contents +subroutine deallocate_MOM_domain(MOM_domain, cursory) + type(MOM_domain_type), pointer :: MOM_domain !< A pointer to the MOM_domain_type being deallocated + logical, optional, intent(in) :: cursory !< If true do not deallocate fields associated + !! with the underlying infrastructure + + if (associated(MOM_domain)) then + call deallocate_domain_contents(MOM_domain, cursory) + deallocate(MOM_domain) + endif + +end subroutine deallocate_MOM_domain + +!> deallocate_domain_contents deallocates memory associated with pointers +!! inside of a MOM_domain_type. +subroutine deallocate_domain_contents(MOM_domain, cursory) + type(MOM_domain_type), intent(inout) :: MOM_domain !< A MOM_domain_type whose contents will be deallocated + logical, optional, intent(in) :: cursory !< If true do not deallocate fields associated + !! with the underlying infrastructure + + logical :: invasive ! If true, deallocate fields associated with the underlying infrastructure + + invasive = .true. ; if (present(cursory)) invasive = .not.cursory + + if (associated(MOM_domain%mpp_domain)) then + if (invasive) call mpp_deallocate_domain(MOM_domain%mpp_domain) + deallocate(MOM_domain%mpp_domain) + endif + if (associated(MOM_domain%mpp_domain_d2)) then + if (invasive) call mpp_deallocate_domain(MOM_domain%mpp_domain_d2) + deallocate(MOM_domain%mpp_domain_d2) + endif + if (associated(MOM_domain%maskmap)) deallocate(MOM_domain%maskmap) + +end subroutine deallocate_domain_contents + +!> MOM_thread_affinity_set returns true if the number of openMP threads have been set to a value greater than 1. +function MOM_thread_affinity_set() + ! Local variables + !$ integer :: ocean_nthreads ! Number of openMP threads + !$ integer :: omp_get_num_threads ! An openMP function that returns the number of threads + logical :: MOM_thread_affinity_set + + MOM_thread_affinity_set = .false. + !$ call fms_affinity_init() + !$OMP PARALLEL + !$OMP MASTER + !$ ocean_nthreads = omp_get_num_threads() + !$OMP END MASTER + !$OMP END PARALLEL + !$ MOM_thread_affinity_set = (ocean_nthreads > 1 ) +end function MOM_thread_affinity_set + +!> set_MOM_thread_affinity sest the number of openMP threads to use with the ocean. +subroutine set_MOM_thread_affinity(ocean_nthreads, ocean_hyper_thread) + integer, intent(in) :: ocean_nthreads !< Number of openMP threads to use for the ocean model + logical, intent(in) :: ocean_hyper_thread !< If true, use hyper threading + + ! Local variables + !$ integer :: omp_get_thread_num, omp_get_num_threads !< These are the results of openMP functions + + !$ call fms_affinity_set('OCEAN', ocean_hyper_thread, ocean_nthreads) + !$ call omp_set_num_threads(ocean_nthreads) + !$ write(6,*) "MOM_domains_mod OMPthreading ", fms_affinity_get(), omp_get_thread_num(), omp_get_num_threads() + !$ flush(6) +end subroutine set_MOM_thread_affinity + +!> This subroutine retrieves the 1-d domains that make up the 2d-domain in a MOM_domain +subroutine get_domain_components_MD(MOM_dom, x_domain, y_domain) + type(MOM_domain_type), intent(in) :: MOM_dom !< The MOM_domain whose contents are being extracted + type(domain1D), optional, intent(inout) :: x_domain !< The 1-d logical x-domain + type(domain1D), optional, intent(inout) :: y_domain !< The 1-d logical y-domain + + call mpp_get_domain_components(MOM_dom%mpp_domain, x_domain, y_domain) +end subroutine get_domain_components_MD + +!> This subroutine retrieves the 1-d domains that make up a 2d-domain +subroutine get_domain_components_d2D(domain, x_domain, y_domain) + type(domain2D), intent(in) :: domain !< The 2D domain whose contents are being extracted + type(domain1D), optional, intent(inout) :: x_domain !< The 1-d logical x-domain + type(domain1D), optional, intent(inout) :: y_domain !< The 1-d logical y-domain + + call mpp_get_domain_components(domain, x_domain, y_domain) +end subroutine get_domain_components_d2D !> clone_MD_to_MD copies one MOM_domain_type into another, while allowing !! some properties of the new type to differ from the original one. @@ -1970,13 +2102,29 @@ end subroutine get_simple_array_j_ind !> Returns the global shape of h-point arrays subroutine get_global_shape(domain, niglobal, njglobal) - type(MOM_domain_type), intent(in) :: domain !< MOM domain + type(MOM_domain_type), intent(in) :: domain !< MOM domain from which to extract information integer, intent(out) :: niglobal !< i-index global size of h-point arrays integer, intent(out) :: njglobal !< j-index global size of h-point arrays niglobal = domain%niglobal njglobal = domain%njglobal - end subroutine get_global_shape +!> Returns arrays of the i- and j- sizes of the h-point computational domains for each +!! element of the grid layout. Any input values in the extent arrays are discarded, so +!! they are effectively intent out despite their declared intent of inout. +subroutine get_layout_extents(Domain, extent_i, extent_j) + type(MOM_domain_type), intent(in) :: domain !< MOM domain from which to extract information + integer, dimension(:), allocatable, intent(inout) :: extent_i !< The number of points in the + !! i-direction in each i-row of the layout + integer, dimension(:), allocatable, intent(inout) :: extent_j !< The number of points in the + !! j-direction in each j-row of the layout + + if (allocated(extent_i)) deallocate(extent_i) + if (allocated(extent_j)) deallocate(extent_j) + allocate(extent_i(domain%layout(1))) ; extent_i(:) = 0 + allocate(extent_j(domain%layout(2))) ; extent_j(:) = 0 + call mpp_get_domain_extents(domain%mpp_domain, extent_i, extent_j) +end subroutine get_layout_extents + end module MOM_domains diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index f755b7f675..c7d7e98e4b 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -3,52 +3,46 @@ module MOM_io ! This file is part of MOM6. See LICENSE.md for the license. - +use MOM_domains, only : MOM_domain_type, AGRID, BGRID_NE, CGRID_NE, get_domain_components +use MOM_domains, only : domain1D, get_simple_array_i_ind, get_simple_array_j_ind +use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING -use MOM_domains, only : MOM_domain_type, AGRID, BGRID_NE, CGRID_NE -use MOM_domains, only : get_simple_array_i_ind, get_simple_array_j_ind use MOM_file_parser, only : log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_io_wrapper, only : MOM_read_data, MOM_read_vector, MOM_write_field, read_axis_data +use MOM_io_wrapper, only : file_exists, field_exists, read_field_chksum +use MOM_io_wrapper, only : open_file, close_file, field_size, fieldtype, get_filename_appendix +use MOM_io_wrapper, only : flush_file, get_file_info, get_file_atts, get_file_fields +use MOM_io_wrapper, only : get_file_times, read_data, axistype, get_axis_data +use MOM_io_wrapper, only : write_field, write_metadata, write_version_number, get_ensemble_id +use MOM_io_wrapper, only : open_namelist_file, check_nml_error, io_infra_init, io_infra_end +use MOM_io_wrapper, only : APPEND_FILE, ASCII_FILE, MULTIPLE, NETCDF_FILE, OVERWRITE_FILE +use MOM_io_wrapper, only : READONLY_FILE, SINGLE_FILE, WRITEONLY_FILE +use MOM_io_wrapper, only : CENTER, CORNER, NORTH_FACE, EAST_FACE use MOM_string_functions, only : lowercase, slasher use MOM_verticalGrid, only : verticalGrid_type -use ensemble_manager_mod, only : get_ensemble_id -use fms_mod, only : write_version_number, open_namelist_file, check_nml_error -use fms_io_mod, only : file_exist, field_size, read_data -use fms_io_mod, only : field_exists=>field_exist, io_infra_end=>fms_io_exit -use fms_io_mod, only : get_filename_appendix=>get_filename_appendix -use mpp_domains_mod, only : domain1d, domain2d, mpp_get_domain_components -use mpp_domains_mod, only : CENTER, CORNER, NORTH_FACE=>NORTH, EAST_FACE=>EAST -use mpp_io_mod, only : open_file => mpp_open, close_file => mpp_close -use mpp_io_mod, only : mpp_write_meta, write_field => mpp_write -use mpp_io_mod, only : mpp_get_atts, mpp_attribute_exist -use mpp_io_mod, only : mpp_get_axes, axistype, get_axis_data=>mpp_get_axis_data -use mpp_io_mod, only : mpp_get_fields, fieldtype, flush_file=>mpp_flush -use mpp_io_mod, only : APPEND_FILE=>MPP_APPEND, ASCII_FILE=>MPP_ASCII -use mpp_io_mod, only : MULTIPLE=>MPP_MULTI, NETCDF_FILE=>MPP_NETCDF -use mpp_io_mod, only : OVERWRITE_FILE=>MPP_OVERWR, READONLY_FILE=>MPP_RDONLY -use mpp_io_mod, only : SINGLE_FILE=>MPP_SINGLE, WRITEONLY_FILE=>MPP_WRONLY -use mpp_io_mod, only : get_file_info=>mpp_get_info, get_file_atts=>mpp_get_atts -use mpp_io_mod, only : get_file_fields=>mpp_get_fields, get_file_times=>mpp_get_times -use mpp_io_mod, only : io_infra_init=>mpp_io_init - use iso_fortran_env, only : stdout_iso=>output_unit, stderr_iso=>error_unit -use netcdf +use netcdf, only : NF90_open, NF90_inquire, NF90_inq_varids, NF90_inquire_variable +use netcdf, only : NF90_Inquire_Dimension, NF90_max_name, NF90_max_var_dims +use netcdf, only : NF90_STRERROR, NF90_NOWRITE, NF90_NOERR implicit none ; private -public :: close_file, create_file, field_exists, field_size, fieldtype, get_filename_appendix +! These interfaces are actually implemented in this file. +public :: create_file, reopen_file, num_timelevels, cmor_long_std, ensembler, MOM_io_init +public :: var_desc, modify_vardesc, query_vardesc +! The following are simple pass throughs of routines from MOM_io_wrapper or other modules +public :: close_file, field_exists, field_size, fieldtype, get_filename_appendix public :: file_exists, flush_file, get_file_info, get_file_atts, get_file_fields public :: get_file_times, open_file, read_axis_data, read_data, read_field_chksum -public :: num_timelevels, MOM_read_data, MOM_read_vector, ensembler -public :: reopen_file, slasher, write_field, write_version_number, MOM_io_init +public :: MOM_read_data, MOM_read_vector, MOM_write_field, get_axis_data +public :: slasher, write_field, write_version_number public :: open_namelist_file, check_nml_error, io_infra_init, io_infra_end +! These are encoding constants. public :: APPEND_FILE, ASCII_FILE, MULTIPLE, NETCDF_FILE, OVERWRITE_FILE public :: READONLY_FILE, SINGLE_FILE, WRITEONLY_FILE public :: CENTER, CORNER, NORTH_FACE, EAST_FACE -public :: var_desc, modify_vardesc, query_vardesc, cmor_long_std -public :: get_axis_data !> Type for describing a variable, typically a tracer type, public :: vardesc @@ -65,27 +59,6 @@ module MOM_io !! convert from intensive to extensive end type vardesc -!> Indicate whether a file exists, perhaps with domain decomposition -interface file_exists - module procedure FMS_file_exists - module procedure MOM_file_exists -end interface - -!> Read a data field from a file -interface MOM_read_data - module procedure MOM_read_data_4d - module procedure MOM_read_data_3d - module procedure MOM_read_data_2d - module procedure MOM_read_data_1d - module procedure MOM_read_data_0d -end interface - -!> Read a pair of data fields representing the two components of a vector from a file -interface MOM_read_vector - module procedure MOM_read_vector_3d - module procedure MOM_read_vector_2d -end interface - integer, public :: stdout = stdout_iso !< standard output unit integer, public :: stderr = stderr_iso !< standard output unit @@ -165,7 +138,7 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit if (one_file) then call open_file(unit, filename, OVERWRITE_FILE, NETCDF_FILE, threading=thread) else - call open_file(unit, filename, OVERWRITE_FILE, NETCDF_FILE, domain=Domain%mpp_domain) + call open_file(unit, filename, OVERWRITE_FILE, NETCDF_FILE, MOM_domain=Domain) endif ! Define the coordinates. @@ -229,39 +202,36 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit if (.not.domain_set) call MOM_error(FATAL, "create_file: "//& "An ocean_grid_type or dyn_horgrid_type is required to create a file with a horizontal coordinate.") - call mpp_get_domain_components(Domain%mpp_domain, x_domain, y_domain) + call get_domain_components(Domain, x_domain, y_domain) endif if ((use_layer .or. use_int) .and. .not.present(GV)) call MOM_error(FATAL, & "create_file: A vertical grid type is required to create a file with a vertical coordinate.") -! Specify all optional arguments to mpp_write_meta: name, units, longname, cartesian, calendar, sense, -! domain, data, min). Otherwise if optional arguments are added to mpp_write_meta the compiler may -! (and in case of GNU does) get confused and crash. if (use_lath) & - call mpp_write_meta(unit, axis_lath, name="lath", units=y_axis_units, longname="Latitude", & - cartesian='Y', domain = y_domain, data=gridLatT(jsg:jeg)) + call write_metadata(unit, axis_lath, name="lath", units=y_axis_units, longname="Latitude", & + cartesian='Y', domain=y_domain, data=gridLatT(jsg:jeg)) if (use_lonh) & - call mpp_write_meta(unit, axis_lonh, name="lonh", units=x_axis_units, longname="Longitude", & - cartesian='X', domain = x_domain, data=gridLonT(isg:ieg)) + call write_metadata(unit, axis_lonh, name="lonh", units=x_axis_units, longname="Longitude", & + cartesian='X', domain=x_domain, data=gridLonT(isg:ieg)) if (use_latq) & - call mpp_write_meta(unit, axis_latq, name="latq", units=y_axis_units, longname="Latitude", & - cartesian='Y', domain = y_domain, data=gridLatB(JsgB:JegB)) + call write_metadata(unit, axis_latq, name="latq", units=y_axis_units, longname="Latitude", & + cartesian='Y', domain=y_domain, data=gridLatB(JsgB:JegB)) if (use_lonq) & - call mpp_write_meta(unit, axis_lonq, name="lonq", units=x_axis_units, longname="Longitude", & - cartesian='X', domain = x_domain, data=gridLonB(IsgB:IegB)) + call write_metadata(unit, axis_lonq, name="lonq", units=x_axis_units, longname="Longitude", & + cartesian='X', domain=x_domain, data=gridLonB(IsgB:IegB)) if (use_layer) & - call mpp_write_meta(unit, axis_layer, name="Layer", units=trim(GV%zAxisUnits), & - longname="Layer "//trim(GV%zAxisLongName), cartesian='Z', & - sense=1, data=GV%sLayer(1:GV%ke)) + call write_metadata(unit, axis_layer, name="Layer", units=trim(GV%zAxisUnits), & + longname="Layer "//trim(GV%zAxisLongName), cartesian='Z', & + sense=1, data=GV%sLayer(1:GV%ke)) if (use_int) & - call mpp_write_meta(unit, axis_int, name="Interface", units=trim(GV%zAxisUnits), & - longname="Interface "//trim(GV%zAxisLongName), cartesian='Z', & - sense=1, data=GV%sInterface(1:GV%ke+1)) + call write_metadata(unit, axis_int, name="Interface", units=trim(GV%zAxisUnits), & + longname="Interface "//trim(GV%zAxisLongName), cartesian='Z', & + sense=1, data=GV%sInterface(1:GV%ke+1)) if (use_time) then ; if (present(timeunit)) then ! Set appropriate units, depending on the value. @@ -279,9 +249,9 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit write(time_units,'(es8.2," s")') timeunit endif - call mpp_write_meta(unit, axis_time, name="Time", units=time_units, longname="Time", cartesian='T') + call write_metadata(unit, axis_time, name="Time", units=time_units, longname="Time", cartesian='T') else - call mpp_write_meta(unit, axis_time, name="Time", units="days", longname="Time",cartesian= 'T') + call write_metadata(unit, axis_time, name="Time", units="days", longname="Time", cartesian= 'T') endif ; endif if (use_periodic) then @@ -290,8 +260,8 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit ! Define a periodic axis with unit labels. allocate(period_val(num_periods)) do k=1,num_periods ; period_val(k) = real(k) ; enddo - call mpp_write_meta(unit, axis_periodic, name="Period", units="nondimensional", & - longname="Periods for cyclical varaiables", cartesian= 't', data=period_val) + call write_metadata(unit, axis_periodic, name="Period", units="nondimensional", & + longname="Periods for cyclical varaiables", cartesian='T', data=period_val) deallocate(period_val) endif @@ -328,14 +298,14 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit call MOM_error(WARNING, "MOM_io create_file: "//trim(vars(k)%name)//& " has unrecognized t_grid "//trim(vars(k)%t_grid)) end select - pack = 1 + pack = 1 if (present(checksums)) then - call mpp_write_meta(unit, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & - vars(k)%longname, pack = pack, checksum=checksums(k,:)) + call write_metadata(unit, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & + vars(k)%longname, pack=pack, checksum=checksums(k,:)) else - call mpp_write_meta(unit, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & - vars(k)%longname, pack = pack) + call write_metadata(unit, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & + vars(k)%longname, pack=pack) endif enddo @@ -407,7 +377,7 @@ subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit if (one_file) then call open_file(unit, filename, APPEND_FILE, NETCDF_FILE, threading=thread) else - call open_file(unit, filename, APPEND_FILE, NETCDF_FILE, domain=Domain%mpp_domain) + call open_file(unit, filename, APPEND_FILE, NETCDF_FILE, MOM_domain=Domain) endif if (unit < 0) return @@ -424,76 +394,21 @@ subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit call MOM_error(FATAL,"MOM_io: "//mesg) endif - if (nvar>0) call mpp_get_fields(unit,fields(1:nvar)) + if (nvar > 0) call get_file_fields(unit, fields(1:nvar)) - ! Check the field names... + ! Check for inconsistent field names... ! do i=1,nvar -! call mpp_get_field_atts(fields(i),name) -! !if (trim(name) /= trim(vars%name) then -! !write (mesg,'("Reopening file ",a," variable ",a," is called ",a,".")',& -! ! filename,vars%name,name) -! !call MOM_error(NOTE,"MOM_io: "//mesg) +! call get_field_atts(fields(i), name) +! !if (trim(name) /= trim(vars%name)) then +! ! write (mesg, '("Reopening file ",a," variable ",a," is called ",a,".")',& +! ! trim(filename), trim(vars%name), trim(name)) +! ! call MOM_error(NOTE, "MOM_io: "//trim(mesg)) +! !endif ! enddo endif end subroutine reopen_file -!> Read the data associated with a named axis in a file -subroutine read_axis_data(filename, axis_name, var) - character(len=*), intent(in) :: filename !< Name of the file to read - character(len=*), intent(in) :: axis_name !< Name of the axis to read - real, dimension(:), intent(out) :: var !< The axis location data - - integer :: i,len,unit, ndim, nvar, natt, ntime - logical :: axis_found - type(axistype), allocatable :: axes(:) - type(axistype) :: time_axis - character(len=32) :: name, units - - call open_file(unit, trim(filename), action=READONLY_FILE, form=NETCDF_FILE, & - threading=MULTIPLE, fileset=SINGLE_FILE) - -!Find the number of variables (nvar) in this file - call get_file_info(unit, ndim, nvar, natt, ntime) -! ------------------------------------------------------------------- -! Allocate space for the number of axes in the data file. -! ------------------------------------------------------------------- - allocate(axes(ndim)) - call mpp_get_axes(unit, axes, time_axis) - - axis_found = .false. - do i = 1, ndim - call get_file_atts(axes(i), name=name, len=len, units=units) - if (name == axis_name) then - axis_found = .true. - call get_axis_data(axes(i),var) - exit - endif - enddo - - if (.not.axis_found) call MOM_error(FATAL, "MOM_io read_axis_data: "//& - "Unable to find axis "//trim(axis_name)//" in file "//trim(filename)) - - deallocate(axes) - -end subroutine read_axis_data - -subroutine read_field_chksum(field, chksum, valid_chksum) - type(fieldtype), intent(in) :: field !< The field whose checksum attribute is to be read. - integer(kind=8), intent(out) :: chksum !< The checksum for the field. - logical, intent(out) :: valid_chksum !< If true, chksum has been successfully read. - ! Local variables - integer(kind=8), dimension(3) :: checksum_file - - checksum_file(:) = -1 - valid_chksum = mpp_attribute_exist(field, "checksum") - if (valid_chksum) then - call mpp_get_atts(field, checksum=checksum_file) - chksum = checksum_file(1) - else - chksum = -1 - endif -end subroutine read_field_chksum !> This function determines how many time levels a variable has. function num_timelevels(filename, varname, min_dims) result(n_time) @@ -518,8 +433,7 @@ function num_timelevels(filename, varname, min_dims) result(n_time) status = NF90_OPEN(filename, NF90_NOWRITE, ncid) if (status /= NF90_NOERR) then call MOM_error(WARNING,"num_timelevels: "//& - " Difficulties opening "//trim(filename)//" - "//& - trim(NF90_STRERROR(status))) + " Difficulties opening "//trim(filename)//" - "//trim(NF90_STRERROR(status))) return endif @@ -570,16 +484,14 @@ function num_timelevels(filename, varname, min_dims) result(n_time) if (.not.found) then call MOM_error(WARNING,"num_timelevels: "//& - " variable "//trim(varname)//" was not found in file "//& - trim(filename)) + " variable "//trim(varname)//" was not found in file "//trim(filename)) return endif status = nf90_inquire_variable(ncid, varid, ndims = ndims) if (status /= NF90_NOERR) then - call MOM_error(WARNING,"num_timelevels: "//& - trim(NF90_STRERROR(status))//" Getting number of dimensions of "//& - trim(varname)//" in "//trim(filename)) + call MOM_error(WARNING,"num_timelevels: "//trim(NF90_STRERROR(status))//& + " Getting number of dimensions of "//trim(varname)//" in "//trim(filename)) return endif @@ -596,9 +508,8 @@ function num_timelevels(filename, varname, min_dims) result(n_time) status = nf90_inquire_variable(ncid, varid, dimids = dimids(1:ndims)) if (status /= NF90_NOERR) then - call MOM_error(WARNING,"num_timelevels: "//& - trim(NF90_STRERROR(status))//" Getting last dimension ID for "//& - trim(varname)//" in "//trim(filename)) + call MOM_error(WARNING,"num_timelevels: "//trim(NF90_STRERROR(status))//& + " Getting last dimension ID for "//trim(varname)//" in "//trim(filename)) return endif @@ -607,8 +518,6 @@ function num_timelevels(filename, varname, min_dims) result(n_time) trim(NF90_STRERROR(status))//" Getting number of time levels of "//& trim(varname)//" in "//trim(filename)) - return - end function num_timelevels @@ -621,7 +530,7 @@ function var_desc(name, units, longname, hor_grid, z_grid, t_grid, & character(len=*), intent(in) :: name !< variable name character(len=*), optional, intent(in) :: units !< variable units character(len=*), optional, intent(in) :: longname !< variable long name - character(len=*), optional, intent(in) :: hor_grid !< variable horizonal staggering + character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1 character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name @@ -662,7 +571,7 @@ subroutine modify_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & character(len=*), optional, intent(in) :: name !< name of variable character(len=*), optional, intent(in) :: units !< units of variable character(len=*), optional, intent(in) :: longname !< long name of variable - character(len=*), optional, intent(in) :: hor_grid !< horizonal staggering of variable + character(len=*), optional, intent(in) :: hor_grid !< horizontal staggering of variable character(len=*), optional, intent(in) :: z_grid !< vertical staggering of variable character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1 character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name @@ -721,8 +630,8 @@ subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & character(len=*), optional, intent(out) :: name !< name of variable character(len=*), optional, intent(out) :: units !< units of variable character(len=*), optional, intent(out) :: longname !< long name of variable - character(len=*), optional, intent(out) :: hor_grid !< horiz staggering of variable - character(len=*), optional, intent(out) :: z_grid !< vert staggering of variable + character(len=*), optional, intent(out) :: hor_grid !< horizontal staggering of variable + character(len=*), optional, intent(out) :: z_grid !< verticle staggering of variable character(len=*), optional, intent(out) :: t_grid !< time description: s, p, or 1 character(len=*), optional, intent(out) :: cmor_field_name !< CMOR name character(len=*), optional, intent(out) :: cmor_units !< CMOR physical dimensions of variable @@ -758,7 +667,6 @@ subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, & end subroutine query_vardesc - !> Copies a string subroutine safe_string_copy(str1, str2, fieldnm, caller) character(len=*), intent(in) :: str1 !< The string being copied @@ -778,7 +686,6 @@ subroutine safe_string_copy(str1, str2, fieldnm, caller) str2 = trim(str1) end subroutine safe_string_copy - !> Returns a name with "%#E" or "%E" replaced with the ensemble member number. function ensembler(name, ens_no_in) result(en_nm) character(len=*), intent(in) :: name !< The name to be modified @@ -836,262 +743,19 @@ function ensembler(name, ens_no_in) result(en_nm) end function ensembler - -!> Returns true if the named file or its domain-decomposed variant exists. -function MOM_file_exists(filename, MOM_Domain) - character(len=*), intent(in) :: filename !< The name of the file being inquired about - type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition - -! This function uses the fms_io function file_exist to determine whether -! a named file (or its decomposed variant) exists. - - logical :: MOM_file_exists - - MOM_file_exists = file_exist(filename, MOM_Domain%mpp_domain) - -end function MOM_file_exists - -!> Returns true if the named file or its domain-decomposed variant exists. -function FMS_file_exists(filename, domain, no_domain) - character(len=*), intent(in) :: filename !< The name of the file being inquired about - type(domain2d), optional, intent(in) :: domain !< The mpp domain2d that describes the decomposition - logical, optional, intent(in) :: no_domain !< This file does not use domain decomposition -! This function uses the fms_io function file_exist to determine whether -! a named file (or its decomposed variant) exists. - - logical :: FMS_file_exists - - FMS_file_exists = file_exist(filename, domain, no_domain) - -end function FMS_file_exists - - -!> This function uses the fms_io function read_data to read a scalar -!! data field named "fieldname" from file "filename". -subroutine MOM_read_data_0d(filename, fieldname, data, timelevel, scale) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, intent(inout) :: data !< The 1-dimensional array into which the data - integer, optional, intent(in) :: timelevel !< The time level in the file to read - real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied - !! by before it is returned. - - call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) - - if (present(scale)) then ; if (scale /= 1.0) then - data = scale*data - endif ; endif - -end subroutine MOM_read_data_0d - -!> This function uses the fms_io function read_data to read a 1-D -!! data field named "fieldname" from file "filename". -subroutine MOM_read_data_1d(filename, fieldname, data, timelevel, scale) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, dimension(:), intent(inout) :: data !< The 1-dimensional array into which the data - integer, optional, intent(in) :: timelevel !< The time level in the file to read - real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied - !! by before they are returned. - - call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) - - if (present(scale)) then ; if (scale /= 1.0) then - data(:) = scale*data(:) - endif ; endif - -end subroutine MOM_read_data_1d - -!> This function uses the fms_io function read_data to read a distributed -!! 2-D data field named "fieldname" from file "filename". Valid values for -!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. -subroutine MOM_read_data_2d(filename, fieldname, data, MOM_Domain, & - timelevel, position, scale) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, dimension(:,:), intent(inout) :: data !< The 2-dimensional array into which the data - !! should be read - type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition - integer, optional, intent(in) :: timelevel !< The time level in the file to read - integer, optional, intent(in) :: position !< A flag indicating where this data is located - real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied - !! by before it is returned. - - integer :: is, ie, js, je - - call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & - timelevel=timelevel, position=position) - - if (present(scale)) then ; if (scale /= 1.0) then - call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) - call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) - data(is:ie,js:je) = scale*data(is:ie,js:je) - endif ; endif - -end subroutine MOM_read_data_2d - -!> This function uses the fms_io function read_data to read a distributed -!! 3-D data field named "fieldname" from file "filename". Valid values for -!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. -subroutine MOM_read_data_3d(filename, fieldname, data, MOM_Domain, & - timelevel, position, scale) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional array into which the data - !! should be read - type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition - integer, optional, intent(in) :: timelevel !< The time level in the file to read - integer, optional, intent(in) :: position !< A flag indicating where this data is located - real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied - !! by before it is returned. - - integer :: is, ie, js, je - - call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & - timelevel=timelevel, position=position) - - if (present(scale)) then ; if (scale /= 1.0) then - call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) - call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) - data(is:ie,js:je,:) = scale*data(is:ie,js:je,:) - endif ; endif - -end subroutine MOM_read_data_3d - -!> This function uses the fms_io function read_data to read a distributed -!! 4-D data field named "fieldname" from file "filename". Valid values for -!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. -subroutine MOM_read_data_4d(filename, fieldname, data, MOM_Domain, & - timelevel, position, scale) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: fieldname !< The variable name of the data in the file - real, dimension(:,:,:,:), intent(inout) :: data !< The 4-dimensional array into which the data - !! should be read - type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition - integer, optional, intent(in) :: timelevel !< The time level in the file to read - integer, optional, intent(in) :: position !< A flag indicating where this data is located - real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied - !! by before it is returned. - - integer :: is, ie, js, je - - call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & - timelevel=timelevel, position=position) - - if (present(scale)) then ; if (scale /= 1.0) then - call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) - call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) - data(is:ie,js:je,:,:) = scale*data(is:ie,js:je,:,:) - endif ; endif - -end subroutine MOM_read_data_4d - - -!> This function uses the fms_io function read_data to read a pair of distributed -!! 2-D data fields with names given by "[uv]_fieldname" from file "filename". Valid values for -!! "stagger" include CGRID_NE, BGRID_NE, and AGRID. -subroutine MOM_read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, & - timelevel, stagger, scalar_pair, scale) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: u_fieldname !< The variable name of the u data in the file - character(len=*), intent(in) :: v_fieldname !< The variable name of the v data in the file - real, dimension(:,:), intent(inout) :: u_data !< The 2 dimensional array into which the - !! u-component of the data should be read - real, dimension(:,:), intent(inout) :: v_data !< The 2 dimensional array into which the - !! v-component of the data should be read - type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition - integer, optional, intent(in) :: timelevel !< The time level in the file to read - integer, optional, intent(in) :: stagger !< A flag indicating where this vector is discretized - logical, optional, intent(in) :: scalar_pair !< If true, a pair of scalars are to be read.cretized - real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied - !! by before they are returned. - integer :: is, ie, js, je - integer :: u_pos, v_pos - - u_pos = EAST_FACE ; v_pos = NORTH_FACE - if (present(stagger)) then - if (stagger == CGRID_NE) then ; u_pos = EAST_FACE ; v_pos = NORTH_FACE - elseif (stagger == BGRID_NE) then ; u_pos = CORNER ; v_pos = CORNER - elseif (stagger == AGRID) then ; u_pos = CENTER ; v_pos = CENTER ; endif - endif - - call read_data(filename, u_fieldname, u_data, MOM_Domain%mpp_domain, & - timelevel=timelevel, position=u_pos) - call read_data(filename, v_fieldname, v_data, MOM_Domain%mpp_domain, & - timelevel=timelevel, position=v_pos) - - if (present(scale)) then ; if (scale /= 1.0) then - call get_simple_array_i_ind(MOM_Domain, size(u_data,1), is, ie) - call get_simple_array_j_ind(MOM_Domain, size(u_data,2), js, je) - u_data(is:ie,js:je) = scale*u_data(is:ie,js:je) - call get_simple_array_i_ind(MOM_Domain, size(v_data,1), is, ie) - call get_simple_array_j_ind(MOM_Domain, size(v_data,2), js, je) - v_data(is:ie,js:je) = scale*v_data(is:ie,js:je) - endif ; endif - -end subroutine MOM_read_vector_2d - - -!> This function uses the fms_io function read_data to read a pair of distributed -!! 3-D data fields with names given by "[uv]_fieldname" from file "filename". Valid values for -!! "stagger" include CGRID_NE, BGRID_NE, and AGRID. -subroutine MOM_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, & - timelevel, stagger, scalar_pair, scale) - character(len=*), intent(in) :: filename !< The name of the file to read - character(len=*), intent(in) :: u_fieldname !< The variable name of the u data in the file - character(len=*), intent(in) :: v_fieldname !< The variable name of the v data in the file - real, dimension(:,:,:), intent(inout) :: u_data !< The 3 dimensional array into which the - !! u-component of the data should be read - real, dimension(:,:,:), intent(inout) :: v_data !< The 3 dimensional array into which the - !! v-component of the data should be read - type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition - integer, optional, intent(in) :: timelevel !< The time level in the file to read - integer, optional, intent(in) :: stagger !< A flag indicating where this vector is discretized - logical, optional, intent(in) :: scalar_pair !< If true, a pair of scalars are to be read.cretized - real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied - !! by before they are returned. - - integer :: is, ie, js, je - integer :: u_pos, v_pos - - u_pos = EAST_FACE ; v_pos = NORTH_FACE - if (present(stagger)) then - if (stagger == CGRID_NE) then ; u_pos = EAST_FACE ; v_pos = NORTH_FACE - elseif (stagger == BGRID_NE) then ; u_pos = CORNER ; v_pos = CORNER - elseif (stagger == AGRID) then ; u_pos = CENTER ; v_pos = CENTER ; endif - endif - - call read_data(filename, u_fieldname, u_data, MOM_Domain%mpp_domain, & - timelevel=timelevel, position=u_pos) - call read_data(filename, v_fieldname, v_data, MOM_Domain%mpp_domain, & - timelevel=timelevel, position=v_pos) - - if (present(scale)) then ; if (scale /= 1.0) then - call get_simple_array_i_ind(MOM_Domain, size(u_data,1), is, ie) - call get_simple_array_j_ind(MOM_Domain, size(u_data,2), js, je) - u_data(is:ie,js:je,:) = scale*u_data(is:ie,js:je,:) - call get_simple_array_i_ind(MOM_Domain, size(v_data,1), is, ie) - call get_simple_array_j_ind(MOM_Domain, size(v_data,2), js, je) - v_data(is:ie,js:je,:) = scale*v_data(is:ie,js:je,:) - endif ; endif - -end subroutine MOM_read_vector_3d - - !> Initialize the MOM_io module subroutine MOM_io_init(param_file) type(param_file_type), intent(in) :: param_file !< structure indicating the open file to !! parse for model parameter values. -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "MOM_io" ! This module's name. call log_version(param_file, mdl, version) end subroutine MOM_io_init - !> \namespace mom_io !! !! This file contains a number of subroutines that manipulate diff --git a/src/framework/MOM_io_wrapper.F90 b/src/framework/MOM_io_wrapper.F90 new file mode 100644 index 0000000000..7437b59db1 --- /dev/null +++ b/src/framework/MOM_io_wrapper.F90 @@ -0,0 +1,532 @@ +!> This module contains a thin inteface to mpp and fms I/O code +module MOM_io_wrapper + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_array_transform, only : allocate_rotated_array, rotate_array +use MOM_domains, only : MOM_domain_type, AGRID, BGRID_NE, CGRID_NE +use MOM_domains, only : get_simple_array_i_ind, get_simple_array_j_ind +use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING + +use ensemble_manager_mod, only : get_ensemble_id +use fms_mod, only : write_version_number, open_namelist_file, check_nml_error +use fms_io_mod, only : file_exist, field_exist, field_size, read_data +use fms_io_mod, only : io_infra_end=>fms_io_exit, get_filename_appendix +use mpp_domains_mod, only : domain2d, CENTER, CORNER, NORTH_FACE=>NORTH, EAST_FACE=>EAST +use mpp_io_mod, only : mpp_open, close_file=>mpp_close +use mpp_io_mod, only : write_metadata=>mpp_write_meta, write_field=>mpp_write +use mpp_io_mod, only : get_field_atts=>mpp_get_atts, mpp_attribute_exist +use mpp_io_mod, only : mpp_get_axes, axistype, get_axis_data=>mpp_get_axis_data +use mpp_io_mod, only : mpp_get_fields, fieldtype, flush_file=>mpp_flush +use mpp_io_mod, only : APPEND_FILE=>MPP_APPEND, ASCII_FILE=>MPP_ASCII +use mpp_io_mod, only : MULTIPLE=>MPP_MULTI, NETCDF_FILE=>MPP_NETCDF +use mpp_io_mod, only : OVERWRITE_FILE=>MPP_OVERWR, READONLY_FILE=>MPP_RDONLY +use mpp_io_mod, only : SINGLE_FILE=>MPP_SINGLE, WRITEONLY_FILE=>MPP_WRONLY +use mpp_io_mod, only : get_file_info=>mpp_get_info, get_file_atts=>mpp_get_atts +use mpp_io_mod, only : get_file_fields=>mpp_get_fields, get_file_times=>mpp_get_times +use mpp_io_mod, only : io_infra_init=>mpp_io_init + +implicit none ; private + +! These interfaces are actually implemented in this file. +public :: MOM_read_data, MOM_read_vector, MOM_write_field, read_axis_data +public :: file_exists, field_exists, read_field_chksum +! The following are simple pass throughs of routines from other modules. +public :: open_file, close_file, field_size, fieldtype, get_filename_appendix +public :: flush_file, get_file_info, get_file_atts, get_file_fields, get_field_atts +public :: get_file_times, read_data, axistype, get_axis_data +public :: write_field, write_metadata, write_version_number, get_ensemble_id +public :: open_namelist_file, check_nml_error, io_infra_init, io_infra_end +! These are encoding constants. +public :: APPEND_FILE, ASCII_FILE, MULTIPLE, NETCDF_FILE, OVERWRITE_FILE +public :: READONLY_FILE, SINGLE_FILE, WRITEONLY_FILE +public :: CENTER, CORNER, NORTH_FACE, EAST_FACE + +!> Indicate whether a file exists, perhaps with domain decomposition +interface file_exists + module procedure FMS_file_exists + module procedure MOM_file_exists +end interface + +!> Read a data field from a file +interface MOM_read_data + module procedure MOM_read_data_4d + module procedure MOM_read_data_3d + module procedure MOM_read_data_2d + module procedure MOM_read_data_1d + module procedure MOM_read_data_0d +end interface + +!> Write a registered field to an output file +interface MOM_write_field + module procedure MOM_write_field_4d + module procedure MOM_write_field_3d + module procedure MOM_write_field_2d + module procedure MOM_write_field_1d + module procedure MOM_write_field_0d +end interface MOM_write_field + +!> Read a pair of data fields representing the two components of a vector from a file +interface MOM_read_vector + module procedure MOM_read_vector_3d + module procedure MOM_read_vector_2d +end interface + +contains + +!> Read the data associated with a named axis in a file +subroutine read_axis_data(filename, axis_name, var) + character(len=*), intent(in) :: filename !< Name of the file to read + character(len=*), intent(in) :: axis_name !< Name of the axis to read + real, dimension(:), intent(out) :: var !< The axis location data + + integer :: i, len, unit, ndim, nvar, natt, ntime + logical :: axis_found + type(axistype), allocatable :: axes(:) + type(axistype) :: time_axis + character(len=32) :: name, units + + call open_file(unit, trim(filename), action=READONLY_FILE, form=NETCDF_FILE, & + threading=MULTIPLE, fileset=SINGLE_FILE) + +!Find the number of variables (nvar) in this file + call get_file_info(unit, ndim, nvar, natt, ntime) +! ------------------------------------------------------------------- +! Allocate space for the number of axes in the data file. +! ------------------------------------------------------------------- + allocate(axes(ndim)) + call mpp_get_axes(unit, axes, time_axis) + + axis_found = .false. + do i = 1, ndim + call get_file_atts(axes(i), name=name, len=len, units=units) + if (name == axis_name) then + axis_found = .true. + call get_axis_data(axes(i), var) + exit + endif + enddo + + if (.not.axis_found) call MOM_error(FATAL, "MOM_io read_axis_data: "//& + "Unable to find axis "//trim(axis_name)//" in file "//trim(filename)) + + deallocate(axes) + +end subroutine read_axis_data + +subroutine read_field_chksum(field, chksum, valid_chksum) + type(fieldtype), intent(in) :: field !< The field whose checksum attribute is to be read. + integer(kind=8), intent(out) :: chksum !< The checksum for the field. + logical, intent(out) :: valid_chksum !< If true, chksum has been successfully read. + ! Local variables + integer(kind=8), dimension(3) :: checksum_file + + checksum_file(:) = -1 + valid_chksum = mpp_attribute_exist(field, "checksum") + if (valid_chksum) then + call get_field_atts(field, checksum=checksum_file) + chksum = checksum_file(1) + else + chksum = -1 + endif +end subroutine read_field_chksum + + +!> Returns true if the named file or its domain-decomposed variant exists. +function MOM_file_exists(filename, MOM_Domain) + character(len=*), intent(in) :: filename !< The name of the file being inquired about + type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + +! This function uses the fms_io function file_exist to determine whether +! a named file (or its decomposed variant) exists. + + logical :: MOM_file_exists + + MOM_file_exists = file_exist(filename, MOM_Domain%mpp_domain) + +end function MOM_file_exists + +!> Returns true if the named file or its domain-decomposed variant exists. +function FMS_file_exists(filename, domain, no_domain) + character(len=*), intent(in) :: filename !< The name of the file being inquired about + type(domain2d), optional, intent(in) :: domain !< The mpp domain2d that describes the decomposition + logical, optional, intent(in) :: no_domain !< This file does not use domain decomposition +! This function uses the fms_io function file_exist to determine whether +! a named file (or its decomposed variant) exists. + + logical :: FMS_file_exists + + FMS_file_exists = file_exist(filename, domain, no_domain) + +end function FMS_file_exists + +!> open_file opens a file for parallel or single-file I/O. +subroutine open_file(unit, file, action, form, threading, fileset, nohdrs, domain, MOM_domain) + integer, intent(out) :: unit !< The I/O unit for the opened file + character(len=*), intent(in) :: file !< The name of the file being opened + integer, optional, intent(in) :: action !< A flag indicating whether the file can be read + !! or written to and how to handle existing files. + integer, optional, intent(in) :: form !< A flag indicating the format of a new file. The + !! default is ASCII_FILE, but NETCDF_FILE is also common. + integer, optional, intent(in) :: threading !< A flag indicating whether one (SINGLE_FILE) + !! or multiple PEs (MULTIPLE) participate in I/O. + !! With the default, the root PE does I/O. + integer, optional, intent(in) :: fileset !< A flag indicating whether multiple PEs doing I/O due + !! to threading=MULTIPLE write to the same file (SINGLE_FILE) + !! or to one file per PE (MULTIPLE, the default). + logical, optional, intent(in) :: nohdrs !< If nohdrs is .TRUE., headers are not written to + !! ASCII files. The default is .false. + type(domain2d), optional, intent(in) :: domain !< A domain2d type that describes the decomposition + type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition + + if (present(MOM_Domain)) then + call mpp_open(unit, file, action=action, form=form, threading=threading, fileset=fileset, & + nohdrs=nohdrs, domain=MOM_Domain%mpp_domain) + else + call mpp_open(unit, file, action=action, form=form, threading=threading, fileset=fileset, & + nohdrs=nohdrs, domain=domain) + endif +end subroutine open_file + +!> Field_exists returns true if the field indicated by field_name is present in the +!! file file_name. If file_name does not exist, it returns false. +function field_exists(filename, field_name, domain, no_domain, MOM_domain) + character(len=*), intent(in) :: filename !< The name of the file being inquired about + character(len=*), intent(in) :: field_name !< The name of the field being sought + type(domain2d), target, optional, intent(in) :: domain !< A domain2d type that describes the decomposition + logical, optional, intent(in) :: no_domain !< This file does not use domain decomposition + type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition + logical :: field_exists !< True if filename exists and field_name is in filename + + if (present(MOM_domain)) then + field_exists = field_exist(filename, field_name, domain=MOM_domain%mpp_domain, no_domain=no_domain) + else + field_exists = field_exist(filename, field_name, domain=domain, no_domain=no_domain) + endif + +end function field_exists + +!> This function uses the fms_io function read_data to read a scalar +!! data field named "fieldname" from file "filename". +subroutine MOM_read_data_0d(filename, fieldname, data, timelevel, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: fieldname !< The variable name of the data in the file + real, intent(inout) :: data !< The 1-dimensional array into which the data + integer, optional, intent(in) :: timelevel !< The time level in the file to read + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied + !! by before it is returned. + + call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) + + if (present(scale)) then ; if (scale /= 1.0) then + data = scale*data + endif ; endif + +end subroutine MOM_read_data_0d + +!> This function uses the fms_io function read_data to read a 1-D +!! data field named "fieldname" from file "filename". +subroutine MOM_read_data_1d(filename, fieldname, data, timelevel, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: fieldname !< The variable name of the data in the file + real, dimension(:), intent(inout) :: data !< The 1-dimensional array into which the data + integer, optional, intent(in) :: timelevel !< The time level in the file to read + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied + !! by before they are returned. + + call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.) + + if (present(scale)) then ; if (scale /= 1.0) then + data(:) = scale*data(:) + endif ; endif + +end subroutine MOM_read_data_1d + +!> This function uses the fms_io function read_data to read a distributed +!! 2-D data field named "fieldname" from file "filename". Valid values for +!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. +subroutine MOM_read_data_2d(filename, fieldname, data, MOM_Domain, & + timelevel, position, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: fieldname !< The variable name of the data in the file + real, dimension(:,:), intent(inout) :: data !< The 2-dimensional array into which the data + !! should be read + type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + integer, optional, intent(in) :: timelevel !< The time level in the file to read + integer, optional, intent(in) :: position !< A flag indicating where this data is located + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied + !! by before it is returned. + + integer :: is, ie, js, je + + call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & + timelevel=timelevel, position=position) + + if (present(scale)) then ; if (scale /= 1.0) then + call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) + data(is:ie,js:je) = scale*data(is:ie,js:je) + endif ; endif + +end subroutine MOM_read_data_2d + +!> This function uses the fms_io function read_data to read a distributed +!! 3-D data field named "fieldname" from file "filename". Valid values for +!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. +subroutine MOM_read_data_3d(filename, fieldname, data, MOM_Domain, & + timelevel, position, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: fieldname !< The variable name of the data in the file + real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional array into which the data + !! should be read + type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + integer, optional, intent(in) :: timelevel !< The time level in the file to read + integer, optional, intent(in) :: position !< A flag indicating where this data is located + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied + !! by before it is returned. + + integer :: is, ie, js, je + + call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & + timelevel=timelevel, position=position) + + if (present(scale)) then ; if (scale /= 1.0) then + call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) + data(is:ie,js:je,:) = scale*data(is:ie,js:je,:) + endif ; endif + +end subroutine MOM_read_data_3d + +!> This function uses the fms_io function read_data to read a distributed +!! 4-D data field named "fieldname" from file "filename". Valid values for +!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE. +subroutine MOM_read_data_4d(filename, fieldname, data, MOM_Domain, & + timelevel, position, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: fieldname !< The variable name of the data in the file + real, dimension(:,:,:,:), intent(inout) :: data !< The 4-dimensional array into which the data + !! should be read + type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + integer, optional, intent(in) :: timelevel !< The time level in the file to read + integer, optional, intent(in) :: position !< A flag indicating where this data is located + real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied + !! by before it is returned. + + integer :: is, ie, js, je + + call read_data(filename, fieldname, data, MOM_Domain%mpp_domain, & + timelevel=timelevel, position=position) + + if (present(scale)) then ; if (scale /= 1.0) then + call get_simple_array_i_ind(MOM_Domain, size(data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(data,2), js, je) + data(is:ie,js:je,:,:) = scale*data(is:ie,js:je,:,:) + endif ; endif + +end subroutine MOM_read_data_4d + + +!> This function uses the fms_io function read_data to read a pair of distributed +!! 2-D data fields with names given by "[uv]_fieldname" from file "filename". Valid values for +!! "stagger" include CGRID_NE, BGRID_NE, and AGRID. +subroutine MOM_read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, & + timelevel, stagger, scalar_pair, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: u_fieldname !< The variable name of the u data in the file + character(len=*), intent(in) :: v_fieldname !< The variable name of the v data in the file + real, dimension(:,:), intent(inout) :: u_data !< The 2 dimensional array into which the + !! u-component of the data should be read + real, dimension(:,:), intent(inout) :: v_data !< The 2 dimensional array into which the + !! v-component of the data should be read + type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + integer, optional, intent(in) :: timelevel !< The time level in the file to read + integer, optional, intent(in) :: stagger !< A flag indicating where this vector is discretized + logical, optional, intent(in) :: scalar_pair !< If true, a pair of scalars are to be read + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied + !! by before they are returned. + integer :: is, ie, js, je + integer :: u_pos, v_pos + + u_pos = EAST_FACE ; v_pos = NORTH_FACE + if (present(stagger)) then + if (stagger == CGRID_NE) then ; u_pos = EAST_FACE ; v_pos = NORTH_FACE + elseif (stagger == BGRID_NE) then ; u_pos = CORNER ; v_pos = CORNER + elseif (stagger == AGRID) then ; u_pos = CENTER ; v_pos = CENTER ; endif + endif + + call read_data(filename, u_fieldname, u_data, MOM_Domain%mpp_domain, & + timelevel=timelevel, position=u_pos) + call read_data(filename, v_fieldname, v_data, MOM_Domain%mpp_domain, & + timelevel=timelevel, position=v_pos) + + if (present(scale)) then ; if (scale /= 1.0) then + call get_simple_array_i_ind(MOM_Domain, size(u_data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(u_data,2), js, je) + u_data(is:ie,js:je) = scale*u_data(is:ie,js:je) + call get_simple_array_i_ind(MOM_Domain, size(v_data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(v_data,2), js, je) + v_data(is:ie,js:je) = scale*v_data(is:ie,js:je) + endif ; endif + +end subroutine MOM_read_vector_2d + +!> This function uses the fms_io function read_data to read a pair of distributed +!! 3-D data fields with names given by "[uv]_fieldname" from file "filename". Valid values for +!! "stagger" include CGRID_NE, BGRID_NE, and AGRID. +subroutine MOM_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, & + timelevel, stagger, scalar_pair, scale) + character(len=*), intent(in) :: filename !< The name of the file to read + character(len=*), intent(in) :: u_fieldname !< The variable name of the u data in the file + character(len=*), intent(in) :: v_fieldname !< The variable name of the v data in the file + real, dimension(:,:,:), intent(inout) :: u_data !< The 3 dimensional array into which the + !! u-component of the data should be read + real, dimension(:,:,:), intent(inout) :: v_data !< The 3 dimensional array into which the + !! v-component of the data should be read + type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition + integer, optional, intent(in) :: timelevel !< The time level in the file to read + integer, optional, intent(in) :: stagger !< A flag indicating where this vector is discretized + logical, optional, intent(in) :: scalar_pair !< If true, a pair of scalars are to be read.cretized + real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied + !! by before they are returned. + + integer :: is, ie, js, je + integer :: u_pos, v_pos + + u_pos = EAST_FACE ; v_pos = NORTH_FACE + if (present(stagger)) then + if (stagger == CGRID_NE) then ; u_pos = EAST_FACE ; v_pos = NORTH_FACE + elseif (stagger == BGRID_NE) then ; u_pos = CORNER ; v_pos = CORNER + elseif (stagger == AGRID) then ; u_pos = CENTER ; v_pos = CENTER ; endif + endif + + call read_data(filename, u_fieldname, u_data, MOM_Domain%mpp_domain, & + timelevel=timelevel, position=u_pos) + call read_data(filename, v_fieldname, v_data, MOM_Domain%mpp_domain, & + timelevel=timelevel, position=v_pos) + + if (present(scale)) then ; if (scale /= 1.0) then + call get_simple_array_i_ind(MOM_Domain, size(u_data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(u_data,2), js, je) + u_data(is:ie,js:je,:) = scale*u_data(is:ie,js:je,:) + call get_simple_array_i_ind(MOM_Domain, size(v_data,1), is, ie) + call get_simple_array_j_ind(MOM_Domain, size(v_data,2), js, je) + v_data(is:ie,js:je,:) = scale*v_data(is:ie,js:je,:) + endif ; endif + +end subroutine MOM_read_vector_3d + + +!> Write a 4d field to an output file, potentially with rotation +subroutine MOM_write_field_4d(io_unit, field_md, MOM_domain, field, tstamp, tile_count, & + fill_value, turns) + integer, intent(in) :: io_unit !< File I/O unit handle + type(fieldtype), intent(in) :: field_md !< Field type with metadata + type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition + real, dimension(:,:,:,:), intent(inout) :: field !< Unrotated field to write + real, optional, intent(in) :: tstamp !< Model timestamp + integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) + real, optional, intent(in) :: fill_value !< Missing data fill value + integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data + + real, allocatable :: field_rot(:,:,:,:) ! A rotated version of field, with the same units + integer :: qturns ! The number of quarter turns through which to rotate field + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call write_field(io_unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & + tile_count=tile_count, default_data=fill_value) + else + call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + call write_field(io_unit, field_md, MOM_domain%mpp_domain, field_rot, tstamp=tstamp, & + tile_count=tile_count, default_data=fill_value) + deallocate(field_rot) + endif +end subroutine MOM_write_field_4d + +!> Write a 3d field to an output file, potentially with rotation +subroutine MOM_write_field_3d(io_unit, field_md, MOM_domain, field, tstamp, tile_count, & + fill_value, turns) + integer, intent(in) :: io_unit !< File I/O unit handle + type(fieldtype), intent(in) :: field_md !< Field type with metadata + type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition + real, dimension(:,:,:), intent(inout) :: field !< Unrotated field to write + real, optional, intent(in) :: tstamp !< Model timestamp + integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) + real, optional, intent(in) :: fill_value !< Missing data fill value + integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data + + real, allocatable :: field_rot(:,:,:) ! A rotated version of field, with the same units + integer :: qturns ! The number of quarter turns through which to rotate field + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call write_field(io_unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & + tile_count=tile_count, default_data=fill_value) + else + call allocate_rotated_array(field, [1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + call write_field(io_unit, field_md, MOM_domain%mpp_domain, field_rot, tstamp=tstamp, & + tile_count=tile_count, default_data=fill_value) + deallocate(field_rot) + endif +end subroutine MOM_write_field_3d + +!> Write a 2d field to an output file, potentially with rotation +subroutine MOM_write_field_2d(io_unit, field_md, MOM_domain, field, tstamp, tile_count, & + fill_value, turns) + integer, intent(in) :: io_unit !< File I/O unit handle + type(fieldtype), intent(in) :: field_md !< Field type with metadata + type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition + real, dimension(:,:), intent(inout) :: field !< Unrotated field to write + real, optional, intent(in) :: tstamp !< Model timestamp + integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) + real, optional, intent(in) :: fill_value !< Missing data fill value + integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data + + real, allocatable :: field_rot(:,:) ! A rotated version of field, with the same units + integer :: qturns ! The number of quarter turns through which to rotate field + + qturns = 0 + if (present(turns)) qturns = modulo(turns, 4) + + if (qturns == 0) then + call write_field(io_unit, field_md, MOM_domain%mpp_domain, field, tstamp=tstamp, & + tile_count=tile_count, default_data=fill_value) + else + call allocate_rotated_array(field, [1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + call write_field(io_unit, field_md, MOM_domain%mpp_domain, field_rot, tstamp=tstamp, & + tile_count=tile_count, default_data=fill_value) + deallocate(field_rot) + endif +end subroutine MOM_write_field_2d + +!> Write a 1d field to an output file +subroutine MOM_write_field_1d(io_unit, field_md, field, tstamp, fill_value) + integer, intent(in) :: io_unit !< File I/O unit handle + type(fieldtype), intent(in) :: field_md !< Field type with metadata + real, dimension(:), intent(inout) :: field !< Field to write + real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: fill_value !< Missing data fill value + + call write_field(io_unit, field_md, field, tstamp=tstamp) +end subroutine MOM_write_field_1d + +!> Write a 0d field to an output file +subroutine MOM_write_field_0d(io_unit, field_md, field, tstamp, fill_value) + integer, intent(in) :: io_unit !< File I/O unit handle + type(fieldtype), intent(in) :: field_md !< Field type with metadata + real, intent(inout) :: field !< Field to write + real, optional, intent(in) :: tstamp !< Model timestamp + real, optional, intent(in) :: fill_value !< Missing data fill value + + call write_field(io_unit, field_md, field, tstamp=tstamp) +end subroutine MOM_write_field_0d + +end module MOM_io_wrapper diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 619aff6f18..73a41c5aa5 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -3,21 +3,20 @@ module MOM_restart ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_checksums, only : chksum => rotated_field_chksum use MOM_domains, only : PE_here, num_PEs use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_string_functions, only : lowercase use MOM_grid, only : ocean_grid_type use MOM_io, only : create_file, fieldtype, file_exists, open_file, close_file -use MOM_io, only : MOM_read_data, read_data, get_filename_appendix, read_field_chksum +use MOM_io, only : MOM_read_data, read_data, MOM_write_field, read_field_chksum use MOM_io, only : get_file_info, get_file_atts, get_file_fields, get_file_times -use MOM_io, only : vardesc, var_desc, query_vardesc, modify_vardesc +use MOM_io, only : vardesc, var_desc, query_vardesc, modify_vardesc, get_filename_appendix use MOM_io, only : MULTIPLE, NETCDF_FILE, READONLY_FILE, SINGLE_FILE use MOM_io, only : CENTER, CORNER, NORTH_FACE, EAST_FACE +use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type, time_type_to_real, real_to_time use MOM_time_manager, only : days_in_month, get_date, set_date -use MOM_transform_FMS, only : chksum => rotated_mpp_chksum -use MOM_transform_FMS, only : write_field => rotated_write_field use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -370,7 +369,7 @@ subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units type(MOM_restart_CS), pointer :: CS !< A pointer to a MOM_restart_CS object (intent in/out) character(len=*), optional, intent(in) :: longname !< variable long name character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: hor_grid !< variable horizonal staggering, 'h' if absent + character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, 'h' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent @@ -397,7 +396,7 @@ subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units type(MOM_restart_CS), pointer :: CS !< A pointer to a MOM_restart_CS object (intent in/out) character(len=*), optional, intent(in) :: longname !< variable long name character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: hor_grid !< variable horizonal staggering, 'h' if absent + character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, 'h' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent @@ -424,7 +423,7 @@ subroutine register_restart_field_2d(f_ptr, name, mandatory, CS, longname, units type(MOM_restart_CS), pointer :: CS !< A pointer to a MOM_restart_CS object (intent in/out) character(len=*), optional, intent(in) :: longname !< variable long name character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: hor_grid !< variable horizonal staggering, 'h' if absent + character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, 'h' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, '1' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent @@ -452,7 +451,7 @@ subroutine register_restart_field_1d(f_ptr, name, mandatory, CS, longname, units type(MOM_restart_CS), pointer :: CS !< A pointer to a MOM_restart_CS object (intent in/out) character(len=*), optional, intent(in) :: longname !< variable long name character(len=*), optional, intent(in) :: units !< variable units - character(len=*), optional, intent(in) :: hor_grid !< variable horizonal staggering, '1' if absent + character(len=*), optional, intent(in) :: hor_grid !< variable horizontal staggering, '1' if absent character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering, 'L' if absent character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1, 's' if absent @@ -503,7 +502,7 @@ function query_initialized_name(name, CS) result(query_initialized) ! This subroutine returns .true. if the field referred to by name has ! initialized from a restart file, and .false. otherwise. - integer :: m,n + integer :: m, n if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "query_initialized: Module must be initialized before it is used.") if (CS%novars > CS%max_fields) call restart_error(CS) @@ -537,7 +536,7 @@ function query_initialized_0d(f_ptr, CS) result(query_initialized) ! This subroutine tests whether the field pointed to by f_ptr has ! been initialized from a restart file. - integer :: m,n + integer :: m, n if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "query_initialized: Module must be initialized before it is used.") if (CS%novars > CS%max_fields) call restart_error(CS) @@ -564,7 +563,7 @@ function query_initialized_1d(f_ptr, CS) result(query_initialized) ! This subroutine tests whether the field pointed to by f_ptr has ! been initialized from a restart file. - integer :: m,n + integer :: m, n if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "query_initialized: Module must be initialized before it is used.") if (CS%novars > CS%max_fields) call restart_error(CS) @@ -592,7 +591,7 @@ function query_initialized_2d(f_ptr, CS) result(query_initialized) ! This subroutine tests whether the field pointed to by f_ptr has ! been initialized from a restart file. - integer :: m,n + integer :: m, n if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "query_initialized: Module must be initialized before it is used.") if (CS%novars > CS%max_fields) call restart_error(CS) @@ -620,7 +619,7 @@ function query_initialized_3d(f_ptr, CS) result(query_initialized) ! This subroutine tests whether the field pointed to by f_ptr has ! been initialized from a restart file. - integer :: m,n + integer :: m, n if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "query_initialized: Module must be initialized before it is used.") if (CS%novars > CS%max_fields) call restart_error(CS) @@ -648,7 +647,7 @@ function query_initialized_4d(f_ptr, CS) result(query_initialized) ! This subroutine tests whether the field pointed to by f_ptr has ! been initialized from a restart file. - integer :: m,n + integer :: m, n if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "query_initialized: Module must be initialized before it is used.") if (CS%novars > CS%max_fields) call restart_error(CS) @@ -677,7 +676,7 @@ function query_initialized_0d_name(f_ptr, name, CS) result(query_initialized) ! This subroutine tests whether the field pointed to by f_ptr or with the ! specified variable name has been initialized from a restart file. - integer :: m,n + integer :: m, n if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "query_initialized: Module must be initialized before it is used.") if (CS%novars > CS%max_fields) call restart_error(CS) @@ -713,7 +712,7 @@ function query_initialized_1d_name(f_ptr, name, CS) result(query_initialized) ! This subroutine tests whether the field pointed to by f_ptr or with the ! specified variable name has been initialized from a restart file. - integer :: m,n + integer :: m, n if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "query_initialized: Module must be initialized before it is used.") if (CS%novars > CS%max_fields) call restart_error(CS) @@ -749,7 +748,7 @@ function query_initialized_2d_name(f_ptr, name, CS) result(query_initialized) ! This subroutine tests whether the field pointed to by f_ptr or with the ! specified variable name has been initialized from a restart file. - integer :: m,n + integer :: m, n if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // & "query_initialized: Module must be initialized before it is used.") if (CS%novars > CS%max_fields) call restart_error(CS) @@ -875,7 +874,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ ! this should be 2 Gb or less. integer :: start_var, next_var ! The starting variables of the ! current and next files. - integer :: unit ! The mpp unit of the open file. + integer :: unit ! The I/O unit of the open file. integer :: m, nz, num_files, var_periods integer :: seconds, days, year, month, hour, minute character(len=8) :: hor_grid, z_grid, t_grid ! Variable grid info. @@ -908,10 +907,10 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ restartname = trim(CS%restartfile) if (present(filename)) restartname = trim(filename) if (PRESENT(time_stamped)) then ; if (time_stamped) then - call get_date(time,year,month,days,hour,minute,seconds) + call get_date(time, year, month, days, hour, minute, seconds) ! Compute the year-day, because I don't like months. - RWH do m=1,month-1 - days = days + days_in_month(set_date(year,m,2,0,0,0)) + days = days + days_in_month(set_date(year, m, 2, 0, 0, 0)) enddo seconds = seconds + 60*minute + 3600*hour if (year <= 9999) then @@ -1030,19 +1029,19 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ do m=start_var,next_var-1 if (associated(CS%var_ptr3d(m)%p)) then - call write_field(unit,fields(m-start_var+1), G%Domain%mpp_domain, & + call MOM_write_field(unit,fields(m-start_var+1), G%Domain, & CS%var_ptr3d(m)%p, restart_time, turns=-turns) elseif (associated(CS%var_ptr2d(m)%p)) then - call write_field(unit,fields(m-start_var+1), G%Domain%mpp_domain, & + call MOM_write_field(unit,fields(m-start_var+1), G%Domain, & CS%var_ptr2d(m)%p, restart_time, turns=-turns) elseif (associated(CS%var_ptr4d(m)%p)) then - call write_field(unit,fields(m-start_var+1), G%Domain%mpp_domain, & + call MOM_write_field(unit,fields(m-start_var+1), G%Domain, & CS%var_ptr4d(m)%p, restart_time, turns=-turns) elseif (associated(CS%var_ptr1d(m)%p)) then - call write_field(unit, fields(m-start_var+1), CS%var_ptr1d(m)%p, & + call MOM_write_field(unit, fields(m-start_var+1), CS%var_ptr1d(m)%p, & restart_time) elseif (associated(CS%var_ptr0d(m)%p)) then - call write_field(unit, fields(m-start_var+1), CS%var_ptr0d(m)%p, & + call MOM_write_field(unit, fields(m-start_var+1), CS%var_ptr0d(m)%p, & restart_time) endif enddo @@ -1087,7 +1086,7 @@ subroutine restore_state(filename, directory, day, G, CS) integer :: sizes(7) integer :: ndim, nvar, natt, ntime, pos - integer :: unit(CS%max_fields) ! The mpp unit of all open files. + integer :: unit(CS%max_fields) ! The I/O units of all open files. character(len=200) :: unit_path(CS%max_fields) ! The file names. logical :: unit_is_global(CS%max_fields) ! True if the file is global. @@ -1364,7 +1363,7 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & type(MOM_restart_CS), pointer :: CS !< The control structure returned by a previous !! call to restart_init. integer, dimension(:), & - optional, intent(out) :: units !< The mpp units of all opened files. + optional, intent(out) :: units !< The I/O units of all opened files. character(len=*), dimension(:), & optional, intent(out) :: file_paths !< The full paths to open files. logical, dimension(:), & @@ -1452,8 +1451,7 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & ! Look for decomposed files using the I/O Layout. fexists = file_exists(filepath, G%Domain) if (fexists .and. (present(units))) & - call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, & - domain=G%Domain%mpp_domain) + call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, MOM_domain=G%Domain) if (fexists .and. present(global_files)) global_files(n) = .false. endif diff --git a/src/framework/MOM_transform_FMS.F90 b/src/framework/MOM_transform_FMS.F90 index 572a9717dc..a4a3f7c2c4 100644 --- a/src/framework/MOM_transform_FMS.F90 +++ b/src/framework/MOM_transform_FMS.F90 @@ -3,41 +3,16 @@ module MOM_transform_FMS -use horiz_interp_mod, only : horiz_interp_type +use MOM_array_transform, only : allocate_rotated_array, rotate_array use MOM_error_handler, only : MOM_error, FATAL -use MOM_io, only : fieldtype, write_field -use mpp_domains_mod, only : domain2D -use mpp_mod, only : mpp_chksum +use horiz_interp_mod, only : horiz_interp_type use time_manager_mod, only : time_type use time_interp_external_mod, only : time_interp_external -use MOM_array_transform, only : allocate_rotated_array, rotate_array +implicit none ; private -implicit none - -private -public rotated_mpp_chksum -public rotated_write_field public rotated_time_interp_external -!> Rotate and compute the FMS (mpp) checksum of a field -interface rotated_mpp_chksum - module procedure rotated_mpp_chksum_real_0d - module procedure rotated_mpp_chksum_real_1d - module procedure rotated_mpp_chksum_real_2d - module procedure rotated_mpp_chksum_real_3d - module procedure rotated_mpp_chksum_real_4d -end interface rotated_mpp_chksum - -!> Rotate and write a registered field to an FMS output file -interface rotated_write_field - module procedure rotated_write_field_real_0d - module procedure rotated_write_field_real_1d - module procedure rotated_write_field_real_2d - module procedure rotated_write_field_real_3d - module procedure rotated_write_field_real_4d -end interface rotated_write_field - !> Read a field based on model time, and rotate to the model domain interface rotated_time_interp_external module procedure rotated_time_interp_external_0d @@ -50,255 +25,6 @@ module MOM_transform_FMS ! NOTE: No transformations are applied to the 0d and 1d field implementations, ! but are provided to maintain compatibility with the FMS interfaces. - -!> Compute the FMS (mpp) checksum of a scalar. -!! This function is provided to support the full FMS mpp_chksum interface. -function rotated_mpp_chksum_real_0d(field, pelist, mask_val, turns) & - result(chksum) - real, intent(in) :: field !> Input scalar - integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum - real, optional, intent(in) :: mask_val !> FMS mask value - integer, optional, intent(in) :: turns !> Number of quarter turns - integer :: chksum !> FMS checksum of scalar - - if (present(turns)) & - call MOM_error(FATAL, "Rotation not supported for 0d fields.") - - chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) -end function rotated_mpp_chksum_real_0d - - -!> Compute the FMS (mpp) checksum of a 1d field. -!! This function is provided to support the full FMS mpp_chksum interface. -function rotated_mpp_chksum_real_1d(field, pelist, mask_val, turns) & - result(chksum) - real, intent(in) :: field(:) !> Input field - integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum - real, optional, intent(in) :: mask_val !> FMS mask value - integer, optional, intent(in) :: turns !> Number of quarter-turns - integer :: chksum !> FMS checksum of field - - if (present(turns)) & - call MOM_error(FATAL, "Rotation not supported for 1d fields.") - - chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) -end function rotated_mpp_chksum_real_1d - - -!> Compute the FMS (mpp) checksum of a rotated 2d field. -function rotated_mpp_chksum_real_2d(field, pelist, mask_val, turns) & - result(chksum) - real, intent(in) :: field(:,:) !> Unrotated input field - integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum - real, optional, intent(in) :: mask_val !> FMS mask value - integer, optional, intent(in) :: turns !> Number of quarter-turns - integer :: chksum !> FMS checksum of field - - real, allocatable :: field_rot(:,:) - integer :: qturns - - qturns = 0 - if (present(turns)) & - qturns = modulo(turns, 4) - - if (qturns == 0) then - chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) - else - call allocate_rotated_array(field, [1,1], qturns, field_rot) - call rotate_array(field, qturns, field_rot) - chksum = mpp_chksum(field_rot, pelist=pelist, mask_val=mask_val) - deallocate(field_rot) - endif -end function rotated_mpp_chksum_real_2d - - -!> Compute the FMS (mpp) checksum of a rotated 3d field. -function rotated_mpp_chksum_real_3d(field, pelist, mask_val, turns) & - result(chksum) - real, intent(in) :: field(:,:,:) !> Unrotated input field - integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum - real, optional, intent(in) :: mask_val !> FMS mask value - integer, optional, intent(in) :: turns !> Number of quarter-turns - integer :: chksum !> FMS checksum of field - - real, allocatable :: field_rot(:,:,:) - integer :: qturns - - qturns = 0 - if (present(turns)) & - qturns = modulo(turns, 4) - - if (qturns == 0) then - chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) - else - call allocate_rotated_array(field, [1,1,1], qturns, field_rot) - call rotate_array(field, qturns, field_rot) - chksum = mpp_chksum(field_rot, pelist=pelist, mask_val=mask_val) - deallocate(field_rot) - endif -end function rotated_mpp_chksum_real_3d - - -!> Compute the FMS (mpp) checksum of a rotated 4d field. -function rotated_mpp_chksum_real_4d(field, pelist, mask_val, turns) & - result(chksum) - real, intent(in) :: field(:,:,:,:) !> Unrotated input field - integer, optional, intent(in) :: pelist(:) !> PE list of ranks to checksum - real, optional, intent(in) :: mask_val !> FMS mask value - integer, optional, intent(in) :: turns !> Number of quarter-turns - integer :: chksum !> FMS checksum of field - - real, allocatable :: field_rot(:,:,:,:) - integer :: qturns - - qturns = 0 - if (present(turns)) & - qturns = modulo(turns, 4) - - if (qturns == 0) then - chksum = mpp_chksum(field, pelist=pelist, mask_val=mask_val) - else - call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot) - call rotate_array(field, qturns, field_rot) - chksum = mpp_chksum(field_rot, pelist=pelist, mask_val=mask_val) - deallocate(field_rot) - endif -end function rotated_mpp_chksum_real_4d - - -! NOTE: In MOM_io, write_field points to mpp_write, which supports a very broad -! range of interfaces. Here, we only support the much more narrow family of -! mpp_write_2ddecomp functions used to write tiled data. - - -!> Write the rotation of a 1d field to an FMS output file -!! This function is provided to support the full FMS write_field interface. -subroutine rotated_write_field_real_0d(io_unit, field_md, field, tstamp, turns) - integer, intent(in) :: io_unit !> File I/O unit handle - type(fieldtype), intent(in) :: field_md !> FMS field metadata - real, intent(inout) :: field !> Unrotated field array - real, optional, intent(in) :: tstamp !> Model timestamp - integer, optional, intent(in) :: turns !> Number of quarter-turns - - if (present(turns)) & - call MOM_error(FATAL, "Rotation not supported for 0d fields.") - - call write_field(io_unit, field_md, field, tstamp=tstamp) -end subroutine rotated_write_field_real_0d - - -!> Write the rotation of a 1d field to an FMS output file -!! This function is provided to support the full FMS write_field interface. -subroutine rotated_write_field_real_1d(io_unit, field_md, field, tstamp, turns) - integer, intent(in) :: io_unit !> File I/O unit handle - type(fieldtype), intent(in) :: field_md !> FMS field metadata - real, intent(inout) :: field(:) !> Unrotated field array - real, optional, intent(in) :: tstamp !> Model timestamp - integer, optional, intent(in) :: turns !> Number of quarter-turns - - if (present(turns)) & - call MOM_error(FATAL, "Rotation not supported for 0d fields.") - - call write_field(io_unit, field_md, field, tstamp=tstamp) -end subroutine rotated_write_field_real_1d - - -!> Write the rotation of a 2d field to an FMS output file -subroutine rotated_write_field_real_2d(io_unit, field_md, domain, field, & - tstamp, tile_count, default_data, turns) - integer, intent(in) :: io_unit !> File I/O unit handle - type(fieldtype), intent(in) :: field_md !> FMS field metadata - type(domain2D), intent(inout) :: domain !> FMS MPP domain - real, intent(inout) :: field(:,:) !> Unrotated field array - real, optional, intent(in) :: tstamp !> Model timestamp - integer, optional, intent(in) :: tile_count !> PEs per tile (default: 1) - real, optional, intent(in) :: default_data !> Default fill value - integer, optional, intent(in) :: turns !> Number of quarter-turns - - real, allocatable :: field_rot(:,:) - integer :: qturns - - qturns = 0 - if (present(turns)) & - qturns = modulo(turns, 4) - - if (qturns == 0) then - call write_field(io_unit, field_md, domain, field, tstamp=tstamp, & - tile_count=tile_count, default_data=default_data) - else - call allocate_rotated_array(field, [1,1], qturns, field_rot) - call rotate_array(field, qturns, field_rot) - call write_field(io_unit, field_md, domain, field_rot, tstamp=tstamp, & - tile_count=tile_count, default_data=default_data) - deallocate(field_rot) - endif -end subroutine rotated_write_field_real_2d - - -!> Write the rotation of a 3d field to an FMS output file -subroutine rotated_write_field_real_3d(io_unit, field_md, domain, field, & - tstamp, tile_count, default_data, turns) - integer, intent(in) :: io_unit !> File I/O unit handle - type(fieldtype), intent(in) :: field_md !> FMS field metadata - type(domain2D), intent(inout) :: domain !> FMS MPP domain - real, intent(inout) :: field(:,:,:) !> Unrotated field array - real, optional, intent(in) :: tstamp !> Model timestamp - integer, optional, intent(in) :: tile_count !> PEs per tile (default: 1) - real, optional, intent(in) :: default_data !> Default fill value - integer, optional, intent(in) :: turns !> Number of quarter-turns - - real, allocatable :: field_rot(:,:,:) - integer :: qturns - - qturns = 0 - if (present(turns)) & - qturns = modulo(turns, 4) - - if (qturns == 0) then - call write_field(io_unit, field_md, domain, field, tstamp=tstamp, & - tile_count=tile_count, default_data=default_data) - else - call allocate_rotated_array(field, [1,1,1], qturns, field_rot) - call rotate_array(field, qturns, field_rot) - call write_field(io_unit, field_md, domain, field_rot, tstamp=tstamp, & - tile_count=tile_count, default_data=default_data) - deallocate(field_rot) - endif -end subroutine rotated_write_field_real_3d - - -!> Write the rotation of a 4d field to an FMS output file -subroutine rotated_write_field_real_4d(io_unit, field_md, domain, field, & - tstamp, tile_count, default_data, turns) - integer, intent(in) :: io_unit !> File I/O unit handle - type(fieldtype), intent(in) :: field_md !> FMS field metadata - type(domain2D), intent(inout) :: domain !> FMS MPP domain - real, intent(inout) :: field(:,:,:,:) !> Unrotated field array - real, optional, intent(in) :: tstamp !> Model timestamp - integer, optional, intent(in) :: tile_count !> PEs per tile (default: 1) - real, optional, intent(in) :: default_data !> Default fill value - integer, optional, intent(in) :: turns !> Number of quarter-turns - - real, allocatable :: field_rot(:,:,:,:) - integer :: qturns - - qturns = 0 - if (present(turns)) & - qturns = modulo(turns, 4) - - if (qturns == 0) then - call write_field(io_unit, field_md, domain, field, tstamp=tstamp, & - tile_count=tile_count, default_data=default_data) - else - call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot) - call rotate_array(field, qturns, field_rot) - call write_field(io_unit, field_md, domain, field_rot, tstamp=tstamp, & - tile_count=tile_count, default_data=default_data) - deallocate(field_rot) - endif -end subroutine rotated_write_field_real_4d - - !> Read a scalar field based on model time !! This function is provided to support the full FMS time_interp_external !! interface. diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index f6c6768a85..b9b42d0cbf 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -15,8 +15,8 @@ module MOM_ice_shelf use MOM_IS_diag_mediator, only : diag_mediator_init, diag_mediator_end, set_diag_mediator_grid use MOM_IS_diag_mediator, only : enable_averages, enable_averaging, disable_averaging use MOM_IS_diag_mediator, only : diag_mediator_infrastructure_init, diag_mediator_close_registration -use MOM_domains, only : MOM_domains_init, clone_MOM_domain -use MOM_domains, only : pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE, CORNER +use MOM_domain_init, only : MOM_domains_init +use MOM_domains, only : clone_MOM_domain, pass_var, pass_vector, TO_ALL, CGRID_NE, BGRID_NE, CORNER use MOM_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid use MOM_dyn_horgrid, only : rescale_dyn_horgrid_bathymetry use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 4526d9e9c7..eee168eefb 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -3,22 +3,19 @@ module MOM_grid_initialize ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_checksums, only : hchksum, Bchksum -use MOM_checksums, only : uvchksum, hchksum_pair, Bchksum_pair -use MOM_domains, only : pass_var, pass_vector, pe_here, root_PE, broadcast -use MOM_domains, only : AGRID, BGRID_NE, CGRID_NE, To_All, Scalar_Pair -use MOM_domains, only : To_North, To_South, To_East, To_West -use MOM_domains, only : MOM_define_domain, MOM_define_IO_domain -use MOM_domains, only : MOM_domain_type -use MOM_dyn_horgrid, only : dyn_horgrid_type, set_derived_dyn_horgrid +use MOM_checksums, only : hchksum, Bchksum, uvchksum, hchksum_pair, Bchksum_pair +use MOM_domains, only : pass_var, pass_vector, pe_here, root_PE, broadcast +use MOM_domains, only : AGRID, BGRID_NE, CGRID_NE, To_All, Scalar_Pair +use MOM_domains, only : To_North, To_South, To_East, To_West +use MOM_domains, only : MOM_define_domain, MOM_define_IO_domain, get_layout_extents +use MOM_domains, only : MOM_domain_type, deallocate_domain_contents +use MOM_dyn_horgrid, only : dyn_horgrid_type, set_derived_dyn_horgrid use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_io, only : MOM_read_data, read_data, slasher, file_exists, stdout -use MOM_io, only : CORNER, NORTH_FACE, EAST_FACE -use MOM_unit_scaling, only : unit_scale_type - -use mpp_domains_mod, only : mpp_get_domain_extents, mpp_deallocate_domain +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_io, only : MOM_read_data, read_data, slasher, file_exists, stdout +use MOM_io, only : CORNER, NORTH_FACE, EAST_FACE +use MOM_unit_scaling, only : unit_scale_type implicit none ; private @@ -192,8 +189,8 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) type(MOM_domain_type) :: SGdom ! Supergrid domain logical :: lon_bug ! If true use an older buggy answer in the tripolar longitude. integer :: i, j, i2, j2 - integer :: npei,npej - integer, dimension(:), allocatable :: exni,exnj + integer, dimension(:), allocatable :: exni ! The extents of the grid for each i-row of the layout + integer, dimension(:), allocatable :: exnj ! The extents of the grid for each j-row of the layout integer :: start(4), nread(4) call callTree_enter("set_grid_metrics_from_mosaic(), MOM_grid_initialize.F90") @@ -224,9 +221,7 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) nj = 2*(G%jec-G%jsc+1) ! j size of supergrid ! Define a domain for the supergrid (SGdom) - npei = G%domain%layout(1) ; npej = G%domain%layout(2) - allocate(exni(npei)) ; allocate(exnj(npej)) - call mpp_get_domain_extents(G%domain%mpp_domain, exni, exnj) + call get_layout_extents(G%domain, exni, exnj) allocate(SGdom%mpp_domain) SGdom%nihalo = 2*G%domain%nihalo+1 SGdom%njhalo = 2*G%domain%njhalo+1 @@ -243,19 +238,18 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) call MOM_define_domain(global_indices, SGdom%layout, SGdom%mpp_domain, & xflags=G%domain%X_FLAGS, yflags=G%domain%Y_FLAGS, & xhalo=SGdom%nihalo, yhalo=SGdom%njhalo, & - xextent=exni,yextent=exnj, & + xextent=exni, yextent=exnj, & symmetry=.true., name="MOM_MOSAIC", maskmap=G%domain%maskmap) else call MOM_define_domain(global_indices, SGdom%layout, SGdom%mpp_domain, & xflags=G%domain%X_FLAGS, yflags=G%domain%Y_FLAGS, & xhalo=SGdom%nihalo, yhalo=SGdom%njhalo, & - xextent=exni,yextent=exnj, & + xextent=exni, yextent=exnj, & symmetry=.true., name="MOM_MOSAIC") endif call MOM_define_IO_domain(SGdom%mpp_domain, SGdom%io_layout) - deallocate(exni) - deallocate(exnj) + deallocate(exni, exnj) ! Read X from the supergrid tmpZ(:,:) = 999. @@ -346,8 +340,7 @@ subroutine set_grid_metrics_from_mosaic(G, param_file, US) ni=SGdom%niglobal nj=SGdom%njglobal - call mpp_deallocate_domain(SGdom%mpp_domain) - deallocate(SGdom%mpp_domain) + call deallocate_domain_contents(SGdom) call pass_vector(dyCu, dxCv, G%Domain, To_All+Scalar_Pair, CGRID_NE) call pass_vector(dxCu, dyCv, G%Domain, To_All+Scalar_Pair, CGRID_NE) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index ec51a045cf..24318954a1 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -13,7 +13,7 @@ module MOM_shared_initialization use MOM_file_parser, only : get_param, log_param, param_file_type, log_version use MOM_io, only : close_file, create_file, fieldtype, file_exists, stdout use MOM_io, only : MOM_read_data, MOM_read_vector, SINGLE_FILE, MULTIPLE -use MOM_io, only : slasher, vardesc, write_field, var_desc +use MOM_io, only : slasher, vardesc, MOM_write_field, var_desc use MOM_string_functions, only : uppercase use MOM_unit_scaling, only : unit_scale_type @@ -1312,60 +1312,60 @@ subroutine write_ocean_geometry_file(G, param_file, directory, geom_file, US) file_threading, dG=G) do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(I,J) = G%geoLatBu(I,J) ; enddo ; enddo - call write_field(unit, fields(1), G%Domain%mpp_domain, out_q) + call MOM_write_field(unit, fields(1), G%Domain, out_q) do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(I,J) = G%geoLonBu(I,J) ; enddo ; enddo - call write_field(unit, fields(2), G%Domain%mpp_domain, out_q) - call write_field(unit, fields(3), G%Domain%mpp_domain, G%geoLatT) - call write_field(unit, fields(4), G%Domain%mpp_domain, G%geoLonT) + call MOM_write_field(unit, fields(2), G%Domain, out_q) + call MOM_write_field(unit, fields(3), G%Domain, G%geoLatT) + call MOM_write_field(unit, fields(4), G%Domain, G%geoLonT) do j=js,je ; do i=is,ie ; out_h(i,j) = Z_to_m_scale*G%bathyT(i,j) ; enddo ; enddo - call write_field(unit, fields(5), G%Domain%mpp_domain, out_h) + call MOM_write_field(unit, fields(5), G%Domain, out_h) do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(i,J) = s_to_T_scale*G%CoriolisBu(I,J) ; enddo ; enddo - call write_field(unit, fields(6), G%Domain%mpp_domain, out_q) + call MOM_write_field(unit, fields(6), G%Domain, out_q) ! I think that all of these copies are holdovers from a much earlier ! ancestor code in which many of the metrics were macros that could have ! had reduced dimensions, and that they are no longer needed in MOM6. -RWH do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = L_to_m_scale*G%dxCv(i,J) ; enddo ; enddo - call write_field(unit, fields(7), G%Domain%mpp_domain, out_v) + call MOM_write_field(unit, fields(7), G%Domain, out_v) do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = L_to_m_scale*G%dyCu(I,j) ; enddo ; enddo - call write_field(unit, fields(8), G%Domain%mpp_domain, out_u) + call MOM_write_field(unit, fields(8), G%Domain, out_u) do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = L_to_m_scale*G%dxCu(I,j) ; enddo ; enddo - call write_field(unit, fields(9), G%Domain%mpp_domain, out_u) + call MOM_write_field(unit, fields(9), G%Domain, out_u) do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = L_to_m_scale*G%dyCv(i,J) ; enddo ; enddo - call write_field(unit, fields(10), G%Domain%mpp_domain, out_v) + call MOM_write_field(unit, fields(10), G%Domain, out_v) do j=js,je ; do i=is,ie ; out_h(i,j) = L_to_m_scale*G%dxT(i,j); enddo ; enddo - call write_field(unit, fields(11), G%Domain%mpp_domain, out_h) + call MOM_write_field(unit, fields(11), G%Domain, out_h) do j=js,je ; do i=is,ie ; out_h(i,j) = L_to_m_scale*G%dyT(i,j) ; enddo ; enddo - call write_field(unit, fields(12), G%Domain%mpp_domain, out_h) + call MOM_write_field(unit, fields(12), G%Domain, out_h) do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(i,J) = L_to_m_scale*G%dxBu(I,J) ; enddo ; enddo - call write_field(unit, fields(13), G%Domain%mpp_domain, out_q) + call MOM_write_field(unit, fields(13), G%Domain, out_q) do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(I,J) = L_to_m_scale*G%dyBu(I,J) ; enddo ; enddo - call write_field(unit, fields(14), G%Domain%mpp_domain, out_q) + call MOM_write_field(unit, fields(14), G%Domain, out_q) - do j=js,je ; do i=is,ie ; out_h(i,j) = G%areaT(i,j) ; enddo ; enddo - call write_field(unit, fields(15), G%Domain%mpp_domain, out_h) - do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(I,J) = G%areaBu(I,J) ; enddo ; enddo - call write_field(unit, fields(16), G%Domain%mpp_domain, out_q) + do j=js,je ; do i=is,ie ; out_h(i,j) = L_to_m_scale**2*G%areaT(i,j) ; enddo ; enddo + call MOM_write_field(unit, fields(15), G%Domain, out_h) + do J=Jsq,Jeq ; do I=Isq,Ieq ; out_q(I,J) = L_to_m_scale**2*G%areaBu(I,J) ; enddo ; enddo + call MOM_write_field(unit, fields(16), G%Domain, out_q) do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = L_to_m_scale*G%dx_Cv(i,J) ; enddo ; enddo - call write_field(unit, fields(17), G%Domain%mpp_domain, out_v) + call MOM_write_field(unit, fields(17), G%Domain, out_v) do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = L_to_m_scale*G%dy_Cu(I,j) ; enddo ; enddo - call write_field(unit, fields(18), G%Domain%mpp_domain, out_u) - call write_field(unit, fields(19), G%Domain%mpp_domain, G%mask2dT) + call MOM_write_field(unit, fields(18), G%Domain, out_u) + call MOM_write_field(unit, fields(19), G%Domain, G%mask2dT) if (G%bathymetry_at_vel) then do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = Z_to_m_scale*G%Dblock_u(I,j) ; enddo ; enddo - call write_field(unit, fields(20), G%Domain%mpp_domain, out_u) + call MOM_write_field(unit, fields(20), G%Domain, out_u) do j=js,je ; do I=Isq,Ieq ; out_u(I,j) = Z_to_m_scale*G%Dopen_u(I,j) ; enddo ; enddo - call write_field(unit, fields(21), G%Domain%mpp_domain, out_u) + call MOM_write_field(unit, fields(21), G%Domain, out_u) do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = Z_to_m_scale*G%Dblock_v(i,J) ; enddo ; enddo - call write_field(unit, fields(22), G%Domain%mpp_domain, out_v) + call MOM_write_field(unit, fields(22), G%Domain, out_v) do J=Jsq,Jeq ; do i=is,ie ; out_v(i,J) = Z_to_m_scale*G%Dopen_v(i,J) ; enddo ; enddo - call write_field(unit, fields(23), G%Domain%mpp_domain, out_v) + call MOM_write_field(unit, fields(23), G%Domain, out_v) endif call close_file(unit) diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index 6064e27726..1f95cb5162 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -536,7 +536,7 @@ subroutine find_in_files(filenames, varname, array, G) do nf=1,size(filenames) if (LEN_TRIM(filenames(nf)) == 0) cycle - if (field_exists(filenames(nf), varname, G%Domain%mpp_domain)) then + if (field_exists(filenames(nf), varname, MOM_domain=G%Domain)) then call MOM_read_data(filenames(nf), varname, array, G%Domain) return endif