MOM_domains.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Describes the decomposed MOM domain and has routines for communications across PEs
6module mom_domains
7
8use mom_coms_infra, only : mom_infra_init, mom_infra_end
9use mom_coms_infra, only : pe_here, root_pe, num_pes, broadcast
10use mom_coms_infra, only : sum_across_pes, min_across_pes, max_across_pes
11use mom_domain_infra, only : mom_domain_type, domain2d, domain1d, group_pass_type
12use mom_domain_infra, only : create_mom_domain, clone_mom_domain, deallocate_mom_domain
13use mom_domain_infra, only : get_domain_extent, get_domain_components, same_domain
14use mom_domain_infra, only : compute_block_extent, get_global_shape
15use mom_domain_infra, only : pass_var, pass_vector, fill_symmetric_edges
16use mom_domain_infra, only : pass_var_start, pass_var_complete
17use mom_domain_infra, only : pass_vector_start, pass_vector_complete
18use mom_domain_infra, only : create_group_pass, do_group_pass
19use mom_domain_infra, only : start_group_pass, complete_group_pass
20use mom_domain_infra, only : rescale_comp_data, global_field, redistribute_array, broadcast_domain
21use mom_domain_infra, only : mom_thread_affinity_set, set_mom_thread_affinity
22use mom_domain_infra, only : agrid, bgrid_ne, cgrid_ne, scalar_pair
23use mom_domain_infra, only : corner, center, north_face, east_face
24use mom_domain_infra, only : to_east, to_west, to_north, to_south, to_all, omit_corners
25use mom_domain_infra, only : compute_extent
26use mom_error_handler, only : mom_error, mom_mesg, note, warning, fatal, is_root_pe
27use mom_file_parser, only : get_param, log_param, log_version, param_file_type
28use mom_io_infra, only : file_exists, read_field, open_ascii_file, close_file, writeonly_file
30use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
32
33implicit none ; private
34
35public :: mom_infra_init, mom_infra_end
36! Domain types and creation and destruction routines
37public :: mom_domain_type, domain2d, domain1d
38public :: mom_domains_init, create_mom_domain, clone_mom_domain, deallocate_mom_domain
39public :: mom_thread_affinity_set, set_mom_thread_affinity
40public :: mom_define_layout
41! Domain query routines
42public :: get_domain_extent, get_domain_components, get_global_shape, same_domain
43public :: pe_here, root_pe, num_pes
44! Blocks are not actively used in MOM6, so this routine could be deprecated.
45public :: compute_block_extent
46! Single call communication routines
47public :: pass_var, pass_vector, fill_symmetric_edges, broadcast
48! Non-blocking communication routines
49public :: pass_var_start, pass_var_complete, pass_vector_start, pass_vector_complete
50! Multi-variable group communication routines and type
51public :: create_group_pass, do_group_pass, group_pass_type, start_group_pass, complete_group_pass
52! Global reduction routines
53public :: sum_across_pes, min_across_pes, max_across_pes
54public :: global_field, redistribute_array, broadcast_domain
55! Simple index-convention-invariant array manipulation routine
56public :: rescale_comp_data
57!> These encoding constants are used to indicate the staggering of scalars and vectors
58public :: agrid, bgrid_ne, cgrid_ne, scalar_pair
59!> These encoding constants are used to indicate the discretization position of a variable
60public :: corner, center, north_face, east_face
61!> These encoding constants indicate communication patterns. In practice they can be added.
62public :: to_east, to_west, to_north, to_south, to_all, omit_corners
63
64contains
65
66!> MOM_domains_init initializes a MOM_domain_type variable, based on the information
67!! read in from a param_file_type, and optionally returns data describing various
68!! properties of the domain type.
69subroutine mom_domains_init(MOM_dom, param_file, symmetric, static_memory, &
70 NIHALO, NJHALO, NIGLOBAL, NJGLOBAL, NIPROC, NJPROC, &
71 min_halo, domain_name, include_name, param_suffix, US, MOM_dom_unmasked)
72 type(mom_domain_type), pointer :: mom_dom !< A pointer to the MOM_domain_type
73 !! being defined here.
74 type(param_file_type), intent(in) :: param_file !< A structure to parse for
75 !! run-time parameters
76 logical, optional, intent(in) :: symmetric !< If present, this specifies
77 !! whether this domain is symmetric, regardless of
78 !! whether the macro SYMMETRIC_MEMORY_ is defined.
79 logical, optional, intent(in) :: static_memory !< If present and true, this
80 !! domain type is set up for static memory and
81 !! error checking of various input values is
82 !! performed against those in the input file.
83 integer, optional, intent(in) :: nihalo !< Default halo sizes, required
84 !! with static memory.
85 integer, optional, intent(in) :: njhalo !< Default halo sizes, required
86 !! with static memory.
87 integer, optional, intent(in) :: niglobal !< Total domain sizes, required
88 !! with static memory.
89 integer, optional, intent(in) :: njglobal !< Total domain sizes, required
90 !! with static memory.
91 integer, optional, intent(in) :: niproc !< Processor counts, required with
92 !! static memory.
93 integer, optional, intent(in) :: njproc !< Processor counts, required with
94 !! static memory.
95 integer, dimension(2), optional, intent(inout) :: min_halo !< If present, this sets the
96 !! minimum halo size for this domain in the i- and j-
97 !! directions, and returns the actual halo size used.
98 character(len=*), optional, intent(in) :: domain_name !< A name for this domain, "MOM"
99 !! if missing.
100 character(len=*), optional, intent(in) :: include_name !< A name for model's include file,
101 !! "MOM_memory.h" if missing.
102 character(len=*), optional, intent(in) :: param_suffix !< A suffix to apply to
103 !! layout-specific parameters.
104 type(unit_scale_type), optional, pointer :: us !< A dimensional unit scaling type
105 type(mom_domain_type), optional, pointer :: mom_dom_unmasked !< Unmasked MOM domain instance.
106 !! Set to null if masking is not enabled.
107
108 ! Local variables
109 integer, dimension(2) :: layout ! The number of logical processors in the i- and j- directions
110 integer, dimension(2) :: auto_layout ! The layout determined by the auto masking routine
111 integer, dimension(2) :: auto_io_layout ! The IO layout determined by the auto masking routine
112 integer, dimension(2) :: layout_unmasked ! A temporary layout for unmasked domain
113 integer, dimension(2) :: io_layout ! The layout of logical processors for input and output
114 !$ integer :: ocean_nthreads ! Number of openMP threads
115 !$ logical :: ocean_omp_hyper_thread ! If true use openMP hyper-threads
116 integer, dimension(2) :: n_global ! The number of i- and j- points in the global computational domain
117 integer, dimension(2) :: n_halo ! The number of i- and j- points in the halos
118 integer :: nihalo_dflt, njhalo_dflt ! The default halo sizes
119 integer :: pes_used ! The number of processors used
120 logical, dimension(2) :: reentrant ! True if the x- and y- directions are periodic.
121 logical :: tripolar_n ! A flag indicating whether there is northern tripolar connectivity
122 logical :: is_static ! If true, static memory is being used for this domain.
123 logical :: is_symmetric ! True if the domain being set up will use symmetric memory.
124 logical :: nonblocking ! If true, nonblocking halo updates will be used.
125 logical :: thin_halos ! If true, If true, optional arguments may be used to specify the
126 ! width of the halos that are updated with each call.
127 logical :: auto_mask_table ! Runtime flag that turns on automatic mask table generator
128 integer :: auto_io_layout_fac ! Used to compute IO layout when auto_mask_table is True.
129 integer :: target_io_pes ! Target number of IO PEs for when auto_mask_table is True.
130 logical :: mask_table_exists ! True if there is a mask table file
131 logical :: is_mom_domain ! True if this domain is being set for MOM, and not another component like SIS2.
132 character(len=128) :: inputdir ! The directory in which to find the diag table
133 character(len=200) :: mask_table ! The file name and later the full path to the diag table
134 character(len=200) :: default_masktable ! Default filename for mask table
135 character(len=64) :: inc_nm ! The name of the memory include file
136 character(len=200) :: mesg ! A string to use for error messages
137
138 integer :: nip_parsed, njp_parsed
139 character(len=8) :: char_xsiz, char_ysiz, char_niglobal, char_njglobal
140 character(len=40) :: nihalo_nm, njhalo_nm, layout_nm, io_layout_nm, masktable_nm
141 character(len=40) :: niproc_nm, njproc_nm
142 character(len=200) :: topo_config
143 integer :: id_clock_auto_mask
144 character(len=:), allocatable :: masktable_desc
145 ! This include declares and sets the variable "version".
146# include "version_variable.h"
147 character(len=40) :: mdl ! This module's name.
148
149 pes_used = num_pes()
150
151 mdl = "MOM_domains"
152
153 is_symmetric = .true. ; if (present(symmetric)) is_symmetric = symmetric
154 if (present(min_halo)) mdl = trim(mdl)//" min_halo"
155
156 inc_nm = "MOM_memory.h" ; if (present(include_name)) inc_nm = trim(include_name)
157
158 nihalo_nm = "NIHALO" ; njhalo_nm = "NJHALO"
159 layout_nm = "LAYOUT" ; io_layout_nm = "IO_LAYOUT" ; masktable_nm = "MASKTABLE"
160 niproc_nm = "NIPROC" ; njproc_nm = "NJPROC"
161 if (present(param_suffix)) then ; if (len(trim(adjustl(param_suffix))) > 0) then
162 nihalo_nm = "NIHALO"//(trim(adjustl(param_suffix)))
163 njhalo_nm = "NJHALO"//(trim(adjustl(param_suffix)))
164 layout_nm = "LAYOUT"//(trim(adjustl(param_suffix)))
165 io_layout_nm = "IO_LAYOUT"//(trim(adjustl(param_suffix)))
166 masktable_nm = "MASKTABLE"//(trim(adjustl(param_suffix)))
167 niproc_nm = "NIPROC"//(trim(adjustl(param_suffix)))
168 njproc_nm = "NJPROC"//(trim(adjustl(param_suffix)))
169 endif ; endif
170
171 is_static = .false. ; if (present(static_memory)) is_static = static_memory
172 if (is_static) then
173 if (.not.present(nihalo)) call mom_error(fatal, "NIHALO must be "// &
174 "present in the call to MOM_domains_init with static memory.")
175 if (.not.present(njhalo)) call mom_error(fatal, "NJHALO must be "// &
176 "present in the call to MOM_domains_init with static memory.")
177 if (.not.present(niglobal)) call mom_error(fatal, "NIGLOBAL must be "// &
178 "present in the call to MOM_domains_init with static memory.")
179 if (.not.present(njglobal)) call mom_error(fatal, "NJGLOBAL must be "// &
180 "present in the call to MOM_domains_init with static memory.")
181 if (.not.present(niproc)) call mom_error(fatal, "NIPROC must be "// &
182 "present in the call to MOM_domains_init with static memory.")
183 if (.not.present(njproc)) call mom_error(fatal, "NJPROC must be "// &
184 "present in the call to MOM_domains_init with static memory.")
185 endif
186
187 ! Read all relevant parameters and write them to the model log.
188 call log_version(param_file, mdl, version, "", log_to_all=.true., layout=.true.)
189 call get_param(param_file, mdl, "REENTRANT_X", reentrant(1), &
190 "If true, the domain is zonally reentrant.", default=.true.)
191 call get_param(param_file, mdl, "REENTRANT_Y", reentrant(2), &
192 "If true, the domain is meridionally reentrant.", &
193 default=.false.)
194 call get_param(param_file, mdl, "TRIPOLAR_N", tripolar_n, &
195 "Use tripolar connectivity at the northern edge of the "//&
196 "domain. With TRIPOLAR_N, NIGLOBAL must be even.", &
197 default=.false.)
198
199# ifndef NOT_SET_AFFINITY
200 !$ if (.not.MOM_thread_affinity_set()) then
201 !$ call get_param(param_file, mdl, "OCEAN_OMP_THREADS", ocean_nthreads, &
202 !$ "The number of OpenMP threads that MOM6 will use.", &
203 !$ default=1, layoutParam=.true.)
204 !$ call get_param(param_file, mdl, "OCEAN_OMP_HYPER_THREAD", ocean_omp_hyper_thread, &
205 !$ "If True, use hyper-threading.", default=.false., layoutParam=.true.)
206 !$ call set_MOM_thread_affinity(ocean_nthreads, ocean_omp_hyper_thread)
207 !$ endif
208# endif
209
210 call log_param(param_file, mdl, "!SYMMETRIC_MEMORY_", is_symmetric, &
211 "If defined, the velocity point data domain includes every face of the "//&
212 "thickness points. In other words, some arrays are larger than others, "//&
213 "depending on where they are on the staggered grid. Also, the starting "//&
214 "index of the velocity-point arrays is usually 0, not 1. "//&
215 "This can only be set at compile time.",&
216 layoutparam=.true.)
217 call get_param(param_file, mdl, "NONBLOCKING_UPDATES", nonblocking, &
218 "If true, non-blocking halo updates may be used.", &
219 default=.false., layoutparam=.true.)
220 call get_param(param_file, mdl, "THIN_HALO_UPDATES", thin_halos, &
221 "If true, optional arguments may be used to specify the width of the "//&
222 "halos that are updated with each call.", &
223 default=.true., layoutparam=.true.)
224
225 nihalo_dflt = 4 ; njhalo_dflt = 4
226 if (present(nihalo)) nihalo_dflt = nihalo
227 if (present(njhalo)) njhalo_dflt = njhalo
228
229 call log_param(param_file, mdl, "!STATIC_MEMORY_", is_static, &
230 "If STATIC_MEMORY_ is defined, the principle variables will have sizes that "//&
231 "are statically determined at compile time. Otherwise the sizes are not "//&
232 "determined until run time. The STATIC option is substantially faster, but "//&
233 "does not allow the PE count to be changed at run time. This can only be "//&
234 "set at compile time.", layoutparam=.true.)
235
236 if (is_static) then
237 call get_param(param_file, mdl, "NIGLOBAL", n_global(1), &
238 "The total number of thickness grid points in the x-direction in the physical "//&
239 "domain. With STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.", &
240 default=niglobal)
241 call get_param(param_file, mdl, "NJGLOBAL", n_global(2), &
242 "The total number of thickness grid points in the y-direction in the physical "//&
243 "domain. With STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.", &
244 default=njglobal)
245 if (n_global(1) /= niglobal) call mom_error(fatal,"MOM_domains_init: " // &
246 "static mismatch for NIGLOBAL_ domain size. Header file does not match input namelist")
247 if (n_global(2) /= njglobal) call mom_error(fatal,"MOM_domains_init: " // &
248 "static mismatch for NJGLOBAL_ domain size. Header file does not match input namelist")
249
250 ! Check the requirement of equal sized compute domains when STATIC_MEMORY_ is used.
251 if ((mod(niglobal, niproc) /= 0) .OR. (mod(njglobal, njproc) /= 0)) then
252 write( char_xsiz, '(I0)' ) niproc
253 write( char_ysiz, '(I0)' ) njproc
254 write( char_niglobal, '(I0)' ) niglobal
255 write( char_njglobal, '(I0)' ) njglobal
256 call mom_error(warning, 'MOM_domains: Processor decomposition (NIPROC_,NJPROC_) = ('//&
257 trim(char_xsiz)//','//trim(char_ysiz)//') does not evenly divide size '//&
258 'set by preprocessor macro ('//trim(char_niglobal)//','//trim(char_njglobal)//').')
259 call mom_error(fatal,'MOM_domains: #undef STATIC_MEMORY_ in '//trim(inc_nm)//' to use '//&
260 'dynamic allocation, or change processor decomposition to evenly divide the domain.')
261 endif
262 else
263 call get_param(param_file, mdl, "NIGLOBAL", n_global(1), &
264 "The total number of thickness grid points in the x-direction in the physical "//&
265 "domain. With STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.", &
266 fail_if_missing=.true.)
267 call get_param(param_file, mdl, "NJGLOBAL", n_global(2), &
268 "The total number of thickness grid points in the y-direction in the physical "//&
269 "domain. With STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.", &
270 fail_if_missing=.true.)
271 endif
272
273 call get_param(param_file, mdl, trim(nihalo_nm), n_halo(1), &
274 "The number of halo points on each side in the x-direction. How this is set "//&
275 "varies with the calling component and static or dynamic memory configuration.", &
276 default=nihalo_dflt)
277 call get_param(param_file, mdl, trim(njhalo_nm), n_halo(2), &
278 "The number of halo points on each side in the y-direction. How this is set "//&
279 "varies with the calling component and static or dynamic memory configuration.", &
280 default=njhalo_dflt)
281 if (present(min_halo)) then
282 n_halo(1) = max(n_halo(1), min_halo(1))
283 min_halo(1) = n_halo(1)
284 n_halo(2) = max(n_halo(2), min_halo(2))
285 min_halo(2) = n_halo(2)
286 ! These are generally used only with static memory, so they are considered layout params.
287 call log_param(param_file, mdl, "!NIHALO min_halo", n_halo(1), layoutparam=.true.)
288 call log_param(param_file, mdl, "!NJHALO min_halo", n_halo(2), layoutparam=.true.)
289 endif
290 if (is_static .and. .not.present(min_halo)) then
291 if (n_halo(1) /= nihalo) call mom_error(fatal,"MOM_domains_init: " // &
292 "static mismatch for "//trim(nihalo_nm)//" domain size")
293 if (n_halo(2) /= njhalo) call mom_error(fatal,"MOM_domains_init: " // &
294 "static mismatch for "//trim(njhalo_nm)//" domain size")
295 endif
296
297 call get_param(param_file, mdl, "INPUTDIR", inputdir, do_not_log=.true., default=".")
298 inputdir = slasher(inputdir)
299
300 is_mom_domain = .true.
301 if (present(domain_name)) then
302 is_mom_domain = (index(domain_name, "MOM") > 1)
303 endif
304
305 if (is_mom_domain) then
306 call get_param(param_file, mdl, "TOPO_CONFIG", topo_config, do_not_log=.true., fail_if_missing=.true.)
307 else ! SIS2 has a default value for TOPO_CONFIG.
308 call get_param(param_file, mdl, "TOPO_CONFIG", topo_config, default="file", do_not_log=.true.)
309 endif
310
311 auto_mask_table = .false.
312 if (.not. present(param_suffix) .and. .not. is_static .and. trim(topo_config) == 'file') then
313 call get_param(param_file, mdl, 'AUTO_MASKTABLE', auto_mask_table, &
314 "Turn on automatic mask table generation to eliminate land blocks.", &
315 default=.false., layoutparam=.true.)
316 endif
317
318 masktable_desc = "A text file to specify n_mask, layout and mask_list. This feature masks out "//&
319 "processors that contain only land points. The first line of mask_table is the "//&
320 "number of regions to be masked out. The second line is the layout of the "//&
321 "model and must be consistent with the actual model layout. The following "//&
322 "(n_mask) lines give the logical positions of the processors that are masked "//&
323 "out. The mask_table can be created by tools like check_mask. The following "//&
324 "example of mask_table masks out 2 processors, (1,2) and (3,6), out of the 24 "//&
325 "in a 4x6 layout: \n 2\n 4,6\n 1,2\n 3,6\n"
326
327 default_masktable = "MOM_mask_table"
328 if (auto_mask_table) default_masktable = "MOM_auto_mask_table"
329
330 call get_param(param_file, mdl, trim(masktable_nm), mask_table, masktable_desc, &
331 default=default_masktable, layoutparam=.true.)
332
333 if (auto_mask_table) then
334 id_clock_auto_mask = cpu_clock_id('(Ocean gen_auto_mask_table)', grain=clock_routine)
335
336 call get_param(param_file, mdl, "TARGET_IO_PES", target_io_pes, &
337 "When AUTO_MASKTABLE is enabled, target number of IO PEs. If the given target number "//&
338 "of IO PEs is not achievable, the target number of IO PEs is set to the nearest smaller "//&
339 "number of PEs that is achievable.", default=1, layoutparam=.true.)
340 if (target_io_pes <= 0) then
341 call mom_error(fatal, 'TARGET_IO_PES must be a nonnegative integer.')
342 endif
343
344 ! Auto-generate a mask file and determine the layout
345 call cpu_clock_begin(id_clock_auto_mask)
346 if (is_root_pe()) then
347 call gen_auto_mask_table(n_global, reentrant, tripolar_n, pes_used, param_file, inputdir, &
348 mask_table, target_io_pes, auto_layout, auto_io_layout, us)
349 endif
350 call broadcast(auto_layout, length=2)
351 call broadcast(auto_io_layout, length=2)
352 call cpu_clock_end(id_clock_auto_mask)
353 endif
354
355 ! First, check the run directory for the mask_table input file.
356 mask_table_exists = file_exists(trim(mask_table))
357 ! If not found, check the input directory
358 if (.not. mask_table_exists) then
359 mask_table = trim(inputdir)//trim(mask_table)
360 mask_table_exists = file_exists(mask_table)
361 endif
362
363 if (is_static) then
364 layout(1) = niproc ; layout(2) = njproc
365 else
366 call get_param(param_file, mdl, trim(layout_nm), layout, &
367 "The processor layout to be used, or 0, 0 to automatically set the layout "//&
368 "based on the number of processors.", defaults=(/0, 0/), do_not_log=.true.)
369 call get_param(param_file, mdl, trim(niproc_nm), nip_parsed, &
370 "The number of processors in the x-direction.", default=-1, do_not_log=.true.)
371 call get_param(param_file, mdl, trim(njproc_nm), njp_parsed, &
372 "The number of processors in the y-direction.", default=-1, do_not_log=.true.)
373 if (nip_parsed > -1) then
374 if ((layout(1) > 0) .and. (layout(1) /= nip_parsed)) &
375 call mom_error(fatal, trim(layout_nm)//" and "//trim(niproc_nm)//" set inconsistently. "//&
376 "Only LAYOUT should be used.")
377 layout(1) = nip_parsed
378 call mom_mesg(trim(niproc_nm)//" used to set "//trim(layout_nm)//" in dynamic mode. "//&
379 "Shift to using "//trim(layout_nm)//" instead.")
380 endif
381 if (njp_parsed > -1) then
382 if ((layout(2) > 0) .and. (layout(2) /= njp_parsed)) &
383 call mom_error(fatal, trim(layout_nm)//" and "//trim(njproc_nm)//" set inconsistently. "//&
384 "Only "//trim(layout_nm)//" should be used.")
385 layout(2) = njp_parsed
386 call mom_mesg(trim(njproc_nm)//" used to set "//trim(layout_nm)//" in dynamic mode. "//&
387 "Shift to using "//trim(layout_nm)//" instead.")
388 endif
389
390 if (auto_mask_table) then
391 if (layout(1) /= 0 .and. layout(1) /= auto_layout(1)) then
392 call mom_error(fatal, "Cannot set LAYOUT or NIPROC when AUTO_MASKTABLE is enabled.")
393 endif
394 if (layout(2) /= 0 .and. layout(2) /= auto_layout(2)) then
395 call mom_error(fatal, "Cannot set LAYOUT or NJPROC when AUTO_MASKTABLE is enabled.")
396 endif
397 layout(:) = auto_layout(:)
398 endif
399
400 if ( (layout(1) == 0) .and. (layout(2) == 0) ) &
401 call mom_define_layout(n_global, pes_used, layout)
402 if ( (layout(1) /= 0) .and. (layout(2) == 0) ) layout(2) = pes_used / layout(1)
403 if ( (layout(1) == 0) .and. (layout(2) /= 0) ) layout(1) = pes_used / layout(2)
404
405 if (layout(1)*layout(2) /= pes_used .and. (.not. mask_table_exists) ) then
406 write(mesg,'("MOM_domains_init: The product of the two components of layout, ", &
407 & I0,", ",I0,", is not the number of PEs used, ",I0,".")') &
408 layout(1), layout(2), pes_used
409 call mom_error(fatal, mesg)
410 endif
411 endif
412 call log_param(param_file, mdl, trim(niproc_nm), layout(1), &
413 "The number of processors in the x-direction. With STATIC_MEMORY_ this "//&
414 "is set in "//trim(inc_nm)//" at compile time.", layoutparam=.true.)
415 call log_param(param_file, mdl, trim(njproc_nm), layout(2), &
416 "The number of processors in the y-direction. With STATIC_MEMORY_ this "//&
417 "is set in "//trim(inc_nm)//" at compile time.", layoutparam=.true.)
418 call log_param(param_file, mdl, trim(layout_nm), layout, &
419 "The processor layout that was actually used.", layoutparam=.true.)
420
421 ! Idiot check that fewer PEs than columns have been requested
422 if (layout(1)*layout(2) > n_global(1)*n_global(2)) then
423 write(mesg,'(a,I0,a,I0,a)') 'You requested to use ', layout(1)*layout(2), &
424 ' PEs but there are only ', n_global(1)*n_global(2), ' columns in the model'
425 call mom_error(fatal, mesg)
426 endif
427
428 if (mask_table_exists) &
429 call mom_error(note, 'MOM_domains_init: reading maskmap information from '//trim(mask_table))
430
431 ! Set up the I/O layout, it will be checked later that it uses an even multiple of the number of
432 ! PEs in each direction.
433 io_layout(:) = (/ 1, 1 /)
434
435 ! Compute a valid IO layout if auto_mask_table is on. Otherwise, read in IO_LAYOUT parameter,
436 if (auto_mask_table) then
437 ! aa: AUTO_IO_LAYOUT_FAC is fragile and should be deprecated/removed in favor of TARGET_IO_PES.
438 call get_param(param_file, mdl, "AUTO_IO_LAYOUT_FAC", auto_io_layout_fac, &
439 "When AUTO_MASKTABLE is enabled, io layout is calculated by performing integer "//&
440 "division of the runtime-determined domain layout with this factor. If the factor "//&
441 "is set to 0 (default), the io layout is set to 1,1. NOTE: TARGET_IO_PES is a more "//&
442 "robust way to set the number of IO PEs when auto masking is turned on.", &
443 default=0, layoutparam=.true.)
444 if (auto_io_layout_fac>1 .and. target_io_pes>1) then
445 call mom_error(fatal, 'AUTO_IO_LAYOUT_FAC and TARGET_IO_PES cannot be set simultaneously.')
446 endif
447 if (target_io_pes>1) then
448 io_layout(1) = auto_io_layout(1)
449 io_layout(2) = auto_io_layout(2)
450 elseif (auto_io_layout_fac>0) then
451 io_layout(1) = max(layout(1)/auto_io_layout_fac, 1)
452 io_layout(2) = max(layout(2)/auto_io_layout_fac, 1)
453 elseif (auto_io_layout_fac<0) then
454 call mom_error(fatal, 'AUTO_IO_LAYOUT_FAC must be a nonnegative integer.')
455 endif
456 call log_param(param_file, mdl, trim(io_layout_nm), io_layout, &
457 "The processor layout to be used, or 0,0 to automatically set the io_layout "//&
458 "to be the same as the layout.", layoutparam=.true.)
459 else
460 call get_param(param_file, mdl, trim(io_layout_nm), io_layout, &
461 "The processor layout to be used, or 0,0 to automatically set the io_layout "//&
462 "to be the same as the layout.", defaults=(/1, 1/), layoutparam=.true.)
463 endif
464
465 ! Create an unmasked domain if requested. This is used for writing out unmasked ocean geometry.
466 if (present(mom_dom_unmasked) .and. mask_table_exists) then
467 call mom_define_layout(n_global, pes_used, layout_unmasked)
468 call create_mom_domain(mom_dom_unmasked, n_global, n_halo, reentrant, tripolar_n, layout_unmasked, &
469 domain_name=domain_name, symmetric=symmetric, thin_halos=thin_halos, &
470 nonblocking=nonblocking)
471 endif
472
473 call create_mom_domain(mom_dom, n_global, n_halo, reentrant, tripolar_n, layout, &
474 io_layout=io_layout, domain_name=domain_name, mask_table=mask_table, &
475 symmetric=symmetric, thin_halos=thin_halos, nonblocking=nonblocking)
476
477end subroutine mom_domains_init
478
479!> Given a global array size and a number of (logical) processors, provide a layout of the
480!! processors in the two directions where the total number of processors is the product of
481!! the two layouts and number of points in the partitioned arrays are as close as possible
482!! to an aspect ratio of 1.
483subroutine mom_define_layout(n_global, ndivs, layout)
484 integer, dimension(2), intent(in) :: n_global !< The total number of gridpoints in 2 directions
485 integer, intent(in) :: ndivs !< The total number of (logical) PEs
486 integer, dimension(2), intent(out) :: layout !< The generated layout of PEs
487
488 ! Local variables
489 integer :: isz, jsz, idiv, jdiv
490
491 ! At present, this algorithm is a copy of mpp_define_layout, but it could perhaps be improved?
492
493 isz = n_global(1) ; jsz = n_global(2)
494 ! First try to divide ndivs to match the domain aspect ratio. If this is not an even
495 ! divisor of ndivs, reduce idiv until a factor is found.
496 idiv = max(nint( sqrt(float(ndivs*isz)/jsz) ), 1)
497 do while( mod(ndivs,idiv) /= 0 )
498 idiv = idiv - 1
499 enddo ! This will terminate at idiv=1 if not before
500 jdiv = ndivs / idiv
501
502 layout = (/ idiv, jdiv /)
503end subroutine mom_define_layout
504
505!> Given a desired number of active npes, generate a layout and mask_table
506subroutine gen_auto_mask_table(n_global, reentrant, tripolar_N, npes, param_file, &
507 inputdir, filename, target_io_pes, layout, io_layout, US)
508 integer, dimension(2), intent(in) :: n_global !< The total number of gridpoints in 2 directions
509 logical, dimension(2), intent(in) :: reentrant !< True if the x- and y- directions are periodic.
510 logical, intent(in) :: tripolar_N !< A flag indicating whether there is n. tripolar connectivity
511 integer, intent(in) :: npes !< The desired number of active PEs.
512 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
513 character(len=128), intent(in) :: inputdir !< INPUTDIR parameter
514 character(len=*), intent(in) :: filename !< Mask table file path (to be auto-generated.)
515 integer, intent(inout) :: target_io_pes !< Target number of IO PEs when auto_mask_table is True.
516 integer, dimension(2), intent(out) :: layout !< The generated layout of PEs (incl. masked blocks)
517 integer, dimension(2), intent(out) :: io_layout !< The generated IO layout based on target_io_pes.
518 type(unit_scale_type), optional, pointer :: US !< A dimensional unit scaling type
519
520 ! Local variables
521 real, dimension(n_global(1), n_global(2)) :: D ! Bathymetric depth (to be read in from TOPO_FILE) [Z ~> m]
522 integer, dimension(:,:), allocatable :: mask ! Cell masks (based on D and MINIMUM_DEPTH)
523 character(len=200) :: topo_filepath, topo_file ! Strings for file/path
524 character(len=200) :: topo_varname ! Variable name in file
525 character(len=200) :: topo_config
526 character(len=40) :: mdl = "gen_auto_mask_table" ! This subroutine's name.
527 integer :: i, j, p, p_up
528 real :: Dmask ! The depth for masking in the same units as D [Z ~> m]
529 real :: min_depth ! The minimum ocean depth in the same units as D [Z ~> m]
530 real :: mask_depth ! The depth shallower than which to mask a point as land. [Z ~> m]
531 real :: glob_ocn_frac ! ratio of ocean points to total number of points [nondim]
532 real :: ar ! layout aspect ratio to check if it is too extreme [nondim]
533 real :: m_to_Z ! A conversion factor from m to height units [Z m-1 ~> 1]
534 integer :: nx, ny ! global domain sizes
535 integer, parameter :: ibuf=2, jbuf=2
536 real, parameter :: r_extreme = 4.0 ! aspect ratio limit (>1) for a layout to be considered [nondim]
537 integer :: num_masked_blocks
538 integer, allocatable :: mask_table(:,:)
539 integer :: max_feasible_p ! max division count that leads to enough land block masking to arrive at target compute PEs
540 real, parameter :: pfrac = 0.01 ! fraction by which to reduce the max_feasible_p if target IO PEs is not achievable.
541 ! (This is to provide a wiggle room for the target IO PEs to be achieved.) [nondim]
542 character(len=200) :: mesg ! A string to use for error messages
543
544 m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
545
546 io_layout = [1, 1]
547
548 ! Read in params necessary for auto-masking
549 call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
550 units="m", default=0.0, scale=m_to_z, do_not_log=.true.)
551 call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, &
552 units="m", default=-9999.0, scale=m_to_z, do_not_log=.true.)
553 call get_param(param_file, mdl, "TOPO_CONFIG", topo_config, default="file", do_not_log=.true.)
554 call get_param(param_file, mdl, "TOPO_FILE", topo_file, do_not_log=.true., default="topog.nc")
555 call get_param(param_file, mdl, "TOPO_VARNAME", topo_varname, do_not_log=.true., default="depth")
556 topo_filepath = trim(inputdir)//trim(topo_file)
557
558 ! Sanity checks
559 if (.not. is_root_pe()) then
560 call mom_error(fatal, 'gen_auto_mask_table should only be called by the root PE.')
561 endif
562 if (trim(topo_config) /= "file") then
563 call mom_error(fatal, 'Auto mask table only works with TOPO_CONFIG="file"')
564 endif
565 if (.not.file_exists(topo_filepath)) then
566 call mom_error(fatal, " gen_auto_mask_table: Unable to open "//trim(topo_filepath))
567 endif
568
569 nx = n_global(1)
570 ny = n_global(2)
571
572 ! Read in bathymetric depth.
573 d(:,:) = -9.0e30 * m_to_z ! Initializing to a very large negative depth (tall mountains) everywhere.
574 call read_field(topo_filepath, trim(topo_varname), d, start=(/1, 1/), nread=n_global, no_domain=.true., &
575 scale=m_to_z)
576
577 allocate(mask(nx+2*ibuf, ny+2*jbuf), source=0)
578
579 ! Determine cell masks
580 dmask = mask_depth
581 if (mask_depth == -9999.0*m_to_z) dmask = min_depth
582 do i=1,nx ; do j=1,ny
583 if (d(i,j) <= dmask) then
584 mask(i+ibuf,j+jbuf) = 0
585 else
586 mask(i+ibuf,j+jbuf) = 1
587 endif
588 enddo ; enddo
589
590 ! fill in buffer cells
591
592 if (reentrant(1)) then ! REENTRANT_X
593 mask(1:ibuf, :) = mask(nx+1:nx+ibuf, :)
594 mask(ibuf+nx+1:nx+2*ibuf, :) = mask(ibuf+1:2*ibuf, :)
595 endif
596
597 if (reentrant(2)) then ! REENTRANT_Y
598 mask(:, 1:jbuf) = mask(:, ny+1:ny+jbuf)
599 mask(:, jbuf+ny+1:ny+2*jbuf) = mask(:, jbuf+1:2*jbuf)
600 endif
601
602 if (tripolar_n) then ! TRIPOLAR_N
603 do i=1,nx+2*ibuf
604 do j=1,jbuf
605 mask(i, jbuf+ny+j) = mask(nx+2*ibuf+1-i, jbuf+ny+1-j)
606 enddo
607 enddo
608 endif
609
610 ! Tripolar Stitch Fix: In cases where masking is asymmetrical across the tripolar stitch, there's a possibility
611 ! that certain unmasked blocks won't be able to obtain grid metrics from the halo points. This occurs when the
612 ! neighboring block on the opposite side of the tripolar stitch is masked. As a consequence, certain metrics like
613 ! dxT and dyT may be calculated through extrapolation (refer to extrapolate_metric), potentially leading to the
614 ! generation of non-positive values. This can result in divide-by-zero errors elsewhere, e.g., in MOM_hor_visc.F90.
615 ! Currently, the safest and most general solution is to prohibit masking along the tripolar stitch:
616 if (tripolar_n) then
617 mask(:, jbuf+ny) = 1
618 endif
619
620 glob_ocn_frac = real(sum(mask(1+ibuf:nx+ibuf, 1+jbuf:ny+jbuf))) / (nx * ny)
621 max_feasible_p = 0
622
623 ! Iteratively check for all possible division counts starting from the upper bound of npes/glob_ocn_frac.
624 ! The first encountered feasible division count is stored in max_feasible_p. If the target_io_pes is not
625 ! achievable with this layout, the iteration continues until max_feasible_p * (1.0 - pfrac) is reached or the
626 ! target_io_pes is satisfiable. If not, the target_io_pes is decremented and the iteration is re-done from
627 ! max_feasible_p to max_feasible_p * (1.0 - pfrac).
628 outer: do i = target_io_pes, 1, -1
629
630 if (max_feasible_p == 0) then ! first iteration
631 p_up = ceiling(npes/glob_ocn_frac)
632 else ! subsequent iterations with reduced target_io_pes
633 p_up = max_feasible_p
634 endif
635
636 do p = p_up, npes, -1
637
638 ! compute the layout for the current division count, p
639 call mom_define_layout(n_global, p, layout)
640
641 ! don't bother checking this p if the aspect ratio is extreme
642 ar = (real(nx)/layout(1)) / (real(ny)/layout(2))
643 if ( ar * r_extreme < 1 .or. r_extreme < ar ) cycle
644
645 ! Get the number of masked_blocks for this particular division count
646 call determine_land_blocks(mask, nx, ny, layout(1), layout(2), ibuf, jbuf, num_masked_blocks)
647
648 ! If we can eliminate enough blocks to reach the target compute npes, check if the target IO PEs can also be
649 ! satisfied with this layout. If so, terminate the whole iteration.
650 if (p-num_masked_blocks <= npes) then ! We can eliminate enough blocks to reach the target compute npes
651 if (max_feasible_p == 0) max_feasible_p = p
652 if (mod(layout(1)*layout(2), i) == 0) then
653 io_layout = auto_determine_io_layout(layout(1), layout(2), i)
654 ! skip if aspect ratio is extreme
655 ar = (real(layout(1))/io_layout(1)) / (real(layout(2))/io_layout(2))
656 if ( ar * r_extreme < 1 .or. r_extreme < ar ) cycle
657 call mom_error(note, "Found the optimum layout for auto-masking. Terminating iteration.")
658 if (i /= target_io_pes) then
659 write(mesg,'(a,i0)') "For compatibility with compute layout, changed number of IO PEs to: ", i
660 call mom_error(note, mesg)
661 endif
662 exit outer
663 endif
664 endif
665
666 ! Do not reduce the division count too much to satisfy the target_io_pes. Instead, re-do the iteration
667 ! with a reduced target_io_pes.
668 if (p <= max_feasible_p * (1.0 - pfrac)) exit ! (inner do loop)
669 enddo
670 enddo outer
671
672 if (num_masked_blocks == 0) then
673 call mom_error(fatal, "Couldn't auto-eliminate any land blocks. Try to increase the number "//&
674 "of MOM6 PEs or set AUTO_MASKTABLE to False.")
675 endif
676
677 ! Call determine_land_blocks once again, this time to retrieve and write out the mask_table.
678 allocate(mask_table(num_masked_blocks,2))
679 call determine_land_blocks(mask, nx, ny, layout(1), layout(2), ibuf, jbuf, num_masked_blocks, mask_table)
680 call write_auto_mask_file(mask_table, layout, npes, filename)
681 deallocate(mask_table)
682 deallocate(mask)
683
684end subroutine gen_auto_mask_table
685
686!> Given a number of domain divisions, compute the max number of land blocks that can be eliminated,
687!! and return the resulting mask table if requested.
688subroutine determine_land_blocks(mask, nx, ny, idiv, jdiv, ibuf, jbuf, num_masked_blocks, mask_table)
689 integer, dimension(:,:), intent(in) :: mask !< cell masks based on depth and MINIMUM_DEPTH
690 integer, intent(in) :: nx !< Total number of gridpoints in x-dir (global)
691 integer, intent(in) :: ny !< Total number of gridpoints in y-dir (global)
692 integer, intent(in) :: idiv !< number of divisions along x-dir
693 integer, intent(in) :: jdiv !< number of divisions along y-dir
694 integer, intent(in) :: ibuf !< number of buffer cells in x-dir.
695 !! (not necessarily the same as NIHALO)
696 integer, intent(in) :: jbuf !< number of buffer cells in y-dir.
697 !! (not necessarily the same as NJHALO)
698 integer, intent(out) :: num_masked_blocks !< the final number of masked blocks
699 integer, intent(out), optional :: mask_table(:,:) !< the resulting array of mask_table
700 ! integer
701 integer, dimension(idiv) :: ibegin !< The starting index of each division along x axis
702 integer, dimension(idiv) :: iend !< The ending index of each division along x axis
703 integer, dimension(jdiv) :: jbegin !< The starting index of each division along y axis
704 integer, dimension(jdiv) :: jend !< The ending index of each division along y axis
705 integer :: i, j, ib, ie, jb, je
706
707 call compute_extent(1, nx, idiv, ibegin, iend)
708 call compute_extent(1, ny, jdiv, jbegin, jend)
709
710 num_masked_blocks = 0
711
712 do i=1,idiv
713 ib = ibegin(i)
714 ie = iend(i) + 2 * ibuf
715 do j=1,jdiv
716 jb = jbegin(j)
717 je = jend(j) + 2 * jbuf
718
719 if (any(mask(ib:ie,jb:je)==1)) cycle
720
721 num_masked_blocks = num_masked_blocks + 1
722
723 if (present(mask_table)) then
724 if ( num_masked_blocks > size(mask_table, dim=1)) then
725 call mom_error(fatal, "The mask_table argument passed to determine_land_blocks() has insufficient size.")
726 endif
727
728 mask_table(num_masked_blocks,1) = i
729 mask_table(num_masked_blocks,2) = j
730 endif
731 enddo
732 enddo
733
734end subroutine determine_land_blocks
735
736!> Write out the auto-generated mask information to a file in the run directory.
737subroutine write_auto_mask_file(mask_table, layout, npes, filename)
738 integer, intent(in) :: mask_table(:,:) !> mask table array to be written out.
739 integer, dimension(2), intent(in) :: layout !> PE layout
740 integer, intent(in) :: npes !> Number of divisions (incl. eliminated ones)
741 character(len=*), intent(in) :: filename !> file name for the mask_table to be written
742 ! local
743 integer :: file_ascii= -1 !< The unit number of the auto-generated mask_file file.
744 integer :: true_num_masked_blocks
745 integer :: p
746
747 ! Eliminate only enough blocks to ensure that the number of active blocks precisely matches the target npes.
748 true_num_masked_blocks = layout(1) * layout(2) - npes
749
750 call open_ascii_file(file_ascii, trim(filename), action=writeonly_file)
751 write(file_ascii, '(I0)') true_num_masked_blocks
752 write(file_ascii, '(I0,",",I0)') layout(1), layout(2)
753 do p = 1, true_num_masked_blocks
754 write(file_ascii, '(I0,",",I0)') mask_table(p,1), mask_table(p,2)
755 enddo
756 call close_file(file_ascii)
757end subroutine write_auto_mask_file
758
759!> Computes the io layout based on the domain layout and a target number of IO PEs.
760function auto_determine_io_layout(idiv, jdiv, nio) result(best_io_layout)
761 integer, intent(in) :: idiv ! The number of compute divisions along the x-axis
762 integer, intent(in) :: jdiv ! The number of compute divisions along the y-axis
763 integer, intent(in) :: nio ! Target number of IO PEs, s.t., (idiv_io * jdiv_io) % nio == 0
764 ! return
765 integer :: io_layout(2) ! Temporary IO layout
766 integer :: best_io_layout(2) ! Best IO layout
767 integer :: f, best_idiv_io, best_jdiv_io
768 real :: ratio_diff, min_ratio_diff
769
770 if (mod(idiv*jdiv, nio) /= 0) then
771 call mom_error(fatal, "The product of the compute layout must be divisible by the target number of IO PEs.")
772 endif
773
774 min_ratio_diff = 1.0e30 ! Large initial value
775 best_io_layout = [1, nio]
776
777 ! Iterate over all factors of nio
778 do f = 1, nio
779 if (mod(nio, f) /= 0) cycle
780 io_layout = [f, nio / f]
781
782 ! Check divisibility constraints
783 if (mod(idiv, io_layout(1)) == 0 .and. mod(jdiv, io_layout(2)) == 0) then
784 ratio_diff = abs(real(io_layout(1)) / real(io_layout(2)) - real(idiv) / real(jdiv))
785
786 ! Update best choice if ratio_diff is smaller
787 if (ratio_diff < min_ratio_diff) then
788 min_ratio_diff = ratio_diff
789 best_io_layout = [io_layout(1), io_layout(2)]
790 endif
791 endif
792 enddo
793
794end function auto_determine_io_layout
795
796end module mom_domains