source: CPL/oasis3/trunk/src/lib/mpp_io/src/mpp_domains_mod.F90 @ 1677

Last change on this file since 1677 was 1677, checked in by aclsce, 12 years ago

Imported oasis3 (tag ipslcm5a) from cvs server to svn server (igcmg project).

File size: 120.7 KB
Line 
1!-----------------------------------------------------------------------
2!   Domain decomposition and domain update for message-passing codes
3!
4! AUTHOR: V. Balaji (vbalaji@noaa.gov)
5!         Princeton University/GFDL
6!
7! MODIFICATIONS: Reiner Vogelsang (reiner@sgi.com)
8!
9! This program is free software; The author agrees that you can
10! redistribute and/or modify this version of the program under the
11! terms of the Lesser GNU General Public License as published
12! by the Free Software Foundation.
13!
14! This program is distributed in the hope that it will be useful,
15! but WITHOUT ANY WARRANTY; without even the implied warranty of
16! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17! Lesser GNU General Public License for more details
18! (http://www.gnu.org/copyleft/lesser.html).
19!-----------------------------------------------------------------------
20!#include <os.h>
21#if defined(_CRAYT3E) || defined(_CRAYT3D) || defined(sgi_mipspro)
22#define SGICRAY_MPP
23#endif
24!shmalloc is used on MPP SGI/Cray systems for shmem
25#if defined(use_libSMA) && defined(SGICRAY_MPP)
26#define use_shmalloc
27#endif
28
29module mpp_domains_mod
30  use mod_kinds_model
31!a generalized domain decomposition package for use with mpp_mod
32!Balaji (vb@gfdl.gov) 15 March 1999
33  use mpp_mod
34#include <os.h>
35  implicit none
36  private
37  character(len=128), private :: version= &
38       '$Id: mpp_domains_mod.F90,v 1.1.1.1 2005/03/23 16:01:11 adm Exp $'
39  character(len=128), private :: tagname= &
40       '$Name: ipslcm5a $'
41  character(len=128), private :: version_update_domains2D, version_global_reduce, version_global_sum, version_global_field
42
43!parameters used to define domains: these are passed to the flags argument of mpp_define_domains
44!  if data domain is to have global extent, set GLOBAL_DATA_DOMAIN
45!  if global domain has periodic boundaries, set CYCLIC_GLOBAL_DOMAIN
46!  sum flags together if more than one of the above conditions is to be met.
47  integer, parameter, private :: GLOBAL=0, CYCLIC=1
48  integer, parameter, private :: WEST=2, EAST=3, SOUTH=4, NORTH=5
49  integer, parameter, private :: SEND=1, RECV=2
50  integer, parameter, public :: GLOBAL_DATA_DOMAIN=2**GLOBAL, CYCLIC_GLOBAL_DOMAIN=2**CYCLIC
51!gridtypes
52  integer, parameter, private :: AGRID=0, BGRID=1, CGRID=2
53  integer, parameter, public :: BGRID_NE=BGRID+2**NORTH+2**EAST
54  integer, parameter, public :: BGRID_SW=BGRID+2**SOUTH+2**WEST
55  integer, parameter, public :: CGRID_NE=CGRID+2**NORTH+2**EAST
56  integer, parameter, public :: CGRID_SW=CGRID+2**SOUTH+2**WEST
57  integer, private :: grid_offset_type=AGRID
58!folds
59  integer, parameter, public :: FOLD_WEST_EDGE = 2**WEST, FOLD_EAST_EDGE = 2**EAST
60  integer, parameter, public :: FOLD_SOUTH_EDGE=2**SOUTH, FOLD_NORTH_EDGE=2**NORTH
61!update
62  integer, parameter, public :: WUPDATE=2**WEST, EUPDATE=2**EAST, SUPDATE=2**SOUTH, NUPDATE=2**NORTH
63  integer, parameter, public :: XUPDATE=WUPDATE+EUPDATE, YUPDATE=SUPDATE+NUPDATE
64!used by mpp_global_sum
65  integer, parameter, public :: BITWISE_EXACT_SUM=1
66
67  type, public :: domain_axis_spec        !type used to specify index limits along an axis of a domain
68     private
69     integer :: begin, end, size, max_size      !start, end of domain axis, size, max size in set
70     logical :: is_global       !TRUE if domain axis extent covers global domain
71  end type domain_axis_spec
72  type, public :: domain1D
73     private
74     type(domain_axis_spec) :: compute, data, global
75     logical :: cyclic
76     type(domain1D), pointer :: list(:)
77     integer :: pe              !PE to which this domain is assigned
78     integer :: pos             !position of this PE within link list, i.e domain%list(pos)%pe = pe
79  end type domain1D
80!domaintypes of higher rank can be constructed from type domain1D
81!typically we only need 1 and 2D, but could need higher (e.g 3D LES)
82!some elements are repeated below if they are needed once per domain, not once per axis
83  type, private :: rectangle
84     integer :: is, ie, js, je
85     logical :: overlap, folded
86  end type rectangle
87  type, public :: domain2D
88     private
89     type(domain1D) :: x
90     type(domain1D) :: y
91     type(domain2D), pointer :: list(:)
92     integer :: pe              !PE to which this domain is assigned
93     integer :: pos             !position of this PE within link list, i.e domain%list(pos)%pe = pe
94     integer :: fold, gridtype
95     logical :: overlap
96     type(rectangle) :: recv_e, recv_se, recv_s, recv_sw, &
97                        recv_w, recv_nw, recv_n, recv_ne
98     type(rectangle) :: send_e, send_se, send_s, send_sw, &
99                        send_w, send_nw, send_n, send_ne
100     logical :: remote_domains_initialized
101     type(rectangle) :: recv_e_off, recv_se_off, recv_s_off, recv_sw_off, &
102                        recv_w_off, recv_nw_off, recv_n_off, recv_ne_off
103     type(rectangle) :: send_e_off, send_se_off, send_s_off, send_sw_off, &
104                        send_w_off, send_nw_off, send_n_off, send_ne_off
105     logical :: remote_off_domains_initialized
106  end type domain2D
107  type(domain1D), public :: NULL_DOMAIN1D
108  type(domain2D), public :: NULL_DOMAIN2D
109
110  integer, private :: pe
111
112  integer, private :: tk
113  logical, private :: verbose=.FALSE., debug=.FALSE., domain_clocks_on=.FALSE.
114  logical, private :: module_is_initialized=.FALSE.
115  integer, parameter, public :: MPP_DOMAIN_TIME=MPP_DEBUG+1
116  integer :: send_clock=0, recv_clock=0, unpk_clock=0, wait_clock=0, pack_clock=0, pack_loop_clock=0
117
118!stack for internal buffers
119!allocated differently if use_shmalloc
120#ifdef use_shmalloc
121  real(DOUBLE_KIND), private :: mpp_domains_stack(1)
122  pointer( ptr_stack, mpp_domains_stack )
123#else
124  real(DOUBLE_KIND), private, allocatable :: mpp_domains_stack(:)
125#endif
126  integer, private :: mpp_domains_stack_size=0, mpp_domains_stack_hwm=0
127
128!used by mpp_define_domains2D_new to transmit data
129  integer :: domain_info_buf(16)
130#ifdef use_shmalloc
131  pointer( ptr_info, domain_info_buf )
132#endif
133
134!public interfaces
135
136! this interface can add halo to a no-halo domain and copy everything else.
137
138  interface mpp_copy_domains
139     module procedure mpp_copy_domains1D
140     module procedure mpp_copy_domains2D
141  end interface
142
143  interface mpp_define_domains
144     module procedure mpp_define_domains1D
145     module procedure mpp_define_domains2D
146  end interface
147
148  interface mpp_update_domains
149     module procedure mpp_update_domain2D_r8_2d
150     module procedure mpp_update_domain2D_r8_3d
151     module procedure mpp_update_domain2D_r8_4d
152     module procedure mpp_update_domain2D_r8_5d
153     module procedure mpp_update_domain2D_r8_2dv
154     module procedure mpp_update_domain2D_r8_3dv
155     module procedure mpp_update_domain2D_r8_4dv
156     module procedure mpp_update_domain2D_r8_5dv
157     module procedure mpp_update_domain2D_c8_2d
158     module procedure mpp_update_domain2D_c8_3d
159     module procedure mpp_update_domain2D_c8_4d
160     module procedure mpp_update_domain2D_c8_5d
161#ifndef no_8byte_integers
162     module procedure mpp_update_domain2D_i8_2d
163     module procedure mpp_update_domain2D_i8_3d
164     module procedure mpp_update_domain2D_i8_4d
165     module procedure mpp_update_domain2D_i8_5d
166     module procedure mpp_update_domain2D_l8_2d
167     module procedure mpp_update_domain2D_l8_3d
168     module procedure mpp_update_domain2D_l8_4d
169     module procedure mpp_update_domain2D_l8_5d
170#endif
171#ifndef no_4byte_reals
172     module procedure mpp_update_domain2D_r4_2d
173     module procedure mpp_update_domain2D_r4_3d
174     module procedure mpp_update_domain2D_r4_4d
175     module procedure mpp_update_domain2D_r4_5d
176#endif
177#ifndef no_4byte_cmplx
178     module procedure mpp_update_domain2D_c4_2d
179     module procedure mpp_update_domain2D_c4_3d
180     module procedure mpp_update_domain2D_c4_4d
181     module procedure mpp_update_domain2D_c4_5d
182#endif
183     module procedure mpp_update_domain2D_i4_2d
184     module procedure mpp_update_domain2D_i4_3d
185     module procedure mpp_update_domain2D_i4_4d
186     module procedure mpp_update_domain2D_i4_5d
187     module procedure mpp_update_domain2D_l4_2d
188     module procedure mpp_update_domain2D_l4_3d
189     module procedure mpp_update_domain2D_l4_4d
190     module procedure mpp_update_domain2D_l4_5d
191#ifndef no_4byte_reals
192     module procedure mpp_update_domain2D_r4_2dv
193     module procedure mpp_update_domain2D_r4_3dv
194     module procedure mpp_update_domain2D_r4_4dv
195     module procedure mpp_update_domain2D_r4_5dv
196#endif
197
198!     module procedure mpp_update_domain1D_r8_2d
199!     module procedure mpp_update_domain1D_r8_3d
200!     module procedure mpp_update_domain1D_r8_4d
201!     module procedure mpp_update_domain1D_r8_5d
202!     module procedure mpp_update_domain1D_c8_2d
203!     module procedure mpp_update_domain1D_c8_3d
204!     module procedure mpp_update_domain1D_c8_4d
205!     module procedure mpp_update_domain1D_c8_5d
206!#ifndef no_8byte_integers
207!     module procedure mpp_update_domain1D_i8_2d
208!     module procedure mpp_update_domain1D_i8_3d
209!     module procedure mpp_update_domain1D_i8_4d
210!     module procedure mpp_update_domain1D_i8_5d
211!     module procedure mpp_update_domain1D_l8_2d
212!     module procedure mpp_update_domain1D_l8_3d
213!     module procedure mpp_update_domain1D_l8_4d
214!     module procedure mpp_update_domain1D_l8_5d
215!#endif
216!     module procedure mpp_update_domain1D_r4_2d
217!     module procedure mpp_update_domain1D_r4_3d
218!     module procedure mpp_update_domain1D_r4_4d
219!     module procedure mpp_update_domain1D_r4_5d
220!     module procedure mpp_update_domain1D_c4_2d
221!     module procedure mpp_update_domain1D_c4_3d
222!     module procedure mpp_update_domain1D_c4_4d
223!     module procedure mpp_update_domain1D_c4_5d
224!     module procedure mpp_update_domain1D_i4_2d
225!     module procedure mpp_update_domain1D_i4_3d
226!     module procedure mpp_update_domain1D_i4_4d
227!     module procedure mpp_update_domain1D_i4_5d
228!     module procedure mpp_update_domain1D_l4_2d
229!     module procedure mpp_update_domain1D_l4_3d
230!     module procedure mpp_update_domain1D_l4_4d
231!     module procedure mpp_update_domain1D_l4_5d
232  end interface
233
234  interface mpp_redistribute
235     module procedure mpp_redistribute_r8_2D
236     module procedure mpp_redistribute_r8_3D
237     module procedure mpp_redistribute_r8_4D
238     module procedure mpp_redistribute_r8_5D
239     module procedure mpp_redistribute_c8_2D
240     module procedure mpp_redistribute_c8_3D
241     module procedure mpp_redistribute_c8_4D
242     module procedure mpp_redistribute_c8_5D
243#ifndef no_8byte_integers
244     module procedure mpp_redistribute_i8_2D
245     module procedure mpp_redistribute_i8_3D
246     module procedure mpp_redistribute_i8_4D
247     module procedure mpp_redistribute_i8_5D
248     module procedure mpp_redistribute_l8_2D
249     module procedure mpp_redistribute_l8_3D
250     module procedure mpp_redistribute_l8_4D
251     module procedure mpp_redistribute_l8_5D
252#endif
253#ifndef no_4byte_reals
254     module procedure mpp_redistribute_r4_2D
255     module procedure mpp_redistribute_r4_3D
256     module procedure mpp_redistribute_r4_4D
257     module procedure mpp_redistribute_r4_5D
258#endif
259#ifndef no_4byte_cmplx
260     module procedure mpp_redistribute_c4_2D
261     module procedure mpp_redistribute_c4_3D
262     module procedure mpp_redistribute_c4_4D
263     module procedure mpp_redistribute_c4_5D
264#endif
265     module procedure mpp_redistribute_i4_2D
266     module procedure mpp_redistribute_i4_3D
267     module procedure mpp_redistribute_i4_4D
268     module procedure mpp_redistribute_i4_5D
269     module procedure mpp_redistribute_l4_2D
270     module procedure mpp_redistribute_l4_3D
271     module procedure mpp_redistribute_l4_4D
272     module procedure mpp_redistribute_l4_5D
273  end interface
274
275  interface mpp_global_field
276     module procedure mpp_global_field2D_r8_2d
277     module procedure mpp_global_field2D_r8_3d
278     module procedure mpp_global_field2D_r8_4d
279     module procedure mpp_global_field2D_r8_5d
280     module procedure mpp_global_field2D_c8_2d
281     module procedure mpp_global_field2D_c8_3d
282     module procedure mpp_global_field2D_c8_4d
283     module procedure mpp_global_field2D_c8_5d
284#ifndef no_8byte_integers
285     module procedure mpp_global_field2D_i8_2d
286     module procedure mpp_global_field2D_i8_3d
287     module procedure mpp_global_field2D_i8_4d
288     module procedure mpp_global_field2D_i8_5d
289     module procedure mpp_global_field2D_l8_2d
290     module procedure mpp_global_field2D_l8_3d
291     module procedure mpp_global_field2D_l8_4d
292     module procedure mpp_global_field2D_l8_5d
293#endif
294#ifndef no_4byte_reals
295     module procedure mpp_global_field2D_r4_2d
296     module procedure mpp_global_field2D_r4_3d
297     module procedure mpp_global_field2D_r4_4d
298     module procedure mpp_global_field2D_r4_5d
299#endif
300#ifndef no_4byte_cmplx
301     module procedure mpp_global_field2D_c4_2d
302     module procedure mpp_global_field2D_c4_3d
303     module procedure mpp_global_field2D_c4_4d
304     module procedure mpp_global_field2D_c4_5d
305#endif
306     module procedure mpp_global_field2D_i4_2d
307     module procedure mpp_global_field2D_i4_3d
308     module procedure mpp_global_field2D_i4_4d
309     module procedure mpp_global_field2D_i4_5d
310     module procedure mpp_global_field2D_l4_2d
311     module procedure mpp_global_field2D_l4_3d
312     module procedure mpp_global_field2D_l4_4d
313     module procedure mpp_global_field2D_l4_5d
314     module procedure mpp_global_field1D_r8_2d
315     module procedure mpp_global_field1D_c8_2d
316#ifndef no_8byte_integers
317     module procedure mpp_global_field1D_i8_2d
318     module procedure mpp_global_field1D_l8_2d
319#endif
320#ifndef no_4byte_reals
321     module procedure mpp_global_field1D_r4_2d
322#endif
323#ifndef no_4byte_cmplx
324     module procedure mpp_global_field1D_c4_2d
325#endif
326     module procedure mpp_global_field1D_i4_2d
327     module procedure mpp_global_field1D_l4_2d
328  end interface
329
330  interface mpp_global_max
331     module procedure mpp_global_max_r8_2d
332     module procedure mpp_global_max_r8_3d
333     module procedure mpp_global_max_r8_4d
334     module procedure mpp_global_max_r8_5d
335#ifndef no_4byte_reals
336     module procedure mpp_global_max_r4_2d
337     module procedure mpp_global_max_r4_3d
338     module procedure mpp_global_max_r4_4d
339     module procedure mpp_global_max_r4_5d
340#endif
341#ifndef no_8byte_integers
342     module procedure mpp_global_max_i8_2d
343     module procedure mpp_global_max_i8_3d
344     module procedure mpp_global_max_i8_4d
345     module procedure mpp_global_max_i8_5d
346#endif
347     module procedure mpp_global_max_i4_2d
348     module procedure mpp_global_max_i4_3d
349     module procedure mpp_global_max_i4_4d
350     module procedure mpp_global_max_i4_5d
351  end interface
352
353  interface mpp_global_min
354     module procedure mpp_global_min_r8_2d
355     module procedure mpp_global_min_r8_3d
356     module procedure mpp_global_min_r8_4d
357     module procedure mpp_global_min_r8_5d
358#ifndef no_4byte_reals
359     module procedure mpp_global_min_r4_2d
360     module procedure mpp_global_min_r4_3d
361     module procedure mpp_global_min_r4_4d
362     module procedure mpp_global_min_r4_5d
363#endif
364#ifndef no_8byte_integers
365     module procedure mpp_global_min_i8_2d
366     module procedure mpp_global_min_i8_3d
367     module procedure mpp_global_min_i8_4d
368     module procedure mpp_global_min_i8_5d
369#endif
370     module procedure mpp_global_min_i4_2d
371     module procedure mpp_global_min_i4_3d
372     module procedure mpp_global_min_i4_4d
373     module procedure mpp_global_min_i4_5d
374  end interface
375
376  interface mpp_global_sum
377     module procedure mpp_global_sum_r8_2d
378     module procedure mpp_global_sum_r8_3d
379     module procedure mpp_global_sum_r8_4d
380     module procedure mpp_global_sum_r8_5d
381#ifndef no_4byte_reals
382     module procedure mpp_global_sum_r4_2d
383     module procedure mpp_global_sum_r4_3d
384     module procedure mpp_global_sum_r4_4d
385     module procedure mpp_global_sum_r4_5d
386#endif
387     module procedure mpp_global_sum_c8_2d
388     module procedure mpp_global_sum_c8_3d
389     module procedure mpp_global_sum_c8_4d
390     module procedure mpp_global_sum_c8_5d
391#ifndef no_4byte_cmplx
392     module procedure mpp_global_sum_c4_2d
393     module procedure mpp_global_sum_c4_3d
394     module procedure mpp_global_sum_c4_4d
395     module procedure mpp_global_sum_c4_5d
396#endif
397#ifndef no_8byte_integers
398     module procedure mpp_global_sum_i8_2d
399     module procedure mpp_global_sum_i8_3d
400     module procedure mpp_global_sum_i8_4d
401     module procedure mpp_global_sum_i8_5d
402#endif
403     module procedure mpp_global_sum_i4_2d
404     module procedure mpp_global_sum_i4_3d
405     module procedure mpp_global_sum_i4_4d
406     module procedure mpp_global_sum_i4_5d
407  end interface
408
409  interface operator(.EQ.)
410     module procedure mpp_domain1D_eq
411     module procedure mpp_domain2D_eq
412  end interface
413
414  interface operator(.NE.)
415     module procedure mpp_domain1D_ne
416     module procedure mpp_domain2D_ne
417  end interface
418
419  interface mpp_get_compute_domain
420     module procedure mpp_get_compute_domain1D
421     module procedure mpp_get_compute_domain2D
422  end interface
423
424  interface mpp_get_compute_domains
425     module procedure mpp_get_compute_domains1D
426     module procedure mpp_get_compute_domains2D
427  end interface
428
429  interface mpp_get_data_domain
430     module procedure mpp_get_data_domain1D
431     module procedure mpp_get_data_domain2D
432  end interface
433
434  interface mpp_get_global_domain
435     module procedure mpp_get_global_domain1D
436     module procedure mpp_get_global_domain2D
437  end interface
438
439  interface mpp_define_layout
440     module procedure mpp_define_layout2D
441  end interface
442
443  interface mpp_get_pelist
444     module procedure mpp_get_pelist1D
445     module procedure mpp_get_pelist2D
446  end interface
447
448  interface mpp_get_layout
449     module procedure mpp_get_layout1D
450     module procedure mpp_get_layout2D
451  end interface
452
453  public :: mpp_broadcast_domain, mpp_define_layout, mpp_define_domains, mpp_domains_init, mpp_domains_set_stack_size, &
454            mpp_domains_exit, mpp_get_compute_domain, mpp_get_compute_domains, mpp_get_data_domain, mpp_get_global_domain, &
455            mpp_get_domain_components, mpp_get_layout, mpp_get_pelist, mpp_redistribute, mpp_update_domains, &
456            mpp_global_field, mpp_global_max, mpp_global_min, mpp_global_sum, operator(.EQ.), operator(.NE.), &
457            mpp_copy_domains
458
459  contains
460
461!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
462!                                                                             !
463!                 MPP_DOMAINS: initialization and termination                 !
464!                                                                             !
465!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
466
467    subroutine mpp_domains_init(flags)
468!initialize domain decomp package
469      integer, intent(in), optional :: flags
470      integer :: l=0
471
472      if( module_is_initialized )return
473      call mpp_init(flags)           !this is a no-op if already initialized
474      pe = mpp_pe()
475      module_is_initialized = .TRUE.
476      if( pe.EQ.mpp_root_pe() )then
477          write( stdlog(),'(/a)' )'MPP_DOMAINS module '//trim(version) &
478               //trim(tagname)
479!          write( stdlog(),'(a)' )trim(version_update_domains2D)
480      end if
481
482      if( PRESENT(flags) )then
483          debug   = flags.EQ.MPP_DEBUG
484          verbose = flags.EQ.MPP_VERBOSE .OR. debug
485          domain_clocks_on = flags.EQ.MPP_DOMAIN_TIME
486      end if
487
488      call mpp_domains_set_stack_size(983040) !default, pretty arbitrary
489#ifdef use_shmalloc
490      call mpp_malloc( ptr_info, 16, l )
491#endif
492
493!NULL_DOMAIN is a domaintype that can be used to initialize to undef
494      NULL_DOMAIN1D%global%begin  = -1; NULL_DOMAIN1D%global%end  = -1; NULL_DOMAIN1D%global%size = 0
495      NULL_DOMAIN1D%data%begin    = -1; NULL_DOMAIN1D%data%end    = -1; NULL_DOMAIN1D%data%size = 0
496      NULL_DOMAIN1D%compute%begin = -1; NULL_DOMAIN1D%compute%end = -1; NULL_DOMAIN1D%compute%size = 0
497      NULL_DOMAIN1D%pe = ANY_PE
498      NULL_DOMAIN2D%x = NULL_DOMAIN1D
499      NULL_DOMAIN2D%y = NULL_DOMAIN1D
500      NULL_DOMAIN2D%pe = ANY_PE
501
502      if( domain_clocks_on )then
503          pack_clock = mpp_clock_id( 'Halo pack' )
504          pack_loop_clock = mpp_clock_id( 'Halo pack loop' )
505          send_clock = mpp_clock_id( 'Halo send' )
506          recv_clock = mpp_clock_id( 'Halo recv' )
507          unpk_clock = mpp_clock_id( 'Halo unpk' )
508          wait_clock = mpp_clock_id( 'Halo wait' )
509      end if
510      return
511    end subroutine mpp_domains_init
512
513    subroutine mpp_domains_set_stack_size(n)
514!set the mpp_domains_stack variable to be at least n LONG words long
515      integer, intent(in) :: n
516      character(len=8) :: text
517
518      if( n.LE.mpp_domains_stack_size )return
519#ifdef use_shmalloc
520      call mpp_malloc( ptr_stack, n, mpp_domains_stack_size )
521#else
522      if( allocated(mpp_domains_stack) )deallocate(mpp_domains_stack)
523      allocate( mpp_domains_stack(n) )
524      mpp_domains_stack_size = n
525#endif
526      write( text,'(i8)' )n
527      if( pe.EQ.mpp_root_pe() )call mpp_error( NOTE, 'MPP_DOMAINS_SET_STACK_SIZE: stack size set to '//text//'.' )
528
529      return
530    end subroutine mpp_domains_set_stack_size
531
532    subroutine mpp_domains_exit()
533!currently does not have much to do, but provides the possibility of re-initialization
534      if( .NOT.module_is_initialized )return
535      call mpp_max(mpp_domains_stack_hwm)
536      if( pe.EQ.mpp_root_pe() )write( stdout(),* )'MPP_DOMAINS_STACK high water mark=', mpp_domains_stack_hwm
537      module_is_initialized = .FALSE.
538      return
539    end subroutine mpp_domains_exit
540
541!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
542!                                                                             !
543!                MPP_DOMAINS: overloaded operators (==, /=)                   !
544!                                                                             !
545!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
546
547    function mpp_domain1D_eq( a, b )
548      logical :: mpp_domain1D_eq
549      type(domain1D), intent(in) :: a, b
550
551      mpp_domain1D_eq = ( a%compute%begin.EQ.b%compute%begin .AND. &
552                          a%compute%end  .EQ.b%compute%end   .AND. &
553                          a%data%begin   .EQ.b%data%begin    .AND. &
554                          a%data%end     .EQ.b%data%end      .AND. & 
555                          a%global%begin .EQ.b%global%begin  .AND. &
556                          a%global%end   .EQ.b%global%end    )
557!compare pelists
558!      if( mpp_domain1D_eq )mpp_domain1D_eq = ASSOCIATED(a%list) .AND. ASSOCIATED(b%list)
559!      if( mpp_domain1D_eq )mpp_domain1D_eq = size(a%list).EQ.size(b%list)
560!      if( mpp_domain1D_eq )mpp_domain1D_eq = ALL(a%list%pe.EQ.b%list%pe)
561     
562      return
563    end function mpp_domain1D_eq
564
565    function mpp_domain1D_ne( a, b )
566      logical :: mpp_domain1D_ne
567      type(domain1D), intent(in) :: a, b
568
569      mpp_domain1D_ne = .NOT. ( a.EQ.b )
570      return
571    end function mpp_domain1D_ne
572
573    function mpp_domain2D_eq( a, b )
574      logical :: mpp_domain2D_eq
575      type(domain2D), intent(in) :: a, b
576
577      mpp_domain2D_eq = a%x.EQ.b%x .AND. a%y.EQ.b%y
578      if( mpp_domain2D_eq .AND. ((a%pe.EQ.ANY_PE).OR.(b%pe.EQ.ANY_PE)) )return !NULL_DOMAIN2D
579!compare pelists
580      if( mpp_domain2D_eq )mpp_domain2D_eq = ASSOCIATED(a%list) .AND. ASSOCIATED(b%list)
581      if( mpp_domain2D_eq )mpp_domain2D_eq = size(a%list).EQ.size(b%list)
582      if( mpp_domain2D_eq )mpp_domain2D_eq = ALL(a%list%pe.EQ.b%list%pe)
583      return
584    end function mpp_domain2D_eq
585
586    function mpp_domain2D_ne( a, b )
587      logical :: mpp_domain2D_ne
588      type(domain2D), intent(in) :: a, b
589
590      mpp_domain2D_ne = .NOT. ( a.EQ.b )
591      return
592    end function mpp_domain2D_ne
593
594
595!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
596!                                                                             !
597!              MPP_COPY_DOMAINS: Copy domain and add halo                     !
598!                                                                             !
599!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
600
601    subroutine mpp_copy_domains1D(domain_in, domain_out, halo)
602! ARGUMENTS
603! domain_in is the input domain
604! domain_out is the returned domain
605! halo (optional) defines halo width (currently the same on both sides)
606
607      type(domain1D), intent(in)    :: domain_in
608      type(domain1D), intent(out)   :: domain_out
609      integer, intent(in), optional :: halo
610
611      integer, dimension(:), allocatable :: pelist, extent
612      integer :: ndivs, global_indices(2) !(/ isg, ieg /)
613      integer :: isg, ieg, halosz
614
615     
616! get the global indices of the input domain
617      call mpp_get_global_domain(domain_in, isg, ieg )
618      global_indices(1) = isg;       global_indices(2) = ieg
619
620! get the number of divisions of domain
621      call mpp_get_layout(domain_in, ndivs)
622
623      allocate(pelist(0:ndivs-1),extent(0:ndivs-1))
624
625! get the pe list of the input domain
626      call mpp_get_pelist( domain_in, pelist)
627
628! get the halo
629      halosz = 0
630      if(present(halo)) halosz = halo
631
632! get the extent
633      call mpp_get_compute_domains(Domain_in, size = extent(0:ndivs-1))     
634
635      call mpp_define_domains( global_indices, ndivs, domain_out, pelist = pelist, &
636                                 halo = halosz, extent = extent )
637
638    end subroutine mpp_copy_domains1D
639
640
641!----------------------------------------------------------------------------------
642
643
644    subroutine mpp_copy_domains2D(domain_in, domain_out, xhalo, yhalo)
645! ARGUMENTS
646! domain_in is the input domain
647! domain_out is the returned domain
648! xhalo, yhalo (optional) defines halo width.
649
650      type(domain2D), intent(in)    :: domain_in
651      type(domain2D), intent(out)   :: domain_out
652      integer, intent(in), optional :: xhalo, yhalo
653
654      integer, dimension(:), allocatable :: pelist, xextent, yextent, xbegin, xend, &
655                                            xsize, ybegin, yend, ysize
656      integer :: ndivx, ndivy, ndivs, npes, global_indices(4), layout(2)
657      integer :: isg, ieg, jsg, jeg, xhalosz, yhalosz, i, j, m, n, pe
658     
659! get the global indices of the input domain
660      call mpp_get_global_domain(domain_in, isg, ieg, jsg, jeg )
661      global_indices(1) = isg;       global_indices(2) = ieg
662      global_indices(3) = jsg;       global_indices(4) = jeg
663
664! get the number of divisions of domain
665      call mpp_get_layout(domain_in, layout)
666      ndivx = layout(1); ndivy = layout(2)
667      ndivs = ndivx * ndivy
668      npes = mpp_npes()
669
670      if( ndivs.NE.npes )call mpp_error( 'mpp_domains_mod', &
671                 'mpp_copy_domains: number of PEs is not consistent with the layout', FATAL )
672      allocate(pelist(0:npes-1), xsize(0:npes-1), ysize(0:npes-1), &
673               xbegin(0:npes-1), xend(0:npes-1),                   &
674               ybegin(0:npes-1), yend(0:npes-1),                   &
675               xextent(0:ndivx), yextent(0:ndivy))
676
677! get the pe list of the input domain
678      call mpp_get_pelist2D( domain_in, pelist)
679
680! get the halo
681      xhalosz = 0; yhalosz = 0
682      if(present(xhalo)) xhalosz = xhalo
683      if(present(yhalo)) yhalosz = yhalo
684
685! get the extent
686      call mpp_get_compute_domains(Domain_in, xbegin, xend, xsize, &
687                                          ybegin, yend, ysize  )
688           
689!  compute the exact x-axis decomposition
690      i = isg
691      m = 0
692      xextent = 0
693      do pe = 0, npes-1
694        if ( xbegin(pe) == i ) then
695          m = m+1
696          xextent (m) = xsize (pe)
697          i = i + xsize(pe)
698          if ( m == size(xextent) .or. i > ieg ) exit
699        endif
700      enddo
701
702!  compute the exact y-axis decomposition
703      j = jsg 
704      n = 0   
705      yextent = 0
706      do pe = 0, npes-1
707        if ( ybegin(pe) == j ) then
708          n = n+1 
709          yextent (n) = ysize (pe)
710          j = j + ysize(pe)
711          if ( n == size(yextent) .or. j > jeg ) exit
712        endif   
713      enddo   
714
715      call mpp_define_domains( global_indices, layout, domain_out, pelist = pelist, &
716                               xhalo = xhalosz, yhalo = yhalosz, xextent = xextent, & 
717                               yextent = yextent)
718
719    end subroutine mpp_copy_domains2D
720
721
722
723!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
724!                                                                             !
725!              MPP_DEFINE_DOMAINS: define layout and decomposition            !
726!                                                                             !
727!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
728
729    subroutine mpp_define_domains1D( global_indices, ndivs, domain, pelist, &
730                                     flags, halo, extent, maskmap,offset)
731!routine to divide global array indices among domains, and assign domains to PEs
732!domain is of type domain1D
733!ARGUMENTS:
734!      global_indices(2)=(isg,ieg) gives the extent of global domain
735!      ndivs is number of divisions of domain: even divisions unless extent is present.
736!      domain is the returned domain1D
737!      pelist (optional) list of PEs to which domains are to be assigned (default 0...npes-1)
738!                 size of pelist must correspond to number of mask=.TRUE. divisions
739!      flags define whether compute and data domains are global (undecomposed) and whether global domain has periodic boundaries
740!      halo (optional) defines halo width (currently the same on both sides)
741!      extent (optional) array defines width of each division (used for non-uniform domain decomp, for e.g load-balancing)
742!      maskmap (optional) a division whose maskmap=.FALSE. is not assigned to any domain
743!  By default we assume decomposition of compute and data domains, non-periodic boundaries, no halo, as close to uniform extents
744!  as the input parameters permit
745      integer, intent(in) :: global_indices(2) !(/ isg, ieg /)
746      integer, intent(in) :: ndivs
747      type(domain1D), intent(inout) :: domain !declared inout so that existing links, if any, can be nullified
748      integer, intent(in), optional :: pelist(0:)
749      integer, intent(in), optional :: flags, halo
750      integer, intent(in), optional :: extent(0:)
751      logical, intent(in), optional :: maskmap(0:)
752      integer, intent(in), optional :: offset(0:)
753
754      logical :: compute_domain_is_global, data_domain_is_global
755      integer :: ndiv, n, isg, ieg, is, ie, i, off, pos, hs, he
756      integer, allocatable :: pes(:)
757      logical, allocatable :: mask(:)
758!rv
759      logical ::blocks
760      integer :: lastie
761!rv
762      integer :: halosz
763!used by symmetry algorithm
764      integer :: imax, ndmax, ndmirror
765      logical :: symmetrize
766!statement functions
767      logical :: even, odd
768      even(n) = (mod(n,2).EQ.0)
769      odd (n) = (mod(n,2).EQ.1)
770     
771      if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS1D: You must first call mpp_domains_init.' )
772!get global indices
773      isg = global_indices(1)
774      ieg = global_indices(2)
775      if( ndivs.GT.ieg-isg+1 )call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS1D: more divisions requested than rows available.' )
776!get the list of PEs on which to assign domains; if pelist is absent use 0..npes-1
777      if( PRESENT(pelist) )then
778          if( .NOT.any(pelist.EQ.pe) )then
779              write( stderr(),* )'pe=', pe, ' pelist=', pelist
780              call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS1D: pe must be in pelist.' )
781          end if
782          allocate( pes(0:size(pelist)-1) )
783          pes(:) = pelist(:)
784      else
785          allocate( pes(0:mpp_npes()-1) )
786          pes(:) = (/ (i,i=0,mpp_npes()-1) /)
787      end if
788
789!rv Block distribution ?
790      blocks=.false.
791      if(present(offset)) blocks=.true.
792      if(blocks) then
793       if((.not.present(maskmap)).and.(.not.present(extent))) &
794         call mpp_error(FATAL, 'MPP_DEFINE_DOMAINS1D: you have to define maskmap and extent!')
795        if(size(offset).ne.ndivs) &
796          call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS1D: offset array size must equal number of domain divisions.' )
797      endif
798!rv
799!get number of real domains: 1 mask domain per PE in pes
800      allocate( mask(0:ndivs-1) )
801      mask = .TRUE.                 !default mask
802      if( PRESENT(maskmap) )then
803          if( size(maskmap).NE.ndivs ) &
804               call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS1D: maskmap array size must equal number of domain divisions.' )
805          mask(:) = maskmap(:)
806      end if
807      if( count(mask).NE.size(pes) ) &
808           call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS1D: number of TRUEs in maskmap array must match PE count.' )
809       
810      if( PRESENT(extent) )then
811          if( size(extent).NE.ndivs ) &
812               call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS1D: extent array size must equal number of domain divisions.' )
813      end if
814!get halosize
815      halosz = 0
816      if( PRESENT(halo) )halosz = halo
817
818!get flags
819      compute_domain_is_global = .FALSE.
820      data_domain_is_global    = .FALSE.
821      domain%cyclic = .FALSE.
822      if( PRESENT(flags) )then
823!NEW: obsolete flag global_compute_domain, since ndivs is non-optional and you cannot have global compute and ndivs.NE.1
824          compute_domain_is_global = ndivs.EQ.1
825!if compute domain is global, data domain must also be
826          data_domain_is_global    = BTEST(flags,GLOBAL) .OR. compute_domain_is_global
827          domain%cyclic  = BTEST(flags,CYCLIC) .AND. halosz.NE.0
828      end if
829
830!set up links list
831      allocate( domain%list(0:ndivs-1) )
832
833!set global domain
834      domain%list(:)%global%begin     = isg
835      domain%list(:)%global%end       = ieg
836      domain%list(:)%global%size      = ieg-isg+1
837      domain%list(:)%global%max_size  = ieg-isg+1
838      domain%list(:)%global%is_global = .TRUE. !always
839
840!get compute domain
841      if( compute_domain_is_global )then
842          domain%list(:)%compute%begin = isg
843          domain%list(:)%compute%end   = ieg
844          domain%list(:)%compute%is_global = .TRUE.
845          domain%list(:)%pe = pes(:)
846          domain%pos = 0
847      else
848          domain%list(:)%compute%is_global = .FALSE.
849          is = isg
850          lastie=-1
851          n = 0
852          do ndiv=0,ndivs-1
853             if( PRESENT(extent) )then
854!rv build the grid by individual block sizes
855               if(blocks) then
856                 if(mask(ndiv)) then
857                   is=isg+offset(ndiv)
858                   ie = is + extent(ndiv) - 1
859!                 print*,pe,':',ndiv,blocks,is,ie,ieg,offset
860                 endif
861!                 if( ndiv.EQ.ndivs-1 .AND. ie.NE.ieg ) &
862!                      call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: extent array limits do not match global domain.' )
863               else
864                 ie = is + extent(ndiv) - 1
865                 if( ndiv.EQ.ndivs-1 .AND. ie.NE.ieg ) &
866                      call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: extent array limits do not match global domain.' )
867               endif
868             else
869!modified for mirror-symmetry
870!original line
871!                 ie = is + CEILING( float(ieg-is+1)/(ndivs-ndiv) ) - 1
872
873!problem of dividing nx points into n domains maintaining symmetry
874!i.e nx=18 n=4 4554 and 5445 are solutions but 4455 is not.
875!this will always work for nx even n even or odd
876!this will always work for nx odd, n odd
877!this will never  work for nx odd, n even: for this case we supersede the mirror calculation
878!                 symmetrize = .NOT. ( mod(ndivs,2).EQ.0 .AND. mod(ieg-isg+1,2).EQ.1 )
879!nx even n odd fails if n>nx/2
880                 symmetrize = ( even(ndivs) .AND. even(ieg-isg+1) ) .OR. &
881                              (  odd(ndivs) .AND.  odd(ieg-isg+1) ) .OR. &
882                              (  odd(ndivs) .AND. even(ieg-isg+1) .AND. ndivs.LT.(ieg-isg+1)/2 )
883
884!mirror domains are stored in the list and retrieved if required.
885                 if( ndiv.EQ.0 )then
886!initialize max points and max domains
887                     imax = ieg
888                     ndmax = ndivs
889                 end if
890!do bottom half of decomposition, going over the midpoint for odd ndivs
891                 if( ndiv.LT.(ndivs-1)/2+1 )then
892!domain is sized by dividing remaining points by remaining domains
893                     ie = is + CEILING( REAL(imax-is+1)/(ndmax-ndiv) ) - 1
894                     ndmirror = (ndivs-1) - ndiv !mirror domain
895                     if( ndmirror.GT.ndiv .AND. symmetrize )then !only for domains over the midpoint
896!mirror extents, the max(,) is to eliminate overlaps
897                         domain%list(ndmirror)%compute%begin = max( isg+ieg-ie, ie+1 )
898                         domain%list(ndmirror)%compute%end   = max( isg+ieg-is, ie+1 )
899                         imax = domain%list(ndmirror)%compute%begin - 1
900                         ndmax = ndmax - 1
901                     end if
902                 else
903                     if( symmetrize )then
904!do top half of decomposition by retrieving saved values
905                         is = domain%list(ndiv)%compute%begin
906                         ie = domain%list(ndiv)%compute%end
907                     else
908                         ie = is + CEILING( REAL(imax-is+1)/(ndmax-ndiv) ) - 1
909                     end if
910                 end if
911             end if
912             if(blocks) then
913               if(mask(ndiv))then
914                 if( ie.LT.is )call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: domain extents must be positive definite.' )
915                 domain%list(ndiv)%compute%begin = is
916                 domain%list(ndiv)%compute%end   = ie
917!                 if(lastie.gt.0.and.is.ne.lastie+1) &
918!                   call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: block extents do not span space completely.' )
919                 lastie=ie
920                 domain%list(ndiv)%pe = pes(n)
921                 if( pe.EQ.pes(n) )domain%pos = ndiv
922                 n = n + 1
923                 is = ie + 1
924               else
925!rv for Arnaud
926                 domain%list(ndiv)%compute%begin = 0
927                 domain%list(ndiv)%compute%end   = 0
928               endif
929               if(ndiv.eq.ndivs-1.and.lastie.gt.ieg) &
930                   call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: block extents exceed max. space.' )
931             else
932             if( ie.LT.is )call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: domain extents must be positive definite.' )
933             domain%list(ndiv)%compute%begin = is
934             domain%list(ndiv)%compute%end   = ie
935             if( ndiv.GT.0 .AND. is.NE.domain%list(ndiv-1)%compute%end+1 ) &
936                  call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: domain extents do not span space completely.' )
937             if( ndiv.EQ.ndivs-1 .AND. domain%list(ndiv)%compute%end.NE.ieg ) &
938                  call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: domain extents do not span space completely.' )
939             if( mask(ndiv) )then
940                 domain%list(ndiv)%pe = pes(n)
941                 if( pe.EQ.pes(n) )domain%pos = ndiv
942                 n = n + 1
943             end if
944             is = ie + 1
945             endif
946          end do
947      end if
948         
949      domain%list(:)%compute%size  = domain%list(:)%compute%end - domain%list(:)%compute%begin + 1
950
951!get data domain
952!data domain is at least equal to compute domain
953      domain%list(:)%data%begin = domain%list(:)%compute%begin
954      domain%list(:)%data%end   = domain%list(:)%compute%end
955      domain%list(:)%data%is_global = .FALSE.
956!apply global flags
957      if( data_domain_is_global )then
958          domain%list(:)%data%begin  = isg
959          domain%list(:)%data%end    = ieg
960          domain%list(:)%data%is_global = .TRUE.
961      end if
962!apply margins
963      domain%list(:)%data%begin = domain%list(:)%data%begin - halosz
964      domain%list(:)%data%end   = domain%list(:)%data%end   + halosz 
965      domain%list(:)%data%size  = domain%list(:)%data%end - domain%list(:)%data%begin + 1
966
967!      domain = domain%list(pos) !load domain from domain%list(pos)
968      domain%compute = domain%list(domain%pos)%compute
969      domain%data = domain%list(domain%pos)%data
970      domain%global = domain%list(domain%pos)%global
971      domain%compute%max_size = MAXVAL( domain%list(:)%compute%size )
972      domain%data%max_size    = MAXVAL( domain%list(:)%data%size )
973      domain%global%max_size  = domain%global%size
974
975!PV786667: the deallocate stmts can be removed when fixed (7.3.1.3m)
976      deallocate( pes, mask )
977      return
978
979      contains
980
981        function if_overlap( hs, he, cs, ce, os, oe )
982!function to look if a portion of the halo [hs,he] lies with in the compute region [cs,ce]
983!if yes, if_overlap returns true, and the overlap region is returned in [os,oe]
984          logical :: if_overlap
985          integer, intent(in) :: hs, he, cs, ce
986          integer, intent(out) :: os, oe
987          os = max(hs,cs)
988          oe = min(he,ce)
989          if( debug )write( stderr(),'(a,7i4)' ) &
990               'MPP_DEFINE_DOMAINS1D: pe, hs, he, cs, ce, os, oe=', pe, hs, he, cs, ce, os, oe
991          if_overlap = oe.GE.os
992          return
993        end function if_overlap
994         
995    end subroutine mpp_define_domains1D
996
997    subroutine mpp_define_domains2D( global_indices, layout, domain, pelist, &
998                                                     xflags, yflags, &
999                                     xhalo, yhalo, xextent, yextent, maskmap,offsetx,offsety, name )
1000!define 2D data and computational domain on global rectilinear cartesian domain (isg:ieg,jsg:jeg) and assign them to PEs
1001      integer, intent(in) :: global_indices(4) !(/ isg, ieg, jsg, jeg /)
1002      integer, intent(in) :: layout(2)
1003      type(domain2D), intent(inout) :: domain
1004      integer, intent(in), optional :: pelist(0:)
1005      integer, intent(in), optional :: xflags, yflags, xhalo, yhalo
1006      integer, intent(in), optional :: xextent(0:), yextent(0:)
1007      logical, intent(in), optional :: maskmap(0:,0:)
1008      integer, intent(in), optional :: offsetx(0:), offsety(0:)
1009      character(len=*), intent(in), optional :: name
1010      integer :: i, j, m, n
1011      integer :: ipos, jpos, pos
1012      integer :: ndivx, ndivy, isg, ieg, jsg, jeg, isd, ied, jsd, jed
1013     
1014      logical, allocatable :: mask(:,:)
1015!rv
1016      logical :: blocks
1017!rv
1018      integer, allocatable :: pes(:), pearray(:,:)
1019      character(len=8) :: text
1020
1021      if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS2D: You must first call mpp_domains_init.' )
1022      ndivx = layout(1); ndivy = layout(2)
1023      isg = global_indices(1); ieg = global_indices(2); jsg = global_indices(3); jeg = global_indices(4)
1024
1025      if( PRESENT(pelist) )then
1026          allocate( pes(0:size(pelist)-1) )
1027          pes = pelist
1028      else
1029          allocate( pes(0:mpp_npes()-1) )
1030          call mpp_get_current_pelist(pes)
1031!          pes = (/ (i,i=0,mpp_npes()-1) /)
1032      end if
1033!rv Block distribution ?
1034      blocks=.false.
1035      if(present(offsetx).and.present(offsety)) blocks=.true.
1036      if(present(offsetx).and.(.not.present(offsety))) then
1037         call mpp_error(FATAL,'MPP_DEFINE_DOMAINS1D: you have to define offsetx AND offsety!')
1038      endif
1039      if((.not.present(offsetx)).and.present(offsety)) then
1040         call mpp_error(FATAL,'MPP_DEFINE_DOMAINS1D: you have to define offsetx AND offsety!')
1041      endif
1042      if(blocks) then
1043       if((.not.present(maskmap)).and.(.not.present(xextent))) &
1044         call mpp_error(FATAL, 'MPP_DEFINE_DOMAINS1D: you have to define maskmap and xextent!')
1045       if((.not.present(maskmap)).and.(.not.present(yextent))) &
1046         call mpp_error(FATAL, 'MPP_DEFINE_DOMAINS1D: you have to define maskmap and yextent!')
1047       if(size(offsetx).ne.ndivx) &
1048         call mpp_error(FATAL,'MPP_DEFINE_DOMAINS2D: offsetx array does not match layout!')
1049       if(size(offsety).ne.ndivy) &
1050         call mpp_error(FATAL,'MPP_DEFINE_DOMAINS2D: offsety array does not match layout!')
1051      endif
1052!rv
1053
1054      allocate( mask(0:ndivx-1,0:ndivy-1) )
1055      mask = .TRUE.
1056      if( PRESENT(maskmap) )then
1057          if( size(maskmap,1).NE.ndivx .OR. size(maskmap,2).NE.ndivy ) &
1058               call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS2D: maskmap array does not match layout.' )
1059          mask(:,:) = maskmap(:,:)
1060      end if
1061!number of unmask domains in layout must equal number of PEs assigned
1062      n = count(mask)
1063      if( n.NE.size(pes) )then
1064          write( text,'(i8)' )n
1065          call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS2D: incorrect number of PEs assigned for this layout and maskmap. Use ' &
1066               //text//' PEs for this domain decomposition.' )
1067      end if
1068
1069!place on PE array; need flag to assign them to j first and then i
1070      allocate( pearray(0:ndivx-1,0:ndivy-1) )
1071      pearray(:,:) = NULL_PE
1072      ipos = NULL_PE; jpos = NULL_PE; pos = NULL_PE
1073      n = 0
1074      do j = 0,ndivy-1
1075         do i = 0,ndivx-1
1076            if( mask(i,j) )then
1077                pearray(i,j) = pes(n)
1078                if( pes(n).EQ.pe )then
1079                    pos = n
1080                    ipos = i
1081                    jpos = j
1082                end if
1083                n = n + 1
1084            end if
1085         end do
1086      end do
1087      if( ipos.EQ.NULL_PE .OR. jpos.EQ.NULL_PE .or. pos.EQ.NULL_PE ) &
1088           call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS2D: pelist must include this PE.' )
1089      if( debug )write( stderr(), * )'pe, ipos, jpos=', pe, ipos, jpos, ' pearray(:,jpos)=', pearray(:,jpos), &
1090                                                                        ' pearray(ipos,:)=', pearray(ipos,:)
1091     
1092      if(blocks) then
1093!do domain decomposition using 1D versions in X and Y
1094      call mpp_define_domains( global_indices(1:2), ndivx, domain%x, &
1095           pack(pearray(:,jpos),mask(:,jpos)), xflags, xhalo, xextent, mask(:,jpos),offset=offsetx )
1096      call mpp_define_domains( global_indices(3:4), ndivy, domain%y, &
1097           pack(pearray(ipos,:),mask(ipos,:)), yflags, yhalo, yextent, mask(ipos,:),offset=offsety )
1098      else
1099!do domain decomposition using 1D versions in X and Y
1100      call mpp_define_domains( global_indices(1:2), ndivx, domain%x, &
1101           pack(pearray(:,jpos),mask(:,jpos)), xflags, xhalo, xextent, mask(:,jpos))
1102      call mpp_define_domains( global_indices(3:4), ndivy, domain%y, &
1103           pack(pearray(ipos,:),mask(ipos,:)), yflags, yhalo, yextent, mask(ipos,:))
1104      endif
1105      if( domain%x%list(domain%x%pos)%pe.NE.domain%y%list(domain%y%pos)%pe ) &
1106           call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS2D: domain%x%list(ipos)%pe.NE.domain%y%list(jpos)%pe.' ) 
1107      domain%pos = pos
1108      domain%pe  = pe
1109
1110!set up fold
1111      domain%fold = 0
1112      if( PRESENT(xflags) )then
1113          if( BTEST(xflags,WEST) )domain%fold = domain%fold + FOLD_WEST_EDGE
1114          if( BTEST(xflags,EAST) )domain%fold = domain%fold + FOLD_EAST_EDGE
1115      end if
1116      if( PRESENT(yflags) )then
1117          if( BTEST(yflags,SOUTH) )domain%fold = domain%fold + FOLD_SOUTH_EDGE
1118          if( BTEST(yflags,NORTH) )domain%fold = domain%fold + FOLD_NORTH_EDGE
1119      end if
1120         
1121      if( BTEST(domain%fold,SOUTH) .OR. BTEST(domain%fold,NORTH) )then
1122          if( domain%y%cyclic )call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: an axis cannot be both folded and cyclic.' )
1123          if( modulo(domain%x%global%size,2).NE.0 ) &
1124               call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: number of points in X must be even when there is a fold in Y.' )
1125!check if folded domain boundaries line up in X: compute domains lining up is a sufficient condition for symmetry
1126          n = ndivx - 1
1127          do i = 0,n/2
1128             if( domain%x%list(i)%compute%size.NE.domain%x%list(n-i)%compute%size ) &
1129                  call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: Folded domain boundaries must line up (mirror-symmetric extents).' )
1130          end do
1131      end if
1132      if( BTEST(domain%fold,WEST) .OR. BTEST(domain%fold,EAST) )then
1133          if( domain%x%cyclic )call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: an axis cannot be both folded and cyclic.' )
1134          if( modulo(domain%y%global%size,2).NE.0 ) &
1135               call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: number of points in Y must be even when there is a fold in X.' )
1136!check if folded domain boundaries line up in Y: compute domains lining up is a sufficient condition for symmetry
1137          n = ndivy - 1
1138          do i = 0,n/2
1139             if( domain%y%list(i)%compute%size.NE.domain%y%list(n-i)%compute%size ) &
1140                  call mpp_error( FATAL, 'MPP_DEFINE_DOMAINS: Folded domain boundaries must line up (mirror-symmetric extents).' )
1141          end do
1142      end if
1143
1144!set up domain%list
1145      if( debug )write( stderr(),'(a,9i4)' )'pe, domain=', pe, domain_info_buf(1:8)
1146      if( pe.EQ.mpp_root_pe() .AND. PRESENT(name) )then
1147          write( stdlog(), '(/a,i3,a,i3)' )trim(name)//' domain decomposition: ', ndivx, ' X', ndivy
1148          write( stdlog(), '(3x,a)' )'pe,   is,  ie,  js,  je,    isd, ied, jsd, jed'
1149      end if
1150      call mpp_sync()
1151      call mpp_get_compute_domain( domain, domain_info_buf(1), domain_info_buf(2), domain_info_buf(3), domain_info_buf(4) )
1152      call mpp_get_data_domain   ( domain, domain_info_buf(5), domain_info_buf(6), domain_info_buf(7), domain_info_buf(8) )
1153      n = size(pes)
1154      allocate( domain%list(0:n-1) ) !this is only used for storage of remote compute and data domain info
1155      do i = 0,n-1
1156         m = mod(pos+i,n)
1157         domain%list(m)%pe = pes(m)
1158         call mpp_transmit( domain_info_buf(1:8), 8, pes(mod(pos+n-i,n)), domain_info_buf(9:16), 8, pes(m) )
1159         domain%list(m)%x%compute%begin = domain_info_buf(9)
1160         domain%list(m)%x%compute%end   = domain_info_buf(10)
1161         domain%list(m)%y%compute%begin = domain_info_buf(11)
1162         domain%list(m)%y%compute%end   = domain_info_buf(12)
1163         domain%list(m)%x%data%begin = domain_info_buf(13)
1164         domain%list(m)%x%data%end   = domain_info_buf(14)
1165         domain%list(m)%y%data%begin = domain_info_buf(15)
1166         domain%list(m)%y%data%end   = domain_info_buf(16)
1167         if( pe.EQ.mpp_root_pe() .AND. PRESENT(name) )write( stdlog(), '(2x,i3,x,4i5,3x,4i5)' )pes(m), domain_info_buf(9:)
1168      end do
1169      call mpp_sync_self(pes)
1170      domain%list(:)%x%compute%size = domain%list(:)%x%compute%end - domain%list(:)%x%compute%begin + 1
1171      domain%list(:)%y%compute%size = domain%list(:)%y%compute%end - domain%list(:)%y%compute%begin + 1
1172      domain%list(:)%x%data%size = domain%list(:)%x%data%end - domain%list(:)%x%data%begin + 1
1173      domain%list(:)%y%data%size = domain%list(:)%y%data%end - domain%list(:)%y%data%begin + 1
1174
1175      domain%remote_domains_initialized = .FALSE.
1176      domain%remote_off_domains_initialized = .FALSE.
1177      call compute_overlaps(domain)
1178!PV786667: the deallocate stmts can be removed when fixed (7.3.1.3m)
1179      deallocate( pes, mask, pearray )
1180         
1181      return
1182    end subroutine mpp_define_domains2D
1183
1184    subroutine mpp_broadcast_domain( domain )
1185!broadcast domain (useful only outside the context of its own pelist)
1186      type(domain2D), intent(inout) :: domain
1187      integer, allocatable :: pes(:)
1188      logical :: native         !true if I'm on the pelist of this domain
1189      integer :: listsize, listpos
1190      integer :: n
1191      integer, dimension(5) :: msg, info         !pe and compute domain of each item in list
1192
1193      if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_BROADCAST_DOMAIN: You must first call mpp_domains_init.' )
1194
1195!get the current pelist
1196      allocate( pes(0:mpp_npes()-1) )
1197      call mpp_get_current_pelist(pes)
1198
1199!am I part of this domain?
1200      native = ASSOCIATED(domain%list)
1201
1202!set local list size
1203      if( native )then
1204          listsize = size(domain%list)
1205      else
1206          listsize = 0
1207      end if
1208      call mpp_max(listsize)
1209
1210      if( .NOT.native )then
1211!initialize domain%list and set null values in message
1212          allocate( domain%list(0:listsize-1) )
1213          domain%pe = NULL_PE
1214          domain%pos = -1
1215          domain%x%compute%begin =  HUGE(1)
1216          domain%x%compute%end   = -HUGE(1)
1217          domain%y%compute%begin =  HUGE(1)
1218          domain%y%compute%end   = -HUGE(1)
1219      end if
1220!initialize values in info
1221      info(1) = domain%pe
1222      call mpp_get_compute_domain( domain, info(2), info(3), info(4), info(5) )
1223
1224!broadcast your info across current pelist and unpack if needed
1225      listpos = 0
1226      do n = 0,mpp_npes()-1
1227         msg = info
1228         if( pe.EQ.pes(n) .AND. debug )write( stderr(),* )'PE ', pe, 'broadcasting msg ', msg
1229         call mpp_broadcast( msg, 5, pes(n) )
1230!no need to unpack message if native
1231!no need to unpack message from non-native PE
1232         if( .NOT.native .AND. msg(1).NE.NULL_PE )then
1233             domain%list(listpos)%pe = msg(1)
1234             domain%list(listpos)%x%compute%begin = msg(2)
1235             domain%list(listpos)%x%compute%end   = msg(3)
1236             domain%list(listpos)%y%compute%begin = msg(4)
1237             domain%list(listpos)%y%compute%end   = msg(5)
1238             listpos = listpos + 1
1239             if( debug )write( stderr(),* )'PE ', pe, 'received domain from PE ', msg(1), 'is,ie,js,je=', msg(2:5)
1240         end if
1241      end do
1242    end subroutine mpp_broadcast_domain
1243
1244    subroutine compute_overlaps( domain )
1245!computes remote domain overlaps
1246!assumes only one in each direction
1247      type(domain2D), intent(inout) :: domain
1248      integer :: i, j, k, m, n
1249      integer :: is, ie, js, je, isc, iec, jsc, jec, isd, ied, jsd, jed, isg, ieg, jsg, jeg, ioff, joff
1250      integer :: list
1251
1252      if( grid_offset_type.EQ.AGRID .AND. domain%remote_domains_initialized     )return
1253      if( grid_offset_type.NE.AGRID .AND. domain%remote_off_domains_initialized )return
1254      domain%gridtype = grid_offset_type
1255      n = size(domain%list)
1256!send
1257      call mpp_get_compute_domain( domain, isc, iec, jsc, jec )
1258      call mpp_get_global_domain ( domain, isg, ieg, jsg, jeg, xsize=ioff, ysize=joff ) !cyclic offsets
1259      domain%list(:)%overlap = .FALSE.
1260      do list = 0,n-1
1261         m = mod( domain%pos+list, n )
1262!to_pe's eastern halo
1263         is = domain%list(m)%x%compute%end+1; ie = domain%list(m)%x%data%end
1264         js = domain%list(m)%y%compute%begin; je = domain%list(m)%y%compute%end
1265         if( is.GT.ieg )then
1266             if( domain%x%cyclic )then !try cyclic offset
1267                 is = is-ioff; ie = ie-ioff
1268             else if( BTEST(domain%fold,EAST) )then
1269                 i=is; is = 2*ieg-ie+1; ie = 2*ieg-i+1
1270                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1271                 if( BTEST(grid_offset_type,EAST) )then
1272                     is = is - 1; ie = ie - 1
1273                 end if
1274             end if
1275         end if
1276         is = max(is,isc); ie = min(ie,iec)
1277         js = max(js,jsc); je = min(je,jec)
1278         if( ie.GE.is .AND. je.GE.js )then
1279             domain%list(m)%overlap = .TRUE.
1280             if( grid_offset_type.NE.AGRID )then
1281                 domain%list(m)%send_w_off%overlap = .TRUE.
1282                 domain%list(m)%send_w_off%is = is
1283                 domain%list(m)%send_w_off%ie = ie
1284                 domain%list(m)%send_w_off%js = js
1285                 domain%list(m)%send_w_off%je = je
1286             else
1287                 domain%list(m)%send_w%overlap = .TRUE.
1288                 domain%list(m)%send_w%is = is
1289                 domain%list(m)%send_w%ie = ie
1290                 domain%list(m)%send_w%js = js
1291                 domain%list(m)%send_w%je = je
1292             end if
1293         else
1294             domain%list(m)%send_w%overlap = .FALSE.
1295             domain%list(m)%send_w_off%overlap = .FALSE.
1296         end if
1297!to_pe's SE halo
1298         is = domain%list(m)%x%compute%end+1; ie = domain%list(m)%x%data%end
1299         js = domain%list(m)%y%data%begin; je = domain%list(m)%y%compute%begin-1
1300         if( is.GT.ieg )then
1301             if( domain%x%cyclic )then !try cyclic offset
1302                 is = is-ioff; ie = ie-ioff
1303             else if( BTEST(domain%fold,EAST) )then
1304                 i=is; is = 2*ieg-ie+1; ie = 2*ieg-i+1
1305                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1306                 if( BTEST(grid_offset_type,EAST) )then
1307                     is = is - 1; ie = ie - 1
1308                 end if
1309             end if
1310         end if
1311         if( jsg.GT.je )then
1312             if( domain%y%cyclic )then !try cyclic offset
1313                 js = js+joff; je = je+joff
1314             else if( BTEST(domain%fold,SOUTH) )then
1315                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1316                 j=js; js = 2*jsg-je-1; je = 2*jsg-j-1
1317                 if( BTEST(grid_offset_type,SOUTH) )then
1318                     js = js + 1; je = je + 1
1319                 end if
1320             end if
1321         end if
1322         is = max(is,isc); ie = min(ie,iec)
1323         js = max(js,jsc); je = min(je,jec)
1324         if( ie.GE.is .AND. je.GE.js )then
1325             domain%list(m)%overlap = .TRUE.
1326             if( grid_offset_type.NE.AGRID )then
1327                 domain%list(m)%send_nw_off%overlap = .TRUE.
1328                 domain%list(m)%send_nw_off%is = is
1329                 domain%list(m)%send_nw_off%ie = ie
1330                 domain%list(m)%send_nw_off%js = js
1331                 domain%list(m)%send_nw_off%je = je
1332             else
1333                 domain%list(m)%send_nw%overlap = .TRUE.
1334                 domain%list(m)%send_nw%is = is
1335                 domain%list(m)%send_nw%ie = ie
1336                 domain%list(m)%send_nw%js = js
1337                 domain%list(m)%send_nw%je = je
1338             end if
1339         else
1340             domain%list(m)%send_nw%overlap = .FALSE.
1341             domain%list(m)%send_nw_off%overlap = .FALSE.
1342         end if
1343!to_pe's southern halo
1344         is = domain%list(m)%x%compute%begin; ie = domain%list(m)%x%compute%end
1345         js = domain%list(m)%y%data%begin; je = domain%list(m)%y%compute%begin-1
1346         if( jsg.GT.je )then
1347             if( domain%y%cyclic )then !try cyclic offset
1348                 js = js+joff; je = je+joff
1349             else if( BTEST(domain%fold,SOUTH) )then
1350                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1351                 j=js; js = 2*jsg-je-1; je = 2*jsg-j-1
1352                 if( BTEST(grid_offset_type,SOUTH) )then
1353                     js = js + 1; je = je + 1
1354                 end if
1355             end if
1356         end if
1357         is = max(is,isc); ie = min(ie,iec)
1358         js = max(js,jsc); je = min(je,jec)
1359         if( ie.GE.is .AND. je.GE.js )then
1360             domain%list(m)%overlap = .TRUE.
1361             if( grid_offset_type.NE.AGRID )then
1362                 domain%list(m)%send_n_off%overlap = .TRUE.
1363                 domain%list(m)%send_n_off%is = is
1364                 domain%list(m)%send_n_off%ie = ie
1365                 domain%list(m)%send_n_off%js = js
1366                 domain%list(m)%send_n_off%je = je
1367             else
1368                 domain%list(m)%send_n%overlap = .TRUE.
1369                 domain%list(m)%send_n%is = is
1370                 domain%list(m)%send_n%ie = ie
1371                 domain%list(m)%send_n%js = js
1372                 domain%list(m)%send_n%je = je
1373             end if
1374         else
1375             domain%list(m)%send_n%overlap = .FALSE.
1376             domain%list(m)%send_n_off%overlap = .FALSE.
1377         end if
1378!to_pe's SW halo
1379         is = domain%list(m)%x%data%begin; ie = domain%list(m)%x%compute%begin-1
1380         js = domain%list(m)%y%data%begin; je = domain%list(m)%y%compute%begin-1
1381         if( isg.GT.ie )then
1382             if( domain%x%cyclic )then !try cyclic offset
1383                 is = is+ioff; ie = ie+ioff
1384             else if( BTEST(domain%fold,WEST) )then
1385                 i=is; is = 2*isg-ie-1; ie = 2*isg-i-1
1386                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1387                 if( BTEST(grid_offset_type,WEST) )then
1388                     is = is + 1; ie = ie + 1
1389                 end if
1390             end if
1391         end if
1392         if( jsg.GT.je )then
1393             if( domain%y%cyclic )then !try cyclic offset
1394                 js = js+joff; je = je+joff
1395             else if( BTEST(domain%fold,SOUTH) )then
1396                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1397                 j=js; js = 2*jsg-je-1; je = 2*jsg-j-1
1398                 if( BTEST(grid_offset_type,SOUTH) )then
1399                     js = js + 1; je = je + 1
1400                 end if
1401             end if
1402         end if
1403         is = max(is,isc); ie = min(ie,iec)
1404         js = max(js,jsc); je = min(je,jec)
1405         if( ie.GE.is .AND. je.GE.js )then
1406             domain%list(m)%overlap = .TRUE.
1407             if( grid_offset_type.NE.AGRID )then
1408                 domain%list(m)%send_ne_off%overlap = .TRUE.
1409                 domain%list(m)%send_ne_off%is = is
1410                 domain%list(m)%send_ne_off%ie = ie
1411                 domain%list(m)%send_ne_off%js = js
1412                 domain%list(m)%send_ne_off%je = je
1413             else
1414                 domain%list(m)%send_ne%overlap = .TRUE.
1415                 domain%list(m)%send_ne%is = is
1416                 domain%list(m)%send_ne%ie = ie
1417                 domain%list(m)%send_ne%js = js
1418                 domain%list(m)%send_ne%je = je
1419             end if
1420         else
1421             domain%list(m)%send_ne%overlap = .FALSE.
1422             domain%list(m)%send_ne_off%overlap = .FALSE.
1423         end if
1424!to_pe's western halo
1425         is = domain%list(m)%x%data%begin; ie = domain%list(m)%x%compute%begin-1
1426         js = domain%list(m)%y%compute%begin; je = domain%list(m)%y%compute%end
1427         if( isg.GT.ie )then
1428             if( domain%x%cyclic )then !try cyclic offset
1429                 is = is+ioff; ie = ie+ioff
1430             else if( BTEST(domain%fold,WEST) )then
1431                 i=is; is = 2*isg-ie-1; ie = 2*isg-i-1
1432                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1433                 if( BTEST(grid_offset_type,WEST) )then
1434                     is = is + 1; ie = ie + 1
1435                 end if
1436             end if
1437         end if
1438         is = max(is,isc); ie = min(ie,iec)
1439         js = max(js,jsc); je = min(je,jec)
1440         if( ie.GE.is .AND. je.GE.js )then
1441             domain%list(m)%overlap = .TRUE.
1442             if( grid_offset_type.NE.AGRID )then
1443                 domain%list(m)%send_e_off%overlap = .TRUE.
1444                 domain%list(m)%send_e_off%is = is
1445                 domain%list(m)%send_e_off%ie = ie
1446                 domain%list(m)%send_e_off%js = js
1447                 domain%list(m)%send_e_off%je = je
1448             else
1449                 domain%list(m)%send_e%overlap = .TRUE.
1450                 domain%list(m)%send_e%is = is
1451                 domain%list(m)%send_e%ie = ie
1452                 domain%list(m)%send_e%js = js
1453                 domain%list(m)%send_e%je = je
1454             end if
1455         else
1456             domain%list(m)%send_e%overlap = .FALSE.
1457             domain%list(m)%send_e_off%overlap = .FALSE.
1458         end if
1459!to_pe's NW halo
1460         is = domain%list(m)%x%data%begin; ie = domain%list(m)%x%compute%begin-1
1461         js = domain%list(m)%y%compute%end+1; je = domain%list(m)%y%data%end
1462         if( isg.GT.ie )then
1463             if( domain%x%cyclic )then !try cyclic offset
1464                 is = is+ioff; ie = ie+ioff
1465             else if( BTEST(domain%fold,WEST) )then
1466                 i=is; is = 2*isg-ie-1; ie = 2*isg-i-1
1467                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1468                 if( BTEST(grid_offset_type,WEST) )then
1469                     is = is + 1; ie = ie + 1
1470                 end if
1471             end if
1472         end if
1473         if( js.GT.jeg )then
1474             if( domain%y%cyclic )then !try cyclic offset
1475                 js = js-joff; je = je-joff
1476             else if( BTEST(domain%fold,NORTH) )then
1477                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1478                 j=js; js = 2*jeg-je+1; je = 2*jeg-j+1
1479                 if( BTEST(grid_offset_type,NORTH) )then
1480                     js = js - 1; je = je - 1
1481                 end if
1482             end if
1483         end if
1484         is = max(is,isc); ie = min(ie,iec)
1485         js = max(js,jsc); je = min(je,jec)
1486         if( ie.GE.is .AND. je.GE.js )then
1487             domain%list(m)%overlap = .TRUE.
1488             if( grid_offset_type.NE.AGRID )then
1489                 domain%list(m)%send_se_off%overlap = .TRUE.
1490                 domain%list(m)%send_se_off%is = is
1491                 domain%list(m)%send_se_off%ie = ie
1492                 domain%list(m)%send_se_off%js = js
1493                 domain%list(m)%send_se_off%je = je
1494             else
1495                 domain%list(m)%send_se%overlap = .TRUE.
1496                 domain%list(m)%send_se%is = is
1497                 domain%list(m)%send_se%ie = ie
1498                 domain%list(m)%send_se%js = js
1499                 domain%list(m)%send_se%je = je
1500             end if
1501         else
1502             domain%list(m)%send_se%overlap = .FALSE.
1503             domain%list(m)%send_se_off%overlap = .FALSE.
1504         end if
1505!to_pe's northern halo
1506         is = domain%list(m)%x%compute%begin; ie = domain%list(m)%x%compute%end
1507         js = domain%list(m)%y%compute%end+1; je = domain%list(m)%y%data%end
1508         if( js.GT.jeg )then
1509             if( domain%y%cyclic )then !try cyclic offset
1510                 js = js-joff; je = je-joff
1511             else if( BTEST(domain%fold,NORTH) )then
1512                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1513                 j=js; js = 2*jeg-je+1; je = 2*jeg-j+1
1514                 if( BTEST(grid_offset_type,NORTH) )then
1515                     js = js - 1; je = je - 1
1516                 end if
1517             end if
1518         end if
1519         is = max(is,isc); ie = min(ie,iec)
1520         js = max(js,jsc); je = min(je,jec)
1521         if( ie.GE.is .AND. je.GE.js )then
1522             domain%list(m)%overlap = .TRUE.
1523             if( grid_offset_type.NE.AGRID )then
1524                 domain%list(m)%send_s_off%overlap = .TRUE.
1525                 domain%list(m)%send_s_off%is = is
1526                 domain%list(m)%send_s_off%ie = ie
1527                 domain%list(m)%send_s_off%js = js
1528                 domain%list(m)%send_s_off%je = je
1529             else
1530                 domain%list(m)%send_s%overlap = .TRUE.
1531                 domain%list(m)%send_s%is = is
1532                 domain%list(m)%send_s%ie = ie
1533                 domain%list(m)%send_s%js = js
1534                 domain%list(m)%send_s%je = je
1535             end if
1536         else
1537             domain%list(m)%send_s%overlap = .FALSE.
1538             domain%list(m)%send_s_off%overlap = .FALSE.
1539         end if
1540!to_pe's NE halo
1541         is = domain%list(m)%x%compute%end+1; ie = domain%list(m)%x%data%end
1542         js = domain%list(m)%y%compute%end+1; je = domain%list(m)%y%data%end
1543         if( is.GT.ieg )then
1544             if( domain%x%cyclic )then !try cyclic offset
1545                 is = is-ioff; ie = ie-ioff
1546             else if( BTEST(domain%fold,EAST) )then
1547                 i=is; is = 2*ieg-ie+1; ie = 2*ieg-i+1
1548                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1549             end if
1550         end if
1551         if( js.GT.jeg )then
1552             if( domain%y%cyclic )then !try cyclic offset
1553                 js = js-joff; je = je-joff
1554             else if( BTEST(domain%fold,NORTH) )then
1555                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1556                 j=js; js = 2*jeg-je+1; je = 2*jeg-j+1
1557                 if( BTEST(grid_offset_type,NORTH) )then
1558                     js = js - 1; je = je - 1
1559                 end if
1560             end if
1561         end if
1562         is = max(is,isc); ie = min(ie,iec)
1563         js = max(js,jsc); je = min(je,jec)
1564         if( ie.GE.is .AND. je.GE.js )then
1565             domain%list(m)%overlap = .TRUE.
1566             if( grid_offset_type.NE.AGRID )then
1567                 domain%list(m)%send_sw_off%overlap = .TRUE.
1568                 domain%list(m)%send_sw_off%is = is
1569                 domain%list(m)%send_sw_off%ie = ie
1570                 domain%list(m)%send_sw_off%js = js
1571                 domain%list(m)%send_sw_off%je = je
1572             else
1573                 domain%list(m)%send_sw%overlap = .TRUE.
1574                 domain%list(m)%send_sw%is = is
1575                 domain%list(m)%send_sw%ie = ie
1576                 domain%list(m)%send_sw%js = js
1577                 domain%list(m)%send_sw%je = je
1578             end if
1579         else
1580             domain%list(m)%send_sw%overlap = .FALSE.
1581             domain%list(m)%send_sw_off%overlap = .FALSE.
1582         end if
1583      end do
1584             
1585!recv
1586      do list = 0,n-1
1587         m = mod( domain%pos+n-list, n )
1588         call mpp_get_compute_domain( domain%list(m), isc, iec, jsc, jec )
1589!recv_e
1590         isd = domain%x%compute%end+1; ied = domain%x%data%end
1591         jsd = domain%y%compute%begin; jed = domain%y%compute%end
1592         is=isc; ie=iec; js=jsc; je=jec
1593         domain%list(m)%recv_e%folded = .FALSE.
1594         if( isd.GT.ieg )then
1595             if( domain%x%cyclic )then !try cyclic offset
1596                 is = is+ioff; ie = ie+ioff
1597             else if( BTEST(domain%fold,EAST) )then
1598                 domain%list(m)%recv_e%folded = .TRUE.
1599                 i=is; is = 2*ieg-ie+1; ie = 2*ieg-i+1
1600                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1601                 if( BTEST(grid_offset_type,EAST) )then
1602                     is = is - 1; ie = ie - 1
1603                 end if
1604             end if
1605         end if
1606         is = max(isd,is); ie = min(ied,ie)
1607         js = max(jsd,js); je = min(jed,je)
1608         if( ie.GE.is .AND. je.GE.js )then
1609             domain%list(m)%overlap = .TRUE.
1610             if( grid_offset_type.NE.AGRID )then
1611                 domain%list(m)%recv_e_off%overlap = .TRUE.
1612                 domain%list(m)%recv_e_off%is = is
1613                 domain%list(m)%recv_e_off%ie = ie
1614                 domain%list(m)%recv_e_off%js = js
1615                 domain%list(m)%recv_e_off%je = je
1616             else
1617                 domain%list(m)%recv_e%overlap = .TRUE.
1618                 domain%list(m)%recv_e%is = is
1619                 domain%list(m)%recv_e%ie = ie
1620                 domain%list(m)%recv_e%js = js
1621                 domain%list(m)%recv_e%je = je
1622             endif
1623         else
1624             domain%list(m)%recv_e%overlap = .FALSE.
1625             domain%list(m)%recv_e_off%overlap = .FALSE.
1626         end if
1627!recv_se
1628         isd = domain%x%compute%end+1; ied = domain%x%data%end
1629         jsd = domain%y%data%begin; jed = domain%y%compute%begin-1
1630         is=isc; ie=iec; js=jsc; je=jec
1631         domain%list(m)%recv_se%folded = .FALSE.
1632         if( jed.LT.jsg )then
1633             if( domain%y%cyclic )then !try cyclic offset
1634                 js = js-joff; je = je-joff
1635             else if( BTEST(domain%fold,SOUTH) )then
1636                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1637                 j=js; js = 2*jsg-je-1; je = 2*jsg-j-1
1638                 domain%list(m)%recv_se%folded = .TRUE.
1639                 if( BTEST(grid_offset_type,SOUTH) )then
1640                     js = js + 1; je = je + 1
1641                 end if
1642             end if
1643         end if
1644         if( isd.GT.ieg )then
1645             if( domain%x%cyclic )then !try cyclic offset
1646                 is = is+ioff; ie = ie+ioff
1647             else if( BTEST(domain%fold,EAST) )then
1648                 i=is; is = 2*ieg-ie+1; ie = 2*ieg-i+1
1649                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1650                 domain%list(m)%recv_se%folded = .TRUE.
1651                 if( BTEST(grid_offset_type,EAST) )then
1652                     is = is - 1; ie = ie - 1
1653                 end if
1654             end if
1655         end if
1656         is = max(isd,is); ie = min(ied,ie)
1657         js = max(jsd,js); je = min(jed,je)
1658         if( ie.GE.is .AND. je.GE.js )then
1659             domain%list(m)%overlap = .TRUE.
1660             if( grid_offset_type.NE.AGRID )then
1661                 domain%list(m)%recv_se_off%overlap = .TRUE.
1662                 domain%list(m)%recv_se_off%is = is
1663                 domain%list(m)%recv_se_off%ie = ie
1664                 domain%list(m)%recv_se_off%js = js
1665                 domain%list(m)%recv_se_off%je = je
1666             else
1667                 domain%list(m)%recv_se%overlap = .TRUE.
1668                 domain%list(m)%recv_se%is = is
1669                 domain%list(m)%recv_se%ie = ie
1670                 domain%list(m)%recv_se%js = js
1671                 domain%list(m)%recv_se%je = je
1672             endif
1673         else
1674             domain%list(m)%recv_se%overlap = .FALSE.
1675             domain%list(m)%recv_se_off%overlap = .FALSE.
1676         end if
1677!recv_s
1678         isd = domain%x%compute%begin; ied = domain%x%compute%end
1679         jsd = domain%y%data%begin; jed = domain%y%compute%begin-1
1680         is=isc; ie=iec; js=jsc; je=jec
1681         domain%list(m)%recv_s%folded = .FALSE.
1682         if( jed.LT.jsg )then
1683             if( domain%y%cyclic )then !try cyclic offset
1684                 js = js-joff; je = je-joff
1685             else if( BTEST(domain%fold,SOUTH) )then
1686                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1687                 j=js; js = 2*jsg-je-1; je = 2*jsg-j-1
1688                 domain%list(m)%recv_s%folded = .TRUE.
1689                 if( BTEST(grid_offset_type,SOUTH) )then
1690                     js = js + 1; je = je + 1
1691                 end if
1692             end if
1693         end if
1694         is = max(isd,is); ie = min(ied,ie)
1695         js = max(jsd,js); je = min(jed,je)
1696         if( ie.GE.is .AND. je.GE.js )then
1697             domain%list(m)%overlap = .TRUE.
1698             if( grid_offset_type.NE.AGRID )then
1699                 domain%list(m)%recv_s_off%overlap = .TRUE.
1700                 domain%list(m)%recv_s_off%is = is
1701                 domain%list(m)%recv_s_off%ie = ie
1702                 domain%list(m)%recv_s_off%js = js
1703                 domain%list(m)%recv_s_off%je = je
1704             else
1705                 domain%list(m)%recv_s%overlap = .TRUE.
1706                 domain%list(m)%recv_s%is = is
1707                 domain%list(m)%recv_s%ie = ie
1708                 domain%list(m)%recv_s%js = js
1709                 domain%list(m)%recv_s%je = je
1710             endif
1711         else
1712             domain%list(m)%recv_s%overlap = .FALSE.
1713             domain%list(m)%recv_s_off%overlap = .FALSE.
1714
1715         end if
1716!recv_sw
1717         isd = domain%x%data%begin; ied = domain%x%compute%begin-1
1718         jsd = domain%y%data%begin; jed = domain%y%compute%begin-1
1719         is=isc; ie=iec; js=jsc; je=jec
1720         domain%list(m)%recv_sw%folded = .FALSE.
1721         if( jed.LT.jsg )then
1722             if( domain%y%cyclic )then !try cyclic offset
1723                 js = js-joff; je = je-joff
1724             else if( BTEST(domain%fold,SOUTH) )then
1725                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1726                 j=js; js = 2*jsg-je-1; je = 2*jsg-j-1
1727                 domain%list(m)%recv_sw%folded = .TRUE.
1728                 if( BTEST(grid_offset_type,SOUTH) )then
1729                     js = js + 1; je = je + 1
1730                 end if
1731             end if
1732         end if
1733         if( ied.LT.isg )then
1734             if( domain%x%cyclic )then !try cyclic offset
1735                 is = is-ioff; ie = ie-ioff
1736             else if( BTEST(domain%fold,WEST) )then
1737                 i=is; is = 2*isg-ie-1; ie = 2*isg-i-1
1738                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1739                 domain%list(m)%recv_sw%folded = .TRUE.
1740                 if( BTEST(grid_offset_type,WEST) )then
1741                     is = is + 1; ie = ie + 1
1742                 end if
1743             end if
1744         end if
1745         is = max(isd,is); ie = min(ied,ie)
1746         js = max(jsd,js); je = min(jed,je)
1747         if( ie.GE.is .AND. je.GE.js )then
1748             domain%list(m)%overlap = .TRUE.
1749             if( grid_offset_type.NE.AGRID )then
1750                 domain%list(m)%recv_sw_off%overlap = .TRUE.
1751                 domain%list(m)%recv_sw_off%is = is
1752                 domain%list(m)%recv_sw_off%ie = ie
1753                 domain%list(m)%recv_sw_off%js = js
1754                 domain%list(m)%recv_sw_off%je = je
1755             else
1756                 domain%list(m)%recv_sw%overlap = .TRUE.
1757                 domain%list(m)%recv_sw%is = is
1758                 domain%list(m)%recv_sw%ie = ie
1759                 domain%list(m)%recv_sw%js = js
1760                 domain%list(m)%recv_sw%je = je
1761             endif
1762         else
1763             domain%list(m)%recv_sw%overlap = .FALSE.
1764             domain%list(m)%recv_sw_off%overlap = .FALSE.
1765         end if
1766!recv_w
1767         isd = domain%x%data%begin; ied = domain%x%compute%begin-1
1768         jsd = domain%y%compute%begin; jed = domain%y%compute%end
1769         is=isc; ie=iec; js=jsc; je=jec
1770         domain%list(m)%recv_w%folded = .FALSE.
1771         if( ied.LT.isg )then
1772             if( domain%x%cyclic )then !try cyclic offset
1773                 is = is-ioff; ie = ie-ioff
1774             else if( BTEST(domain%fold,WEST) )then
1775                 i=is; is = 2*isg-ie-1; ie = 2*isg-i-1
1776                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1777                 domain%list(m)%recv_w%folded = .TRUE.
1778                 if( BTEST(grid_offset_type,WEST) )then
1779                     is = is + 1; ie = ie + 1
1780                 end if
1781             end if
1782         end if
1783         is = max(isd,is); ie = min(ied,ie)
1784         js = max(jsd,js); je = min(jed,je)
1785         if( ie.GE.is .AND. je.GE.js )then
1786             domain%list(m)%overlap = .TRUE.
1787             if( grid_offset_type.NE.AGRID )then
1788                 domain%list(m)%recv_w_off%overlap = .TRUE.
1789                 domain%list(m)%recv_w_off%is = is
1790                 domain%list(m)%recv_w_off%ie = ie
1791                 domain%list(m)%recv_w_off%js = js
1792                 domain%list(m)%recv_w_off%je = je
1793             else
1794                 domain%list(m)%recv_w%overlap = .TRUE.
1795                 domain%list(m)%recv_w%is = is
1796                 domain%list(m)%recv_w%ie = ie
1797                 domain%list(m)%recv_w%js = js
1798                 domain%list(m)%recv_w%je = je
1799             endif
1800         else
1801             domain%list(m)%recv_w%overlap = .FALSE.
1802             domain%list(m)%recv_w_off%overlap = .FALSE.
1803         end if
1804!recv_nw
1805         isd = domain%x%data%begin; ied = domain%x%compute%begin-1
1806         jsd = domain%y%compute%end+1; jed = domain%y%data%end
1807         is=isc; ie=iec; js=jsc; je=jec
1808         domain%list(m)%recv_nw%folded = .FALSE.
1809         if( jsd.GT.jeg )then
1810             if( domain%y%cyclic )then !try cyclic offset
1811                 js = js+joff; je = je+joff
1812             else if( BTEST(domain%fold,NORTH) )then
1813                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1814                 j=js; js = 2*jeg-je+1; je = 2*jeg-j+1
1815                 domain%list(m)%recv_nw%folded = .TRUE.
1816                 if( BTEST(grid_offset_type,NORTH) )then
1817                     js = js - 1; je = je - 1
1818                 end if
1819             end if
1820         end if
1821         if( ied.LT.isg )then
1822             if( domain%x%cyclic )then !try cyclic offset
1823                 is = is-ioff; ie = ie-ioff
1824             else if( BTEST(domain%fold,WEST) )then
1825                 i=is; is = 2*isg-ie-1; ie = 2*isg-i-1
1826                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1827                 domain%list(m)%recv_nw%folded = .TRUE.
1828                 if( BTEST(grid_offset_type,WEST) )then
1829                     is = is + 1; ie = ie + 1
1830                 end if
1831             end if
1832         end if
1833         is = max(isd,is); ie = min(ied,ie)
1834         js = max(jsd,js); je = min(jed,je)
1835         if( ie.GE.is .AND. je.GE.js )then
1836             domain%list(m)%overlap = .TRUE.
1837             if( grid_offset_type.NE.AGRID )then
1838                 domain%list(m)%recv_nw_off%overlap = .TRUE.
1839                 domain%list(m)%recv_nw_off%is = is
1840                 domain%list(m)%recv_nw_off%ie = ie
1841                 domain%list(m)%recv_nw_off%js = js
1842                 domain%list(m)%recv_nw_off%je = je
1843             else
1844                 domain%list(m)%recv_nw%overlap = .TRUE.
1845                 domain%list(m)%recv_nw%is = is
1846                 domain%list(m)%recv_nw%ie = ie
1847                 domain%list(m)%recv_nw%js = js
1848                 domain%list(m)%recv_nw%je = je
1849             endif
1850         else
1851             domain%list(m)%recv_nw%overlap = .FALSE.
1852             domain%list(m)%recv_nw_off%overlap = .FALSE.
1853         end if
1854!recv_n
1855         isd = domain%x%compute%begin; ied = domain%x%compute%end
1856         jsd = domain%y%compute%end+1; jed = domain%y%data%end
1857         is=isc; ie=iec; js=jsc; je=jec
1858         domain%list(m)%recv_n%folded = .FALSE.
1859         if( jsd.GT.jeg )then
1860             if( domain%y%cyclic )then !try cyclic offset
1861                 js = js+joff; je = je+joff
1862             else if( BTEST(domain%fold,NORTH) )then
1863                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1864                 j=js; js = 2*jeg-je+1; je = 2*jeg-j+1
1865                 domain%list(m)%recv_n%folded = .TRUE.
1866                 if( BTEST(grid_offset_type,NORTH) )then
1867                     js = js - 1; je = je - 1
1868                 end if
1869             end if
1870         end if
1871         is = max(isd,is); ie = min(ied,ie)
1872         js = max(jsd,js); je = min(jed,je)
1873         if( ie.GE.is .AND. je.GE.js )then
1874             domain%list(m)%overlap = .TRUE.
1875             if( grid_offset_type.NE.AGRID )then
1876                 domain%list(m)%recv_n_off%overlap = .TRUE.
1877                 domain%list(m)%recv_n_off%is = is
1878                 domain%list(m)%recv_n_off%ie = ie
1879                 domain%list(m)%recv_n_off%js = js
1880                 domain%list(m)%recv_n_off%je = je
1881             else
1882                 domain%list(m)%recv_n%overlap = .TRUE.
1883                 domain%list(m)%recv_n%is = is
1884                 domain%list(m)%recv_n%ie = ie
1885                 domain%list(m)%recv_n%js = js
1886                 domain%list(m)%recv_n%je = je
1887             end if
1888         else
1889             domain%list(m)%recv_n%overlap = .FALSE.
1890             domain%list(m)%recv_n_off%overlap = .FALSE.
1891         end if
1892!recv_ne
1893         isd = domain%x%compute%end+1; ied = domain%x%data%end
1894         jsd = domain%y%compute%end+1; jed = domain%y%data%end
1895         is=isc; ie=iec; js=jsc; je=jec
1896         domain%list(m)%recv_ne%folded = .FALSE.
1897         if( jsd.GT.jeg )then
1898             if( domain%y%cyclic )then !try cyclic offset
1899                 js = js+joff; je = je+joff
1900             else if( BTEST(domain%fold,NORTH) )then
1901                 i=is; is = isg+ieg-ie; ie = isg+ieg-i
1902                 j=js; js = 2*jeg-je+1; je = 2*jeg-j+1
1903                 domain%list(m)%recv_ne%folded = .TRUE.
1904                 if( BTEST(grid_offset_type,NORTH) )then
1905                     js = js - 1; je = je - 1
1906                 end if
1907             end if
1908         end if
1909         if( isd.GT.ieg )then
1910             if( domain%x%cyclic )then !try cyclic offset
1911                 is = is+ioff; ie = ie+ioff
1912             else if( BTEST(domain%fold,EAST) )then
1913                 i=is; is = 2*ieg-ie+1; ie = 2*ieg-i+1
1914                 j=js; js = jsg+jeg-je; je = jsg+jeg-j
1915                 domain%list(m)%recv_ne%folded = .TRUE.
1916                 if( BTEST(grid_offset_type,EAST) )then
1917                     is = is - 1; ie = ie - 1
1918                 end if
1919             end if
1920         end if
1921         is = max(isd,is); ie = min(ied,ie)
1922         js = max(jsd,js); je = min(jed,je)
1923         if( ie.GE.is .AND. je.GE.js )then
1924             domain%list(m)%overlap = .TRUE.
1925             if( grid_offset_type.NE.AGRID )then
1926                 domain%list(m)%recv_ne_off%overlap = .TRUE.
1927                 domain%list(m)%recv_ne_off%is = is
1928                 domain%list(m)%recv_ne_off%ie = ie
1929                 domain%list(m)%recv_ne_off%js = js
1930                 domain%list(m)%recv_ne_off%je = je
1931             else
1932                 domain%list(m)%recv_ne%overlap = .TRUE.
1933                 domain%list(m)%recv_ne%is = is
1934                 domain%list(m)%recv_ne%ie = ie
1935                 domain%list(m)%recv_ne%js = js
1936                 domain%list(m)%recv_ne%je = je
1937             end if
1938         else
1939             domain%list(m)%recv_ne%overlap = .FALSE.
1940             domain%list(m)%recv_ne_off%overlap = .FALSE.
1941         end if
1942      end do
1943      if( grid_offset_type.EQ.AGRID )domain%remote_domains_initialized = .TRUE.
1944      if( grid_offset_type.NE.AGRID )domain%remote_off_domains_initialized = .TRUE.
1945      return
1946    end subroutine compute_overlaps
1947
1948    subroutine mpp_define_layout2D( global_indices, ndivs, layout )
1949      integer, intent(in) :: global_indices(4) !(/ isg, ieg, jsg, jeg /)
1950      integer, intent(in) :: ndivs !number of divisions to divide global domain
1951      integer, intent(out) :: layout(2)
1952
1953      integer :: isg, ieg, jsg, jeg, isz, jsz, idiv, jdiv
1954
1955      isg = global_indices(1)
1956      ieg = global_indices(2)
1957      jsg = global_indices(3)
1958      jeg = global_indices(4)
1959
1960      isz = ieg - isg + 1
1961      jsz = jeg - jsg + 1
1962!first try to divide ndivs in the domain aspect ratio: if imperfect aspect, reduce idiv till it divides ndivs
1963      idiv = nint( sqrt(float(ndivs*isz)/jsz) )
1964      idiv = max(idiv,1) !for isz=1 line above can give 0
1965      do while( mod(ndivs,idiv).NE.0 )
1966         idiv = idiv - 1
1967      end do                 !will terminate at idiv=1 if not before
1968      jdiv = ndivs/idiv
1969
1970      layout = (/ idiv, jdiv /)
1971      return
1972    end subroutine mpp_define_layout2D
1973
1974!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1975!                                                                             !
1976!     MPP_GET and SET routiness: retrieve various components of domains       !
1977!                                                                             !
1978!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1979
1980    subroutine mpp_get_compute_domain1D( domain, begin, end, size, max_size, is_global )
1981      type(domain1D), intent(in) :: domain
1982      integer, intent(out), optional :: begin, end, size, max_size
1983      logical, intent(out), optional :: is_global
1984
1985      if( PRESENT(begin)     )begin     = domain%compute%begin
1986      if( PRESENT(end)       )end       = domain%compute%end
1987      if( PRESENT(size)      )size      = domain%compute%size
1988      if( PRESENT(max_size)  )max_size  = domain%compute%max_size
1989      if( PRESENT(is_global) )is_global = domain%compute%is_global
1990      return
1991    end subroutine mpp_get_compute_domain1D
1992
1993    subroutine mpp_get_data_domain1D( domain, begin, end, size, max_size, is_global )
1994      type(domain1D), intent(in) :: domain
1995      integer, intent(out), optional :: begin, end, size, max_size
1996      logical, intent(out), optional :: is_global
1997
1998      if( PRESENT(begin)     )begin     = domain%data%begin
1999      if( PRESENT(end)       )end       = domain%data%end
2000      if( PRESENT(size)      )size      = domain%data%size
2001      if( PRESENT(max_size)  )max_size  = domain%data%max_size
2002      if( PRESENT(is_global) )is_global = domain%data%is_global
2003      return
2004    end subroutine mpp_get_data_domain1D
2005
2006    subroutine mpp_get_global_domain1D( domain, begin, end, size, max_size )
2007      type(domain1D), intent(in) :: domain
2008      integer, intent(out), optional :: begin, end, size, max_size
2009
2010      if( PRESENT(begin)    )begin    = domain%global%begin
2011      if( PRESENT(end)      )end      = domain%global%end
2012      if( PRESENT(size)     )size     = domain%global%size
2013      if( PRESENT(max_size) )max_size = domain%global%max_size
2014      return
2015    end subroutine mpp_get_global_domain1D
2016
2017    subroutine mpp_get_compute_domain2D( domain, xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size, &
2018         x_is_global, y_is_global )
2019      type(domain2D), intent(in) :: domain
2020      integer, intent(out), optional :: xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size
2021      logical, intent(out), optional :: x_is_global, y_is_global
2022      call mpp_get_compute_domain( domain%x, xbegin, xend, xsize, xmax_size, x_is_global )
2023      call mpp_get_compute_domain( domain%y, ybegin, yend, ysize, ymax_size, y_is_global )
2024      return
2025    end subroutine mpp_get_compute_domain2D
2026
2027    subroutine mpp_get_data_domain2D( domain, xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size, &
2028         x_is_global, y_is_global )
2029      type(domain2D), intent(in) :: domain
2030      integer, intent(out), optional :: xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size
2031      logical, intent(out), optional :: x_is_global, y_is_global
2032      call mpp_get_data_domain( domain%x, xbegin, xend, xsize, xmax_size, x_is_global )
2033      call mpp_get_data_domain( domain%y, ybegin, yend, ysize, ymax_size, y_is_global )
2034      return
2035    end subroutine mpp_get_data_domain2D
2036
2037    subroutine mpp_get_global_domain2D( domain, xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size )
2038      type(domain2D), intent(in) :: domain
2039      integer, intent(out), optional :: xbegin, xend, ybegin, yend, xsize, xmax_size, ysize, ymax_size
2040      call mpp_get_global_domain( domain%x, xbegin, xend, xsize, xmax_size )
2041      call mpp_get_global_domain( domain%y, ybegin, yend, ysize, ymax_size )
2042      return
2043    end subroutine mpp_get_global_domain2D
2044
2045    subroutine mpp_get_domain_components( domain, x, y )
2046      type(domain2D), intent(in) :: domain
2047      type(domain1D), intent(out), optional :: x, y
2048      if( PRESENT(x) )x = domain%x
2049      if( PRESENT(y) )y = domain%y
2050      return
2051    end subroutine mpp_get_domain_components
2052
2053    subroutine mpp_get_compute_domains1D( domain, begin, end, size )
2054      type(domain1D), intent(in) :: domain
2055      integer, intent(out), optional, dimension(:) :: begin, end, size 
2056
2057      if( .NOT.module_is_initialized ) &
2058           call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: must first call mpp_domains_init.' )
2059!we use shape instead of size for error checks because size is used as an argument
2060      if( PRESENT(begin) )then
2061          if( any(shape(begin).NE.shape(domain%list)) ) &
2062               call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: begin array size does not match domain.' )
2063          begin(:) = domain%list(:)%compute%begin
2064      end if
2065      if( PRESENT(end) )then
2066          if( any(shape(end).NE.shape(domain%list)) ) &
2067               call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: end array size does not match domain.' )
2068          end(:) = domain%list(:)%compute%end
2069      end if
2070      if( PRESENT(size) )then
2071          if( any(shape(size).NE.shape(domain%list)) ) &
2072               call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: size array size does not match domain.' )
2073          size(:) = domain%list(:)%compute%size
2074      end if
2075      return
2076    end subroutine mpp_get_compute_domains1D
2077
2078    subroutine mpp_get_compute_domains2D( domain, xbegin, xend, xsize, ybegin, yend, ysize )
2079      type(domain2D), intent(in) :: domain
2080      integer, intent(out), optional, dimension(:) :: xbegin, xend, xsize, ybegin, yend, ysize
2081
2082      if( .NOT.module_is_initialized ) &
2083           call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: must first call mpp_domains_init.' )
2084
2085      if( PRESENT(xbegin) )then
2086          if( size(xbegin).NE.size(domain%list) ) &
2087               call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: xbegin array size does not match domain.' )
2088          xbegin(:) = domain%list(:)%x%compute%begin
2089      end if
2090      if( PRESENT(xend) )then
2091          if( size(xend).NE.size(domain%list) ) &
2092               call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: xend array size does not match domain.' )
2093          xend(:) = domain%list(:)%x%compute%end
2094      end if
2095      if( PRESENT(xsize) )then
2096          if( size(xsize).NE.size(domain%list) ) &
2097               call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: xsize array size does not match domain.' )
2098          xsize(:) = domain%list(:)%x%compute%size
2099      end if
2100      if( PRESENT(ybegin) )then
2101          if( size(ybegin).NE.size(domain%list) ) &
2102               call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: ybegin array size does not match domain.' )
2103          ybegin(:) = domain%list(:)%y%compute%begin
2104      end if
2105      if( PRESENT(yend) )then
2106          if( size(yend).NE.size(domain%list) ) &
2107               call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: yend array size does not match domain.' )
2108          yend(:) = domain%list(:)%y%compute%end
2109      end if
2110      if( PRESENT(ysize) )then
2111          if( size(ysize).NE.size(domain%list) ) &
2112               call mpp_error( FATAL, 'MPP_GET_COMPUTE_DOMAINS: ysize array size does not match domain.' )
2113          ysize(:) = domain%list(:)%y%compute%size
2114      end if
2115      return
2116    end subroutine mpp_get_compute_domains2D
2117
2118    subroutine mpp_get_pelist1D( domain, pelist, pos )
2119      type(domain1D), intent(in) :: domain
2120      integer, intent(out) :: pelist(:)
2121      integer, intent(out), optional :: pos
2122      integer :: ndivs
2123
2124      if( .NOT.module_is_initialized ) &
2125           call mpp_error( FATAL, 'MPP_GET_PELIST: must first call mpp_domains_init.' )
2126      ndivs = size(domain%list)
2127     
2128      if( size(pelist).NE.ndivs ) &
2129           call mpp_error( FATAL, 'MPP_GET_PELIST: pelist array size does not match domain.' )
2130
2131      pelist(:) = domain%list(0:ndivs-1)%pe
2132      if( PRESENT(pos) )pos = domain%pos
2133      return
2134    end subroutine mpp_get_pelist1D
2135
2136    subroutine mpp_get_pelist2D( domain, pelist, pos )
2137      type(domain2D), intent(in) :: domain
2138      integer, intent(out) :: pelist(:)
2139      integer, intent(out), optional :: pos
2140
2141      if( .NOT.module_is_initialized ) &
2142           call mpp_error( FATAL, 'MPP_GET_PELIST: must first call mpp_domains_init.' )
2143      if( size(pelist).NE.size(domain%list) ) &
2144           call mpp_error( FATAL, 'MPP_GET_PELIST: pelist array size does not match domain.' )
2145
2146      pelist(:) = domain%list(:)%pe
2147      if( PRESENT(pos) )pos = domain%pos
2148      return
2149    end subroutine mpp_get_pelist2D
2150
2151    subroutine mpp_get_layout1D( domain, layout )
2152      type(domain1D), intent(in) :: domain
2153      integer, intent(out) :: layout
2154     
2155      if( .NOT.module_is_initialized ) &
2156           call mpp_error( FATAL, 'MPP_GET_LAYOUT: must first call mpp_domains_init.' )
2157
2158      layout = size(domain%list)
2159      return
2160    end subroutine mpp_get_layout1D
2161
2162    subroutine mpp_get_layout2D( domain, layout )
2163      type(domain2D), intent(in) :: domain
2164      integer, intent(out) :: layout(2)
2165     
2166      if( .NOT.module_is_initialized ) &
2167           call mpp_error( FATAL, 'MPP_GET_LAYOUT: must first call mpp_domains_init.' )
2168
2169      layout(1) = size(domain%x%list)
2170      layout(2) = size(domain%y%list)
2171      return
2172    end subroutine mpp_get_layout2D
2173
2174!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2175!                                                                             !
2176!              MPP_UPDATE_DOMAINS: fill halos for 2D decomposition            !
2177!                                                                             !
2178!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2179
2180#define VECTOR_FIELD_
2181#define MPP_TYPE_ real(DOUBLE_KIND)
2182#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain2D_r8_2D
2183#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain2D_r8_3D
2184#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain2D_r8_4D
2185#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain2D_r8_5D
2186
2187#ifdef  VECTOR_FIELD_
2188#define MPP_UPDATE_DOMAINS_2D_V_ mpp_update_domain2D_r8_2Dv
2189#define MPP_UPDATE_DOMAINS_3D_V_ mpp_update_domain2D_r8_3Dv
2190#define MPP_UPDATE_DOMAINS_4D_V_ mpp_update_domain2D_r8_4Dv
2191#define MPP_UPDATE_DOMAINS_5D_V_ mpp_update_domain2D_r8_5Dv
2192#endif
2193#define MPP_REDISTRIBUTE_2D_ mpp_redistribute_r8_2D
2194#define MPP_REDISTRIBUTE_3D_ mpp_redistribute_r8_3D
2195#define MPP_REDISTRIBUTE_4D_ mpp_redistribute_r8_4D
2196#define MPP_REDISTRIBUTE_5D_ mpp_redistribute_r8_5D
2197#include <mpp_update_domains2D.h>
2198#undef  VECTOR_FIELD_
2199
2200#define MPP_TYPE_ complex(DOUBLE_KIND)
2201#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain2D_c8_2D
2202#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain2D_c8_3D
2203#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain2D_c8_4D
2204#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain2D_c8_5D
2205#define MPP_REDISTRIBUTE_2D_ mpp_redistribute_c8_2D
2206#define MPP_REDISTRIBUTE_3D_ mpp_redistribute_c8_3D
2207#define MPP_REDISTRIBUTE_4D_ mpp_redistribute_c8_4D
2208#define MPP_REDISTRIBUTE_5D_ mpp_redistribute_c8_5D
2209#include <mpp_update_domains2D.h>
2210
2211#ifndef no_8byte_integers
2212#define MPP_TYPE_ integer(LONG_KIND)
2213#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain2D_i8_2D
2214#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain2D_i8_3D
2215#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain2D_i8_4D
2216#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain2D_i8_5D
2217#define MPP_REDISTRIBUTE_2D_ mpp_redistribute_i8_2D
2218#define MPP_REDISTRIBUTE_3D_ mpp_redistribute_i8_3D
2219#define MPP_REDISTRIBUTE_4D_ mpp_redistribute_i8_4D
2220#define MPP_REDISTRIBUTE_5D_ mpp_redistribute_i8_5D
2221#include <mpp_update_domains2D.h>
2222
2223#define MPP_TYPE_ logical(LONG_KIND)
2224#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain2D_l8_2D
2225#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain2D_l8_3D
2226#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain2D_l8_4D
2227#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain2D_l8_5D
2228#define MPP_REDISTRIBUTE_2D_ mpp_redistribute_l8_2D
2229#define MPP_REDISTRIBUTE_3D_ mpp_redistribute_l8_3D
2230#define MPP_REDISTRIBUTE_4D_ mpp_redistribute_l8_4D
2231#define MPP_REDISTRIBUTE_5D_ mpp_redistribute_l8_5D
2232#include <mpp_update_domains2D.h>
2233#endif
2234
2235#ifndef no_4byte_reals
2236#define VECTOR_FIELD_
2237#define MPP_TYPE_ real(FLOAT_KIND)
2238#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain2D_r4_2D
2239#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain2D_r4_3D
2240#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain2D_r4_4D
2241#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain2D_r4_5D
2242#ifdef  VECTOR_FIELD_
2243#define MPP_UPDATE_DOMAINS_2D_V_ mpp_update_domain2D_r4_2Dv
2244#define MPP_UPDATE_DOMAINS_3D_V_ mpp_update_domain2D_r4_3Dv
2245#define MPP_UPDATE_DOMAINS_4D_V_ mpp_update_domain2D_r4_4Dv
2246#define MPP_UPDATE_DOMAINS_5D_V_ mpp_update_domain2D_r4_5Dv
2247#endif
2248#define MPP_REDISTRIBUTE_2D_ mpp_redistribute_r4_2D
2249#define MPP_REDISTRIBUTE_3D_ mpp_redistribute_r4_3D
2250#define MPP_REDISTRIBUTE_4D_ mpp_redistribute_r4_4D
2251#define MPP_REDISTRIBUTE_5D_ mpp_redistribute_r4_5D
2252#include <mpp_update_domains2D.h>
2253#undef  VECTOR_FIELD_
2254#endif
2255
2256#ifndef no_4byte_cmplx
2257#define MPP_TYPE_ complex(FLOAT_KIND)
2258#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain2D_c4_2D
2259#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain2D_c4_3D
2260#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain2D_c4_4D
2261#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain2D_c4_5D
2262#define MPP_REDISTRIBUTE_2D_ mpp_redistribute_c4_2D
2263#define MPP_REDISTRIBUTE_3D_ mpp_redistribute_c4_3D
2264#define MPP_REDISTRIBUTE_4D_ mpp_redistribute_c4_4D
2265#define MPP_REDISTRIBUTE_5D_ mpp_redistribute_c4_5D
2266#include <mpp_update_domains2D.h>
2267#endif
2268
2269#define MPP_TYPE_ integer(INT_KIND)
2270#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain2D_i4_2D
2271#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain2D_i4_3D
2272#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain2D_i4_4D
2273#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain2D_i4_5D
2274#define MPP_REDISTRIBUTE_2D_ mpp_redistribute_i4_2D
2275#define MPP_REDISTRIBUTE_3D_ mpp_redistribute_i4_3D
2276#define MPP_REDISTRIBUTE_4D_ mpp_redistribute_i4_4D
2277#define MPP_REDISTRIBUTE_5D_ mpp_redistribute_i4_5D
2278#include <mpp_update_domains2D.h>
2279
2280#define MPP_TYPE_ logical(INT_KIND)
2281#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain2D_l4_2D
2282#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain2D_l4_3D
2283#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain2D_l4_4D
2284#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain2D_l4_5D
2285#define MPP_REDISTRIBUTE_2D_ mpp_redistribute_l4_2D
2286#define MPP_REDISTRIBUTE_3D_ mpp_redistribute_l4_3D
2287#define MPP_REDISTRIBUTE_4D_ mpp_redistribute_l4_4D
2288#define MPP_REDISTRIBUTE_5D_ mpp_redistribute_l4_5D
2289#include <mpp_update_domains2D.h>
2290
2291!#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain1D_r8_2D
2292!#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain1D_r8_3D
2293!#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain1D_r8_4D
2294!#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain1D_r8_5D
2295!#define MPP_TYPE_ real(DOUBLE_KIND)
2296!#include <mpp_update_domains1D.h>
2297!
2298!#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain1D_c8_2D
2299!#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain1D_c8_3D
2300!#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain1D_c8_4D
2301!#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain1D_c8_5D
2302!#define MPP_TYPE_ complex(DOUBLE_KIND)
2303!#include <mpp_update_domains1D.h>
2304!
2305!#ifndef no_8byte_integers
2306!#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain1D_i8_2D
2307!#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain1D_i8_3D
2308!#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain1D_i8_4D
2309!#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain1D_i8_5D
2310!#define MPP_TYPE_ integer(LONG_KIND)
2311!#include <mpp_update_domains1D.h>
2312!
2313!#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain1D_l8_2D
2314!#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain1D_l8_3D
2315!#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain1D_l8_4D
2316!#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain1D_l8_5D
2317!#define MPP_TYPE_ logical(LONG_KIND)
2318!#include <mpp_update_domains1D.h>
2319!#endif
2320!
2321!#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain1D_r4_2D
2322!#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain1D_r4_3D
2323!#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain1D_r4_4D
2324!#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain1D_r4_5D
2325!#define MPP_TYPE_ real(FLOAT_KIND)
2326!#include <mpp_update_domains1D.h>
2327!
2328!#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain1D_c4_2D
2329!#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain1D_c4_3D
2330!#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain1D_c4_4D
2331!#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain1D_c4_5D
2332!#define MPP_TYPE_ complex(FLOAT_KIND)
2333!#include <mpp_update_domains1D.h>
2334!
2335!#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain1D_i4_2D
2336!#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain1D_i4_3D
2337!#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain1D_i4_4D
2338!#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain1D_i4_5D
2339!#define MPP_TYPE_ integer(INT_KIND)
2340!#include <mpp_update_domains1D.h>
2341!
2342!#define MPP_UPDATE_DOMAINS_2D_ mpp_update_domain1D_l4_2D
2343!#define MPP_UPDATE_DOMAINS_3D_ mpp_update_domain1D_l4_3D
2344!#define MPP_UPDATE_DOMAINS_4D_ mpp_update_domain1D_l4_4D
2345!#define MPP_UPDATE_DOMAINS_5D_ mpp_update_domain1D_l4_5D
2346!#define MPP_TYPE_ logical(INT_KIND)
2347!#include <mpp_update_domains1D.h>
2348
2349!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2350!                                                                             !
2351!              MPP_GLOBAL_REDUCE: get global max/min of field                 !
2352!                                                                             !
2353!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2354
2355#define MPP_GLOBAL_REDUCE_2D_ mpp_global_max_r8_2d
2356#define MPP_GLOBAL_REDUCE_3D_ mpp_global_max_r8_3d
2357#define MPP_GLOBAL_REDUCE_4D_ mpp_global_max_r8_4d
2358#define MPP_GLOBAL_REDUCE_5D_ mpp_global_max_r8_5d
2359#define MPP_TYPE_ real(DOUBLE_KIND)
2360#define REDUCE_VAL_ maxval
2361#define REDUCE_LOC_ maxloc
2362#define MPP_REDUCE_ mpp_max
2363#include <mpp_global_reduce.h>
2364
2365#define MPP_GLOBAL_REDUCE_2D_ mpp_global_min_r8_2d
2366#define MPP_GLOBAL_REDUCE_3D_ mpp_global_min_r8_3d
2367#define MPP_GLOBAL_REDUCE_4D_ mpp_global_min_r8_4d
2368#define MPP_GLOBAL_REDUCE_5D_ mpp_global_min_r8_5d
2369#define MPP_TYPE_ real(DOUBLE_KIND)
2370#define REDUCE_VAL_ minval
2371#define REDUCE_LOC_ minloc
2372#define MPP_REDUCE_ mpp_min
2373#include <mpp_global_reduce.h>
2374
2375#ifndef no_4byte_reals
2376#define MPP_GLOBAL_REDUCE_2D_ mpp_global_max_r4_2d
2377#define MPP_GLOBAL_REDUCE_3D_ mpp_global_max_r4_3d
2378#define MPP_GLOBAL_REDUCE_4D_ mpp_global_max_r4_4d
2379#define MPP_GLOBAL_REDUCE_5D_ mpp_global_max_r4_5d
2380#define MPP_TYPE_ real(FLOAT_KIND)
2381#define REDUCE_VAL_ maxval
2382#define REDUCE_LOC_ maxloc
2383#define MPP_REDUCE_ mpp_max
2384#include <mpp_global_reduce.h>
2385
2386#define MPP_GLOBAL_REDUCE_2D_ mpp_global_min_r4_2d
2387#define MPP_GLOBAL_REDUCE_3D_ mpp_global_min_r4_3d
2388#define MPP_GLOBAL_REDUCE_4D_ mpp_global_min_r4_4d
2389#define MPP_GLOBAL_REDUCE_5D_ mpp_global_min_r4_5d
2390#define MPP_TYPE_ real(FLOAT_KIND)
2391#define REDUCE_VAL_ minval
2392#define REDUCE_LOC_ minloc
2393#define MPP_REDUCE_ mpp_min
2394#include <mpp_global_reduce.h>
2395#endif
2396
2397#ifndef no_8byte_integers
2398#define MPP_GLOBAL_REDUCE_2D_ mpp_global_max_i8_2d
2399#define MPP_GLOBAL_REDUCE_3D_ mpp_global_max_i8_3d
2400#define MPP_GLOBAL_REDUCE_4D_ mpp_global_max_i8_4d
2401#define MPP_GLOBAL_REDUCE_5D_ mpp_global_max_i8_5d
2402#define MPP_TYPE_ integer(LONG_KIND)
2403#define REDUCE_VAL_ maxval
2404#define REDUCE_LOC_ maxloc
2405#define MPP_REDUCE_ mpp_max
2406#include <mpp_global_reduce.h>
2407
2408#define MPP_GLOBAL_REDUCE_2D_ mpp_global_min_i8_2d
2409#define MPP_GLOBAL_REDUCE_3D_ mpp_global_min_i8_3d
2410#define MPP_GLOBAL_REDUCE_4D_ mpp_global_min_i8_4d
2411#define MPP_GLOBAL_REDUCE_5D_ mpp_global_min_i8_5d
2412#define MPP_TYPE_ integer(LONG_KIND)
2413#define REDUCE_VAL_ minval
2414#define REDUCE_LOC_ minloc
2415#define MPP_REDUCE_ mpp_min
2416#include <mpp_global_reduce.h>
2417#endif
2418
2419#define MPP_GLOBAL_REDUCE_2D_ mpp_global_max_i4_2d
2420#define MPP_GLOBAL_REDUCE_3D_ mpp_global_max_i4_3d
2421#define MPP_GLOBAL_REDUCE_4D_ mpp_global_max_i4_4d
2422#define MPP_GLOBAL_REDUCE_5D_ mpp_global_max_i4_5d
2423#define MPP_TYPE_ integer(INT_KIND)
2424#define REDUCE_VAL_ maxval
2425#define REDUCE_LOC_ maxloc
2426#define MPP_REDUCE_ mpp_max
2427#include <mpp_global_reduce.h>
2428
2429#define MPP_GLOBAL_REDUCE_2D_ mpp_global_min_i4_2d
2430#define MPP_GLOBAL_REDUCE_3D_ mpp_global_min_i4_3d
2431#define MPP_GLOBAL_REDUCE_4D_ mpp_global_min_i4_4d
2432#define MPP_GLOBAL_REDUCE_5D_ mpp_global_min_i4_5d
2433#define MPP_TYPE_ integer(INT_KIND)
2434#define REDUCE_VAL_ minval
2435#define REDUCE_LOC_ minloc
2436#define MPP_REDUCE_ mpp_min
2437#include <mpp_global_reduce.h>
2438
2439!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2440!                                                                             !
2441!                   MPP_GLOBAL_SUM: global sum of field                       !
2442!                                                                             !
2443!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2444
2445#define MPP_GLOBAL_SUM_ mpp_global_sum_r8_2d
2446#define MPP_EXTRA_INDICES_
2447#define MPP_TYPE_ real(DOUBLE_KIND)
2448#include <mpp_global_sum.h>
2449
2450#define MPP_GLOBAL_SUM_ mpp_global_sum_r8_3d
2451#define MPP_EXTRA_INDICES_ ,:
2452#define MPP_TYPE_ real(DOUBLE_KIND)
2453#include <mpp_global_sum.h>
2454
2455#define MPP_GLOBAL_SUM_ mpp_global_sum_r8_4d
2456#define MPP_EXTRA_INDICES_ ,:,:
2457#define MPP_TYPE_ real(DOUBLE_KIND)
2458#include <mpp_global_sum.h>
2459
2460#define MPP_GLOBAL_SUM_ mpp_global_sum_r8_5d
2461#define MPP_EXTRA_INDICES_ ,:,:,:
2462#define MPP_TYPE_ real(DOUBLE_KIND)
2463#include <mpp_global_sum.h>
2464
2465#ifndef no_4byte_reals
2466#define MPP_GLOBAL_SUM_ mpp_global_sum_r4_2d
2467#define MPP_EXTRA_INDICES_
2468#define MPP_TYPE_ real(FLOAT_KIND)
2469#include <mpp_global_sum.h>
2470
2471#define MPP_GLOBAL_SUM_ mpp_global_sum_r4_3d
2472#define MPP_EXTRA_INDICES_ ,:
2473#define MPP_TYPE_ real(FLOAT_KIND)
2474#include <mpp_global_sum.h>
2475
2476#define MPP_GLOBAL_SUM_ mpp_global_sum_r4_4d
2477#define MPP_EXTRA_INDICES_ ,:,:
2478#define MPP_TYPE_ real(FLOAT_KIND)
2479#include <mpp_global_sum.h>
2480
2481#define MPP_GLOBAL_SUM_ mpp_global_sum_r4_5d
2482#define MPP_EXTRA_INDICES_ ,:,:,:
2483#define MPP_TYPE_ real(FLOAT_KIND)
2484#include <mpp_global_sum.h>
2485#endif
2486
2487#define MPP_GLOBAL_SUM_ mpp_global_sum_c8_2d
2488#define MPP_EXTRA_INDICES_
2489#define MPP_TYPE_ complex(DOUBLE_KIND)
2490#include <mpp_global_sum.h>
2491
2492#define MPP_GLOBAL_SUM_ mpp_global_sum_c8_3d
2493#define MPP_EXTRA_INDICES_ ,:
2494#define MPP_TYPE_ complex(DOUBLE_KIND)
2495#include <mpp_global_sum.h>
2496
2497#define MPP_GLOBAL_SUM_ mpp_global_sum_c8_4d
2498#define MPP_EXTRA_INDICES_ ,:,:
2499#define MPP_TYPE_ complex(DOUBLE_KIND)
2500#include <mpp_global_sum.h>
2501
2502#define MPP_GLOBAL_SUM_ mpp_global_sum_c8_5d
2503#define MPP_EXTRA_INDICES_ ,:,:,:
2504#define MPP_TYPE_ complex(DOUBLE_KIND)
2505#include <mpp_global_sum.h>
2506
2507#ifndef no_4byte_cmplx
2508#define MPP_GLOBAL_SUM_ mpp_global_sum_c4_2d
2509#define MPP_EXTRA_INDICES_
2510#define MPP_TYPE_ complex(FLOAT_KIND)
2511#include <mpp_global_sum.h>
2512
2513#define MPP_GLOBAL_SUM_ mpp_global_sum_c4_3d
2514#define MPP_EXTRA_INDICES_ ,:
2515#define MPP_TYPE_ complex(FLOAT_KIND)
2516#include <mpp_global_sum.h>
2517
2518#define MPP_GLOBAL_SUM_ mpp_global_sum_c4_4d
2519#define MPP_EXTRA_INDICES_ ,:,:
2520#define MPP_TYPE_ complex(FLOAT_KIND)
2521#include <mpp_global_sum.h>
2522
2523#define MPP_GLOBAL_SUM_ mpp_global_sum_c4_5d
2524#define MPP_EXTRA_INDICES_ ,:,:,:
2525#define MPP_TYPE_ complex(FLOAT_KIND)
2526#include <mpp_global_sum.h>
2527#endif
2528
2529#ifndef no_8byte_integers
2530#define MPP_GLOBAL_SUM_ mpp_global_sum_i8_2d
2531#define MPP_EXTRA_INDICES_
2532#define MPP_TYPE_ integer(LONG_KIND)
2533#include <mpp_global_sum.h>
2534
2535#define MPP_GLOBAL_SUM_ mpp_global_sum_i8_3d
2536#define MPP_EXTRA_INDICES_ ,:
2537#define MPP_TYPE_ integer(LONG_KIND)
2538#include <mpp_global_sum.h>
2539
2540#define MPP_GLOBAL_SUM_ mpp_global_sum_i8_4d
2541#define MPP_EXTRA_INDICES_ ,:,:
2542#define MPP_TYPE_ integer(LONG_KIND)
2543#include <mpp_global_sum.h>
2544
2545#define MPP_GLOBAL_SUM_ mpp_global_sum_i8_5d
2546#define MPP_EXTRA_INDICES_ ,:,:,:
2547#define MPP_TYPE_ integer(LONG_KIND)
2548#include <mpp_global_sum.h>
2549#endif
2550
2551#define MPP_GLOBAL_SUM_ mpp_global_sum_i4_2d
2552#define MPP_EXTRA_INDICES_
2553#define MPP_TYPE_ integer(INT_KIND)
2554#include <mpp_global_sum.h>
2555
2556#define MPP_GLOBAL_SUM_ mpp_global_sum_i4_3d
2557#define MPP_EXTRA_INDICES_ ,:
2558#define MPP_TYPE_ integer(INT_KIND)
2559#include <mpp_global_sum.h>
2560
2561#define MPP_GLOBAL_SUM_ mpp_global_sum_i4_4d
2562#define MPP_EXTRA_INDICES_ ,:,:
2563#define MPP_TYPE_ integer(INT_KIND)
2564#include <mpp_global_sum.h>
2565
2566#define MPP_GLOBAL_SUM_ mpp_global_sum_i4_5d
2567#define MPP_EXTRA_INDICES_ ,:,:,:
2568#define MPP_TYPE_ integer(INT_KIND)
2569#include <mpp_global_sum.h>
2570   
2571!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2572!                                                                             !
2573!              MPP_GLOBAL_FIELD: get global field from domain field           !
2574!                                                                             !
2575!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2576
2577#define MPP_GLOBAL_FIELD_2D_ mpp_global_field2D_r8_2d
2578#define MPP_GLOBAL_FIELD_3D_ mpp_global_field2D_r8_3d
2579#define MPP_GLOBAL_FIELD_4D_ mpp_global_field2D_r8_4d
2580#define MPP_GLOBAL_FIELD_5D_ mpp_global_field2D_r8_5d
2581#define MPP_GLOBAL1D_FIELD_2D_ mpp_global_field1D_r8_2d
2582#define MPP_TYPE_ real(DOUBLE_KIND)
2583#include <mpp_global_field.h>
2584
2585#define MPP_GLOBAL_FIELD_2D_ mpp_global_field2D_c8_2d
2586#define MPP_GLOBAL_FIELD_3D_ mpp_global_field2D_c8_3d
2587#define MPP_GLOBAL_FIELD_4D_ mpp_global_field2D_c8_4d
2588#define MPP_GLOBAL_FIELD_5D_ mpp_global_field2D_c8_5d
2589#define MPP_GLOBAL1D_FIELD_2D_ mpp_global_field1D_c8_2d
2590#define MPP_TYPE_ complex(DOUBLE_KIND)
2591#include <mpp_global_field.h>
2592
2593#ifndef no_8byte_integers
2594#define MPP_GLOBAL_FIELD_2D_ mpp_global_field2D_i8_2d
2595#define MPP_GLOBAL_FIELD_3D_ mpp_global_field2D_i8_3d
2596#define MPP_GLOBAL_FIELD_4D_ mpp_global_field2D_i8_4d
2597#define MPP_GLOBAL_FIELD_5D_ mpp_global_field2D_i8_5d
2598#define MPP_GLOBAL1D_FIELD_2D_ mpp_global_field1D_i8_2d
2599#define MPP_TYPE_ integer(LONG_KIND)
2600#include <mpp_global_field.h>
2601
2602#define MPP_GLOBAL_FIELD_2D_ mpp_global_field2D_l8_2d
2603#define MPP_GLOBAL_FIELD_3D_ mpp_global_field2D_l8_3d
2604#define MPP_GLOBAL_FIELD_4D_ mpp_global_field2D_l8_4d
2605#define MPP_GLOBAL_FIELD_5D_ mpp_global_field2D_l8_5d
2606#define MPP_GLOBAL1D_FIELD_2D_ mpp_global_field1D_l8_2d
2607#define MPP_TYPE_ logical(LONG_KIND)
2608#include <mpp_global_field.h>
2609#endif
2610
2611#ifndef no_4byte_reals
2612#define MPP_GLOBAL_FIELD_2D_ mpp_global_field2D_r4_2d
2613#define MPP_GLOBAL_FIELD_3D_ mpp_global_field2D_r4_3d
2614#define MPP_GLOBAL_FIELD_4D_ mpp_global_field2D_r4_4d
2615#define MPP_GLOBAL_FIELD_5D_ mpp_global_field2D_r4_5d
2616#define MPP_GLOBAL1D_FIELD_2D_ mpp_global_field1D_r4_2d
2617#define MPP_TYPE_ real(FLOAT_KIND)
2618#include <mpp_global_field.h>
2619#endif
2620
2621#ifndef no_4byte_cmplx
2622#define MPP_GLOBAL_FIELD_2D_ mpp_global_field2D_c4_2d
2623#define MPP_GLOBAL_FIELD_3D_ mpp_global_field2D_c4_3d
2624#define MPP_GLOBAL_FIELD_4D_ mpp_global_field2D_c4_4d
2625#define MPP_GLOBAL_FIELD_5D_ mpp_global_field2D_c4_5d
2626#define MPP_GLOBAL1D_FIELD_2D_ mpp_global_field1D_c4_2d
2627#define MPP_TYPE_ complex(FLOAT_KIND)
2628#include <mpp_global_field.h>
2629#endif
2630
2631#define MPP_GLOBAL_FIELD_2D_ mpp_global_field2D_i4_2d
2632#define MPP_GLOBAL_FIELD_3D_ mpp_global_field2D_i4_3d
2633#define MPP_GLOBAL_FIELD_4D_ mpp_global_field2D_i4_4d
2634#define MPP_GLOBAL_FIELD_5D_ mpp_global_field2D_i4_5d
2635#define MPP_GLOBAL1D_FIELD_2D_ mpp_global_field1D_i4_2d
2636#define MPP_TYPE_ integer(INT_KIND)
2637#include <mpp_global_field.h>
2638
2639#define MPP_GLOBAL_FIELD_2D_ mpp_global_field2D_l4_2d
2640#define MPP_GLOBAL_FIELD_3D_ mpp_global_field2D_l4_3d
2641#define MPP_GLOBAL_FIELD_4D_ mpp_global_field2D_l4_4d
2642#define MPP_GLOBAL_FIELD_5D_ mpp_global_field2D_l4_5d
2643#define MPP_GLOBAL1D_FIELD_2D_ mpp_global_field1D_l4_2d
2644#define MPP_TYPE_ logical(INT_KIND)
2645#include <mpp_global_field.h>
2646
2647end module mpp_domains_mod
2648
2649#ifdef test_mpp_domains
2650program mpp_domains_test
2651  use mpp_mod
2652  use mpp_domains_mod
2653  implicit none
2654  integer :: pe, npes
2655  integer :: nx=128, ny=128, nz=40, halo=2, stackmax=32768
2656  real, dimension(:,:,:), allocatable :: global, gcheck
2657  integer :: unit=7
2658  logical :: debug=.FALSE., opened
2659  namelist / mpp_domains_nml / nx, ny, nz, halo, stackmax, debug
2660  integer :: i, j, k
2661  integer :: layout(2)
2662  integer :: is, ie, js, je, isd, ied, jsd, jed
2663  integer :: id
2664
2665  call mpp_init()
2666
2667  call mpp_set_warn_level(FATAL)
2668!possibly open a file called mpp_domains.nml
2669  do
2670     inquire( unit=unit, opened=opened )
2671     if( .NOT.opened )exit
2672     unit = unit + 1
2673     if( unit.EQ.100 )call mpp_error( FATAL, 'Unable to locate unit number.' )
2674  end do
2675  open( unit=unit, status='OLD', file='mpp_domains.nml', err=10 )
2676  read( unit,mpp_domains_nml )
2677  close(unit)
267810 continue
2679
2680  pe = mpp_pe()
2681  npes = mpp_npes()
2682
2683  if( debug )then
2684      call mpp_domains_init(MPP_DEBUG)
2685  else
2686      call mpp_domains_init(MPP_DOMAIN_TIME)
2687  end if
2688  call mpp_domains_set_stack_size(stackmax)
2689
2690  if( pe.EQ.mpp_root_pe() )print '(a,5i4)', 'npes, nx, ny, nz, halo=', npes, nx, ny, nz, halo
2691
2692  allocate( global(1-halo:nx+halo,1-halo:ny+halo,nz) )
2693  allocate( gcheck(nx,ny,nz) )
2694
2695!fill in global array: with k.iiijjj
2696  do k = 1,nz
2697     do j = 1,ny
2698        do i = 1,nx
2699           global(i,j,k) = k + i*1e-3 + j*1e-6
2700        end do
2701     end do
2702  end do
2703
2704  call test_halo_update( 'Simple' ) !includes global field, global sum tests
2705  call test_halo_update( 'Cyclic' )
2706  call test_halo_update( 'Folded' ) !includes vector field test
2707
2708  call test_redistribute( 'Complete pelist' )
2709  call test_redistribute( 'Overlap  pelist' )
2710  call test_redistribute( 'Disjoint pelist' )
2711 
2712  call mpp_domains_exit()
2713  call mpp_exit()
2714
2715contains
2716
2717  subroutine test_redistribute( type )
2718!test redistribute between two domains
2719    character(len=*), intent(in) :: type
2720    type(domain2D) :: domainx, domainy
2721    real, allocatable, dimension(:,:,:) :: x, y
2722    integer, allocatable :: pelist(:)
2723    integer :: pemax
2724
2725    pemax = npes/2              !the partial pelist will run from 0...pemax
2726
2727!select pelists
2728    select case(type)
2729    case( 'Complete pelist' )
2730!both pelists run from 0...npes-1
2731        allocate( pelist(0:npes-1) )
2732        pelist = (/ (i,i=0,npes-1) /)
2733        call mpp_declare_pelist( pelist )
2734    case( 'Overlap  pelist' )
2735!one pelist from 0...pemax, other from 0...npes-1
2736        allocate( pelist(0:pemax) )
2737        pelist = (/ (i,i=0,pemax) /)
2738        call mpp_declare_pelist( pelist )
2739    case( 'Disjoint pelist' )
2740!one pelist from 0...pemax, other from pemax+1...npes-1
2741        if( pemax+1.GE.npes )return
2742        allocate( pelist(0:pemax) )
2743        pelist = (/ (i,i=0,pemax) /)
2744        call mpp_declare_pelist( pelist )
2745        call mpp_declare_pelist( (/ (i,i=pemax+1,npes-1) /))
2746    case default
2747        call mpp_error( FATAL, 'TEST_REDISTRIBUTE: no such test: '//type )
2748    end select
2749
2750!set up x and y arrays
2751    select case(type)
2752    case( 'Complete pelist' )
2753!set up x array
2754        call mpp_define_layout( (/1,nx,1,ny/), npes, layout )
2755        call mpp_define_domains( (/1,nx,1,ny/), layout, domainx, name=type )
2756        call mpp_get_compute_domain( domainx, is,  ie,  js,  je  )
2757        call mpp_get_data_domain   ( domainx, isd, ied, jsd, jed )
2758        allocate( x(isd:ied,jsd:jed,nz) )
2759        x = 0.
2760        x(is:ie,js:je,:) = global(is:ie,js:je,:)
2761!set up y array
2762        call mpp_define_domains( (/1,nx,1,ny/), (/npes,1/), domainy, name=type )
2763        call mpp_get_compute_domain( domainy, is,  ie,  js,  je  )
2764        call mpp_get_data_domain   ( domainy, isd, ied, jsd, jed )
2765        allocate( y(isd:ied,jsd:jed,nz) )
2766        y = 0.
2767    case( 'Overlap  pelist' )
2768!one pelist from 0...pemax, other from 0...npes-1
2769!set up x array
2770        call mpp_define_layout( (/1,nx,1,ny/), npes, layout )
2771        call mpp_define_domains( (/1,nx,1,ny/), layout, domainx, name=type )
2772        call mpp_get_compute_domain( domainx, is,  ie,  js,  je  )
2773        call mpp_get_data_domain   ( domainx, isd, ied, jsd, jed )
2774        allocate( x(isd:ied,jsd:jed,nz) )
2775        x = 0.
2776        x(is:ie,js:je,:) = global(is:ie,js:je,:)
2777!set up y array
2778        if( ANY(pelist.EQ.pe) )then
2779            call mpp_set_current_pelist(pelist)
2780            call mpp_define_layout( (/1,nx,1,ny/), mpp_npes(), layout )
2781            call mpp_define_domains( (/1,nx,1,ny/), layout, domainy, name=type )
2782            call mpp_get_compute_domain( domainy, is,  ie,  js,  je  )
2783            call mpp_get_data_domain   ( domainy, isd, ied, jsd, jed )
2784            allocate( y(isd:ied,jsd:jed,nz) )
2785            y = 0.
2786        end if
2787    case( 'Disjoint pelist' )
2788!one pelist from 0...pemax, other from pemax+1...npes-1
2789
2790!set up y array
2791        if( ANY(pelist.EQ.pe) )then
2792            call mpp_set_current_pelist(pelist)
2793            call mpp_define_layout( (/1,nx,1,ny/), mpp_npes(), layout )
2794            call mpp_define_domains( (/1,nx,1,ny/), layout, domainy, name=type )
2795            call mpp_get_compute_domain( domainy, is,  ie,  js,  je  )
2796            call mpp_get_data_domain   ( domainy, isd, ied, jsd, jed )
2797            allocate( y(isd:ied,jsd:jed,nz) )
2798            y = 0.
2799        else
2800!set up x array
2801            call mpp_set_current_pelist( (/ (i,i=pemax+1,npes-1) /) )
2802            call mpp_define_layout( (/1,nx,1,ny/), mpp_npes(), layout )
2803            call mpp_define_domains( (/1,nx,1,ny/), layout, domainx, name=type )
2804            call mpp_get_compute_domain( domainx, is,  ie,  js,  je  )
2805            call mpp_get_data_domain   ( domainx, isd, ied, jsd, jed )
2806            allocate( x(isd:ied,jsd:jed,nz) )
2807            x = 0.
2808            x(is:ie,js:je,:) = global(is:ie,js:je,:)
2809         end if
2810    end select
2811
2812!go global and redistribute
2813    call mpp_set_current_pelist()
2814    call mpp_broadcast_domain(domainx)
2815    call mpp_broadcast_domain(domainy)
2816
2817    id = mpp_clock_id( type, flags=MPP_CLOCK_SYNC+MPP_CLOCK_DETAILED )
2818    call mpp_clock_begin(id)
2819    call mpp_redistribute( domainx, x, domainy, y )
2820    call mpp_clock_end  (id)
2821
2822!check answers on pelist
2823    if( ANY(pelist.EQ.pe) )then
2824        call mpp_set_current_pelist(pelist)
2825        call mpp_global_field( domainy, y, gcheck )
2826        call compare_checksums( global(1:nx,1:ny,:), gcheck, type )
2827    end if
2828
2829    call mpp_set_current_pelist()
2830    call mpp_sync()
2831
2832    return
2833  end subroutine test_redistribute
2834
2835  subroutine test_halo_update( type )
2836    character(len=*), intent(in) :: type
2837    real, allocatable, dimension(:,:,:) :: x, y
2838    type(domain2D) :: domain
2839    real :: lsum, gsum
2840   
2841    call mpp_define_layout( (/1,nx,1,ny/), npes, layout )
2842    select case(type)
2843    case( 'Simple' )
2844        call mpp_define_domains( (/1,nx,1,ny/), layout, domain, xhalo=halo, yhalo=halo, name=type )
2845        global(1-halo:0,    :,:) = 0
2846        global(nx+1:nx+halo,:,:) = 0
2847        global(:,    1-halo:0,:) = 0
2848        global(:,ny+1:ny+halo,:) = 0
2849    case( 'Cyclic' )
2850        call mpp_define_domains( (/1,nx,1,ny/), layout, domain, xhalo=halo, yhalo=halo, &
2851             xflags=CYCLIC_GLOBAL_DOMAIN, yflags=CYCLIC_GLOBAL_DOMAIN, name=type )
2852        global(1-halo:0,    1:ny,:) = global(nx-halo+1:nx,1:ny,:)
2853        global(nx+1:nx+halo,1:ny,:) = global(1:halo,      1:ny,:)
2854        global(1-halo:nx+halo,    1-halo:0,:) = global(1-halo:nx+halo,ny-halo+1:ny,:)
2855        global(1-halo:nx+halo,ny+1:ny+halo,:) = global(1-halo:nx+halo,1:halo,      :)
2856    case( 'Folded' )
2857        call mpp_define_domains( (/1,nx,1,ny/), layout, domain, xhalo=halo, yhalo=halo, &
2858             xflags=CYCLIC_GLOBAL_DOMAIN, yflags=FOLD_NORTH_EDGE, name=type )
2859        global(1-halo:0,    1:ny,:) = global(nx-halo+1:nx,1:ny,:)
2860        global(nx+1:nx+halo,1:ny,:) = global(1:halo,      1:ny,:)
2861        global(1-halo:nx+halo,ny+1:ny+halo,:) = global(nx+halo:1-halo:-1,ny:ny-halo+1:-1,:)
2862        global(1-halo:nx+halo,1-halo:0,:) = 0
2863    case default
2864        call mpp_error( FATAL, 'TEST_MPP_DOMAINS: no such test: '//type )
2865    end select
2866
2867!set up x array
2868    call mpp_get_compute_domain( domain, is,  ie,  js,  je  )
2869    call mpp_get_data_domain   ( domain, isd, ied, jsd, jed )
2870    allocate( x(isd:ied,jsd:jed,nz) )
2871    x = 0.
2872    x(is:ie,js:je,:) = global(is:ie,js:je,:)
2873
2874!partial update
2875    id = mpp_clock_id( type//' partial', flags=MPP_CLOCK_SYNC+MPP_CLOCK_DETAILED )
2876    call mpp_clock_begin(id)
2877    call mpp_update_domains( x, domain, NUPDATE+EUPDATE )
2878    call mpp_clock_end  (id)
2879    call compare_checksums( x(is:ied,js:jed,:), global(is:ied,js:jed,:), type//' partial' )
2880
2881!full update
2882    id = mpp_clock_id( type, flags=MPP_CLOCK_SYNC+MPP_CLOCK_DETAILED )
2883    call mpp_clock_begin(id)
2884    call mpp_update_domains( x, domain )
2885    call mpp_clock_end  (id)
2886    call compare_checksums( x, global(isd:ied,jsd:jed,:), type )
2887
2888    select case(type)           !extra tests
2889    case( 'Simple' )
2890
2891!test mpp_global_field
2892        id = mpp_clock_id( type//' global field', flags=MPP_CLOCK_SYNC+MPP_CLOCK_DETAILED )
2893        call mpp_clock_begin(id)
2894        call mpp_global_field( domain, x, gcheck )
2895        call mpp_clock_end  (id)
2896!compare checksums between global and x arrays
2897        call compare_checksums( global(1:nx,1:ny,:), gcheck, 'mpp_global_field' )
2898
2899!test mpp_global_sum
2900        gsum = sum( global(1:nx,1:ny,:) )
2901        id = mpp_clock_id( type//' sum', flags=MPP_CLOCK_SYNC+MPP_CLOCK_DETAILED )
2902        call mpp_clock_begin(id)
2903        lsum = mpp_global_sum( domain, x )
2904        call mpp_clock_end  (id)
2905        if( pe.EQ.mpp_root_pe() )print '(a,2es15.8,a,es12.4)', 'Fast sum=', lsum, gsum
2906!test exact mpp_global_sum
2907        id = mpp_clock_id( type//' exact sum', flags=MPP_CLOCK_SYNC+MPP_CLOCK_DETAILED )
2908        call mpp_clock_begin(id)
2909        lsum = mpp_global_sum( domain, x, BITWISE_EXACT_SUM )
2910        call mpp_clock_end  (id)
2911        if( pe.EQ.mpp_root_pe() )print '(a,2es15.8,a,es12.4)', 'Bitwise-exact sum=', lsum, gsum
2912    case( 'Folded' )!test vector update: folded, with sign flip at fold
2913!fill in folded north edge, cyclic east and west edge
2914        global(1-halo:0,    1:ny,:) = global(nx-halo+1:nx,1:ny,:)
2915        global(nx+1:nx+halo,1:ny,:) = global(1:halo,      1:ny,:)
2916        global(1-halo:nx+halo-1,ny+1:ny+halo,:) = -global(nx+halo-1:1-halo:-1,ny-1:ny-halo:-1,:)
2917        global(nx+halo,ny+1:ny+halo,:) = -global(nx-halo,ny-1:ny-halo:-1,:)
2918        global(1-halo:nx+halo,1-halo:0,:) = 0
2919
2920        x = 0.
2921        x(is:ie,js:je,:) = global(is:ie,js:je,:)
2922!set up y array
2923        allocate( y(isd:ied,jsd:jed,nz) )
2924        y = x
2925        id = mpp_clock_id( type//' vector', flags=MPP_CLOCK_SYNC+MPP_CLOCK_DETAILED )
2926        call mpp_clock_begin(id)
2927        call mpp_update_domains( x, y, domain, gridtype=BGRID_NE )
2928        call mpp_clock_end  (id)
2929!redundant points must be equal and opposite
2930        global(nx/2,ny,:) = 0.  !pole points must have 0 velocity
2931        global(nx  ,ny,:) = 0.  !pole points must have 0 velocity
2932        global(nx/2+1:nx-1, ny,:) = -global(nx/2-1:1:-1, ny,:)
2933        global(1-halo:0,    ny,:) = -global(nx-halo+1:nx,ny,:)
2934        global(nx+1:nx+halo,ny,:) = -global(1:halo,      ny,:)
2935        call compare_checksums( x, global(isd:ied,jsd:jed,:), type//' X' )
2936        call compare_checksums( y, global(isd:ied,jsd:jed,:), type//' Y' )
2937    end select
2938    return
2939  end subroutine test_halo_update
2940
2941  subroutine compare_checksums( a, b, string )
2942    real, intent(in), dimension(:,:,:) :: a, b
2943    character(len=*), intent(in) :: string
2944    integer :: i, j
2945
2946    call mpp_sync()
2947    i = mpp_chksum( a, (/pe/) )
2948    j = mpp_chksum( b, (/pe/) )
2949    if( i.EQ.j )then
2950        if( pe.EQ.mpp_root_pe() )call mpp_error( NOTE, trim(string)//': OK.' )
2951    else
2952        call mpp_error( FATAL, trim(string)//': chksums are not OK.' )
2953    end if
2954  end subroutine compare_checksums
2955end program mpp_domains_test
2956#endif
Note: See TracBrowser for help on using the repository browser.