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

Last change on this file since 3432 was 3432, checked in by trackstand2, 12 years ago
Merge branch 'ksection_partition'
File size: 70.2 KB

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

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

Download in other formats: