New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

mapcomm_mod.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/mapcomm_mod.F90 @ 3849

Last change on this file since 3849 was 3849, checked in by trackstand2, 11 years ago
Merge branch 'partitioner'
File size: 77.4 KB

Line
1	MODULE mapcomm_mod
2	USE in_out_manager, ONLY: numout, lwp
3	USE par_oce, ONLY: jpiglo, jpjglo, jpreci, jprecj, jpi, jpk
4	USE dom_oce, ONLY: nlei, nlej, nldi, nldj, nlci, nlcj, narea, &
5	nleit, nlejt, nldit, nldjt, nlcit, nlcjt, &
6	nimpp, nimppt
7	IMPLICIT NONE
8
9	#define ARPDEBUG
10
11	PRIVATE
12
13	! No. of processes being partitioned over. Set in partition_mod() and
14	! used in both that module and this.
15	INTEGER, SAVE :: nprocp
16
17	! Process ids of ensemble member processes in a linear list.
18	INTEGER, ALLOCATABLE, DIMENSION(:) :: procid
19
20	! Information held by the ensemble leader about all processes.
21	! pnactive Number of active points in each sub-domain.
22	! pielb Lower (west) longitude bound index.
23	! pieub Upper (east) longitude bound index.
24	! piesub Number of longitude gridpoints.
25	! pilbext True if the lower longitude boundary is external.
26	! piubext True if the upper longitude boundary is external.
27	! pjelb Lower (south) latitude bound index.
28	! pjeub Upper (north) latitude bound index.
29	! pjesub Number of latitude gridpoints.
30	! pjlbext True if the lower latitude boundary is external.
31	! pjubext True if the upper latitude boundary is external.
32
33	INTEGER, ALLOCATABLE, DIMENSION(:) :: pnactive, &
34	pielb, pieub, piesub, pjelb, pjeub, pjesub
35	LOGICAL, ALLOCATABLE, DIMENSION(:) :: pilbext, piubext, pjlbext, pjubext
36
37	! Communications lists (one for sending, one for receiving)
38
39	! nsend Number of messages to be sent.
40	! dirsend Direction.
41	! destination Destination process id.
42	! isrcsend X coordinate of source data.
43	! jsrcsend Y coordinate of source data.
44	! idessend X coordinate of destination.
45	! jdessend Y coordinate of destination.
46	! nxsend Size in X of data to be sent.
47	! nysend Size in Y of data to be sent.
48
49	! nrecv Number of messages to be received.
50	! dirrecv Direction.
51	! source Source process id.
52	! idesrecv X coordinate of destination.
53	! jdesrecv Y coordinate of destination.
54	! nxrecv Size in X of data to be received.
55	! nyrecv Size in Y of data to be received.
56
57	INTEGER, PARAMETER :: MaxComm=16
58	INTEGER, SAVE, DIMENSION(MaxComm) :: dirsend,destination, &
59	dirrecv,source
60	! INTEGER, SAVE, DIMENSION(MaxComm,jpreci ) :: isrcsend,jsrcsend &
61	INTEGER, SAVE, DIMENSION(MaxComm) :: isrcsend,jsrcsend, &
62	isrcrecv,jsrcrecv, &
63	idessend,jdessend, &
64	nxsend,nysend,nzsend, &
65	idesrecv,jdesrecv, &
66	nxrecv,nyrecv,nzrecv
67	INTEGER, SAVE :: nsend,nrecv
68
69	! SMP 22 Sep 2009
70	! Alternative, run-length encoded communications lists
71	! omitting permanently dry points.
72	! Of these, isrcrecp, jsrcrecvp
73	! are set up but not currently used,
74	! and could be eliminated.
75
76	! Maximum number of patches a single halo communication can be broken
77	! into when trimming dry points in msg_trim()
78	INTEGER, PARAMETER :: MaxPatch=8
79	INTEGER, SAVE, DIMENSION(MaxPatch,MaxComm,jpreci) :: &
80	isrcsendp, jsrcsendp,&
81	nxsendp, nysendp, nzsendp, &
82	isrcrecvp, jsrcrecvp,&
83	idesrecvp, jdesrecvp,&
84	nxrecvp, nyrecvp, nzrecvp
85	INTEGER, SAVE, DIMENSION(MaxComm,jpreci) :: npatchsend, npatchrecv
86	! Total number of points in each message
87	INTEGER, SAVE, DIMENSION(MaxComm,jpreci) :: nsendp, nsendp2d, nrecvp, nrecvp2d
88
89	! Process dependent partitioning information.
90
91	! ielb Lower (west) longitude bound index.
92	! ieub Upper (east) longitude bound index.
93	! iesub Number of longitude gridpoints.
94	! ilbext True if the lower longitude boundary is external.
95	! iubext True if the upper longitude boundary is external.
96	! jelb Lower (south) latitude bound index.
97	! jeub Upper (north) latitude bound index.
98	! jesub Number of latitude gridpoints.
99	! jlbext True if the lower latitude boundary is external.
100	! jubext True if the upper latitude boundary is external.
101
102	INTEGER, SAVE :: ielb, ieub, iesub
103	INTEGER, SAVE :: jelb, jeub, jesub
104	LOGICAL, SAVE :: ilbext, iubext, jlbext, jubext
105
106	! Global definitions for parallel execution.
107
108	! Direction flags for communications.
109	! Listed so that opposite directions are given values maximally spaced
110	INTEGER, PARAMETER :: NONE=0 &
111	,Iplus=1 &
112	,Iminus=2 &
113	,Jplus=3 &
114	,Jminus=4 &
115	,IplusJplus=5 &
116	,IminusJminus=6 &
117	,IplusJminus=7 &
118	,IminusJplus=8 &
119	,MaxCommDir=8
120	! Array to hold direction flags for looking up the
121	! direction that is opposite to one we have
122	INTEGER, DIMENSION(MaxCommDir) :: opp_dirn
123
124	! Set up indices indicating the north-south and east-west
125	! attributes of the eight basic communication directions:
126	! four edges: W, E, S, N;
127	! four corners: SW, SE, NW, NE.
128	INTEGER, PARAMETER, DIMENSION(MaxCommDir) :: &
129	west = (/ 1, 0, 0, 0, 1, 0, 1, 0 /) &
130	,east = (/ 0, 1, 0, 0, 0, 1, 0, 1 /) &
131	,south = (/ 0, 0, 1, 0, 1, 0, 0, 1 /) &
132	,north = (/ 0, 0, 0, 1, 0, 1, 1, 0 /)
133	! 1 2 3 4 5 6 7 8
134	! W E S N SW NE NW SE
135
136	! cyclic_bc True if a cyclic boundary condition is to be applied
137	! Set using the value of jperio.
138	LOGICAL :: cyclic_bc
139
140	! Stores whether a domain's boundaries have been trimmed as
141	! trimmed(boundary, PE) where boundary is one of {n,e,s,w}idx
142	! for the Northern, Eastern, Southern or Western boundary, respectively.
143	! Allocated in finish_partition(), deallocated in...ARPDBG
144	LOGICAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: trimmed
145	INTEGER, PARAMETER :: nidx = 1, eidx= 2, sidx = 3, widx = 4
146
147	! Value representing land in the mask used for partitioning and message
148	! trimming
149	INTEGER, PARAMETER :: LAND = 0
150
151	! Rather than trim to a point immediately next to a wet point, we
152	! back off nextra points. If we don't do this then the sea-ice
153	! computation goes wrong because it does use values over the land
154	! that immediately border the ocean.
155	INTEGER, SAVE :: nextra
156
157	! Public routines
158	PUBLIC :: mapcomms, iprocmap, set_num_subdomains
159
160	! Public variables
161	PUBLIC :: MaxComm,nsend,nrecv,nxsend,nysend,destination,dirrecv, &
162	dirsend,isrcsend,jsrcsend,idesrecv, jdesrecv, &
163	nxrecv, nyrecv, source, cyclic_bc, idessend, jdessend
164
165	PUBLIC :: nsendp,nsendp2d,nrecvp,nrecvp2d,npatchsend,npatchrecv, &
166	nxsendp,nysendp,nzsendp,nxrecvp,nyrecvp,nzrecvp, &
167	idesrecvp,jdesrecvp,isrcsendp,jsrcsendp, nprocp
168
169	PUBLIC :: ielb, ieub, pielb, pjelb, pieub, pjeub, &
170	iesub, jesub, jeub, ilbext, iubext, jubext, jlbext, pnactive,&
171	piesub, pjesub, jelb, pilbext, pjlbext, pjubext, piubext
172
173	PUBLIC :: NONE &
174	,Iplus &
175	,Iminus &
176	,Jplus &
177	,Jminus &
178	,IplusJplus &
179	,IminusJminus &
180	,IplusJminus &
181	,IminusJplus &
182	,MaxCommDir
183
184	PUBLIC :: opp_dirn, land
185
186	PUBLIC :: trimmed, nidx, eidx, sidx, widx, nextra
187
188	! Switch for trimming dry points from halo swaps
189	LOGICAL, PARAMETER :: msgtrim = .TRUE.
190
191	! Switch for trimming points below ocean floor from halo swaps
192	! Defaults to true unless set via NEMO_MSGTRIM_Z environment var.
193	LOGICAL, PUBLIC, SAVE :: msgtrim_z
194
195	CONTAINS
196
197	SUBROUTINE set_num_subdomains(npes)
198	IMPLICIT none
199	!!------------------------------------------------------------------
200	!! Set the number of subdomains to partition the domain into.
201	!!------------------------------------------------------------------
202	INTEGER, INTENT(in) :: npes
203
204	nprocp = npes
205
206	END SUBROUTINE set_num_subdomains
207
208
209	SUBROUTINE mapcomms ( depth, ibotlvl, nx, ny, jperio, ierr )
210	!!------------------------------------------------------------------
211	! Maps out the communications requirements for the partitioned
212	! domain, adding communications descriptions to the list.
213	!
214	! Mike Ashworth, CLRC Daresbury Laboratory, July 1999
215	!!------------------------------------------------------------------
216	IMPLICIT NONE
217
218	! Subroutine arguments.
219	INTEGER, INTENT(in) :: nx, ny
220	INTEGER, INTENT(in) :: depth(nx,ny) ! Global mask: 0 for land, 1 for ocean
221	INTEGER, INTENT(in) :: ibotlvl(nx,ny)! Last vert level above sea floor
222	INTEGER, INTENT(in) :: jperio ! Periodicity of the mesh
223	INTEGER, INTENT(out):: ierr
224
225	! Local variables.
226	INTEGER :: i, i1, i2, ihalo, iproc, iprocc, iprocx, &
227	iprocy, j, j1, j2, nadd, naddmaxr, naddmaxs
228	INTEGER :: ldiff0, ldiff1 ! Local vars for coping with wrapping of coords
229	INTEGER :: imax, imin ! Max/min value of i that a halo strip can run
230	! to/from (used to avoid including E/W halo cells
231	! in N/S transfers)
232	INTEGER :: ielb_iproc, ieub_iproc ! Lower and upper bounds for proc
233	! iproc corrected for E & W halos
234	INTEGER :: ielb_no_halo, ieub_no_halo ! Lower and upper bounds for local
235	! domain corrected for E & W halos
236	INTEGER, DIMENSION(jpreci) :: idesr, jdesr, idess, jdess &
237	, isrcr, jsrcr, isrcs, jsrcs &
238	, nxr, nyr, nxs, nys
239	LOGICAL :: addcorner
240
241	! Clear the error code.
242	ierr = 0
243
244	! Store whether we have cyclic east-west boundary condition
245	! (See line 49 in domcfg.F90 for the definitions of the jperio values.)
246	IF(jperio == 1 .OR. jperio == 4 .OR. jperio == 6)THEN
247	cyclic_bc = .TRUE.
248	ELSE
249	cyclic_bc = .FALSE.
250	END IF
251
252	ALLOCATE(procid(nprocp), Stat=ierr)
253	IF (ierr > 0) THEN
254	WRITE (numout,*) 'ERROR: mapcomms: Allocate failed for iproc'
255	RETURN
256	END IF
257	! Create ordered list of process ids
258	DO i=1,nprocp,1
259	procid(i) = i-1
260	END DO
261
262	! Populate the look-up table of opposite directions
263	opp_dirn(Iplus ) = Iminus
264	opp_dirn(Iminus ) = Iplus
265	opp_dirn(Jplus ) = Jminus
266	opp_dirn(Jminus ) = Jplus
267	opp_dirn(IplusJplus ) = IminusJminus
268	opp_dirn(IminusJminus) = IplusJplus
269	opp_dirn(IplusJminus ) = IminusJplus
270	opp_dirn(IminusJplus ) = IplusJminus
271
272	dirsend = -999
273	destination = -999
274	isrcsend = -999
275	jsrcsend = -999
276	idessend = -999
277	jdessend = -999
278	nxsend = -999
279	nysend = -999
280	nzsend = -999
281	dirrecv = -999
282	source = -999
283	idesrecv = -999
284	jdesrecv = -999
285	nxrecv = -999
286	nyrecv = -999
287	nzrecv = -999
288
289	! For each of the eight communication directions on a 2d grid of
290	! processes, add a communication to the list.
291
292	! Note: the parameters Iplus, Iminus, etc. refer to array subscript
293	! references in the code like (I+1), (I-1), etc. which represent
294	! a requirement for communications in the OPPOSITE direction. So
295	! Iplus associates with sending data in the minus I direction,
296	! Iminus associates with sending data in the plus I direction, etc.
297
298	nsend = 0
299	nrecv = 0
300
301	! =================================================================
302	! Looking at the border where we will
303	! send data in the minus I direction (Iplus) and
304	! receive data that has been sent in the plus I direction (Iminus).
305	! =================================================================
306
307	! Start from the lower bound of the sub-domain, and carry on looking
308	! for communications with neighbours until we have reached
309	! the upper bound.
310
311	j1 = jelb
312	DO WHILE (j1.LE.jeub)
313
314	! Look for the process which owns the neighbouring point in the
315	! minus I direction.
316
317	iproc = iprocmap(ielb-1,j1)
318	IF ( iproc.GT.0 ) THEN
319
320	! Find where in the j direction the common border between these
321	! sub-domains ends.
322
323	j2 = MIN(jeub,pjeub(iproc))
324
325	#if defined ARPDEBUG
326	IF ( lwp ) THEN
327	WRITE (*,FMT="(I3,': strip along y at: ',I4,',',I4)") narea-1,j1,j2
328	ENDIF
329	#endif
330
331	! ( {i,I}nternal cells, {h,H}alo cells )
332	!
333	! PE=iproc PE=myself
334	!
335	! nleit(iproc) data nldi
336	! s
337	! \| <-------- \|
338	! r
339	! -------->
340	! ----------------------- ---------------------
341	! \| In-1 \| In \| H1 \| H2 \| \| h2 \| h1 \| i1 \| i2 \|
342	! ----------------------- ---------------------
343	! In -> h1
344	! In-1 -> h2
345	! H1 <- i1
346	! H2 <- i2
347
348	! Construct the rest of the data describing the zone,
349	! convert to local indices and extend to multiple halo widths.
350
351	! NEMO sets nldi=1 for nbondi==-1 but then calculates the source
352	! of data to send as follows...
353	! Since this is not == nldi, we don't use nldi here as one
354	! might have expected
355	isrcs(:) = jpreci + 1
356
357	DO ihalo=1,jpreci
358	! Source for the receive must be within internal domain of the
359	! (sending) PE, iproc
360	isrcr(ihalo) = nlcit(iproc) - jpreci - ihalo + 1 ! nleit(iproc)-ihalo+1
361	idesr(ihalo) = ihalo ! Halo goes from 1..jpreci
362	nxr(ihalo) = ihalo
363	nxs(ihalo) = ihalo
364	ENDDO
365
366	! MIN below allows for fact that NEMO sets nlci==nlei at the E
367	! boundary of the global domain when using cyclic bc's
368	idess(:) = MIN(nleit(iproc)+1,nlcit(iproc)) ! Send _to_ E halo column of iproc
369	! Source for a send is always within internal region
370	jsrcs(:) = j1-jelb+nldj !Add nldj 'cos jsrcs is local address in domain
371	jdess(:) = j1-pjelb(iproc)+nldjt(iproc) ! ditto
372	jdesr(:) = jsrcs(:)
373	jsrcr(:) = jdess(:)
374	nyr(:) = j2-j1+1
375	nys(:) = nyr(:)
376
377	! Examine whether corner points should be added to the start.
378
379	naddmaxr = 0
380	naddmaxs = 0
381	! ARPDBG - why loop below when naddmaxs and naddmaxr are scalars?
382	DO ihalo=1,jpreci,1
383
384	! Send corner data while we have data to send
385	! and while there is a point that depends on it.
386
387	IF ( j1-ihalo.GE.jelb .AND. &
388	iprocmap(ielb-ihalo,j1).GT.0 ) THEN
389	naddmaxs = ihalo
390	ENDIF
391
392	! Receive corner data only when we are at the corner
393	! and while the sending sub-domain is the same as for the edge.
394
395	IF ( j1.EQ.jelb .AND. &
396	iprocmap(ielb-ihalo,j1-ihalo).EQ.iproc ) THEN
397	naddmaxr = ihalo
398	ENDIF
399	ENDDO
400
401	! Add the extra points.
402
403	DO ihalo=1,jpreci
404	nadd = MIN(ihalo,naddmaxs)
405	jdess(ihalo) = jdess(ihalo) - nadd
406	jsrcs(ihalo) = jsrcs(ihalo) - nadd
407	nys(ihalo) = nys(ihalo)+nadd
408	#if defined ARPDEBUG
409	IF ( nadd.GT.0 ) THEN
410	WRITE (*,"(I3,': Adding starting corner to send for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
411	END IF
412	#endif
413	nadd = MIN(ihalo,naddmaxr)
414	jdesr(ihalo) = jdesr(ihalo) - nadd
415	jsrcr(ihalo) = jsrcr(ihalo) - nadd
416	nyr(ihalo) = nyr(ihalo)+nadd
417	#if defined ARPDEBUG
418	IF ( nadd.GT.0 ) THEN
419	WRITE (*,"(I3,': Adding starting corner to receive for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
420	ENDIF
421	#endif
422	ENDDO ! Loop over 'overlap' points in i direction
423
424	! Examine whether corner points should be added to the end.
425
426	naddmaxr = 0
427	naddmaxs = 0
428	DO ihalo=1,jpreci,1
429
430	! Send corner data while we have data to send
431	! and while there is a point that depends on it.
432
433	IF ( j2+ihalo.LE.jeub .AND. &
434	iprocmap(ielb-ihalo,j2).GT.0 ) THEN
435	naddmaxs = ihalo
436	ENDIF
437
438	! Receive corner data only when we are at the corner
439	! and while the sending sub-domain is the same as for the edge.
440
441	IF ( j2.EQ.jeub .AND. &
442	iprocmap(ielb-ihalo,j2+ihalo).EQ.iproc ) THEN
443	naddmaxr = ihalo
444	ENDIF
445	ENDDO
446
447	! Add the extra points.
448
449	DO ihalo=1,jpreci,1
450	nadd = MIN(ihalo,naddmaxs)
451	nys(ihalo) = nys(ihalo)+nadd
452	#if defined ARPDEBUG
453	IF ( nadd.GT.0 ) THEN
454	WRITE (*,"(I3,': Adding starting corner to send for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
455	ENDIF
456	#endif
457	nadd = MIN(ihalo,naddmaxr)
458	nyr(ihalo) = nyr(ihalo)+nadd
459	#if defined ARPDEBUG
460	IF ( nadd.GT.0 ) THEN
461	WRITE (*,"(I3,': Adding starting corner to receive for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
462	ENDIF
463	#endif
464	ENDDO
465
466	! Add a send and a receive to the lists for this section
467	! of border.
468
469	CALL addsend (nsend,Iplus,procid(iproc), &
470	isrcs,jsrcs,idess,jdess, &
471	nxs,nys,depth,ibotlvl,ierr)
472	IF ( ierr.NE.0 ) RETURN
473
474	! This is a receive for data _from_ the neighbouring domain
475	! - it is NOT the corresponding receive for the above send.
476	CALL addrecv (nrecv,Iminus,procid(iproc), &
477	isrcr,jsrcr,idesr,jdesr, &
478	nxr,nyr,depth,ibotlvl,ierr)
479	IF ( ierr.NE.0 ) RETURN
480	#if defined ARPDEBUG
481	IF ( lwp ) THEN
482	WRITE (*,'(a,7i6)') 'Adding receive ',procid(iproc) &
483	,isrcr(1),jsrcr(1),idesr(1),jdesr(1),nxr(1),nyr(1)
484	! WRITE (numout,'(21x,6i6)') &
485	! isrcr(2),jsrcr(2),idesr(2),jdesr(2),nxr(2),nyr(2)
486	ENDIF
487	#endif
488
489	! Move the start point to one beyond this strip.
490
491	j1 = j2+1
492
493	ELSE
494
495	! No process found, continue searching up the boundary.
496
497	j1 = j1+1
498	ENDIF
499	ENDDO
500
501	! =================================================================
502	! Looking at the border where we will
503	! send data in the plus I direction (Iminus) and
504	! receive data that has been sent in the minus I direction (Iplus).
505	! =================================================================
506
507	! Start from the lower bound of the sub-domain, and carry on looking
508	! for communications with neighbours until we have reached
509	! the upper bound.
510
511	j1 = jelb
512	DO WHILE (j1.LE.jeub)
513
514	! Look for the process which owns the neighbouring point in the
515	! plus I direction.
516
517	iproc = iprocmap(ieub+1,j1)
518	IF ( iproc.GT.0 ) THEN
519
520	! Find where in the j direction the common border between these
521	! sub-domains ends.
522
523	j2 = MIN(jeub,pjeub(iproc))
524
525	#if defined ARPDEBUG
526	WRITE (*,FMT="(I3,': ARPDBG strip for plus I is ',I3,',',I3)") &
527	narea-1,j1,j2
528	#endif
529
530	! ( {i,I}nternal cells, {h,H}alo cells )
531	!
532	! PE=myself PE=iproc
533	!
534	! nlei data nldit(iproc)
535	! s
536	! \| --------> \|
537	! r
538	! <--------
539	! ----------------------- ---------------------
540	! \| In-1 \| In \| H1 \| H2 \| \| h2 \| h1 \| i1 \| i2 \|
541	! ----------------------- ---------------------
542	! In -> h1
543	! In-1 -> h2
544	! H1 <- i1
545	! H2 <- i2
546
547	! Construct the rest of the data describing the zone,
548	! convert to local indexes and extend to multiple halo widths.
549
550	isrcr(:) = 1 + jpreci ! nldit(iproc) ARPDBG because NEMO sets nldi
551	! to unity if nbondi == -1 (W boundary of
552	! global domain)
553	DO ihalo=1,jpreci
554	! NEMO sets nlei = nlci if nbondi == 1 (E boundary of
555	! global domain). Normally, nlci = jpi = jpreci + iesub + jpreci
556	isrcs(ihalo) = jpi - jpreci - ihalo + 1 ! Outermost halo -> innermost col.
557	idess(ihalo) = ihalo ! Halo runs from 1..jpreci
558	nxr(ihalo) = ihalo
559	nxs(ihalo) = ihalo
560	ENDDO
561	idesr(:) = MIN(nlei + 1, nlci) ! Allow for case when on boundary
562	! of global domain and thus have
563	! no (explicit) halo
564	! Source for a send is within local internal domain
565	jsrcs(:) = j1-jelb+nldj
566	! Destination for a send is within halo on remote domain
567	jdess(:) = j1-pjelb(iproc)+nldjt(iproc)
568
569	jdesr(:) = jsrcs(:)
570	jsrcr(:) = jdess(:)
571	nyr(:) = j2-j1+1
572	nys(:) = nyr(:)
573	! Examine whether corner points should be added to the START.
574
575	naddmaxr = 0
576	naddmaxs = 0
577	DO ihalo=1,jpreci,1
578
579	! Send corner data while we have data to send
580	! and while there is a point that depends on it.
581
582	IF ( j1-ihalo.GE.jelb .AND. &
583	iprocmap(ieub+ihalo,j1).GT.0 ) THEN
584	naddmaxs = ihalo
585	ENDIF
586
587	! Receive corner data only when we are at the corner
588	! and while the sending sub-domain is the same as for the edge.
589
590	IF ( j1.EQ.jelb .AND. &
591	iprocmap(ieub+ihalo,j1-ihalo).EQ.iproc ) THEN
592	naddmaxr = ihalo
593	ENDIF
594	ENDDO
595
596	! Add the extra points.
597	DO ihalo=1,jpreci,1
598	nadd = MIN(ihalo,naddmaxs)
599	jdess(ihalo) = jdess(ihalo) - nadd
600	jsrcs(ihalo) = jsrcs(ihalo) - nadd
601	nys(ihalo) = nys(ihalo)+nadd
602	nadd = MIN(ihalo,naddmaxr)
603	jdesr(ihalo) = jdesr(ihalo) - nadd
604	jsrcr(ihalo) = jsrcr(ihalo) - nadd
605	nyr(ihalo) = nyr(ihalo)+nadd
606	ENDDO
607
608	! Examine whether corner points should be added to the end.
609
610	naddmaxr = 0
611	naddmaxs = 0
612	DO ihalo=1,jpreci,1
613
614	! Send corner data while we have data to send
615	! and while there is a point that depends on it.
616
617	IF ( j2+ihalo.LE.jeub .AND. &
618	iprocmap(ieub+ihalo,j2).GT.0 ) THEN
619	naddmaxs = ihalo
620	ENDIF
621
622	! Receive corner data only when we are at the corner
623	! and while the sending sub-domain is the same as for the edge.
624
625	IF ( j2.EQ.jeub .AND. &
626	iprocmap(ieub+ihalo,j2+ihalo).EQ.iproc ) THEN
627	naddmaxr = ihalo
628	ENDIF
629	ENDDO
630
631	! Add the extra points.
632
633	DO ihalo=1,jpreci,1
634	nadd = MIN(ihalo,naddmaxs)
635	nys(ihalo) = nys(ihalo)+nadd
636	#if defined ARPDEBUG
637	IF ( nadd.GT.0 .AND. lwp ) THEN
638	WRITE (,) 'Adding starting corner to send' &
639	,' for halo ',ihalo,' with ',nadd,' points'
640	ENDIF
641	#endif
642	nadd = MIN(ihalo,naddmaxr)
643	nyr(ihalo) = nyr(ihalo)+nadd
644	#if defined ARPDEBUG
645	IF ( nadd.GT.0 .AND. lwp ) THEN
646	WRITE (,) 'Adding starting corner to receive' &
647	,' for halo ',ihalo,' with ',nadd,' points'
648	ENDIF
649	#endif
650	ENDDO
651
652	! Add a send and a receive to the lists for this section
653	! of border.
654
655	CALL addsend (nsend,Iminus,procid(iproc), &
656	isrcs,jsrcs,idess,jdess,nxs,nys,&
657	depth,ibotlvl,ierr)
658	IF ( ierr.NE.0 ) RETURN
659	#if defined ARPDEBUG
660	IF ( lwp ) THEN
661	WRITE (*,'(a,7i6)') 'Adding send -1111 ',procid(iproc) &
662	,isrcs(1),jsrcs(1),idess(1),jdess(1),nxs(1),nys(1)
663	ENDIF
664	#endif
665
666	CALL addrecv (nrecv,Iplus,procid(iproc), &
667	isrcr,jsrcr,idesr,jdesr,nxr,nyr, &
668	depth,ibotlvl,ierr)
669	IF ( ierr.NE.0 ) RETURN
670
671	! Move the start point to one beyond this strip.
672
673	j1 = j2+1
674
675	ELSE
676
677	! No process found, continue searching up the boundary.
678
679	j1 = j1+1
680	ENDIF
681	ENDDO
682
683	! =================================================================
684	! Looking at the border where we will
685	! send data in the minus J direction (Jplus) and
686	! receive data that has been sent in the plus J direction (Jminus).
687	! =================================================================
688
689	! Ensure we don't include any values from the model domain W and E
690	! halos if we have cyclic b.c.'s
691	ielb_no_halo = ielb
692	ieub_no_halo = ieub
693	IF(cyclic_bc)THEN
694	! West
695	IF((.NOT. trimmed(widx,narea)) .AND. &
696	pilbext(narea)) ielb_no_halo = ielb_no_halo + jpreci
697	! East
698	IF((.NOT. trimmed(eidx,narea)) .AND. &
699	piubext(narea)) ieub_no_halo = ieub_no_halo - jpreci
700	END IF
701
702	! Start from the lower bound of the sub-domain (in global coords),
703	! and carry on looking
704	! for communications with neighbours until we have reached
705	! the upper bound.
706
707	imin = ielb_no_halo
708	imax = ieub_no_halo
709
710	i1 = imin
711	DO WHILE (i1.LE.imax)
712
713	! Look for the process which owns the neighbouring point in the
714	! minus J direction.
715
716	iproc = iprocmap(i1,jelb-1)
717	IF ( iproc.GT.0 ) THEN
718
719	! Ensure we don't include halos from the global domain borders if
720	! we have cyclic bc's.
721	! ielb_iproc = pielb(iproc)
722	ieub_iproc = pieub(iproc)
723	IF(cyclic_bc)THEN
724	! IF(pilbext(iproc))ielb_iproc = pielb(iproc)+jpreci
725	IF( (.NOT. trimmed(eidx,iproc)) .AND. &
726	piubext(iproc)) ieub_iproc = pieub(iproc)-jpreci
727	END IF
728
729	! Find where in the i direction the common border between these
730	! sub-domains ends.
731
732	i2 = MIN(imax,ieub_iproc)
733
734	#if defined ARPDEBUG
735	WRITE (*,FMT="(I3,': ARPDBG strip for minus J is ',I3,',',I3)") &
736	narea-1,i1,i2
737	#endif
738
739	! \| \| \| \| \|\|
740	! \| \| \| \| \|\|
741	! \| \| \| \| \|\|
742	! ------------------------------------------------
743	! \|\| \| \| \|
744	! \|\| \| \| \|
745	! \|\| \| \| \|
746
747
748	! Construct the rest of the data describing the zone,
749	! convert to local indexes and extend to multiple halo widths.
750	! Convert to local coords:
751	! Dist into zone = ipos - start + 1
752	! Pos in zone in local = (start of internal region) + (dist into zone) - 1
753	! Convert from global i1 to local i in current domain
754	! if i1 == nimpp then we must start at i=1 (because nimpp is absolute position
755	! of starting edge of domain including overlap regions)
756	isrcs(:) = i1 - nimpp + 1
757	! Convert from global i1 to local i in the destination domain
758	idess(:) = i1- nimppt(iproc) + 1
759	idesr(:) = isrcs(:)
760	isrcr(:) = idess(:)
761	nxr(:) = i2-i1+1
762	nxs(:) = nxr(:)
763
764	jsrcs(:) = nldj ! First row of 'internal' region of domain
765	DO ihalo=1,jprecj,1
766	! Source for a receive must be within internal region
767	jsrcr(ihalo) = nlejt(iproc)-ihalo+1
768	jdesr(ihalo) = ihalo ! Halo runs from 1..jprecj
769	nyr(ihalo) = ihalo
770	nys(ihalo) = ihalo
771	ENDDO
772	! Destination for a send must be a halo region
773	! nlcjt(iproc) is always in a halo region. Not sure what
774	! happens when halo wider than 1.
775	jdess(:) = nlcjt(iproc)
776
777	! Examine whether corner points should be added to the START.
778
779	naddmaxr = 0
780	naddmaxs = 0
781	DO ihalo=1,jprecj,1
782
783	! Send corner data while we have data to send
784	! and while there is a point that depends on it.
785
786	IF ( i1-ihalo.GE.imin .AND. &
787	iprocmap(i1-ihalo,jelb-ihalo).GT.0 ) THEN
788	naddmaxs = ihalo
789	ENDIF
790
791	! Receive corner data only when we are at the corner
792	! and while the sending sub-domain is the same as for the edge.
793
794	IF ( i1.EQ.imin .AND. &
795	iprocmap(i1-ihalo,jelb-ihalo).EQ.iproc ) THEN
796	naddmaxr = ihalo
797	ENDIF
798	ENDDO
799
800	! Add the extra points.
801
802	DO ihalo=1,jprecj,1
803	nadd = MIN(ihalo,naddmaxs)
804	idess(ihalo) = idess(ihalo)-nadd
805	isrcs(ihalo) = isrcs(ihalo)-nadd
806	nxs(ihalo) = nxs(ihalo)+nadd
807	#if defined ARPDEBUG
808	IF ( nadd.GT.0 ) THEN
809	WRITE (*,"(I3,': Adding starting corner to send for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
810	ENDIF
811	#endif
812	nadd = MIN(ihalo,naddmaxr)
813	idesr(ihalo) = idesr(ihalo)-nadd
814	isrcr(ihalo) = isrcr(ihalo)-nadd
815	nxr(ihalo) = nxr(ihalo)+nadd
816
817	#if defined ARPDEBUG
818	IF ( nadd.GT.0 ) THEN
819	WRITE (*,"(I3,': Adding starting corner to receive for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
820	ENDIF
821	#endif
822	ENDDO
823
824	! Examine whether corner points should be added to the END.
825
826	naddmaxr = 0
827	naddmaxs = 0
828	DO ihalo=1,jprecj,1
829
830	! Send corner data while we have data to send
831	! and while there is a point that depends on it.
832
833	IF ( i2+ihalo.LE.imax .AND. &
834	iprocmap(i2,jelb-ihalo).GT.0 ) THEN
835	naddmaxs = ihalo
836	ENDIF
837
838	! Receive corner data only when we are at the corner
839	! and while the sending sub-domain is the same as for the edge.
840
841	IF ( i2.EQ.imax .AND. &
842	iprocmap(i2+ihalo,jelb-ihalo).EQ.iproc ) THEN
843	naddmaxr = ihalo
844	ENDIF
845	ENDDO
846
847	! Add the extra points.
848
849	DO ihalo=1,jprecj,1
850	nadd = MIN(ihalo,naddmaxs)
851	nxs(ihalo) = nxs(ihalo)+nadd
852	#if defined ARPDEBUG
853	IF ( nadd.GT.0 ) THEN
854	WRITE (*,"(I3,': Adding starting corner to send for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
855	ENDIF
856	#endif
857	nadd = MIN(ihalo,naddmaxr)
858	nxr(ihalo) = nxr(ihalo)+nadd
859	#if defined ARPDEBUG
860	IF ( nadd.GT.0 ) THEN
861	WRITE (*,"(I3,': Adding starting corner to receive for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
862	ENDIF
863	#endif
864	ENDDO
865
866	! Add a send and a receive to the lists for this section
867	! of border.
868
869	CALL addsend (nsend,Jplus,procid(iproc), &
870	isrcs,jsrcs,idess,jdess,nxs,nys, &
871	depth,ibotlvl,ierr)
872	IF ( ierr.NE.0 ) RETURN
873
874	CALL addrecv (nrecv,Jminus,procid(iproc), &
875	isrcr,jsrcr,idesr,jdesr,nxr,nyr, &
876	depth,ibotlvl,ierr)
877	IF ( ierr.NE.0 ) RETURN
878
879	! Move the start point to one beyond this strip.
880
881	i1 = i2+1
882
883	ELSE
884
885	! No process found, continue searching up the boundary.
886
887	i1 = i1+1
888	ENDIF
889	ENDDO
890
891	! =================================================================
892	! Looking at the border where we will
893	! send data in the plus J direction (Jminus) and
894	! receive data that has been sent in the minus J direction (Jplus).
895	! =================================================================
896
897	! Start from the lower bound of the sub-domain, and carry on looking
898	! for communications with neighbours until we have reached
899	! the upper bound.
900
901	imin = ielb_no_halo
902	imax = ieub_no_halo
903	i1 = imin
904
905	DO WHILE (i1.LE.imax)
906
907	! Look for the process which owns the neighbouring point in the
908	! plus J direction.
909
910	iproc = iprocmap(i1,jeub+1)
911	IF ( iproc.GT.0 ) THEN
912	! Ensure we don't include halos from the global borders if we
913	! have cyclic b.c.'s
914	! ielb_iproc = pielb(iproc)
915	ieub_iproc = pieub(iproc)
916	IF(cyclic_bc)THEN
917	! IF(pilbext(iproc))ielb_iproc = pielb(iproc)+jpreci
918	IF((.NOT. trimmed(eidx,iproc)) .AND. &
919	piubext(iproc))ieub_iproc = pieub(iproc)-jpreci
920	END IF
921
922	! Find where in the i direction the common border between these
923	! sub-domains ends.
924
925	i2 = MIN(imax, ieub_iproc)
926
927	#if defined ARPDEBUG
928	WRITE (*,FMT="(I3,': ARPDBG strip for plus J is ',I3,',',I3)") &
929	narea-1,i1,i2
930	#endif
931
932	! Construct the rest of the data describing the zone,
933	! convert to local indexes and extend to multiple halo widths.
934
935	isrcs(:) = i1 - nimpp + 1
936	idess(:) = i1 - nimppt(iproc) + 1
937	idesr(:) = isrcs(:)
938	isrcr(:) = idess(:)
939	nxr(:) = i2-i1+1
940	nxs(:) = nxr(:)
941
942	! Source for a receive must be within an internal region
943	! nldj incorporates whether or not lower halo exists
944	jsrcr(:) = nldjt(iproc)
945
946	DO ihalo=1,jprecj,1
947	jsrcs(ihalo) = nlej-ihalo+1 ! innermost row -> outermost halo
948	jdess(ihalo) = ihalo ! Halo runs from 1..jprecj
949	nyr(ihalo) = ihalo
950	nys(ihalo) = ihalo
951	ENDDO
952	jdesr(:) = nlcj
953
954	! Examine whether corner points should be added to the START.
955
956	naddmaxr = 0
957	naddmaxs = 0
958	DO ihalo=1,jprecj,1
959
960	! Send corner data while we have data to send
961	! and while there is a point that depends on it.
962
963	IF ( i1-ihalo.GE.imin .AND. iprocmap(i1,jeub+ihalo).GT.0 ) THEN
964	naddmaxs = ihalo
965	ENDIF
966
967	! Receive corner data only when we are at the corner
968	! and while the sending sub-domain is the same as for the edge.
969
970	IF ( i1.EQ.imin .AND. &
971	iprocmap(i1-ihalo,jeub+ihalo).EQ.iproc ) THEN
972	naddmaxr = ihalo
973	ENDIF
974	ENDDO
975
976	! Add the extra points.
977
978	DO ihalo=1,jprecj,1
979	nadd = MIN(ihalo,naddmaxs)
980	idess(ihalo) = idess(ihalo) -nadd
981	isrcs(ihalo) = isrcs(ihalo) -nadd
982	nxs(ihalo) = nxs(ihalo)+nadd
983	#if defined ARPDEBUG
984	IF ( nadd.GT.0 ) THEN
985	WRITE (*,"(I3,': Adding starting corner to send for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
986	ENDIF
987	#endif
988	nadd = MIN(ihalo,naddmaxr)
989	idesr(ihalo) = idesr(ihalo) - nadd
990	isrcr(ihalo) = isrcr(ihalo) - nadd
991	nxr(ihalo) = nxr(ihalo)+nadd
992
993	#if defined ARPDEBUG
994	IF ( nadd.GT.0 ) THEN
995	WRITE (*,"(I3,': Adding starting corner to receive for halo ',I2,' with ',I3,' points')") narea-1,ihalo, nadd
996	ENDIF
997	#endif
998	ENDDO
999
1000	! Examine whether corner points should be added to the END.
1001
1002	naddmaxr = 0
1003	naddmaxs = 0
1004	DO ihalo=1,jprecj,1
1005
1006	! Send corner data while we have data to send
1007	! and while there is a point that depends on it.
1008
1009	IF ( i2+ihalo.LE.imax .AND. &
1010	iprocmap(i2,jeub+ihalo).GT.0 ) THEN
1011	naddmaxs = ihalo
1012	ENDIF
1013
1014	! Receive corner data only when we are at the corner
1015	! and while the sending sub-domain is the same as for the edge.
1016
1017	IF ( i2.EQ.imax .AND. & ! Are we at the corner?
1018	iprocmap(i2+ihalo,jeub+ihalo).EQ.iproc ) THEN
1019	naddmaxr = ihalo
1020	ENDIF
1021	ENDDO
1022
1023	! Add the extra points.
1024
1025	DO ihalo=1,jprecj,1
1026	nadd = MIN(ihalo,naddmaxs)
1027	nxs(ihalo) = nxs(ihalo)+nadd
1028	#if defined ARPDEBUG
1029	IF ( nadd.GT.0 .AND. lwp ) THEN
1030	WRITE (,) narea-1,': Adding starting corner to send' &
1031	,' for halo ',ihalo,' with ',nadd,' points'
1032	ENDIF
1033	#endif
1034	nadd = MIN(ihalo,naddmaxr)
1035	nxr(ihalo) = nxr(ihalo)+nadd
1036	#if defined ARPDEBUG
1037	IF ( nadd.GT.0 .AND. lwp ) THEN
1038	WRITE (,) narea-1,': Adding starting corner to receive' &
1039	,' for halo ',ihalo,' with ',nadd,' points'
1040	ENDIF
1041	#endif
1042	ENDDO
1043
1044	! Add a send and a receive to the lists for this section
1045	! of border.
1046
1047	CALL addsend (nsend,Jminus,procid(iproc), &
1048	isrcs,jsrcs,idess,jdess,nxs,nys, &
1049	depth,ibotlvl,ierr)
1050	IF ( ierr.NE.0 ) RETURN
1051
1052	CALL addrecv (nrecv,Jplus,procid(iproc), &
1053	isrcr,jsrcr,idesr,jdesr,nxr,nyr, &
1054	depth,ibotlvl,ierr)
1055	IF ( ierr.NE.0 ) RETURN
1056
1057	! Move the start point to one beyond this strip.
1058
1059	i1 = i2+1
1060
1061	ELSE ! iproc < 0
1062
1063	! No process found, continue searching up the boundary.
1064
1065	i1 = i1+1
1066	ENDIF
1067	ENDDO
1068
1069	! =================================================================
1070	! Corner points are sent with the edge data where possible.
1071	! Check for cases where corner data resides on a processor which
1072	! is not participating in edge communications and set up extra
1073	! diagonal communications for these cases.
1074	! =================================================================
1075
1076	! Loop over the four corner directions
1077	! i = 1 2 3 4 5 6 7 8
1078	! W E S N SW NE NW SE
1079
1080
1081	DO i=5,8
1082
1083	! At first assume there is to be no corner communication
1084
1085	addcorner = .FALSE.
1086
1087	! i1 is to be x-coord just OUTSIDE our domain
1088	! i2 is to be x-coord just WITHIN our domain
1089	! All the following complexity is to allow for fact that the first
1090	! and last columns of the global domain are actually halos when
1091	! we have cyclic E/W boundary conditions.
1092	IF( (iubext .OR. ilbext) .AND. cyclic_bc) THEN
1093	i1 = ielb
1094	i2 = ielb
1095	IF(ilbext .AND. (.NOT. trimmed(widx,narea)) )THEN
1096	i2 = i2+west(i) ! If on W boundary with cyclic bc's, ielb _is_ the halo column
1097	! so add 1 to move inside domain
1098	ELSE
1099	i1 = i1-west(i)
1100	END IF
1101	IF(iubext .AND. (.NOT. trimmed(eidx,narea)) )THEN
1102	! If upper bound is on domain boundary then iesub already
1103	! includes the halo column
1104	i1 = i1+east(i)*(iesub-1)
1105	i2 = i2+east(i)*(iesub-2)
1106	ELSE
1107	i1 = i1+east(i)*iesub
1108	i2 = i2+east(i)*(iesub-1)
1109	END IF
1110
1111	ELSE
1112	i1 = ielb-west(i)+east(i)*iesub
1113	i2 = ielb+east(i)*(iesub-1)
1114	END IF
1115
1116	! For a NW corner:
1117	! \|
1118	! iproc \| iprocy
1119	! ________\|______
1120	! \|
1121	! iprocx \| Me
1122	! \|
1123
1124	! x coord just OUTSIDE our domain but y INSIDE
1125	iprocx = iprocmap(i1, jelb+north(i)*(jesub-1))
1126	! x coord just INSIDE our domain but y OUTSIDE
1127	iprocy = iprocmap(i2, jelb-south(i)+north(i)*jesub)
1128
1129	iprocc = 0
1130
1131	#if defined ARPDEBUG
1132	WRITE(*,FMT="(I3,': ARPDBG: i, i1 (outside), i2 (inside), iprocx, iprocy = ',5I4)") &
1133	narea-1, i,i1,i2,iprocx,iprocy
1134	#endif
1135
1136	! Loop over all required halo widths
1137
1138	DO ihalo=jpreci,1,-1
1139
1140	! Look at the processor in the corner at width ihalo from the
1141	! corner of the sub-domain. We want both x and y to be just
1142	! outside our domain. i1 is already just outside our domain
1143	! so we subtract one from ihalo below:
1144	iproc = iprocmap( i1 - west(i)(ihalo-1)+ east(i)(ihalo-1) &
1145	,south(i)(jelb-ihalo)+north(i)(jeub+ihalo))
1146	! iproc = iprocmap( west(i)(ielb-ihalo)+ east(i)(ieub+ihalo) &
1147	! ,south(i)(jelb-ihalo)+north(i)(jeub+ihalo))
1148
1149	! If the corner processor is different from those to X and Y
1150	! we will need a corner communication.
1151
1152	IF ( iproc.GT.0 .AND. iprocx.GT.0 .AND. iprocy.GT.0 .AND. &
1153	iproc.NE.iprocx .AND. iproc.NE.iprocy ) THEN
1154
1155	! Ensure we don't include halos from the global borders if we
1156	! have cyclic E/W boundaries.
1157	ielb_iproc = pielb(iproc)
1158	ieub_iproc = pieub(iproc)
1159	IF( cyclic_bc )THEN
1160	IF(pilbext(iproc))ielb_iproc = pielb(iproc)+ihalo
1161	IF(piubext(iproc))ieub_iproc = pieub(iproc)-ihalo
1162	END IF
1163
1164	#if defined ARPDEBUG
1165	WRITE (*,FMT="(I3,': adding corner as ',I3,' differs from ',2I3)")&
1166	narea-1, iproc,iprocx,iprocy
1167	#endif
1168	! If the furthest corner point needs a communication,
1169	! we will need them all.
1170
1171	IF ( ihalo.EQ.jpreci ) THEN
1172	iprocc = iproc
1173
1174	! Ensure we don't include halos from the global borders if we
1175	! have cyclic E/W boundaries.
1176	ielb_iproc = pielb(iprocc)
1177	ieub_iproc = pieub(iprocc)
1178	IF( cyclic_bc )THEN
1179	IF(pilbext(iprocc))ielb_iproc = pielb(iprocc)+jpreci
1180	IF(piubext(iprocc))ieub_iproc = pieub(iprocc)-jpreci
1181	END IF
1182	! Set the flag to add everything to the communications list.
1183
1184	addcorner = .TRUE.
1185
1186	ENDIF
1187
1188	! Set the parameters for the communication.
1189	ldiff0 = ielb_iproc - ieub_no_halo
1190	ldiff1 = ielb_no_halo - ieub_iproc
1191	! Allow for wrap-around if necessary
1192	IF(cyclic_bc)THEN
1193	IF(ldiff0 < 1)THEN
1194	!ARPDBG -2 for consistency with procmap
1195	ldiff0 = ldiff0 + (jpiglo - 2)
1196	END IF
1197	IF(ldiff1 < 1)THEN
1198	!ARPDBG -2 for consistency with procmap
1199	ldiff1 = ldiff1 + (jpiglo - 2)
1200	END IF
1201	END IF
1202	nxs (ihalo) = ihalo - east(i)*(ldiff0-1) &
1203	- west(i)*(ldiff1-1)
1204	! Have no cyclic b.c.'s in N/S direction so probably don't need
1205	! the following checks on ldiff{0,1}
1206	ldiff0 = pjelb(iprocc) - jeub
1207	IF(ldiff0 < 1) ldiff0 = ldiff0 + jpjglo
1208	ldiff1 = jelb - pjeub(iprocc)
1209	IF(ldiff1 < 1) ldiff1 = ldiff1 + jpjglo
1210	nys (ihalo) = ihalo - north(i)*(ldiff0-1) &
1211	- south(i)*(ldiff1-1)
1212
1213	! Source for a send must be within the internal region of
1214	! the LOCAL domain
1215	isrcs(ihalo) = east(i) *(iesub-nxs(ihalo)) + nldi
1216	jsrcs(ihalo) = north(i)*(jesub-nys(ihalo)) + nldj
1217	IF( cyclic_bc )THEN
1218	IF( ilbext )THEN
1219	! nldi is still within halo for domains on W edge of
1220	! global domain
1221	isrcs(ihalo) = isrcs(ihalo) + west(i)
1222	ELSE IF( iubext )THEN
1223	! Final column is actually halo for domains on E edge of
1224	! global domain
1225	isrcs(ihalo) = isrcs(ihalo) - east(i)
1226	END IF
1227	END IF
1228
1229	! Destination for a send must be within a halo region on the
1230	! REMOTE domain
1231	! MAX and MIN below to allow for edge domains that do not have
1232	! explicit halo
1233	idess(ihalo) = west(i)*(MIN(nleit(iprocc)+ihalo, nlcit(iprocc)) &
1234	- nxs(ihalo)+1) &
1235	+ east(i)*MAX(nldit(iprocc)-ihalo, 1)
1236
1237	! MAX and MIN below to allow for edge domains that do not have
1238	! explicit halo
1239	jdess(ihalo) = south(i)*(MIN(nlejt(iprocc)+ihalo,nlcjt(iprocc))-nys(ihalo)+1) &
1240	+ north(i)*MAX(nldjt(iprocc) - ihalo,1)
1241
1242	! Source for a receive must be in an internal region of the REMOTE domain
1243	isrcr(ihalo) = west(i)*(piesub(iprocc)-nxs(ihalo)) + nldit(iprocc)
1244	IF( cyclic_bc )THEN
1245
1246	! This _could_ be a corner exchange wrapped around by the cyclic
1247	! boundary conditions:
1248	!
1249	! \|\|------\| \|\|
1250	! \|\| \| \| \|\|
1251	! \|\|a_____\|__ _ _ \| \|\|
1252	! \|\| -----------\|\|
1253	! \|\| \| a\|\|
1254	! \|\| \|________\|\|
1255
1256	IF(pilbext(iprocc))THEN
1257	! nldi is still within halo for domains on E edge of
1258	! global domain
1259	isrcr(ihalo) = isrcr(ihalo) + east(i)
1260	ELSE IF(piubext(iprocc))THEN
1261	! Final column is actually halo for domains on W edge of
1262	! global domain
1263	isrcr(ihalo) = isrcr(ihalo) - west(i)
1264	END IF
1265	END IF
1266	jsrcr(ihalo) = south(i)*(pjesub(iprocc)-nys(ihalo)) + nldjt(iprocc)
1267
1268	! Destination for a receive must be in a halo region (on LOCAL
1269	! domain) and therefore:
1270	! 1 <= {i,j}desr <= jprec{i,j} or >= nle{i,j} + 1
1271	idesr(ihalo) = east(i)*(MIN(nlci,nlei+ihalo)-nxs(ihalo)+1) & !ARPDBG s/iesub/nlei/
1272	+ west(i)*MAX(1,nldi-ihalo) !ARPDBG incl. nldi-
1273
1274	jdesr(ihalo) = north(i)*(MIN(nlcj,nlej+ihalo)-nys(ihalo) + 1) &
1275	+ south(i)*MAX(1,nldj - ihalo)
1276
1277	ELSE
1278
1279	#if defined ARPDEBUG
1280	IF ( iprocc.GT.0 ) THEN
1281	WRITE (*,FMT="(I3,': skipping corner for halo ',I3,' PE ',I3)")&
1282	narea-1,ihalo,iprocc-1
1283	ENDIF
1284	#endif
1285
1286	! No communication for this halo width.
1287	! Clear all the parameters
1288	! in case there are comms for other halo widths.
1289
1290	isrcs(ihalo) = 0
1291	jsrcs(ihalo) = 0
1292	idess(ihalo) = 0
1293	jdess(ihalo) = 0
1294	isrcr(ihalo) = 0
1295	jsrcr(ihalo) = 0
1296	idesr(ihalo) = 0
1297	jdesr(ihalo) = 0
1298	nxs (ihalo) = 0
1299	nys (ihalo) = 0
1300
1301	ENDIF
1302
1303	ENDDO
1304
1305	! The size of the received data is always the same as
1306	! that of the sent data.
1307
1308	nxr(:) = nxs(:)
1309	nyr(:) = nys(:)
1310
1311	! Add the data to the communications lists.
1312
1313	IF ( addcorner ) THEN
1314	#if defined ARPDEBUG
1315	WRITE (*,FMT="(I3,': ARPDBG adding corner send to ',I4,', dir = ',I1)") &
1316	narea-1, procid(iprocc),i
1317	#endif
1318	CALL addsend (nsend,i,procid(iprocc), &
1319	isrcs,jsrcs,idess,jdess,nxs,nys, &
1320	depth,ibotlvl,ierr)
1321	IF ( ierr.NE.0 ) RETURN
1322
1323	! Manually reverse the direction indicator for the receive.
1324	j = opp_dirn(i)
1325
1326	#if defined ARPDEBUG
1327	WRITE (*,FMT="(I3,': ARPDBG adding corner recv. from ',I4,', old dir = ',I1,' new dir = ',I1)") &
1328	narea-1, procid(iprocc),i, j
1329	#endif
1330	CALL addrecv (nrecv,j,procid(iprocc), &
1331	isrcr,jsrcr,idesr,jdesr,nxr,nyr, &
1332	depth,ibotlvl,ierr)
1333	IF ( ierr.NE.0 ) RETURN
1334
1335	ENDIF
1336
1337	ENDDO
1338
1339	DEALLOCATE(procid)
1340
1341	END SUBROUTINE mapcomms
1342
1343
1344	FUNCTION iprocmap ( ia, ja )
1345	!!------------------------------------------------------------------
1346	! Returns the process number (1...N) of the process whose sub-domain
1347	! contains the point with global coordinates (i,j).
1348	! If no process owns the point, returns zero.
1349
1350	! i int input global x-coordinate
1351	! j int input global y-coordinate
1352
1353	! Mike Ashworth, CLRC Daresbury Laboratory, July 1999
1354	! Andrew Porter, STFC Daresbury Laboratory, May 2008
1355	!!------------------------------------------------------------------
1356	IMPLICIT NONE
1357
1358	! Function arguments.
1359	INTEGER :: iprocmap
1360	INTEGER, INTENT(in) :: ia, ja
1361	! Local variables.
1362	INTEGER :: iproc, i, j, iwidth
1363
1364	iprocmap = 0
1365
1366	! Make sure we don't change variable values in calling
1367	! routine...
1368	i = ia
1369	j = ja
1370	IF(cyclic_bc)THEN
1371	! Allow for fact that first and last columns in global domain
1372	! are actually halos
1373	iwidth = jpiglo - 2*jpreci
1374	IF(i >= jpiglo) i = ia - iwidth
1375	IF(i <= 1 ) i = ia + iwidth
1376	! No cyclic BCs in North/South direction
1377	!IF(j > jpjglo) j = ja - jpjglo
1378	!IF(j < 1 ) j = ja + jpjglo
1379	END IF
1380
1381	! Search the processes for the one owning point (i,j).
1382
1383	DO iproc=1,nprocp
1384	IF ( pielb(iproc).LE.i .AND. i.LE.pieub(iproc) .AND. &
1385	pjelb(iproc).LE.j .AND. j.LE.pjeub(iproc) ) THEN
1386	iprocmap = iproc
1387	EXIT
1388	ENDIF
1389	ENDDO
1390
1391	! ARP - for debugging only
1392	!!$ IF(iprocmap == 0)THEN
1393	!!$ WRITE(*,"('iprocmap: failed to find owner PE for (',I3,1x,I3,')')") ia, ja
1394	!!$ WRITE(,) 'PE domains are [xmin:xmax][ymin:ymax]:'
1395	!!$ DO iproc=1,nprocp,1
1396	!!$ WRITE(*,"(I3,': [',I3,':',I3,'][',I3,':',I3,']')") &
1397	!!$ iproc, pielb(iproc), pieub(iproc), pjelb(iproc), pjeub(iproc)
1398	!!$ END DO
1399	!!$ END IF
1400
1401	END FUNCTION iprocmap
1402
1403	SUBROUTINE addsend ( icomm, dir, proc, isrc, jsrc, &
1404	ides, jdes, nx, ny, depth, ibotlvl, ierr )
1405	!!------------------------------------------------------------------
1406	! Adds a send communication specified by the parameters dir through
1407	! to ny to the send communication list at the next position.
1408	! icomm points to the last entry and is incremented and returned
1409	! if successful.
1410	!
1411	! icomm int in/out Location in comms list.
1412	! dir int input Direction.
1413	! proc int input Process id.
1414	! isrc int input X coordinate of source data.
1415	! jsrc int input Y coordinate of source data.
1416	! ides int input X coordinate of destination data.
1417	! jdes int input Y coordinate of destination data.
1418	! nx int input Size in X of data to be sent.
1419	! ny int input Size in Y of data to be sent.
1420	! depth input Global mask, 0 for land, 1 for wet
1421	! ierr int output Error flag.
1422	!
1423	! Mike Ashworth, CLRC Daresbury Laboratory, March 1999
1424	! Stephen Pickles, STFC Daresbury Laboratory
1425	! - Aug 2009: added depth argument and message clipping
1426	! - Sep 2009: added run-length encoding
1427	!!------------------------------------------------------------------
1428	IMPLICIT NONE
1429
1430	! Subroutine arguments.
1431	INTEGER, INTENT(inout) :: icomm
1432	INTEGER, INTENT( in ) :: dir, proc
1433	! Global mask: 0 for land, 1 for ocean
1434	INTEGER, DIMENSION(:,:), INTENT( in ) :: depth
1435	INTEGER, DIMENSION(:,:), INTENT( in ) :: ibotlvl
1436	INTEGER, INTENT( out ) :: ierr
1437	INTEGER, DIMENSION(jpreci), INTENT( in ) :: isrc, jsrc, &
1438	ides, jdes, nx, ny
1439	! Values of corresponding input arguments after clipping
1440	INTEGER, DIMENSION(jpreci) :: cisrc,cjsrc,cides,cjdes,cnx,cny,cnz
1441	! Run-length encoded versions corresponding to above
1442	INTEGER, DIMENSION(MaxPatch,jpreci) :: risrc,rjsrc,rides,rjdes,rnx,rny,rnz
1443	! Number of patches in run-length encoded message
1444	INTEGER, DIMENSION(jpreci) :: npatches
1445	INTEGER :: ihalo, ipatch
1446	INTEGER :: nsendp_untrimmedz ! How many pts we'd be sending without
1447	! trimming in z direction
1448	! Whether there is still a message after clipping
1449	LOGICAL :: something_left
1450
1451	! Clear the error flag.
1452	ierr = 0
1453
1454	! Return if the process id is not set.
1455
1456	IF ( proc.LT.0 ) THEN
1457	RETURN
1458	ENDIF
1459
1460	! Can the message be clipped ?
1461
1462	CALL clip_msg(depth, ibotlvl, &
1463	isrc, jsrc, ides, jdes, nx, ny, &
1464	cisrc,cjsrc,cides,cjdes,cnx,cny,cnz, &
1465	risrc,rjsrc,rides,rjdes,rnx,rny,rnz, &
1466	npatches, something_left)
1467
1468	IF (something_left) THEN
1469
1470	icomm = icomm+1
1471
1472	! Check that the comms list has space for another entry.
1473
1474	IF ( icomm.GT.MaxComm ) THEN
1475	IF ( lwp ) THEN
1476	WRITE (numout,*) 'ERROR: Number of separate ', &
1477	'communications exceeds maximum of ',MaxComm
1478	ENDIF
1479	ierr = -12
1480	RETURN
1481	ENDIF
1482
1483	! Add the data into the comms list at the new location.
1484
1485	dirsend(icomm) = dir
1486	destination(icomm) = proc
1487	isrcsend(icomm) = cisrc(1)
1488	jsrcsend(icomm) = cjsrc(1)
1489	idessend(icomm) = cides(1)
1490	jdessend(icomm) = cjdes(1)
1491	nxsend(icomm) = cnx(1)
1492	nysend(icomm) = cny(1)
1493	IF(msgtrim_z)THEN
1494	nzsend(icomm) = cnz(1)
1495	ELSE
1496	nzsend(icomm) = jpk
1497	END IF
1498
1499	! Zero count of untrimmed pts to send
1500	nsendp_untrimmedz = 0
1501
1502	! Also set up the comms lists encoded as the start points and
1503	! lengths of the contiguous runs of wet points.
1504	DO ihalo=1,jpreci
1505
1506	nsendp2d(icomm,ihalo) = 0
1507	nsendp(icomm,ihalo) = 0
1508	npatchsend(icomm,ihalo) = npatches(ihalo)
1509
1510	DO ipatch=1,npatches(ihalo)
1511
1512	isrcsendp(ipatch,icomm,ihalo) = risrc(ipatch,ihalo)
1513	jsrcsendp(ipatch,icomm,ihalo) = rjsrc(ipatch,ihalo)
1514
1515	nxsendp(ipatch,icomm,ihalo) = rnx(ipatch,ihalo)
1516	nysendp(ipatch,icomm,ihalo) = rny(ipatch,ihalo)
1517	IF(msgtrim_z)THEN
1518	nzsendp(ipatch,icomm,ihalo)= rnz(ipatch,ihalo)
1519	ELSE
1520	nzsendp(ipatch,icomm,ihalo) = jpk
1521	END IF
1522
1523	! Sum the no. of points to be sent over all
1524	! patches for both 2D-array halos and 3D-array halos
1525	nsendp2d(icomm,ihalo) = nsendp2d(icomm,ihalo) + &
1526	nxsendp(ipatch,icomm,ihalo)* &
1527	nysendp(ipatch,icomm,ihalo)
1528	nsendp(icomm,ihalo) = nsendp(icomm,ihalo) + &
1529	nxsendp(ipatch,icomm,ihalo)* &
1530	nysendp(ipatch,icomm,ihalo)* &
1531	nzsendp(ipatch,icomm,ihalo)
1532	IF(msgtrim_z)THEN
1533	nsendp_untrimmedz = nsendp_untrimmedz + &
1534	nxsendp(ipatch,icomm,ihalo)* &
1535	nysendp(ipatch,icomm,ihalo)* &
1536	jpk
1537	END IF
1538	END DO
1539	END DO
1540
1541	#if defined ARPDEBUG
1542	WRITE (*,FMT="(I4,': ARPDBG adding SEND:')") narea-1
1543	WRITE (*,FMT="(I4,': ARPDBG: icomm = ',I2)") narea-1,icomm
1544	WRITE (*,FMT="(I4,': ARPDBG: dir = ',I2)") narea-1,dirsend(icomm)
1545	WRITE (*,FMT="(I4,': ARPDBG: proc = ',I4)") narea-1,destination(icomm)
1546	WRITE (*,FMT="(I4,': ARPDBG: isrc = ',I4)") narea-1,isrcsend(icomm)
1547	WRITE (*,FMT="(I4,': ARPDBG: jsrc = ',I4)") narea-1,jsrcsend(icomm)
1548	WRITE (*,FMT="(I4,': ARPDBG: ides = ',I4)") narea-1,idessend(icomm)
1549	WRITE (*,FMT="(I4,': ARPDBG: jdes = ',I4)") narea-1,jdessend(icomm)
1550	WRITE (*,FMT="(I4,': ARPDBG: nx = ',I4)") narea-1,nxsend(icomm)
1551	WRITE (*,FMT="(I4,': ARPDBG: ny = ',I4)") narea-1,nysend(icomm)
1552	WRITE (*,FMT="(I4,': ARPDBG: nz = ',I4)") narea-1,nzsend(icomm)
1553	WRITE (*,FMT="(I4,': ARPDBG:npatch = ',I3)") narea-1,npatches(1)
1554
1555	DO ipatch=1,npatches(1)
1556	WRITE (*,FMT="(I4,': ARPDBG: patch ',I2,': isrc = ',I4)") &
1557	narea-1,ipatch,isrcsendp(ipatch,icomm,1)
1558	WRITE (*,FMT="(I4,': ARPDBG: patch ',I2,': jsrc = ',I4)") &
1559	narea-1,ipatch,jsrcsendp(ipatch,icomm,1)
1560	WRITE (*,FMT="(I4,': ARPDBG: patch ',I2,': nx = ',I4)") &
1561	narea-1,ipatch,nxsendp(ipatch,icomm,1)
1562	WRITE (*,FMT="(I4,': ARPDBG: patch ',I2,': ny = ',I4)") &
1563	narea-1,ipatch,nysendp(ipatch,icomm,1)
1564	WRITE (*,FMT="(I4,': ARPDBG: patch ',I2,': nz = ',I4)") &
1565	narea-1,ipatch,nzsendp(ipatch,icomm,1)
1566	END DO
1567
1568	WRITE (*,FMT="(I4,': ARPDBG:nsendp = ',I4)") narea-1,nsendp(icomm,1)
1569	IF(msgtrim_z)THEN
1570	WRITE (*,FMT="(I4,': ARPDBG:nsendp WITHOUT z trim = ',I4)") &
1571	narea-1,nsendp_untrimmedz
1572	END IF
1573	WRITE (*,FMT="(I4,': ARPDBG SEND ends')") narea-1
1574	#endif
1575
1576	END IF ! something left after message trimming
1577
1578	END SUBROUTINE addsend
1579
1580	SUBROUTINE addrecv ( icomm, dir, proc, isrc, jsrc, &
1581	ides, jdes, nx, ny, depth, ibotlvl, ierr )
1582	!!------------------------------------------------------------------
1583	! Adds a recv communication specified by the parameters dir through
1584	! to ny to the recv communication list at the next position.
1585	! icomm points to the last entry and is incremented and returned
1586	! if successful.
1587	!
1588	! icomm int in/out Location in comms list.
1589	! dir int input Direction.
1590	! proc int input Process id.
1591	! isrc int input X coordinate of source data.
1592	! jsrc int input Y coordinate of source data.
1593	! ides int input X coordinate of dest. data.
1594	! jdes int input Y coordinate of dest. data.
1595	! nx int input Size in X of data to be sent.
1596	! ny int input Size in Y of data to be sent.
1597	! ierr int output Error flag.
1598	!
1599	! Mike Ashworth, CLRC Daresbury Laboratory, March 1999
1600	!!------------------------------------------------------------------
1601	IMPLICIT NONE
1602
1603	! Subroutine arguments.
1604	INTEGER, INTENT(inout) :: icomm
1605	INTEGER, INTENT( in ) :: dir, proc
1606	INTEGER, INTENT(out) :: ierr
1607	INTEGER, DIMENSION(:,:), INTENT( in ) :: depth
1608	INTEGER, DIMENSION(:,:), INTENT( in ) :: ibotlvl
1609	INTEGER, DIMENSION(jpreci) :: isrc, jsrc, ides, jdes, nx, ny
1610
1611	! Local variables.
1612
1613	! Values of corresponding input arguments after clipping
1614	INTEGER, DIMENSION(jpreci) :: cisrc,cjsrc,cides,cjdes,cnx,cny,cnz
1615	! Run-length encoded versions corresponding to above
1616	INTEGER, dimension(MaxPatch,jpreci) :: risrc,rjsrc,rides,rjdes,rnx,rny,rnz
1617	! Number of patches in run-length encoded message
1618	INTEGER, DIMENSION(jpreci) :: npatches
1619	INTEGER :: ihalo, ipatch
1620	! Whether there is still a message after clipping
1621	LOGICAL :: something_left
1622
1623	! Clear the error flag.
1624
1625	ierr = 0
1626
1627	! Return if the process id is not set.
1628
1629	IF ( proc.LT.0 ) THEN
1630	RETURN
1631	ENDIF
1632
1633	! Can the message be clipped ?
1634
1635	CALL clip_msg(depth, ibotlvl, &
1636	ides, jdes, isrc, jsrc, nx, ny, &
1637	cides,cjdes,cisrc,cjsrc,cnx,cny,cnz, &
1638	rides,rjdes,risrc,rjsrc,rnx,rny,rnz, &
1639	npatches, something_left)
1640
1641	IF (something_left) THEN
1642
1643	icomm = icomm+1
1644
1645	! Check that the comms list has space for another entry.
1646
1647	IF ( icomm.GT.MaxComm ) THEN
1648	IF ( lwp ) THEN
1649	WRITE (numout,*) 'ERROR: Number of separate ', &
1650	'communications exceeds maximum of ',MaxComm
1651	ENDIF
1652	ierr = -12
1653	RETURN
1654	ENDIF
1655
1656	! Add the data into the comms list at the new location.
1657
1658	dirrecv(icomm) = dir
1659	source(icomm) = proc
1660	isrcrecv(icomm) = cisrc(1)
1661	jsrcrecv(icomm) = cjsrc(1)
1662	idesrecv(icomm) = cides(1)
1663	jdesrecv(icomm) = cjdes(1)
1664
1665	nxrecv(icomm) = cnx(1)
1666	nyrecv(icomm) = cny(1)
1667	IF(msgtrim_z)THEN
1668	nzrecv(icomm) = cnz(1)
1669	ELSE
1670	nzrecv(icomm) = jpk
1671	END IF
1672
1673	DO ihalo=1,jpreci
1674
1675	nrecvp2d(icomm,ihalo) = 0
1676	nrecvp(icomm,ihalo) = 0
1677	npatchrecv(icomm,ihalo) = npatches(ihalo)
1678
1679	DO ipatch=1,npatches(ihalo)
1680	isrcrecvp(ipatch,icomm,ihalo) = risrc(ipatch,ihalo)
1681	jsrcrecvp(ipatch,icomm,ihalo) = rjsrc(ipatch,ihalo)
1682	idesrecvp(ipatch,icomm,ihalo) = rides(ipatch,ihalo)
1683	jdesrecvp(ipatch,icomm,ihalo) = rjdes(ipatch,ihalo)
1684	nxrecvp(ipatch,icomm,ihalo) = rnx(ipatch,ihalo)
1685	nyrecvp(ipatch,icomm,ihalo) = rny(ipatch,ihalo)
1686	IF(msgtrim_z)THEN
1687	nzrecvp(ipatch,icomm,ihalo) = rnz(ipatch,ihalo)
1688	ELSE
1689	nzrecvp(ipatch,icomm,ihalo) = jpk
1690	END IF
1691
1692	! Sum the no. of points to be received over all
1693	! patches
1694	nrecvp2d(icomm,ihalo) = nrecvp2d(icomm,ihalo) + &
1695	nxrecvp(ipatch,icomm,ihalo)* &
1696	nyrecvp(ipatch,icomm,ihalo)
1697
1698	nrecvp(icomm,ihalo) = nrecvp(icomm,ihalo) + &
1699	nxrecvp(ipatch,icomm,ihalo)* &
1700	nyrecvp(ipatch,icomm,ihalo)* &
1701	nzrecvp(ipatch,icomm,ihalo)
1702	END DO
1703	END DO
1704
1705	#if defined ARPDEBUG
1706	WRITE (*,FMT="(I3,': ARPDBG adding RECV:')") narea-1
1707	WRITE (*,FMT="(I3,': ARPDBG: icomm = ',I2)") narea-1,icomm
1708	WRITE (*,FMT="(I3,': ARPDBG: dir = ',I2)") narea-1,dir
1709	WRITE (*,FMT="(I3,': ARPDBG: proc = ',I4)") narea-1,proc
1710	WRITE (*,FMT="(I3,': ARPDBG: isrc = ',I4)") narea-1,isrcrecv(icomm)
1711	WRITE (*,FMT="(I3,': ARPDBG: jsrc = ',I4)") narea-1,jsrcrecv(icomm)
1712	WRITE (*,FMT="(I3,': ARPDBG: ides = ',I4)") narea-1,idesrecv(icomm)
1713	WRITE (*,FMT="(I3,': ARPDBG: jdes = ',I4)") narea-1,jdesrecv(icomm)
1714	WRITE (*,FMT="(I3,': ARPDBG: nx = ',I4)") narea-1,nxrecv(icomm)
1715	WRITE (*,FMT="(I3,': ARPDBG: ny = ',I4)") narea-1,nyrecv(icomm)
1716	WRITE (*,FMT="(I3,': ARPDBG: nz = ',I4)") narea-1,nzrecv(icomm)
1717	WRITE (*,FMT="(I3,': ARPDBG:npatch = ',I3)") narea-1,npatches(1)
1718	DO ipatch=1,npatches(1)
1719	WRITE (*,FMT="(I3,': ARPDBG: patch ',I2,': isrc = ',I4)") &
1720	narea-1,ipatch,isrcrecvp(ipatch,icomm,1)
1721	WRITE (*,FMT="(I3,': ARPDBG: patch ',I2,': jsrc = ',I4)") &
1722	narea-1,ipatch,jsrcrecvp(ipatch,icomm,1)
1723	WRITE (*,FMT="(I3,': ARPDBG: patch ',I2,': ides = ',I4)") &
1724	narea-1,ipatch,idesrecvp(ipatch,icomm,1)
1725	WRITE (*,FMT="(I3,': ARPDBG: patch ',I2,': jdes = ',I4)") &
1726	narea-1,ipatch,jdesrecvp(ipatch,icomm,1)
1727	WRITE (*,FMT="(I3,': ARPDBG: patch ',I2,': nx = ',I4)") &
1728	narea-1,ipatch,nxrecvp(ipatch,icomm,1)
1729	WRITE (*,FMT="(I3,': ARPDBG: patch ',I2,': ny = ',I4)") &
1730	narea-1,ipatch,nyrecvp(ipatch,icomm,1)
1731	WRITE (*,FMT="(I3,': ARPDBG: patch ',I2,': nz = ',I4)") &
1732	narea-1,ipatch,nzrecvp(ipatch,icomm,1)
1733	END DO
1734	WRITE (*,FMT="(I3,': ARPDBG:nrecvp = ',I4)") narea-1,nrecvp(icomm,1)
1735	WRITE (*,FMT="(I3,': ARPDBG RECV ends')") narea-1
1736	#endif
1737
1738	END IF ! something left
1739
1740	END SUBROUTINE addrecv
1741
1742
1743	SUBROUTINE clip_msg(depth, ibotlvl, &
1744	iloc, jloc, irem, jrem, nx, ny, &
1745	ciloc,cjloc,cirem,cjrem,cnx,cny,cnz, &
1746	riloc,rjloc,rirem,rjrem,rnx,rny,rnz, &
1747	npatches, something_left)
1748	!!------------------------------------------------------------------
1749	!
1750	! Clip any exterior rows or columns that are permanently dry
1751	! from the message. Also remove any vertical levels that are
1752	! beneath the ocean floor.
1753	!
1754	! depth int input Land/sea mask - global coords
1755	! ibotlvl int input Index of the last vertical level
1756	! above sea floor
1757	! iloc int input local X coordinate of data start
1758	! jloc int input local Y coordinate of data start
1759	! irem int input remote X coordinate of data
1760	! jrem int input remote Y coordinate of data
1761	! nx int input Size in X of data to be sent
1762	! ny int input Size in Y of data to be sent
1763	! ciloc int output As iloc, after clipping
1764	! cjloc int output As jloc, after clipping
1765	! cirem int output As irem, after clipping
1766	! cjrem int output As jrem, after clipping
1767	! cnx int output As nx, after clipping
1768	! cny int output As ny, after clipping
1769	!
1770	! The run-length encoded versions split a message into one
1771	! or more patches, leaving out permanently dry rows/columns
1772	!
1773	! riloc int output As iloc, run-length encoded
1774	! rjloc int output As jloc, run-length encoded
1775	! rirem int output As irem, run-length encoded
1776	! rjrem int output As jrem, run-length encoded
1777	! rnx int output As nx, run-length encoded
1778	! rny int output As ny, run-length encoded
1779	! rnz int output Max depth (level) of this patch
1780	! npatches int output Number of patches
1781	!
1782	! something_left logical output
1783	!
1784	! Stephen Pickles, STFC Daresbury Laboratory, August 2009
1785	! - Written
1786	! Stephen Pickles, STFC Daresbury Laboratory, September 2009
1787	! - Added run-length encoding
1788	! Andrew Porter, STFC Daresbury Laboratory, January 2013
1789	! - Added trimming of levels below sea floor
1790	!!------------------------------------------------------------------
1791	USE dom_oce, ONLY: nimpp, njmpp
1792	IMPLICIT none
1793	! Subroutine arguments.
1794	INTEGER, DIMENSION(:,:), INTENT(in) :: depth ! Global mask (0 dry, 1 wet)
1795	INTEGER, DIMENSION(:,:), INTENT(in) :: ibotlvl ! Bottom level of ocean
1796	INTEGER, DIMENSION(jpreci) :: iloc, jloc, irem, jrem, nx, ny
1797	INTEGER, DIMENSION(jpreci) :: ciloc,cjloc,cirem,cjrem,cnx,cny,cnz
1798	INTEGER, DIMENSION(MaxPatch,jpreci) :: riloc,rjloc,rirem,rjrem,rnx,rny,rnz
1799	INTEGER, DIMENSION(jpreci), INTENT(out) :: npatches
1800	LOGICAL, INTENT(out) :: something_left
1801	! Local variables.
1802	INTEGER :: h, i, j, patch
1803	LOGICAL :: all_dry
1804
1805	! i, j, k limits of the halo patch, in local co-ordinates
1806	! These are set from input arguments, then updated as we trim
1807	INTEGER :: ilo, ihi, jlo, jhi
1808
1809	ciloc(:) = iloc(:)
1810	cjloc(:) = jloc(:)
1811	cirem(:) = irem(:)
1812	cjrem(:) = jrem(:)
1813	cnx(:) = nx(:)
1814	cny(:) = ny(:)
1815	cnz(:) = jpk
1816	riloc(1,:) = iloc(:)
1817	rjloc(1,:) = jloc(:)
1818	rirem(1,:) = irem(:)
1819	rjrem(1,:) = jrem(:)
1820	rnx(1,:) = nx(:)
1821	rny(1,:) = ny(:)
1822	rnz(:,:) = jpk
1823	npatches(:) = 1
1824	something_left = .TRUE.
1825
1826	IF (.NOT. msgtrim) RETURN
1827
1828	something_left = .FALSE.
1829
1830	! Loop over halo widths
1831	haloes: DO h=1, jpreci
1832
1833	ilo = iloc(h)
1834	ihi = iloc(h) + cnx(h) - 1
1835	jlo = jloc(h)
1836	jhi = jloc(h) + cny(h) - 1
1837
1838	! Can any points along the left (low i) edge be trimmed?
1839	left_edge: DO i=ilo, ihi
1840	DO j=jlo, jhi
1841	! depth is global mask, i and j are local coords
1842	! ARPDBG - not sure that nextra needed below?
1843	!IF (depth(i+nimpp-1+nextra,j+njmpp-1) .NE. LAND) EXIT left_edge
1844	IF (depth(i+nimpp-1,j+njmpp-1) .NE. LAND) EXIT left_edge
1845	END DO
1846	ciloc(h) = ciloc(h) + 1
1847	cirem(h) = cirem(h) + 1
1848	cnx(h) = cnx(h) - 1
1849	END DO left_edge
1850
1851	IF (cnx(h) .LE. 0) THEN
1852	cnx(h) = 0
1853	cny(h) = 0
1854	cnz(h) = 0
1855	ciloc(h) = iloc(h)
1856	npatches(h) = 0
1857	rnx(1,h) = 0
1858	rny(1,h) = 0
1859	riloc(1,h) = iloc(h)
1860	CYCLE haloes
1861	END IF
1862	ilo = ciloc(h)
1863
1864	! We now know that the trimmed patch must contain at least
1865	! one seapoint for this halo width.
1866	something_left = .TRUE.
1867
1868	! Can any points along the right (high i) edge be trimmed?
1869	right_edge: DO i=ihi, ilo, -1
1870	DO j=jlo, jhi
1871	! IF (depth(i+nimpp-1-nextra,j+njmpp-1) .ne. land) exit right_edge
1872	IF (depth(i+nimpp-1,j+njmpp-1) .ne. land) exit right_edge
1873	END DO
1874	cnx(h) = cnx(h) - 1
1875	END DO right_edge
1876	ihi = ilo + cnx(h) - 1
1877
1878	! Can any points along the bottom (low j) edge be trimmed?
1879	bottom_edge: DO j=jlo, jhi
1880	DO i=ilo, ihi
1881	IF (depth(i+nimpp-1,j+njmpp-1) .ne. land) exit bottom_edge
1882	END do
1883	cjloc(h) = cjloc(h) + 1
1884	cjrem(h) = cjrem(h) + 1
1885	cny(h) = cny(h) - 1
1886	END DO bottom_edge
1887	jlo = cjloc(h)
1888
1889	! Can any points along the top (high j) edge be trimmed?
1890	top_edge: DO j=jhi, jlo, -1
1891	DO i=ilo, ihi
1892	IF (depth(i+nimpp-1,j+njmpp-1) .ne. land) exit top_edge
1893	END do
1894	cny(h) = cny(h) - 1
1895	END DO top_edge
1896	jhi = jlo + cny(h) - 1
1897
1898	! Exterior dry rows and columns have now been clipped from
1899	! the message. We can do something about interior
1900	! dry rows/columns by splitting a message into patches
1901	! of wet points.
1902
1903	! Start first patch
1904	patch = 1
1905	npatches(h) = patch
1906	riloc(1,h) = ilo
1907	rjloc(1,h) = jlo
1908	rirem(1,h) = cirem(h)
1909	rjrem(1,h) = cjrem(h)
1910	rnx(1,h) = 0
1911	rny(1,h) = 0
1912
1913	IF (cnx(h) .GE. cny(h)) THEN
1914
1915	! Longer in x dimension!
1916	i = ilo
1917	rny(1,h) = cny(h)
1918
1919	make_patches_x: DO WHILE (patch .lt. MaxPatch)
1920
1921	IF(i == ihi)THEN
1922
1923	rnx(patch,h) = 1
1924	ELSE
1925
1926	add_sea_cols: DO WHILE (i .lt. ihi)
1927	! Check this strip in y to see whether all points are dry
1928	! !all_dry = .TRUE.
1929	!IF( ANY( depth(ilo+nimpp-1:ihi+nimpp-2,j+njmpp-1) .NE. LAND ) )all_dry = .FALSE.
1930	!IF( ALL( depth(i+nimpp-1,jlo+njmpp-1:jhi+njmpp-1) == LAND ) )EXIT add_sea_cols
1931
1932	all_dry = .TRUE.
1933	DO j=jlo, jhi
1934	IF (depth(i+nimpp-1,j+njmpp-1) .NE. land) THEN
1935	all_dry = .FALSE.
1936	END IF
1937	END DO
1938	IF (all_dry) EXIT add_sea_cols
1939
1940	rnx(patch,h) = rnx(patch,h) + 1
1941	i = i+1
1942	END DO add_sea_cols
1943	END IF
1944
1945	! This patch is now finished.
1946
1947	! Store max depth of ocean bottom in this patch. riloc holds the starting
1948	! point of current patch in local coords.
1949	! riloc(patch,h) + nimpp - 1 is same point in global coords
1950	! End of patch is then at <start> + <length> - 1
1951	rnz(patch,h) = MAXVAL(ibotlvl(riloc(patch,h)+nimpp-1: &
1952	riloc(patch,h)+rnx(patch,h)+nimpp-2, &
1953	jlo+njmpp-1:jhi+njmpp-1) )
1954
1955	! Skip land cols before starting the next patch.
1956
1957	skip_land_cols: DO WHILE (i .lt. ihi)
1958	DO j=jlo, jhi
1959	IF (depth(i+nimpp-1,j+njmpp-1) .NE. land) THEN
1960	EXIT skip_land_cols
1961	END IF
1962	END DO
1963	i = i+1
1964	END DO skip_land_cols
1965
1966	! No more wet points?
1967	IF (i .GE. ihi) EXIT make_patches_x
1968
1969	! Start next patch
1970	patch = patch + 1
1971	npatches(h) = patch
1972	riloc(patch,h) = i
1973	rjloc(patch,h) = jlo
1974	rirem(patch,h) = cirem(h)+i-ilo
1975	rjrem(patch,h) = cjrem(h)
1976	rnx(patch,h) = 0
1977	rny(patch,h) = cny(h)
1978
1979	END DO make_patches_x
1980
1981	! Finish the last patch
1982	rnx(npatches(h),h) = ihi - riloc(npatches(h),h) + 1
1983	rnz(npatches(h),h) = MAXVAL(ibotlvl(riloc(npatches(h),h)+nimpp-1: &
1984	riloc(npatches(h),h)+rnx(npatches(h),h)+nimpp-2, &
1985	jlo+njmpp-1:jhi+njmpp-1) )
1986
1987	ELSE
1988
1989	! Longer in y dimension!
1990	j = jlo
1991	rnx(1,h) = cnx(h)
1992
1993	make_patches_y: DO WHILE (patch .lt. MaxPatch)
1994
1995	add_sea_rows: DO WHILE (j .lt. jhi)
1996
1997	! IF( ALL( depth(ilo+nimpp-1:ihi+nimpp-1, &
1998	! j+njmpp-1) == LAND ) )EXIT add_sea_rows
1999
2000	all_dry = .TRUE.
2001	DO i=ilo, ihi
2002	if (depth(i+nimpp-1,j+njmpp-1) .ne. land) then
2003	all_dry = .FALSE.
2004	end if
2005	END DO
2006	IF (all_dry) EXIT add_sea_rows
2007
2008	rny(patch,h) = rny(patch,h) + 1
2009	j = j+1
2010	END DO add_sea_rows
2011
2012	! This patch is now finished.
2013
2014	! Store max depth of ocean bottom in this patch
2015	rnz(patch,h) = MAXVAL(ibotlvl(ilo+nimpp-1:ihi+nimpp-1, &
2016	rjloc(patch,h)+njmpp-1: &
2017	rjloc(patch,h)+rny(patch,h)+njmpp-2) )
2018
2019	! Skip land rows before starting the next patch.
2020
2021	skip_land_rows: DO WHILE (j .lt. jhi)
2022	DO i=ilo, ihi
2023	IF (depth(i+nimpp-1,j+njmpp-1) .NE. LAND) THEN
2024	EXIT skip_land_rows
2025	END IF
2026	END DO
2027	j = j+1
2028	END DO skip_land_rows
2029
2030	! No more wet points?
2031	IF (j .ge. jhi) EXIT make_patches_y
2032
2033	! Start next patch
2034	patch = patch + 1
2035	npatches(h) = patch
2036	riloc(patch,h) = ilo
2037	rjloc(patch,h) = j
2038	rirem(patch,h) = cirem(h)
2039	rjrem(patch,h) = cjrem(h)+j-jlo
2040	rnx(patch,h) = cnx(h)
2041	rny(patch,h) = 0
2042
2043	END DO make_patches_y
2044
2045	! Finish the last patch
2046	rny(npatches(h),h) = jhi - rjloc(npatches(h),h) + 1
2047	rnz(npatches(h),h) = MAXVAL(ibotlvl(ilo+nimpp-1:ihi+nimpp-1, &
2048	rjloc(npatches(h),h)+njmpp-1: &
2049	rjloc(npatches(h),h)+rny(npatches(h),h)+njmpp-2) )
2050
2051	END IF
2052
2053	! Max depth for whole message is the maximum of the maximum depth of each
2054	! patch.
2055	cnz(h) = MAXVAL(rnz(:,h))
2056
2057	END DO haloes
2058
2059	END SUBROUTINE clip_msg
2060
2061	END MODULE mapcomm_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats: