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 @ 3837

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

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

Download in other formats: