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=64) :: inc_nm ! The name of the memory include file
135 character(len=200) :: mesg ! A string to use for error messages
136
137 integer :: nip_parsed, njp_parsed
138 character(len=8) :: char_xsiz, char_ysiz, char_niglobal, char_njglobal
139 character(len=40) :: nihalo_nm, njhalo_nm, layout_nm, io_layout_nm, masktable_nm
140 character(len=40) :: niproc_nm, njproc_nm
141 character(len=200) :: topo_config
142 integer :: id_clock_auto_mask
143 character(len=:), allocatable :: masktable_desc
144 character(len=:), allocatable :: auto_mask_table_fname ! Auto-generated mask table file name
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 if (auto_mask_table) then
328 id_clock_auto_mask = cpu_clock_id('(Ocean gen_auto_mask_table)', grain=clock_routine)
329 auto_mask_table_fname = "MOM_auto_mask_table"
330
331 call get_param(param_file, mdl, "TARGET_IO_PES", target_io_pes, &
332 "When AUTO_MASKTABLE is enabled, target number of IO PEs. If the given target number "//&
333 "of IO PEs is not achievable, the target number of IO PEs is set to the nearest smaller "//&
334 "number of PEs that is achievable.", default=1, layoutparam=.true.)
335 if (target_io_pes <= 0) then
336 call mom_error(fatal, 'TARGET_IO_PES must be a nonnegative integer.')
337 endif
338
339 ! Auto-generate a mask file and determine the layout
340 call cpu_clock_begin(id_clock_auto_mask)
341 if (is_root_pe()) then
342 call gen_auto_mask_table(n_global, reentrant, tripolar_n, pes_used, param_file, inputdir, &
343 auto_mask_table_fname, target_io_pes, auto_layout, auto_io_layout, us)
344 endif
345 call broadcast(auto_layout, length=2)
346 call broadcast(auto_io_layout, length=2)
347 call cpu_clock_end(id_clock_auto_mask)
348
349 mask_table = auto_mask_table_fname
350 call log_param(param_file, mdl, trim(masktable_nm), mask_table, masktable_desc, &
351 default="MOM_mask_table", layoutparam=.true.)
352 else
353 call get_param(param_file, mdl, trim(masktable_nm), mask_table, masktable_desc, &
354 default="MOM_mask_table", layoutparam=.true.)
355 endif
356
357 ! First, check the run directory for the mask_table input file.
358 mask_table_exists = file_exists(trim(mask_table))
359 ! If not found, check the input directory
360 if (.not. mask_table_exists) then
361 mask_table = trim(inputdir)//trim(mask_table)
362 mask_table_exists = file_exists(mask_table)
363 endif
364
365 if (is_static) then
366 layout(1) = niproc ; layout(2) = njproc
367 else
368 call get_param(param_file, mdl, trim(layout_nm), layout, &
369 "The processor layout to be used, or 0, 0 to automatically set the layout "//&
370 "based on the number of processors.", defaults=(/0, 0/), do_not_log=.true.)
371 call get_param(param_file, mdl, trim(niproc_nm), nip_parsed, &
372 "The number of processors in the x-direction.", default=-1, do_not_log=.true.)
373 call get_param(param_file, mdl, trim(njproc_nm), njp_parsed, &
374 "The number of processors in the y-direction.", default=-1, do_not_log=.true.)
375 if (nip_parsed > -1) then
376 if ((layout(1) > 0) .and. (layout(1) /= nip_parsed)) &
377 call mom_error(fatal, trim(layout_nm)//" and "//trim(niproc_nm)//" set inconsistently. "//&
378 "Only LAYOUT should be used.")
379 layout(1) = nip_parsed
380 call mom_mesg(trim(niproc_nm)//" used to set "//trim(layout_nm)//" in dynamic mode. "//&
381 "Shift to using "//trim(layout_nm)//" instead.")
382 endif
383 if (njp_parsed > -1) then
384 if ((layout(2) > 0) .and. (layout(2) /= njp_parsed)) &
385 call mom_error(fatal, trim(layout_nm)//" and "//trim(njproc_nm)//" set inconsistently. "//&
386 "Only "//trim(layout_nm)//" should be used.")
387 layout(2) = njp_parsed
388 call mom_mesg(trim(njproc_nm)//" used to set "//trim(layout_nm)//" in dynamic mode. "//&
389 "Shift to using "//trim(layout_nm)//" instead.")
390 endif
391
392 if (auto_mask_table) then
393 if (layout(1) /= 0 .and. layout(1) /= auto_layout(1)) then
394 call mom_error(fatal, "Cannot set LAYOUT or NIPROC when AUTO_MASKTABLE is enabled.")
395 endif
396 if (layout(2) /= 0 .and. layout(2) /= auto_layout(2)) then
397 call mom_error(fatal, "Cannot set LAYOUT or NJPROC when AUTO_MASKTABLE is enabled.")
398 endif
399 layout(:) = auto_layout(:)
400 endif
401
402 if ( (layout(1) == 0) .and. (layout(2) == 0) ) &
403 call mom_define_layout(n_global, pes_used, layout)
404 if ( (layout(1) /= 0) .and. (layout(2) == 0) ) layout(2) = pes_used / layout(1)
405 if ( (layout(1) == 0) .and. (layout(2) /= 0) ) layout(1) = pes_used / layout(2)
406
407 if (layout(1)*layout(2) /= pes_used .and. (.not. mask_table_exists) ) then
408 write(mesg,'("MOM_domains_init: The product of the two components of layout, ", &
409 & I0,", ",I0,", is not the number of PEs used, ",I0,".")') &
410 layout(1), layout(2), pes_used
411 call mom_error(fatal, mesg)
412 endif
413 endif
414 call log_param(param_file, mdl, trim(niproc_nm), layout(1), &
415 "The number of processors in the x-direction. With STATIC_MEMORY_ this "//&
416 "is set in "//trim(inc_nm)//" at compile time.", layoutparam=.true.)
417 call log_param(param_file, mdl, trim(njproc_nm), layout(2), &
418 "The number of processors in the y-direction. With STATIC_MEMORY_ this "//&
419 "is set in "//trim(inc_nm)//" at compile time.", layoutparam=.true.)
420 call log_param(param_file, mdl, trim(layout_nm), layout, &
421 "The processor layout that was actually used.", layoutparam=.true.)
422
423 ! Idiot check that fewer PEs than columns have been requested
424 if (layout(1)*layout(2) > n_global(1)*n_global(2)) then
425 write(mesg,'(a,I0,a,I0,a)') 'You requested to use ', layout(1)*layout(2), &
426 ' PEs but there are only ', n_global(1)*n_global(2), ' columns in the model'
427 call mom_error(fatal, mesg)
428 endif
429
430 if (mask_table_exists) &
431 call mom_error(note, 'MOM_domains_init: reading maskmap information from '//trim(mask_table))
432
433 ! Set up the I/O layout, it will be checked later that it uses an even multiple of the number of
434 ! PEs in each direction.
435 io_layout(:) = (/ 1, 1 /)
436
437 ! Compute a valid IO layout if auto_mask_table is on. Otherwise, read in IO_LAYOUT parameter,
438 if (auto_mask_table) then
439 ! aa: AUTO_IO_LAYOUT_FAC is fragile and should be deprecated/removed in favor of TARGET_IO_PES.
440 call get_param(param_file, mdl, "AUTO_IO_LAYOUT_FAC", auto_io_layout_fac, &
441 "When AUTO_MASKTABLE is enabled, io layout is calculated by performing integer "//&
442 "division of the runtime-determined domain layout with this factor. If the factor "//&
443 "is set to 0 (default), the io layout is set to 1,1. NOTE: TARGET_IO_PES is a more "//&
444 "robust way to set the number of IO PEs when auto masking is turned on.", &
445 default=0, layoutparam=.true.)
446 if (auto_io_layout_fac>1 .and. target_io_pes>1) then
447 call mom_error(fatal, 'AUTO_IO_LAYOUT_FAC and TARGET_IO_PES cannot be set simultaneously.')
448 endif
449 if (target_io_pes>1) then
450 io_layout(1) = auto_io_layout(1)
451 io_layout(2) = auto_io_layout(2)
452 elseif (auto_io_layout_fac>0) then
453 io_layout(1) = max(layout(1)/auto_io_layout_fac, 1)
454 io_layout(2) = max(layout(2)/auto_io_layout_fac, 1)
455 elseif (auto_io_layout_fac<0) then
456 call mom_error(fatal, 'AUTO_IO_LAYOUT_FAC must be a nonnegative integer.')
457 endif
458 call log_param(param_file, mdl, trim(io_layout_nm), io_layout, &
459 "The processor layout to be used, or 0,0 to automatically set the io_layout "//&
460 "to be the same as the layout.", layoutparam=.true.)
461 else
462 call get_param(param_file, mdl, trim(io_layout_nm), io_layout, &
463 "The processor layout to be used, or 0,0 to automatically set the io_layout "//&
464 "to be the same as the layout.", defaults=(/1, 1/), layoutparam=.true.)
465 endif
466
467 ! Create an unmasked domain if requested. This is used for writing out unmasked ocean geometry.
468 if (present(mom_dom_unmasked) .and. mask_table_exists) then
469 call mom_define_layout(n_global, pes_used, layout_unmasked)
470 call create_mom_domain(mom_dom_unmasked, n_global, n_halo, reentrant, tripolar_n, layout_unmasked, &
471 domain_name=domain_name, symmetric=symmetric, thin_halos=thin_halos, &
472 nonblocking=nonblocking)
473 endif
474
475 call create_mom_domain(mom_dom, n_global, n_halo, reentrant, tripolar_n, layout, &
476 io_layout=io_layout, domain_name=domain_name, mask_table=mask_table, &
477 symmetric=symmetric, thin_halos=thin_halos, nonblocking=nonblocking)
478
479end subroutine mom_domains_init
480
481!> Given a global array size and a number of (logical) processors, provide a layout of the
482!! processors in the two directions where the total number of processors is the product of
483!! the two layouts and number of points in the partitioned arrays are as close as possible
484!! to an aspect ratio of 1.
485subroutine mom_define_layout(n_global, ndivs, layout)
486 integer, dimension(2), intent(in) :: n_global !< The total number of gridpoints in 2 directions
487 integer, intent(in) :: ndivs !< The total number of (logical) PEs
488 integer, dimension(2), intent(out) :: layout !< The generated layout of PEs
489
490 ! Local variables
491 integer :: isz, jsz, idiv, jdiv
492
493 ! At present, this algorithm is a copy of mpp_define_layout, but it could perhaps be improved?
494
495 isz = n_global(1) ; jsz = n_global(2)
496 ! First try to divide ndivs to match the domain aspect ratio. If this is not an even
497 ! divisor of ndivs, reduce idiv until a factor is found.
498 idiv = max(nint( sqrt(float(ndivs*isz)/jsz) ), 1)
499 do while( mod(ndivs,idiv) /= 0 )
500 idiv = idiv - 1
501 enddo ! This will terminate at idiv=1 if not before
502 jdiv = ndivs / idiv
503
504 layout = (/ idiv, jdiv /)
505end subroutine mom_define_layout
506
507!> Given a desired number of active npes, generate a layout and mask_table
508subroutine gen_auto_mask_table(n_global, reentrant, tripolar_N, npes, param_file, &
509 inputdir, filename, target_io_pes, layout, io_layout, US)
510 integer, dimension(2), intent(in) :: n_global !< The total number of gridpoints in 2 directions
511 logical, dimension(2), intent(in) :: reentrant !< True if the x- and y- directions are periodic.
512 logical, intent(in) :: tripolar_N !< A flag indicating whether there is n. tripolar connectivity
513 integer, intent(in) :: npes !< The desired number of active PEs.
514 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
515 character(len=128), intent(in) :: inputdir !< INPUTDIR parameter
516 character(len=:), allocatable, intent(in) :: filename !< Mask table file path (to be auto-generated.)
517 integer, intent(inout) :: target_io_pes !< Target number of IO PEs when auto_mask_table is True.
518 integer, dimension(2), intent(out) :: layout !< The generated layout of PEs (incl. masked blocks)
519 integer, dimension(2), intent(out) :: io_layout !< The generated IO layout based on target_io_pes.
520 type(unit_scale_type), optional, pointer :: US !< A dimensional unit scaling type
521
522 ! Local variables
523 real, dimension(n_global(1), n_global(2)) :: D ! Bathymetric depth (to be read in from TOPO_FILE) [Z ~> m]
524 integer, dimension(:,:), allocatable :: mask ! Cell masks (based on D and MINIMUM_DEPTH)
525 character(len=200) :: topo_filepath, topo_file ! Strings for file/path
526 character(len=200) :: topo_varname ! Variable name in file
527 character(len=200) :: topo_config
528 character(len=40) :: mdl = "gen_auto_mask_table" ! This subroutine's name.
529 integer :: i, j, p, p_up
530 real :: Dmask ! The depth for masking in the same units as D [Z ~> m]
531 real :: min_depth ! The minimum ocean depth in the same units as D [Z ~> m]
532 real :: mask_depth ! The depth shallower than which to mask a point as land. [Z ~> m]
533 real :: glob_ocn_frac ! ratio of ocean points to total number of points [nondim]
534 real :: ar ! layout aspect ratio to check if it is too extreme [nondim]
535 real :: m_to_Z ! A conversion factor from m to height units [Z m-1 ~> 1]
536 integer :: nx, ny ! global domain sizes
537 integer, parameter :: ibuf=2, jbuf=2
538 real, parameter :: r_extreme = 4.0 ! aspect ratio limit (>1) for a layout to be considered [nondim]
539 integer :: num_masked_blocks
540 integer, allocatable :: mask_table(:,:)
541 integer :: max_feasible_p ! max division count that leads to enough land block masking to arrive at target compute PEs
542 real, parameter :: pfrac = 0.01 ! fraction by which to reduce the max_feasible_p if target IO PEs is not achievable.
543 ! (This is to provide a wiggle room for the target IO PEs to be achieved.) [nondim]
544 character(len=200) :: mesg ! A string to use for error messages
545
546 m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
547
548 io_layout = [1, 1]
549
550 ! Read in params necessary for auto-masking
551 call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
552 units="m", default=0.0, scale=m_to_z, do_not_log=.true.)
553 call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, &
554 units="m", default=-9999.0, scale=m_to_z, do_not_log=.true.)
555 call get_param(param_file, mdl, "TOPO_CONFIG", topo_config, default="file", do_not_log=.true.)
556 call get_param(param_file, mdl, "TOPO_FILE", topo_file, do_not_log=.true., default="topog.nc")
557 call get_param(param_file, mdl, "TOPO_VARNAME", topo_varname, do_not_log=.true., default="depth")
558 topo_filepath = trim(inputdir)//trim(topo_file)
559
560 ! Sanity checks
561 if (.not. is_root_pe()) then
562 call mom_error(fatal, 'gen_auto_mask_table should only be called by the root PE.')
563 endif
564 if (trim(topo_config) /= "file") then
565 call mom_error(fatal, 'Auto mask table only works with TOPO_CONFIG="file"')
566 endif
567 if (.not.file_exists(topo_filepath)) then
568 call mom_error(fatal, " gen_auto_mask_table: Unable to open "//trim(topo_filepath))
569 endif
570
571 nx = n_global(1)
572 ny = n_global(2)
573
574 ! Read in bathymetric depth.
575 d(:,:) = -9.0e30 * m_to_z ! Initializing to a very large negative depth (tall mountains) everywhere.
576 call read_field(topo_filepath, trim(topo_varname), d, start=(/1, 1/), nread=n_global, no_domain=.true., &
577 scale=m_to_z)
578
579 allocate(mask(nx+2*ibuf, ny+2*jbuf), source=0)
580
581 ! Determine cell masks
582 dmask = mask_depth
583 if (mask_depth == -9999.0*m_to_z) dmask = min_depth
584 do i=1,nx ; do j=1,ny
585 if (d(i,j) <= dmask) then
586 mask(i+ibuf,j+jbuf) = 0
587 else
588 mask(i+ibuf,j+jbuf) = 1
589 endif
590 enddo ; enddo
591
592 ! fill in buffer cells
593
594 if (reentrant(1)) then ! REENTRANT_X
595 mask(1:ibuf, :) = mask(nx+1:nx+ibuf, :)
596 mask(ibuf+nx+1:nx+2*ibuf, :) = mask(ibuf+1:2*ibuf, :)
597 endif
598
599 if (reentrant(2)) then ! REENTRANT_Y
600 mask(:, 1:jbuf) = mask(:, ny+1:ny+jbuf)
601 mask(:, jbuf+ny+1:ny+2*jbuf) = mask(:, jbuf+1:2*jbuf)
602 endif
603
604 if (tripolar_n) then ! TRIPOLAR_N
605 do i=1,nx+2*ibuf
606 do j=1,jbuf
607 mask(i, jbuf+ny+j) = mask(nx+2*ibuf+1-i, jbuf+ny+1-j)
608 enddo
609 enddo
610 endif
611
612 ! Tripolar Stitch Fix: In cases where masking is asymmetrical across the tripolar stitch, there's a possibility
613 ! that certain unmasked blocks won't be able to obtain grid metrics from the halo points. This occurs when the
614 ! neighboring block on the opposite side of the tripolar stitch is masked. As a consequence, certain metrics like
615 ! dxT and dyT may be calculated through extrapolation (refer to extrapolate_metric), potentially leading to the
616 ! generation of non-positive values. This can result in divide-by-zero errors elsewhere, e.g., in MOM_hor_visc.F90.
617 ! Currently, the safest and most general solution is to prohibit masking along the tripolar stitch:
618 if (tripolar_n) then
619 mask(:, jbuf+ny) = 1
620 endif
621
622 glob_ocn_frac = real(sum(mask(1+ibuf:nx+ibuf, 1+jbuf:ny+jbuf))) / (nx * ny)
623 max_feasible_p = 0
624
625 ! Iteratively check for all possible division counts starting from the upper bound of npes/glob_ocn_frac.
626 ! The first encountered feasible division count is stored in max_feasible_p. If the target_io_pes is not
627 ! achievable with this layout, the iteration continues until max_feasible_p * (1.0 - pfrac) is reached or the
628 ! target_io_pes is satisfiable. If not, the target_io_pes is decremented and the iteration is re-done from
629 ! max_feasible_p to max_feasible_p * (1.0 - pfrac).
630 outer: do i = target_io_pes, 1, -1
631
632 if (max_feasible_p == 0) then ! first iteration
633 p_up = ceiling(npes/glob_ocn_frac)
634 else ! subsequent iterations with reduced target_io_pes
635 p_up = max_feasible_p
636 endif
637
638 do p = p_up, npes, -1
639
640 ! compute the layout for the current division count, p
641 call mom_define_layout(n_global, p, layout)
642
643 ! don't bother checking this p if the aspect ratio is extreme
644 ar = (real(nx)/layout(1)) / (real(ny)/layout(2))
645 if ( ar * r_extreme < 1 .or. r_extreme < ar ) cycle
646
647 ! Get the number of masked_blocks for this particular division count
648 call determine_land_blocks(mask, nx, ny, layout(1), layout(2), ibuf, jbuf, num_masked_blocks)
649
650 ! If we can eliminate enough blocks to reach the target compute npes, check if the target IO PEs can also be
651 ! satisfied with this layout. If so, terminate the whole iteration.
652 if (p-num_masked_blocks <= npes) then ! We can eliminate enough blocks to reach the target compute npes
653 if (max_feasible_p == 0) max_feasible_p = p
654 if (mod(layout(1)*layout(2), i) == 0) then
655 io_layout = auto_determine_io_layout(layout(1), layout(2), i)
656 ! skip if aspect ratio is extreme
657 ar = (real(layout(1))/io_layout(1)) / (real(layout(2))/io_layout(2))
658 if ( ar * r_extreme < 1 .or. r_extreme < ar ) cycle
659 call mom_error(note, "Found the optimum layout for auto-masking. Terminating iteration.")
660 if (i /= target_io_pes) then
661 write(mesg,'(a,i0)') "For compatibility with compute layout, changed number of IO PEs to: ", i
662 call mom_error(note, mesg)
663 endif
664 exit outer
665 endif
666 endif
667
668 ! Do not reduce the division count too much to satisfy the target_io_pes. Instead, re-do the iteration
669 ! with a reduced target_io_pes.
670 if (p <= max_feasible_p * (1.0 - pfrac)) exit ! (inner do loop)
671 enddo
672 enddo outer
673
674 if (num_masked_blocks == 0) then
675 call mom_error(fatal, "Couldn't auto-eliminate any land blocks. Try to increase the number "//&
676 "of MOM6 PEs or set AUTO_MASKTABLE to False.")
677 endif
678
679 ! Call determine_land_blocks once again, this time to retrieve and write out the mask_table.
680 allocate(mask_table(num_masked_blocks,2))
681 call determine_land_blocks(mask, nx, ny, layout(1), layout(2), ibuf, jbuf, num_masked_blocks, mask_table)
682 call write_auto_mask_file(mask_table, layout, npes, filename)
683 deallocate(mask_table)
684 deallocate(mask)
685
686end subroutine gen_auto_mask_table
687
688!> Given a number of domain divisions, compute the max number of land blocks that can be eliminated,
689!! and return the resulting mask table if requested.
690subroutine determine_land_blocks(mask, nx, ny, idiv, jdiv, ibuf, jbuf, num_masked_blocks, mask_table)
691 integer, dimension(:,:), intent(in) :: mask !< cell masks based on depth and MINIMUM_DEPTH
692 integer, intent(in) :: nx !< Total number of gridpoints in x-dir (global)
693 integer, intent(in) :: ny !< Total number of gridpoints in y-dir (global)
694 integer, intent(in) :: idiv !< number of divisions along x-dir
695 integer, intent(in) :: jdiv !< number of divisions along y-dir
696 integer, intent(in) :: ibuf !< number of buffer cells in x-dir.
697 !! (not necessarily the same as NIHALO)
698 integer, intent(in) :: jbuf !< number of buffer cells in y-dir.
699 !! (not necessarily the same as NJHALO)
700 integer, intent(out) :: num_masked_blocks !< the final number of masked blocks
701 integer, intent(out), optional :: mask_table(:,:) !< the resulting array of mask_table
702 ! integer
703 integer, dimension(idiv) :: ibegin !< The starting index of each division along x axis
704 integer, dimension(idiv) :: iend !< The ending index of each division along x axis
705 integer, dimension(jdiv) :: jbegin !< The starting index of each division along y axis
706 integer, dimension(jdiv) :: jend !< The ending index of each division along y axis
707 integer :: i, j, ib, ie, jb,je
708
709 call compute_extent(1, nx, idiv, ibegin, iend)
710 call compute_extent(1, ny, jdiv, jbegin, jend)
711
712 num_masked_blocks = 0
713
714 do i=1,idiv
715 ib = ibegin(i)
716 ie = iend(i) + 2 * ibuf
717 do j=1,jdiv
718 jb = jbegin(j)
719 je = jend(j) + 2 * jbuf
720
721 if (any(mask(ib:ie,jb:je)==1)) cycle
722
723 num_masked_blocks = num_masked_blocks + 1
724
725 if (present(mask_table)) then
726 if ( num_masked_blocks > size(mask_table, dim=1)) then
727 call mom_error(fatal, "The mask_table argument passed to determine_land_blocks() has insufficient size.")
728 endif
729
730 mask_table(num_masked_blocks,1) = i
731 mask_table(num_masked_blocks,2) = j
732 endif
733 enddo
734 enddo
735
736end subroutine determine_land_blocks
737
738!> Write out the auto-generated mask information to a file in the run directory.
739subroutine write_auto_mask_file(mask_table, layout, npes, filename)
740 integer, intent(in) :: mask_table(:,:) !> mask table array to be written out.
741 integer, dimension(2), intent(in) :: layout !> PE layout
742 integer, intent(in) :: npes !> Number of divisions (incl. eliminated ones)
743 character(len=:), allocatable, intent(in) :: filename !> file name for the mask_table to be written
744 ! local
745 integer :: file_ascii= -1 !< The unit number of the auto-generated mask_file file.
746 integer :: true_num_masked_blocks
747 integer :: p
748
749 ! Eliminate only enough blocks to ensure that the number of active blocks precisely matches the target npes.
750 true_num_masked_blocks = layout(1) * layout(2) - npes
751
752 call open_ascii_file(file_ascii, trim(filename), action=writeonly_file)
753 write(file_ascii, '(I0)') true_num_masked_blocks
754 write(file_ascii, '(I0,",",I0)') layout(1), layout(2)
755 do p = 1, true_num_masked_blocks
756 write(file_ascii, '(I0,",",I0)') mask_table(p,1), mask_table(p,2)
757 enddo
758 call close_file(file_ascii)
759end subroutine write_auto_mask_file
760
761!> Computes the io layout based on the domain layout and a target number of IO PEs.
762function auto_determine_io_layout(idiv, jdiv, nio) result(best_io_layout)
763 integer, intent(in) :: idiv ! The number of compute divisions along the x-axis
764 integer, intent(in) :: jdiv ! The number of compute divisions along the y-axis
765 integer, intent(in) :: nio ! Target number of IO PEs, s.t., (idiv_io * jdiv_io) % nio == 0
766 ! return
767 integer :: io_layout(2) ! Temporary IO layout
768 integer :: best_io_layout(2) ! Best IO layout
769 integer :: f, best_idiv_io, best_jdiv_io
770 real :: ratio_diff, min_ratio_diff
771
772 if (mod(idiv*jdiv, nio) /= 0) then
773 call mom_error(fatal, "The product of the compute layout must be divisible by the target number of IO PEs.")
774 endif
775
776 min_ratio_diff = 1.0e30 ! Large initial value
777 best_io_layout = [1, nio]
778
779 ! Iterate over all factors of nio
780 do f = 1, nio
781 if (mod(nio, f) /= 0) cycle
782 io_layout = [f, nio / f]
783
784 ! Check divisibility constraints
785 if (mod(idiv, io_layout(1)) == 0 .and. mod(jdiv, io_layout(2)) == 0) then
786 ratio_diff = abs(real(io_layout(1)) / real(io_layout(2)) - real(idiv) / real(jdiv))
787
788 ! Update best choice if ratio_diff is smaller
789 if (ratio_diff < min_ratio_diff) then
790 min_ratio_diff = ratio_diff
791 best_io_layout = [io_layout(1), io_layout(2)]
792 endif
793 endif
794 enddo
795
796end function auto_determine_io_layout
797
798end module mom_domains