From ee7dd32c8534154f0cacb5e113873eb31613b363 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 8 Jan 2021 05:18:49 -0500 Subject: [PATCH 01/11] +Add create_MOM_domain and MOM_domain_init.F90 Added the new file MOM_domain_init.F90, with a copy of MOM_domains_init, to facilitate the separation of the framework code into FMS-specific and FMS-independent directories, and use this new module in MOM.F90 and MOM_ice_shelf.F90. There is also a copy of MOM_domains_init still in MOM_domains.F90, so other modules, like SIS2, will continue to work. Some of the previous contents of MOM_domains_init have been transferred into the new publicly visible routines create_MOM_domain, MOM_thread_affinity_set and set_MOM_thread_affinity. Also removed several module use statements for MOM_domains_init that are not needed. All answers and output files are bitwise identical, but there are new public interfaces. --- src/core/MOM.F90 | 4 +- src/core/MOM_dynamics_split_RK2.F90 | 1 - src/core/MOM_dynamics_unsplit.F90 | 5 +- src/core/MOM_dynamics_unsplit_RK2.F90 | 5 +- src/framework/MOM_domains.F90 | 359 +++++++++++++++----------- src/ice_shelf/MOM_ice_shelf.F90 | 4 +- 6 files changed, 223 insertions(+), 155 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index f11ce42407..9ed4959abd 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 e7c5a71930..30bf23819a 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/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 46cc9c526a..c71ec6b848 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -12,7 +12,7 @@ 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 @@ -36,7 +36,8 @@ 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 +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 @@ -47,6 +48,7 @@ module MOM_domains 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 :: MOM_thread_affinity_set, set_MOM_thread_affinity public :: get_simple_array_i_ind, get_simple_array_j_ind public :: domain2D @@ -1169,7 +1171,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 +1185,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 +1201,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 +1214,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 +1283,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 +1313,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 +1335,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 +1452,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 +1476,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 - if (io_layout(2) == 0) io_layout(2) = layout(2) - if (io_layout(1) == 0) io_layout(1) = layout(1) +!> 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 + + 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%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%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 +1609,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 +1632,44 @@ 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 + + +!> 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 !> 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. diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 5829e49ed3..a5dd92b640 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 From 7091a09c1c0f77e86efdaf08b928d210400dc6dd Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 8 Jan 2021 10:42:56 -0500 Subject: [PATCH 02/11] +Add MOM_write_field Moved rotated_write_field from MOM_transform_FMS.F90 to MOM_io.F90 and renamed it to MOM_write_field, with the mpp_domain argument replaced with a MOM_domain argument. Also changed the calls in save_restart to reflect these changes. All answers and output files are identical. --- src/framework/MOM_io.F90 | 144 ++++++++++++++++++++++++--- src/framework/MOM_restart.F90 | 45 +++++---- src/framework/MOM_transform_FMS.F90 | 148 +--------------------------- 3 files changed, 156 insertions(+), 181 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 529c725274..eda00bfda0 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -3,13 +3,13 @@ module MOM_io ! This file is part of MOM6. See LICENSE.md for the license. - -use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING +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_dyn_horgrid, only : dyn_horgrid_type +use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING 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_string_functions, only : lowercase, slasher use MOM_verticalGrid, only : verticalGrid_type @@ -41,7 +41,7 @@ module MOM_io public :: close_file, create_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 :: num_timelevels, MOM_read_data, MOM_read_vector, MOM_write_field, ensembler public :: reopen_file, slasher, write_field, write_version_number, MOM_io_init public :: open_namelist_file, check_nml_error, io_infra_init, io_infra_end public :: APPEND_FILE, ASCII_FILE, MULTIPLE, NETCDF_FILE, OVERWRITE_FILE @@ -80,6 +80,15 @@ module MOM_io 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 @@ -621,7 +630,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 +671,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 +730,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 @@ -1002,7 +1011,7 @@ subroutine MOM_read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data 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 + 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 @@ -1078,13 +1087,126 @@ subroutine MOM_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data 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 + + !> 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) diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index d9206f5bef..73c5cc94e1 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -9,15 +9,14 @@ module MOM_restart 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_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) @@ -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 diff --git a/src/framework/MOM_transform_FMS.F90 b/src/framework/MOM_transform_FMS.F90 index 572a9717dc..a6c28717b3 100644 --- a/src/framework/MOM_transform_FMS.F90 +++ b/src/framework/MOM_transform_FMS.F90 @@ -5,19 +5,15 @@ module MOM_transform_FMS use horiz_interp_mod, only : horiz_interp_type 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 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 +implicit none ; private -private public rotated_mpp_chksum -public rotated_write_field public rotated_time_interp_external !> Rotate and compute the FMS (mpp) checksum of a field @@ -29,15 +25,6 @@ module MOM_transform_FMS 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 @@ -166,139 +153,6 @@ function rotated_mpp_chksum_real_4d(field, pelist, mask_val, turns) & 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. From 502eb301ee8012796ff8f73b79e7d814f20441a9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 9 Jan 2021 03:38:26 -0500 Subject: [PATCH 03/11] (*)Call MOM_write_field in write_ocean_geometry_file Use calls to MOM_write_field in write_ocean_geometry_file, instead of calls write_field. Also added missing dimensional rescaling factors to the two area fields Ah and Aq, which will cause the output fields in the ocean_geometry.nc files to change, but they are now invariant to the choice of dimensional rescaling, whereas previously they had not been. All answers and output are otherwise bitwise identical, and even ocean_geometry.nc is identical if L_RESCALE_POWER = 0. --- .../MOM_shared_initialization.F90 | 52 +++++++++---------- 1 file changed, 26 insertions(+), 26 deletions(-) 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) From adb8ec46bebdfb8d43cddb186520d28c6fe17313 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 9 Jan 2021 05:23:51 -0500 Subject: [PATCH 04/11] +Add deallocate_MOM_domain and get_layout_extents Added the new routines deallocate_MOM_domain, deallocate_domain_contents and get_layout_extents to standardize the clean-up of memory associated with MOM_domains, provide an interface for obtaining information about the global grid decomposition and limit the dependencies on mpp functions to calls that go through the MOM framework directory. All answers are bitwise identical, although there are new public interfaces. --- src/core/MOM_grid.F90 | 7 +-- src/framework/MOM_domains.F90 | 58 ++++++++++++++++++++-- src/initialization/MOM_grid_initialize.F90 | 43 +++++++--------- 3 files changed, 77 insertions(+), 31 deletions(-) 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/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index c71ec6b848..dc1f8ff867 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -17,6 +17,7 @@ module MOM_domains 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 : 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 @@ -37,6 +38,7 @@ module MOM_domains public :: MOM_domains_init, MOM_infra_init, MOM_infra_end, get_domain_extent, get_domain_extent_dsamp2 public :: create_MOM_domain, clone_MOM_domain +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 @@ -47,7 +49,7 @@ 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 @@ -1639,6 +1641,42 @@ subroutine create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar, lay 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() @@ -2041,13 +2079,27 @@ 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 + integer, dimension(:), allocatable, intent(inout) :: extent_j + + 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/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) From d8806f48e4e81402c524d6cc47838f9ba9181c44 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 10 Jan 2021 06:07:20 -0500 Subject: [PATCH 05/11] +Rename rotated_mpp_chksum to rotated_field_chksum Renamed rotated_mpp_chksum to rotated_field_chksum and moved the routines wrapped by this overloaded interface from MOM_transform_FMS.F90 to MOM_checksums.F90. Also provided access to mpp_chksum as field_chksum via MOM_coms.F90. Both of these are steps to clean up the MOM6 framework code and reduce the direct use of mpp routines in the rest of the MOM6 code. All answers are bitwise identical, but there are effectively new interfaces, and one existing interface was renamed. --- src/diagnostics/MOM_sum_output.F90 | 7 +- src/framework/MOM_checksums.F90 | 139 +++++++++++- src/framework/MOM_coms.F90 | 4 +- src/framework/MOM_domain_init.F90 | 330 ++++++++++++++++++++++++++++ src/framework/MOM_restart.F90 | 10 +- src/framework/MOM_transform_FMS.F90 | 132 +---------- 6 files changed, 474 insertions(+), 148 deletions(-) create mode 100644 src/framework/MOM_domain_init.F90 diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 03204e4322..8b24463af9 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_restart.F90 b/src/framework/MOM_restart.F90 index 73c5cc94e1..8034e19048 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -3,10 +3,10 @@ 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, MOM_write_field, read_field_chksum @@ -14,9 +14,9 @@ module MOM_restart 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_verticalGrid, only : verticalGrid_type implicit none ; private @@ -874,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. @@ -1086,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. @@ -1363,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(:), & diff --git a/src/framework/MOM_transform_FMS.F90 b/src/framework/MOM_transform_FMS.F90 index a6c28717b3..a4a3f7c2c4 100644 --- a/src/framework/MOM_transform_FMS.F90 +++ b/src/framework/MOM_transform_FMS.F90 @@ -3,28 +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 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 -public rotated_mpp_chksum 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 - !> 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 @@ -37,122 +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 - - !> Read a scalar field based on model time !! This function is provided to support the full FMS time_interp_external !! interface. From 81a6ff8ee8f197709fc8357271d907f4e759f7d1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 10 Jan 2021 07:04:33 -0500 Subject: [PATCH 06/11] +Add explicit interface for field_exists to MOM_io Added an explicit interface for field_exists to MOM_io.F90, and added a new optional argument of a MOM_domain_type to field_exists. All answers are bitwise identical, and all previous calls still work exactly as before, but there is a new optional argument. --- src/framework/MOM_io.F90 | 24 +++++++++++++++---- .../lateral/MOM_tidal_forcing.F90 | 2 +- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index eda00bfda0..4e5cd621c4 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -15,9 +15,8 @@ module MOM_io 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 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 : 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 @@ -862,7 +861,7 @@ 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 + 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 @@ -874,6 +873,23 @@ function FMS_file_exists(filename, domain, no_domain) end function FMS_file_exists +!> 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". 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 From 37ef28bc0c32a4d20b92737b76fad2a2d9feadd4 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 10 Jan 2021 11:19:21 -0500 Subject: [PATCH 07/11] +Add get_domain_components to MOM_domains.F90 Added the new interface get_domain_components to MOM_domains.F90 to return the 1-d domains that are make up a 2-d domain, with overloaded variants working on MOM_domain_type or domain2D arguments. The MOM_domains module also now provides access to the domain2D type. All answers are bitwise identical. --- src/framework/MOM_domains.F90 | 31 +++++++++++++++++++++++++++---- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index dc1f8ff867..cc33b238f4 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -15,7 +15,7 @@ module MOM_domains 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 @@ -37,7 +37,7 @@ module MOM_domains implicit none ; private public :: MOM_domains_init, MOM_infra_init, MOM_infra_end, get_domain_extent, get_domain_extent_dsamp2 -public :: create_MOM_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 @@ -52,7 +52,7 @@ module MOM_domains 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 @@ -104,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. @@ -1709,6 +1714,24 @@ subroutine set_MOM_thread_affinity(ocean_nthreads, ocean_hyper_thread) !$ 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. subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, & From 51720fb3357866ec5e6c09c1bee67759f7c7accb Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 10 Jan 2021 11:21:41 -0500 Subject: [PATCH 08/11] +Split MOM_io_wrapper.F90 out from MOM_io.F90 Created the new module MOM_io_wrapper from the contents of MOM_io, so that there are two separate modules, one of which use MOM-specific calls and structures, and another that directly calls or wraps the mpp and FMS I/O interfaces. All of the previous interfaces that were accessible via MOM_io are still being served from the this module, so no changes are needed outside of these two modules. All answers are bitwise identical. --- src/framework/MOM_io.F90 | 596 ++++--------------------------- src/framework/MOM_io_wrapper.F90 | 504 ++++++++++++++++++++++++++ 2 files changed, 565 insertions(+), 535 deletions(-) create mode 100644 src/framework/MOM_io_wrapper.F90 diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 4e5cd621c4..cbede5e3eb 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -3,51 +3,46 @@ module MOM_io ! 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_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_file_parser, only : log_version, param_file_type use MOM_grid, only : ocean_grid_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_exist, field_size, read_data -use fms_io_mod, only : io_infra_end=>fms_io_exit, 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, MOM_write_field, 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 @@ -64,36 +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 - -!> 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 - integer, public :: stdout = stdout_iso !< standard output unit integer, public :: stderr = stderr_iso !< standard output unit @@ -237,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. @@ -287,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 @@ -298,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 @@ -336,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 @@ -432,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) @@ -526,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 @@ -578,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 @@ -604,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 @@ -615,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 @@ -766,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 @@ -786,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 @@ -844,378 +743,6 @@ 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 - -!> 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 - - !> 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 @@ -1229,7 +756,6 @@ subroutine MOM_io_init(param_file) 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..ed4405dbd9 --- /dev/null +++ b/src/framework/MOM_io_wrapper.F90 @@ -0,0 +1,504 @@ +!> 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 : open_file=>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 + +!> 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 From 45d29a99caa03b803974f9986c25c765224374d3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 10 Jan 2021 18:26:24 -0500 Subject: [PATCH 09/11] +Added an explicit interface to open_file Added an explicit interface for open_file to MOM_io_wrapper.F90. This version only includes the optional arguments that are actually used in the MOM6 or SIS2 code, but adds a new optional MOM_domain_type argument. The optional arguments to mpp_open that are not being carried over to the new MOM interface pertain to archaic file formats and will never be used in MOM6. All answers are bitwise identical. --- src/framework/MOM_io.F90 | 4 ++-- src/framework/MOM_io_wrapper.F90 | 30 +++++++++++++++++++++++++++++- src/framework/MOM_restart.F90 | 3 +-- 3 files changed, 32 insertions(+), 5 deletions(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index cbede5e3eb..2cdfcae4b9 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -138,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. @@ -377,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 diff --git a/src/framework/MOM_io_wrapper.F90 b/src/framework/MOM_io_wrapper.F90 index ed4405dbd9..7437b59db1 100644 --- a/src/framework/MOM_io_wrapper.F90 +++ b/src/framework/MOM_io_wrapper.F90 @@ -13,7 +13,7 @@ module MOM_io_wrapper 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 : open_file=>mpp_open, close_file=>mpp_close +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 @@ -160,6 +160,34 @@ function FMS_file_exists(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) diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index 8034e19048..e376d9917b 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -1451,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 From e1ca9a905334880ff720111b15c069f9ad70841d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 10 Jan 2021 21:32:15 -0500 Subject: [PATCH 10/11] Doxygen comments for get_layout_extents arguments Added doxygen comments describing two arguments to get_layout_extents. All answers are bitwise identical. --- src/framework/MOM_domains.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index cc33b238f4..65e7337964 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -2115,8 +2115,10 @@ end subroutine get_global_shape !! 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 - integer, dimension(:), allocatable, intent(inout) :: extent_j + 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) From 77d0cbe9cb948bbe9a7630e26803521ed090d19c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 13 Jan 2021 15:06:48 -0500 Subject: [PATCH 11/11] Corrected a bug using IO_LAYOUT in place of LAYOUT Corrected a bug that was incorrectly using the IO_LAYOUT to set LAYOUT when the IO_LAYOUT uses non-default values greater than 1. All answers are bitwise identical in any cases that worked. --- src/framework/MOM_domains.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 65e7337964..56ac0b3ccf 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -1571,7 +1571,7 @@ subroutine create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar, lay if (io_layout(1) == 0) then MOM_dom%io_layout(1) = layout(1) elseif (io_layout(1) > 1) then - MOM_dom%layout(1) = io_layout(1) + 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) @@ -1582,7 +1582,7 @@ subroutine create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar, lay if (io_layout(2) == 0) then MOM_dom%io_layout(2) = layout(2) elseif (io_layout(2) > 1) then - MOM_dom%layout(2) = io_layout(2) + 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)