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.

partition_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/partition_mod.F90 @ 4467

Last change on this file since 4467 was 4467, checked in by trackstand2, 10 years ago
Changes to make global depth mask the same as that generated locally later. Increase nextra from 2 to 3 to get 192 PE case to work.
File size: 110.7 KB

Line
1	MODULE partition_mod
2	USE par_oce, ONLY: jpni, jpnj, jpnij, jpi, jpj, jpim1, jpjm1, jpij, &
3	jpreci, jprecj, jpk, jpkm1, jperio, jpiglo, jpjglo
4	USE dom_oce, ONLY: ln_zco, nbondi, nbondj, nidom, npolj, &
5	nlci, nlcj, & ! i- & j-dimss of the local subdomain
6	nldi, nlei, & ! first and last indoor i- and j-indexes
7	nldj, nlej, & !
8	nlcit, nlcjt,& !
9	nldit, nldjt,& ! (Unity-indexed) arrays storing above
10	nleit, nlejt,& ! values for each domain.
11	nimpp,njmpp, & ! i- & j-indices for mpp-subdom. left bottom
12	nimppt,njmppt,& ! Unity-indexed arrays storing above
13	! values for each domain.
14	nperio, & ! Local periodicity
15	nwidthmax, & ! Width of widest northern domain
16	narea, & ! ID of local area (= rank + 1)
17	mbkmax ! Deepest level above ocean floor
18	#if defined key_mpp_mpi
19	USE lib_mpp, ONLY: mppsize, mppsync, mpi_comm_opa, &
20	MAX_FACTORS, xfactors, yfactors, nn_pttrim, &
21	nn_cpnode, nn_readpart
22	#endif
23	USE lib_mpp, ONLY: ctl_stop, ctl_warn
24	USE in_out_manager, ONLY: numout, lwp
25	USE mapcomm_mod, ONLY: ielb, ieub, mapcomms, pielb, pjelb, pieub, pjeub,&
26	iesub, jesub, jeub, ilbext, iubext, jubext, &
27	jlbext, pnactive, piesub, pjesub, jelb, pilbext, &
28	piubext, pjlbext, pjubext, nextra, &
29	nprocp ! No. of PEs to partition over
30	USE iom, ONLY: wp, jpdom_unknown, iom_open, iom_get, iom_close
31	IMPLICIT NONE
32	PRIVATE
33
34	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:,:) :: imask ! Mask used for partitioning
35	! (1 for ocean, 0 for land)
36	! set in nemogcm.F90
37	! Holds the bottom level of the ocean at each grid point - used for
38	! trimming halos in z direction
39	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:,:), TARGET :: ibotlevel
40
41	! Parameters for the cost function used when evaluating different
42	! possible domain partitions
43
44	! Mnemonics:
45	! p = process (i.e. core)
46	! n = node
47	! l = length (number of points)
48	! c = count (number of messages)
49	! i = internal (intra-node, or on-node)
50	! x = external (inter-node, or off-node)
51
52	INTEGER, PARAMETER :: pv_index_overall = 1
53	INTEGER, PARAMETER :: pv_index_wet = 2
54	INTEGER, PARAMETER :: pv_index_dry = 3
55	INTEGER, PARAMETER :: pv_index_pli = 4
56	INTEGER, PARAMETER :: pv_index_plx = 5
57	INTEGER, PARAMETER :: pv_index_pci = 6
58	INTEGER, PARAMETER :: pv_index_pcx = 7
59	INTEGER, PARAMETER :: pv_index_nli = 8
60	INTEGER, PARAMETER :: pv_index_nlx = 9
61	INTEGER, PARAMETER :: pv_index_nci = 10
62	INTEGER, PARAMETER :: pv_index_ncx = 11
63	INTEGER, PARAMETER :: pv_index_tli = 12
64	INTEGER, PARAMETER :: pv_index_tlx = 13
65	INTEGER, PARAMETER :: pv_index_tci = 14
66	INTEGER, PARAMETER :: pv_index_tcx = 15
67
68	INTEGER, PARAMETER :: pv_index_pcomms = 16
69	INTEGER, PARAMETER :: pv_index_ncomms = 17
70	INTEGER, PARAMETER :: pv_index_tcomms = 18
71
72	INTEGER, PARAMETER :: pv_num_scores = 18
73
74	REAL(wp),PARAMETER :: pv_awful = 1.0e20
75
76	#define PARTIT_DEBUG
77
78	PUBLIC imask, ibotlevel, smooth_global_bathy, global_bot_level, partition_mask_alloc
79	PUBLIC mpp_init3, partition_rk, partition_mca_rk, write_partition_map
80	PUBLIC read_partition, write_partition
81
82	CONTAINS
83
84	SUBROUTINE partition_mask_alloc(xsize, ysize, ierr)
85	!!------------------------------------------------------------------
86	!! * ROUTINE partition_mask_alloc *
87	!!
88	!! Called from nemogcm to allocate the masks that are members of
89	!! this module
90	!!
91	!!------------------------------------------------------------------
92	INTEGER, INTENT(in) :: xsize, ysize
93	INTEGER, INTENT(out):: ierr
94
95	ALLOCATE(imask(xsize,ysize), ibotlevel(xsize,ysize), Stat=ierr)
96
97	END SUBROUTINE partition_mask_alloc
98
99
100	SUBROUTINE mpp_init3()
101	!!------------------------------------------------------------------
102	!! * ROUTINE mpp_init3 *
103	!!
104	!! History:
105	!! Code adapted from POLCOMS code 17/01/2008 A. Porter
106	!!
107	!!------------------------------------------------------------------
108	#if defined key_mpp_mpi
109	!$AGRIF_DO_NOT_TREAT
110	USE mpi ! For better interface checking
111	#endif
112	USE exchtestmod, ONLY : mpp_test_comms
113	IMPLICIT NONE
114	! Local vars
115	INTEGER :: inum ! temporary logical unit
116	INTEGER :: iproc, jproc ! Loop index
117	INTEGER :: ierr ! Error flag
118	INTEGER :: ji, jj
119	CHARACTER(LEN=4) :: intStr
120
121	! Also set original NEMO arrays as they store internal limits of
122	! each domain in local coordinates
123	nldit(narea) = nldi
124	nleit(narea) = nlei
125	nlcit(narea) = nlci
126	nimppt(narea) = nimpp
127	nldjt(narea) = nldj
128	nlejt(narea) = nlej
129	nlcjt(narea) = nlcj
130	njmppt(narea) = njmpp
131	! Make sure all PEs have all these values
132	! ARPDBG - wrap this MPI in lib_mpp ??
133	#if defined key_mpp_mpi
134	CALL MPI_ALLGATHER(nldi,1,MPI_INTEGER, &
135	nldit,1,MPI_INTEGER,mpi_comm_opa,ierr)
136	CALL MPI_ALLGATHER(nlei,1,MPI_INTEGER, &
137	nleit,1,MPI_INTEGER,mpi_comm_opa,ierr)
138	CALL MPI_ALLGATHER(nlci,1,MPI_INTEGER, &
139	nlcit,1,MPI_INTEGER,mpi_comm_opa,ierr)
140	CALL MPI_ALLGATHER(nimpp,1,MPI_INTEGER, &
141	nimppt,1,MPI_INTEGER,mpi_comm_opa,ierr)
142	CALL MPI_ALLGATHER(nldj,1,MPI_INTEGER, &
143	nldjt,1,MPI_INTEGER,mpi_comm_opa,ierr)
144	CALL MPI_ALLGATHER(nlej,1,MPI_INTEGER, &
145	nlejt,1,MPI_INTEGER,mpi_comm_opa,ierr)
146	CALL MPI_ALLGATHER(nlcj,1,MPI_INTEGER, &
147	nlcjt,1,MPI_INTEGER,mpi_comm_opa,ierr)
148	CALL MPI_ALLGATHER(njmpp,1,MPI_INTEGER, &
149	njmppt,1,MPI_INTEGER,mpi_comm_opa,ierr)
150	CALL MPI_ALLGATHER(iesub,1,MPI_INTEGER, &
151	piesub,1,MPI_INTEGER,mpi_comm_opa,ierr)
152	CALL MPI_ALLGATHER(ielb,1,MPI_INTEGER, &
153	pielb,1,MPI_INTEGER,mpi_comm_opa,ierr)
154	CALL MPI_ALLGATHER(ieub,1,MPI_INTEGER, &
155	pieub,1,MPI_INTEGER,mpi_comm_opa,ierr)
156	CALL MPI_ALLGATHER(jesub,1,MPI_INTEGER, &
157	pjesub,1,MPI_INTEGER,mpi_comm_opa,ierr)
158	CALL MPI_ALLGATHER(jelb,1,MPI_INTEGER, &
159	pjelb,1,MPI_INTEGER,mpi_comm_opa,ierr)
160	CALL MPI_ALLGATHER(jeub,1,MPI_INTEGER, &
161	pjeub,1,MPI_INTEGER,mpi_comm_opa,ierr)
162	#endif
163
164	IF(lwp)THEN
165	! Write out domains in postscript
166
167	OPEN(UNIT=997, FILE="domain_decomp.ps", &
168	STATUS='REPLACE', ACTION='WRITE', IOSTAT=iproc)
169
170	IF(iproc .EQ. 0)THEN ! Check file opened OK
171
172	! Header of ps file
173	WRITE (997,FMT='("%!PS-Adobe-1.0")')
174	WRITE (997,FMT='("/Helvetica findfont %load the font dictionary")')
175	WRITE (997,FMT='("12 scalefont %scale to 12pt")')
176	WRITE (997,FMT='("setfont %make this the current font")')
177	WRITE (997,FMT='("/u { 2 mul } def %set up a scaling factor")')
178
179	! Put green circles at dry points
180	WRITE (997,FMT='("% Filled green circles at dry points")')
181	WRITE (997,FMT='("0.1 setlinewidth")') ! Thin green line
182	WRITE (997,FMT='("0 1 0 setrgbcolor")')
183	WRITE (997,FMT='("newpath")')
184	DO jj = 1,jpjglo,1
185	DO ji = 1,jpiglo,1
186	IF(imask(ji,jj) /= 1)THEN
187	WRITE (997,FMT='(I3," u ",I3," u 1 u 0 360 arc closepath")') ji, jj
188	! Use 'fill' instead of 'stroke' to get filled circles
189	WRITE (997,FMT='("fill")')
190	END IF
191	END DO
192	END DO
193
194	! Draw the outline of the global domain in red
195	WRITE (997,FMT='("% Draw the outline of the global domain in red")')
196	WRITE (997,FMT='("3.0 setlinewidth")') ! Fat line of 3mm for global dom.
197	WRITE (997,FMT='("1 0 0 setrgbcolor")')
198	WRITE (997,FMT='("newpath")')
199	WRITE (997,FMT='("% Cursor initialization")')
200	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') 1,1
201	WRITE (997,FMT='("% Drawing the rectangle")')
202	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo,1
203	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo,jpjglo
204	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') 1,jpjglo
205	WRITE (997,FMT='("closepath stroke")')
206
207	! Now draw the outline of each individual domain, nprocp should have been
208	! set at the very beginning of the recursive partioning process.
209	WRITE (997,FMT='("0.7 setlinewidth")') ! Back to default width
210	WRITE (997,FMT='("0 0 0 setrgbcolor")')! and colour
211	DO iproc=1,nprocp
212	WRITE (997,FMT='("newpath")')
213	WRITE (997,FMT='("% Cursor initialization")')
214	WRITE (997,FMT='(I3," u ",1x,I3," u moveto %BL Corner")') pielb(iproc),pjelb(iproc)
215	WRITE (997,FMT='("% Drawing the rectangle")')
216	WRITE (997,FMT='(I3," u ",1x,I3," u lineto %BR Corner")') pieub(iproc),pjelb(iproc)
217	WRITE (997,FMT='(I3," u ",1x,I3," u lineto %TR Corner")') pieub(iproc),pjeub(iproc)
218	WRITE (997,FMT='(I3," u ",1x,I3," u lineto %TL Corner")') pielb(iproc),pjeub(iproc)
219	WRITE (997,FMT='("closepath stroke")')
220	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') &
221	INT(0.5*(pieub(iproc)+pielb(iproc))), &
222	INT(0.5*(pjeub(iproc)+pjelb(iproc)-1))
223	! Write processor label, left justified
224	WRITE (intStr, FMT='(I4)') iproc-1
225	! Use postscipt operator 'stringwidth' to get the approximate width of the
226	! string containing the PE id and then adjust position so it is centred
227	! about point we've just moveto'd. This doesn't attempt to deal with
228	! vertical centering.
229	WRITE (997,FMT='("(",(A),") dup stringwidth pop 2 div neg 0 rmoveto show")') TRIM(ADJUSTL(intStr))
230	END DO
231	WRITE (997,FMT='("showpage")')
232	WRITE (997,FMT='("%%EOF")')
233	CLOSE(997)
234
235	END IF ! File opened OK
236	END IF ! lwp
237
238	! Test for overlaps of domains (there shouldn't be any!)
239	DO iproc=1, nprocp,1
240	DO jproc=iproc+1, nprocp, 1
241	IF( ( ( (pielb(iproc) .LE. pieub(jproc)) .AND. &
242	(pielb(iproc) .GE. pielb(jproc)) ) &
243	.OR. &
244	( (pieub(iproc) .LE. pieub(jproc)) .AND. &
245	(pieub(iproc) .GE. pielb(jproc)) ) ) &
246	.AND. &
247	( ( (pjelb(iproc) .LE. pjeub(jproc)) .AND. &
248	(pjelb(iproc) .GE. pjelb(jproc)) ) &
249	.OR. &
250	( (pjeub(iproc) .LE. pjeub(jproc)) .AND. &
251	(pjeub(iproc) .GE. pjelb(jproc)) ) &
252	) )THEN
253	WRITE(*,"('ERROR: domains ',I3,' and ',I3,' overlap!')") &
254	iproc, jproc
255	WRITE(*,"(I3,': [',I3,':',I3,'][',I3,':',I3,']')") &
256	iproc, pielb(iproc), pieub(iproc), &
257	pjelb(iproc), pjeub(iproc)
258	WRITE(*,"(I3,': [',I3,':',I3,'][',I3,':',I3,']')") &
259	jproc, pielb(jproc), pieub(jproc), pjelb(jproc), pjeub(jproc)
260
261	END IF
262	END DO
263	END DO
264
265	! Map out the communications for the partitioned domain.
266	CALL mapcomms (imask, ibotlevel, jpiglo, jpjglo, jperio, ierr)
267	IF ( ierr.NE.0 ) THEN
268	IF ( lwp ) WRITE(numout,*) 'Communications mapping failed : ',ierr
269	RETURN
270	ENDIF
271
272	! Prepare mpp north fold
273	#if defined key_mpp_mpi
274	! This invokes the version of the routine contained in this module
275	! and not the original in lib_mpp.F90
276	CALL mpp_ini_north()
277	#endif
278
279	! From mppini_2.h90:
280	! Defined npolj, either 0, 3 , 4 , 5 , 6
281	! In this case the important thing is that npolj /= 0
282	! Because if we go through these line it is because jpni >1 and thus
283	! we must use lbcnorthmpp, which tests only npolj =0 or npolj /= 0
284	npolj = 0
285	IF( jperio == 3 .OR. jperio == 4 ) THEN
286
287	IF( pjubext(narea) ) npolj = 3 ! Top row?
288	IF( pjubext(narea) .AND. piubext(narea) ) npolj = 4 ! Top right? ARPDBG almost certainly wrong!
289	ENDIF
290	IF( jperio == 5 .OR. jperio == 6 ) THEN
291	IF( pjubext(narea) ) npolj = 5 ! Top left?
292	IF( pjubext(narea) .AND. piubext(narea) ) npolj = 6 ! Top right? ARPDBG almost certainly wrong!
293	ENDIF
294
295	! Prepare NetCDF output file (if necessary)
296	CALL mpp_init_ioipsl()
297
298	! Write this partition to file in the format that the code can
299	! read
300	CALL write_partition()
301
302	! Initialise mbkmax because it's needed in any halo swap calls but we didn't have
303	! jpi and jpj until the partitioning was done.
304	DO jj=1,jpj,1
305	DO ji=1,jpi,1
306	mbkmax(ji,jj) = ibotlevel(MIN(jpiglo,ji+nimpp-1), &
307	MIN(jpjglo,jj+njmpp-1))
308	END DO
309	END DO
310
311	! ARPDBG - test comms setup
312	CALL mpp_test_comms(imask, ibotlevel)
313
314	! Free array holding mask used for partitioning
315	DEALLOCATE(imask)
316
317	END SUBROUTINE mpp_init3
318
319	# if defined key_dimgout
320	!!----------------------------------------------------------------------
321	!! 'key_dimgout' NO use of NetCDF files
322	!!----------------------------------------------------------------------
323	SUBROUTINE mpp_init_ioipsl ! Dummy routine
324	END SUBROUTINE mpp_init_ioipsl
325	# else
326	SUBROUTINE mpp_init_ioipsl
327	!!----------------------------------------------------------------------
328	!! * ROUTINE mpp_init_ioipsl *
329	!!
330	!! ** Purpose :
331	!!
332	!! ** Method :
333	!!
334	!! History :
335	!! 9.0 ! 04-03 (G. Madec) MPP-IOIPSL
336	!!----------------------------------------------------------------------
337	USE ioipsl
338	IMPLICIT NONE
339	!! Local declarations
340	INTEGER, DIMENSION(2) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
341	!!----------------------------------------------------------------------
342
343	! The domain is decomposed in 2D only along i- or/and j- direction
344	! So we need at the most only 1D arrays with 2 elements
345	iglo(1) = jpiglo
346	iglo(2) = jpjglo
347	iloc(1) = iesub
348	iloc(2) = jesub
349	iabsf(1) = ielb
350	iabsf(2) = jelb
351	iabsl(:) = iabsf(:) + iloc(:) - 1
352	ihals(1) = jpreci
353	ihals(2) = jprecj
354	ihale(1) = jpreci
355	ihale(2) = jprecj
356	idid(1) = 1
357	idid(2) = 2
358
359	IF(lwp) THEN
360	WRITE(numout,*)
361	WRITE(numout,*) 'partmod: mpp_init_ioipsl : iloc = ', iloc (1), iloc (2)
362	WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~ iabsf = ', iabsf(1), iabsf(2)
363	WRITE(numout,*) ' ihals = ', ihals(1), ihals(2)
364	WRITE(numout,*) ' ihale = ', ihale(1), ihale(2)
365	ENDIF
366
367	#if defined key_mpp_mpi
368	CALL flio_dom_set ( mppsize, narea-1, idid, iglo, iloc, iabsf, iabsl, &
369	ihals, ihale, 'BOX', nidom)
370	#endif
371
372	END SUBROUTINE mpp_init_ioipsl
373	# endif
374
375	SUBROUTINE partition_rk ( mask, istart, istop, jstart, jstop, ierr )
376	!!------------------------------------------------------------------
377	! Partition the domain of nx x ny points
378	! according to the land/sea mask array
379	! using a recursive k-section algorithm,
380	! into nprocx x nprocy sub-domains.
381	!
382	! mask real input Land/sea mask.
383	! gnx int input Size of the full domain.
384	! gny int input
385	! ierr int output Error flag.
386	!!------------------------------------------------------------------
387
388	USE iom, ONLY: jpiglo, jpjglo, wp
389	USE par_oce, ONLY: jpni, jpnj
390	#if defined key_mpp_mpi
391	USE lib_mpp, ONLY: mppsize
392	#endif
393	IMPLICIT NONE
394
395	! Subroutine arguments.
396	INTEGER, INTENT(out) :: ierr
397	INTEGER, INTENT(in) :: istart, istop, jstart, jstop
398	INTEGER, INTENT(in) :: mask(:,:)
399	! Local variables
400	#if defined key_mpp_mpi
401	INTEGER, DIMENSION(MAX_FACTORS) :: fx,fy
402	#endif
403	INTEGER :: f,gnactive &
404	,i,ifax,ifin,ifx,ify,ilb,iproc,ist,isub,isub_old &
405	,isub_new,iub &
406	,j,jfin,jlb,jst,jub,line &
407	,nfax,nfx,nfy,ngone,nsub_old,nsub_new,ntarget,ntry
408	LOGICAL :: partx
409
410	! Clear the error flag.
411	ierr = 0
412
413	#if defined key_mpp_mpi
414
415	#if defined PARTIT_DEBUG
416	IF(lwp)WRITE(,) 'ARPDBG partition_rk: jpn{i,j} = ',jpni,jpnj
417	#endif
418	! Factorise the nprocx and nprocy parameters.
419	CALL factor (fx,MAX_FACTORS,nfx,jpni,ierr)
420	IF ( lwp .AND. ierr.NE.0 ) THEN
421	WRITE (numout,*) 'partition_rk: factorisation in x failed ',ierr
422	RETURN
423	ENDIF
424	CALL factor (fy,MAX_FACTORS,nfy,jpnj,ierr)
425	IF ( lwp .AND. ierr.NE.0 ) THEN
426	WRITE (numout,*) 'partition_rk: factorisation in y failed ',ierr
427	RETURN
428	ENDIF
429
430	#if defined PARTIT_DEBUG
431	IF(lwp)THEN
432	WRITE(,) 'ARPDBG partition_rk: nf{x,y} = ',nfx,nfy
433	WRITE(,) 'ARPDBG partition_rk: fx = ',fx(:nfx)
434	WRITE(,) 'ARPDBG partition_rk: fy = ',fx(:nfy)
435	END IF
436	#endif
437
438	CALL partition_rk_core(mask, jpiglo, jpjglo, &
439	MAX_FACTORS, fx, nfx, fy, nfy, ierr)
440
441	CALL finish_partition()
442
443	#endif
444
445	END SUBROUTINE partition_rk
446
447
448	SUBROUTINE partition_mca_rk(mask, istart, istop, jstart, jstop, ierr)
449	#if defined key_mpp_mpi
450	USE mpi
451	USE lib_mpp, ONLY: mppsize, mpi_comm_opa, &
452	nxfactors, nyfactors, xfactors, yfactors
453	#endif
454	USE lib_mpp, ONLY: ctl_stop
455	USE dom_oce, ONLY: narea
456	IMPLICIT NONE
457	!!------------------------------------------------------------------
458	!! Multi-Core Aware recursive partitioning of the domain. As for
459	!! partition_rk but choose the partion
460	!! so as to minimize off-node MPI communication
461	!!
462	!! Original by Stephen Pickles for POLCOMS code.
463	!! Implementation in NEMO by Andrew Porter, 26/01/2012
464	!!------------------------------------------------------------------
465	! Subroutine arguments.
466	INTEGER, INTENT(in) :: istart, istop, jstart, jstop
467	INTEGER, INTENT(in) :: mask(:,:)
468	INTEGER, INTENT(out) :: ierr
469	! Local variables
470	INTEGER :: ii
471	#if defined key_mpp_mpi
472	INTEGER, DIMENSION(MAX_FACTORS) :: fx, fy, factors
473	INTEGER, DIMENSION(MAX_FACTORS) :: df, multiplicity
474	#endif
475	INTEGER :: nfx, nfy, nfactors, ndf, nperms
476	INTEGER :: check_nprocx, check_nprocy, check_nprocp
477	INTEGER :: iperm
478	CHARACTER(LEN=256) :: fstr
479	INTEGER :: myinst ! MPI process ID of this process
480	INTEGER :: nprocx, nprocy
481
482	! A high score is bad
483	REAL(wp), DIMENSION(pv_num_scores) :: score
484	REAL(wp) :: best_score
485	INTEGER :: best_perm
486	REAL(wp), DIMENSION(2,pv_num_scores) :: best, gbest, wrst, gwrst
487
488	#if defined key_mpp_mpi
489
490	! NEMO only has narea public and not the actual PE rank so
491	! set that here
492	myinst = narea - 1
493
494	! Factorise the total number of MPI processes that we have
495	CALL factor (factors,MAX_FACTORS,nfactors,nprocp,ierr)
496
497	IF ( lwp ) THEN
498	IF ( ierr.NE.0 ) THEN
499	WRITE (numout,*) 'partition_mca_rk: factorisation failed ',ierr
500	RETURN
501	ELSE
502	WRITE (numout,*) 'partition_mca_rk: factors are: ', factors(:nfactors)
503	ENDIF
504	ENDIF
505
506	CALL calc_perms( nfactors, factors, &
507	ndf, df, multiplicity, &
508	nperms )
509
510	DO ii=1,pv_num_scores
511	best(1,ii) = pv_awful
512	best(2,ii) = -1.0_wp
513	END DO
514	DO ii=1,pv_num_scores
515	wrst(1,ii) = 0.0_wp
516	wrst(2,ii) = -1.0_wp
517	END DO
518
519	IF (lwp) THEN
520	WRITE(numout,"('% Partn',2X,10(A4,2X),4(A9,1X),A7)") &
521	'Wet', 'Dry', &
522	'pli', 'plx', 'pci', 'pcx', &
523	'nlx', 'ncx', 'tlx', 'tcx', &
524	'P comms', 'N comms', 'T comms', 'Overall', &
525	'Factors'
526	END IF
527
528	perm_loop: DO iperm=myinst, nperms-1, nprocp
529
530	CALL get_perm_factors( iperm, nfactors, ndf, df, multiplicity, &
531	fx, nfx, fy, nfy, &
532	nprocx, nprocy, .FALSE. )
533
534	CALL partition_rk_core(mask, jpiglo, jpjglo, &
535	MAX_FACTORS, fx, nfx, fy, nfy, ierr)
536
537	IF (ierr .NE. 0) CYCLE perm_loop
538	CALL finish_partition()
539
540	! Compute the cost function for this partition
541	!
542	CALL eval_partition(jpiglo, jpjglo, mask, score)
543	CALL factor_string(fstr,nfx,fx,nfy,fy)
544
545	WRITE (6,'(''%'',I6,1X,10(I5,1X),3(F9.2,1X),E12.4,1x,(A))') &
546	iperm, &
547	INT(score(pv_index_wet)), &
548	INT(score(pv_index_dry)), &
549	INT(score(pv_index_pli)), &
550	INT(score(pv_index_plx)), &
551	INT(score(pv_index_pci)), &
552	INT(score(pv_index_pcx)), &
553	INT(score(pv_index_nlx)), &
554	INT(score(pv_index_ncx)), &
555	INT(score(pv_index_tlx)), &
556	INT(score(pv_index_tcx)), &
557	score(pv_index_pcomms), &
558	score(pv_index_ncomms), &
559	score(pv_index_tcomms), &
560	score(pv_index_overall), &
561	TRIM(fstr)
562
563	DO ii=1,pv_num_scores
564	IF (score(ii) .LT. best(1,ii)) THEN
565	best(1,ii) = score(ii)
566	best(2,ii) = iperm
567	END IF
568	IF (score(ii) .GT. wrst(1,ii)) THEN
569	wrst(1,ii) = score(ii)
570	wrst(2,ii) = iperm
571	END IF
572	END DO
573
574	END DO perm_loop
575
576	! Now choose the "best" partition
577
578	#if defined key_mpp_mpi
579	CALL MPI_ALLREDUCE(best, gbest, pv_num_scores, &
580	MPI_2DOUBLE_PRECISION, &
581	MPI_MINLOC, mpi_comm_opa, ierr)
582	CALL MPI_ALLREDUCE(wrst, gwrst, pv_num_scores, &
583	MPI_2DOUBLE_PRECISION, &
584	MPI_MAXLOC, mpi_comm_opa, ierr)
585	#else
586	CALL ctl_stop('STOP', 'partition_mca_rk: this version requires MPI')
587	#endif
588	best_score = gbest(1,pv_index_overall)
589	best_perm = gbest(2,pv_index_overall)
590
591	IF ( lwp ) THEN
592	WRITE (numout,'(A32,A20,A20)') &
593	' ',' best score / perm ',' worst score / perm'
594	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') 'overall: ', &
595	gbest(1,pv_index_overall), gbest(2,pv_index_overall), &
596	gwrst(1,pv_index_overall), gwrst(2,pv_index_overall)
597	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') 'wet points: ', &
598	gbest(1,pv_index_wet), gbest(2,pv_index_wet), &
599	gwrst(1,pv_index_wet), gwrst(2,pv_index_wet)
600	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') 'dry points: ', &
601	gbest(1,pv_index_dry), gbest(2,pv_index_dry), &
602	gwrst(1,pv_index_dry), gwrst(2,pv_index_dry)
603	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
604	'proc max on-node wet points: ', &
605	gbest(1,pv_index_pli), gbest(2,pv_index_pli), &
606	gwrst(1,pv_index_pli), gwrst(2,pv_index_pli)
607	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
608	'proc max off-node wet points: ', &
609	gbest(1,pv_index_plx), gbest(2,pv_index_plx), &
610	gwrst(1,pv_index_plx), gwrst(2,pv_index_plx)
611	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
612	'proc max on-node messages: ', &
613	gbest(1,pv_index_pci), gbest(2,pv_index_pci), &
614	gwrst(1,pv_index_pci), gwrst(2,pv_index_pci)
615	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
616	'proc max off-node messages: ', &
617	gbest(1,pv_index_pcx), gbest(2,pv_index_pcx), &
618	gwrst(1,pv_index_pcx), gwrst(2,pv_index_pcx)
619	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
620	'node max off-node wet points: ', &
621	gbest(1,pv_index_nlx), gbest(2,pv_index_nlx), &
622	gwrst(1,pv_index_nlx), gwrst(2,pv_index_nlx)
623	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
624	'node max off-node messages: ', &
625	gbest(1,pv_index_ncx), gbest(2,pv_index_ncx), &
626	gwrst(1,pv_index_ncx), gwrst(2,pv_index_ncx)
627	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
628	'total off-node wet points: ', &
629	gbest(1,pv_index_tlx), gbest(2,pv_index_tlx), &
630	gwrst(1,pv_index_tlx), gwrst(2,pv_index_tlx)
631	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
632	'per core communications cost: ', &
633	gbest(1,pv_index_pcomms), gbest(2,pv_index_pcomms), &
634	gwrst(1,pv_index_pcomms), gwrst(2,pv_index_pcomms)
635	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
636	'per node communications cost: ', &
637	gbest(1,pv_index_ncomms), gbest(2,pv_index_ncomms), &
638	gwrst(1,pv_index_ncomms), gwrst(2,pv_index_ncomms)
639	WRITE (numout,'(A32,F12.3,F8.0,E12.3,F8.0)') &
640	'network communications cost: ', &
641	gbest(1,pv_index_tcomms), gbest(2,pv_index_tcomms), &
642	gwrst(1,pv_index_tcomms), gwrst(2,pv_index_tcomms)
643
644	WRITE (numout,"('partition_mca_rk: overall best perm is ',I6,', score=',F12.3)") &
645	best_perm, best_score
646	END IF
647
648	! Use the same partition on all processes
649
650	! If a particular factorisation has been forced by
651	! the nn_{x,y}factors fields in the nammpp section of the namelist
652	! then use that one instead
653
654	IF ((nxfactors + nyfactors) > 0) THEN
655	check_nprocx = 1
656	check_nprocy = 1
657	DO ii=1,nxfactors
658	check_nprocx = check_nprocx*xfactors(ii)
659	END DO
660	DO ii=1,nyfactors
661	check_nprocy = check_nprocy*yfactors(ii)
662	END DO
663	check_nprocp = check_nprocx*check_nprocy
664	IF (check_nprocp .EQ. nprocp) THEN
665	nprocx = check_nprocx
666	nprocy = check_nprocy
667	nfx = nxfactors
668	nfy = nyfactors
669	fx(1:nfx) = xfactors(1:nfx)
670	fy(1:nfy) = yfactors(1:nfy)
671	ELSE
672	CALL ctl_stop('STOP', 'part_mca_rk: '// &
673	'requested factorisation does not match nprocp')
674	END IF
675	ELSE
676	! Use the best factorisation found above
677	IF (best_perm < 0.0_wp) THEN
678	IF (lwp) THEN
679	WRITE (numout,*) 'partition_mca_rk: no feasible partition found'
680	END IF
681	ierr = 1
682	RETURN
683	END IF
684	CALL get_perm_factors(best_perm, nfactors, ndf, df, multiplicity, &
685	fx, nfx, fy, nfy, &
686	nprocx, nprocy, lwp )
687	END IF
688
689	! Set corresponding NEMO variables for PE grid, even though it is now
690	! rather irregular
691	jpni = nprocx
692	jpnj = nprocy
693	jpnij = jpni*jpnj
694
695	IF (lwp) THEN
696	WRITE (numout,'(A39)',advance='no') &
697	'partition_mca_rk: partitioning with factors '
698	CALL print_factors(numout, nfx, fx, nfy, fy)
699	END IF
700
701
702	CALL partition_rk_core ( mask, jpiglo, jpjglo, &
703	MAX_FACTORS, &
704	fx, nfx, fy, nfy, &
705	ierr )
706
707	IF (ierr .NE. 0) THEN
708	IF (lwp) THEN
709	WRITE (numout,*) 'partition_mca_rk: failed to apply selection partition'
710	END IF
711	RETURN
712	END IF
713	CALL finish_partition()
714
715	IF (lwp) THEN
716	CALL eval_partition(jpiglo, jpjglo, mask, score)
717	CALL factor_string(fstr,nfx,fx,nfy,fy)
718	WRITE(numout,'(10(A7,1X),4(A9,1X),A7)') &
719	'Wet', 'Dry', &
720	'pli', 'plx', 'pci', 'pcx', &
721	'nlx', 'ncx', 'tlx', 'tcx', &
722	'P comms', 'N comms', 'T comms', 'Overall', &
723	'Factors'
724	WRITE (6,'(10(F7.0,1X),4(F9.2,1X),A50)') &
725	score(pv_index_wet), &
726	score(pv_index_dry), &
727	score(pv_index_pli), &
728	score(pv_index_plx), &
729	score(pv_index_pci), &
730	score(pv_index_pcx), &
731	score(pv_index_nlx), &
732	score(pv_index_ncx), &
733	score(pv_index_tlx), &
734	score(pv_index_tcx), &
735	score(pv_index_pcomms), &
736	score(pv_index_ncomms), &
737	score(pv_index_tcomms), &
738	score(pv_index_overall), &
739	fstr
740	END IF
741
742	#endif
743
744	END SUBROUTINE partition_mca_rk
745
746
747	SUBROUTINE partition_rk_core( mask, nx, ny, maxfax, fx, nfx, fy, nfy, &
748	ierr )
749	#if defined key_mpp_mpi
750	USE lib_mpp, ONLY: mppsize
751	#endif
752	IMPLICIT NONE
753	!!------------------------------------------------------------------
754	!!------------------------------------------------------------------
755	INTEGER, INTENT(in) :: nx, ny
756	INTEGER, INTENT(in) :: mask(nx,ny)
757	INTEGER, INTENT(in) :: maxfax, nfx, nfy
758	INTEGER, INTENT(in) :: fx(maxfax), fy(maxfax)
759	INTEGER, INTENT(out) :: ierr
760	! Local variables
761	INTEGER :: istart, jstart, istop, jstop
762	INTEGER :: f,gnactive
763	INTEGER :: i,ifax,nfax,ifin,ifx,ify,ilb,iproc
764	INTEGER :: ist,isub,isub_old,isub_new,iub
765	INTEGER :: j,jfin,jlb,jst,jub,line
766	INTEGER :: ngone,nsub_old,nsub_new,ntarget,ntry
767	LOGICAL :: partx
768
769	! Zero the error flag
770	ierr = 0
771
772	! Count the significant (non-zero) factors.
773	nfax = nfx+nfy
774	#if defined PARTIT_DEBUG
775	IF ( lwp ) THEN
776	WRITE (numout,"('jpni = ',I3,1x,'nfx = ',I3,/'Factorized into: ',(I3,1x))") &
777	jpni, nfx, fx(:nfx)
778	WRITE (numout,"('jpnj = ',I3,1x,'nfy = ',I3,/'Factorized into: ',(I3,1x))") &
779	jpnj, nfy, fy(:nfy)
780	WRITE (numout,"('Partitioning a total of ',I3,' times')") nfax
781	ENDIF
782	#endif
783
784	! Define the full domain as the one and only sub-domain.
785	istart = 1
786	istop = nx
787	jstart = 1
788	jstop = ny
789
790	nsub_old = 1
791	pielb(nsub_old) = istart
792	pjelb(nsub_old) = jstart
793	pieub(nsub_old) = istop
794	pjeub(nsub_old) = jstop
795	! Count the number of active points.
796	gnactive = 0
797	DO j=jstart,jstop,1
798	DO i=istart,istop,1
799	IF ( mask(i,j) == 1 ) gnactive = gnactive+1
800	ENDDO
801	ENDDO
802	#if defined PARTIT_DEBUG
803	IF ( lwp )WRITE (numout,*) 'Total wet points ',gnactive
804	#endif
805	pnactive(nsub_old) = gnactive
806
807	! Partition for each factor.
808	DO ifax=1,nfax
809	IF ( ifax.EQ.1 ) THEN
810	! Start by partitioning the dimension with more factors.
811	partx = nfx.GE.nfy
812	ifx = 0
813	ify = 0
814	ELSE
815	! Try to toggle the partitioning direction.
816	partx = .NOT.partx
817	IF ( partx ) THEN
818	! If we have run out of factors in x, switch to y.
819	partx = .NOT. ifx+1.GT.nfx
820	ELSE
821	! If we have run out of factors in y, switch to x.
822	partx = ify+1.GT.nfy
823	ENDIF
824	ENDIF
825	#if defined PARTIT_DEBUG
826	IF ( lwp ) THEN
827	WRITE (numout,*)
828	WRITE (numout,*) '###########################################'
829	WRITE (numout,*)
830	WRITE (numout,*) 'Partition ',ifax,'partx = ',partx
831	ENDIF
832	#endif
833	IF ( partx ) THEN
834	! ============================================================
835	! Partition in x.
836	! ============================================================
837	ifx = ifx+1
838	f = fx(ifx)
839	nsub_new = nsub_old*f
840	#if defined PARTIT_DEBUG
841	IF ( lwp ) THEN
842	WRITE (numout,*) 'Partitioning in x from ',nsub_old &
843	,' to ',nsub_new
844	ENDIF
845	#endif
846	DO isub_old=1,nprocp,nprocp/nsub_old
847	! Check that there are sufficient points to factorise.
848	IF ( pieub(isub_old)-pielb(isub_old)+1.LT.f ) THEN
849	WRITE (numout,*) &
850	'partition_rk: insufficient points to partition'
851	ierr = 1
852	RETURN
853	ENDIF
854
855	! Set the target number of points in the new sub-domains.
856	ntarget = NINT(REAL(pnactive(isub_old))/REAL(f))
857	ist = pielb(isub_old)
858	iub = pieub(isub_old)
859	jlb = pjelb(isub_old)
860	jub = pjeub(isub_old)
861	#if defined PARTIT_DEBUG
862	IF ( lwp ) THEN
863	! WRITE (numout,*)
864	WRITE (numout,"(/'=======================================')")
865	! WRITE (numout,*)
866	WRITE (numout,"(/'Partitioning sub-domain ',I3,': (',I3,':',I3,')(',I3,':',I3,')')") &
867	isub_old,ist,iub,jlb,jub
868	WRITE (numout,"('Old partition has ',I8,' points')") pnactive(isub_old)
869	WRITE (numout,"('Target is ',I8,' points')") ntarget
870	ENDIF
871	#endif
872	! Create the new sub-domains.
873	ngone = 0
874	DO isub=1,f,1
875	isub_new = isub_old + (isub-1)*nprocp/nsub_new
876	#if defined PARTIT_DEBUG
877	IF ( lwp ) THEN
878	WRITE (numout,*)
879	WRITE (numout,*) 'Creating new sub-domain ',isub_new
880	ENDIF
881	#endif
882	! The last new sub-domain takes the remaining data.
883	IF ( isub.EQ.f ) THEN
884	ifin = iub
885	ELSE
886	! Scan lines in x counting active grid points until we
887	! exceed the target.
888	ntry = 0
889	ifin = -1
890	scanrows: DO i=ist,iub,1
891	! Count up the contribution from the next line.
892	line = 0
893	DO j=jlb,jub,1
894	IF ( mask(i,j) == 1 ) line = line+1
895	ENDDO
896
897	! Check against the target.
898	IF ( ntry+line.GT.ntarget ) THEN
899	! Is it nearer the target to include this line or not ?
900	IF ( ntry+line-ntarget.GT.ntarget-ntry ) THEN
901	! Nearer start than end. Finish at previous line.
902	ifin = i-1
903	ELSE
904	! Nearer end than start. Include this line.
905	ifin = i
906	ntry = ntry + line
907	ENDIF
908	#if defined PARTIT_DEBUG
909	IF ( lwp ) THEN
910	WRITE (numout,*) 'Reached target at ',ifin &
911	,' ntry = ',ntry
912	ENDIF
913	#endif
914	EXIT scanrows
915	ENDIF
916	! Add in the line count to the running total.
917	ntry = ntry + line
918	ENDDO scanrows
919	IF ( ifin.LT.0 ) ifin = iub
920	ENDIF
921
922	! Set up the parameters of the new sub-domain.
923	pnactive(isub_new) = 0
924	DO j=jlb,jub
925	DO i=ist,ifin
926	IF ( mask(i,j) == 1 ) THEN
927	pnactive(isub_new) = pnactive(isub_new)+1
928	ENDIF
929	ENDDO
930	ENDDO
931	pielb(isub_new) = ist
932	pieub(isub_new) = ifin
933	pjelb(isub_new) = jlb
934	pjeub(isub_new) = jub
935	#if defined PARTIT_DEBUG
936	IF ( lwp ) THEN
937	WRITE (numout,*) 'New subdomain is ',ist,ifin,jlb,jub &
938	,' with ',pnactive(isub_new),' points'
939	ENDIF
940	#endif
941	! Reset the target based on the actual number of points
942	! remaining, split between the partitions remaining.
943	ngone = ngone+ntry
944	#if defined PARTIT_DEBUG
945	IF ( lwp ) THEN
946	! write (numout,*) 'Target reset to ',ntarget
947	ENDIF
948	#endif
949	! Start the next search at the next point.
950	ist = ifin+1
951	ENDDO
952	ENDDO
953
954	ELSE
955
956	! ============================================================
957	! Partition in y.
958	! ============================================================
959	ify = ify+1
960	f = fy(ify)
961	nsub_new = nsub_old*f
962	#if defined PARTIT_DEBUG
963	IF ( lwp ) THEN
964	WRITE (numout,*) 'Partitioning in y from ',nsub_old &
965	,' to ',nsub_new
966	ENDIF
967	#endif
968	DO isub_old=1,nprocp,nprocp/nsub_old
969	! Check that there are sufficient points to factorise.
970	IF ( pjeub(isub_old)-pjelb(isub_old)+1.LT.f ) THEN
971	WRITE (numout,*) &
972	'partition_rk: insufficient points to partition'
973	ierr = 1
974	RETURN
975	ENDIF
976	! Set the target number of points in the new sub-domains.
977	ntarget = NINT(REAL(pnactive(isub_old))/REAL(f))
978	ilb = pielb(isub_old)
979	iub = pieub(isub_old)
980	jst = pjelb(isub_old)
981	jub = pjeub(isub_old)
982	#if defined PARTIT_DEBUG
983	IF ( lwp ) THEN
984	WRITE (numout,*)
985	WRITE (numout,*) '======================================='
986	WRITE (numout,*)
987	WRITE (numout,*) 'Partitioning sub-domain ',isub_old,' : ' &
988	,ilb,iub,jst,jub
989	WRITE (numout,*) 'Old partition has ',pnactive(isub_old) &
990	,' points'
991	WRITE (numout,*) 'Target is ',ntarget
992	ENDIF
993	#endif
994	! Create the new sub-domains.
995	ngone = 0
996	DO isub=1,f
997	isub_new = isub_old + (isub-1)*nprocp/nsub_new
998	#if defined PARTIT_DEBUG
999	IF ( lwp ) THEN
1000	WRITE (numout,*)
1001	WRITE (numout,*) 'Creating new sub-domain ',isub_new
1002	ENDIF
1003	#endif
1004	! The last new sub-domain takes the remaining data.
1005	IF ( isub.EQ.f ) THEN
1006	jfin = jub
1007	ELSE
1008	! Scan lines in y counting active grid points until we
1009	! exceed the target.
1010	ntry = 0
1011	jfin = -1
1012	scancols: DO j=jst,jub
1013	! Count up the contribution from the next line.
1014	line = 0
1015	DO i=ilb,iub
1016	IF ( mask(i,j) == 1 ) line = line+1
1017	ENDDO
1018	#if defined PARTIT_DEBUG
1019	IF ( lwp ) THEN
1020	!dbg write (numout,*) 'Line ',j,' has ',line
1021	ENDIF
1022	#endif
1023	! Check against the target.
1024	IF ( ntry+line.GT.ntarget ) THEN
1025	! Is it nearer the target to include this line or not ?
1026	IF ( ntry+line-ntarget.GT.ntarget-ntry ) THEN
1027	! Nearer start than end. Finish at previous line.
1028	jfin = j-1
1029	ELSE
1030	! Nearer end than start. Include this line.
1031	jfin = j
1032	ntry = ntry + line
1033	ENDIF
1034	#if defined PARTIT_DEBUG
1035	IF ( lwp ) THEN
1036	WRITE (numout,*) 'Reached target at ',jfin &
1037	,' ntry = ',ntry
1038	ENDIF
1039	#endif
1040	EXIT scancols
1041	ENDIF
1042	! Add in the line count to the running total.
1043	ntry = ntry + line
1044	ENDDO scancols
1045	IF ( jfin.LT.0 ) jfin = jub
1046	ENDIF
1047	! Set up the parameters of the new sub-domain.
1048	pnactive(isub_new) = 0
1049	DO j=jst,jfin
1050	DO i=ilb,iub
1051	IF ( mask(i,j) == 1 ) THEN
1052	pnactive(isub_new) = pnactive(isub_new)+1
1053	ENDIF
1054	ENDDO
1055	ENDDO
1056	pielb(isub_new) = ilb
1057	pieub(isub_new) = iub
1058	pjelb(isub_new) = jst
1059	pjeub(isub_new) = jfin
1060	#if defined PARTIT_DEBUG
1061	IF ( lwp ) THEN
1062	WRITE (numout,*) 'New subdomain is ',ilb,iub,jst,jfin &
1063	,' with ',pnactive(isub_new),' points'
1064	ENDIF
1065	#endif
1066	! Reset the target based on the actual number of points
1067	! remaining, split between the partitions remaining.
1068	ngone = ngone+ntry
1069	#if defined PARTIT_DEBUG
1070	IF(lwp)WRITE (numout,*) 'Target reset to ',ntarget
1071	#endif
1072	! Start the next search at the next point.
1073	jst = jfin+1
1074	ENDDO
1075	ENDDO
1076	ENDIF
1077	! Update the number of sub-domains ready for the next loop.
1078	nsub_old = nsub_new
1079	ENDDO
1080	#if defined PARTIT_DEBUG
1081	IF ( lwp ) THEN
1082	WRITE (numout,*)
1083	WRITE (numout,*) '========================================='
1084	WRITE (numout,*)
1085	ENDIF
1086	#endif
1087	! Set the size of each sub-domain.
1088	DO iproc=1,nprocp
1089	piesub(iproc) = pieub(iproc)-pielb(iproc)+1
1090	pjesub(iproc) = pjeub(iproc)-pjelb(iproc)+1
1091	#if defined PARTIT_DEBUG
1092	IF(lwp)THEN
1093	WRITE (numout,*) 'Domain ',iproc,' has ',pnactive(iproc), ' ocean points'
1094	WRITE(*,"(I3,': [',I3,':',I3,'][',I3,':',I3,']')") &
1095	iproc, pielb(iproc), pieub(iproc), pjelb(iproc), pjeub(iproc)
1096
1097	END IF
1098	#endif
1099	ENDDO
1100
1101	END SUBROUTINE partition_rk_core
1102
1103
1104	SUBROUTINE factor ( ifax, maxfax, nfax, n, ierr )
1105
1106	!!------------------------------------------------------------------
1107	! Subroutine to return the prime factors of n.
1108	! nfax factors are returned in array ifax which is of maximum
1109	! dimension maxfax.
1110	!!------------------------------------------------------------------
1111	IMPLICIT NONE
1112	! Subroutine arguments
1113	INTEGER ierr, n, nfax, maxfax
1114	INTEGER ifax(maxfax)
1115	! Local variables.
1116	INTEGER i, ifac, l, nu
1117	INTEGER lfax(20)
1118	! lfax contains the set of allowed factors.
1119	DATA (lfax(i),i=1,20) / 71,67,59,53,47,43,41,37,31,29 &
1120	,23,19,17,13,11, 7, 5, 3, 2, 1 /
1121	! Clear the error flag.
1122	ierr = 0
1123	! Find the factors of n.
1124	IF ( n.EQ.1 ) THEN
1125	nfax = 0
1126	GOTO 20
1127	ENDIF
1128	! nu holds the unfactorised part of the number.
1129	! nfax holds the number of factors found.
1130	! l points to the allowed factor list.
1131	! ifac holds the current factor.
1132	nu = n
1133	nfax = 0
1134	l = 1
1135	ifac = lfax(l)
1136	! Label 10 is the start of the factor search loop.
1137	10 CONTINUE
1138	! Test whether the factor will divide.
1139	IF ( MOD(nu,ifac).EQ.0 ) THEN
1140	! Add the factor to the list.
1141	nfax = nfax+1
1142	IF ( nfax.GT.maxfax ) THEN
1143	ierr = 6
1144	WRITE (,) &
1145	'FACTOR: insufficient space in factor array ',nfax
1146	RETURN
1147	ENDIF
1148	ifax(nfax) = ifac
1149	! Divide out the factor from the remaining number.
1150	! If unity remains, the number has been
1151	! successfully factorized.
1152	nu = nu/ifac
1153	IF ( nu.EQ.1 ) GO TO 20
1154	! Loop to try the factor again.
1155	GOTO 10
1156	ENDIF
1157	! Move on to the next factor in the list.
1158	l = l+1
1159	ifac = lfax(l)
1160	! Unless the end of the factor list has been reached, loop.
1161	IF ( ifac.GT.1 ) go to 10
1162	! There is a factor in n which is not in the lfac list.
1163	! Add the remaining number onto the end of the list.
1164	nfax = nfax+1
1165	IF ( nfax.GT.maxfax ) THEN
1166	ierr = 6
1167	WRITE (,) 'FACTOR: insufficient space in factor array ',nfax
1168	RETURN
1169	ENDIF
1170	ifax(nfax) = nu
1171	! Label 20 is the exit point from the factor search loop.
1172	20 CONTINUE
1173
1174	END SUBROUTINE factor
1175
1176	!#define TRIM_DEBUG
1177
1178	SUBROUTINE part_trim ( depth, isTrimmed, ierr )
1179	!!------------------------------------------------------------------
1180	!!
1181	!! Examines all the sub-domains and trims off boundary rows and
1182	!! columns which are all land.
1183	!!
1184	!! depth real input Depth map.
1185	!! isTrimmed logical output Whether N,E,S,W edge
1186	!! of domain has been
1187	!! trimmed
1188	!! ierr int output Error flag.
1189	!!
1190	!! Mike Ashworth, CLRC Daresbury Laboratory, July 1999
1191	!!------------------------------------------------------------------
1192	USE par_oce, ONLY: jpreci, jprecj
1193	USE iom, ONLY: jpiglo, jpjglo
1194	IMPLICIT NONE
1195
1196	! Subroutine arguments.
1197	INTEGER, INTENT(in) :: depth(jpiglo,jpjglo)
1198	LOGICAL, DIMENSION(:,:), INTENT(out) :: isTrimmed
1199	INTEGER, INTENT(out) :: ierr
1200	! Local variables.
1201	INTEGER :: i, iproc, j, newbound
1202
1203	! Clear the error code.
1204	ierr = 0
1205
1206	! Clear all flags
1207	! N E S W
1208	! 1 2 3 4
1209	isTrimmed(1:4,1:nprocp) = .FALSE.
1210
1211	! Examine each sub-domain in turn.
1212
1213	subdomain_loop: DO iproc=1,nprocp
1214
1215	! Do not trim if there are no active points at all.
1216	! Otherwise we will trim away the whole sub-domain and we
1217	! will be in big trouble.
1218
1219	! Look at the low-i columns (Western edge of domain)
1220
1221	newbound = pielb(iproc)
1222	lowi: DO i=pielb(iproc),pieub(iproc)
1223
1224	! Scan this column for wet points
1225
1226	DO j=MAX(1,pjelb(iproc)-jprecj),MIN(jpjglo,pjeub(iproc)+jprecj)
1227
1228	IF ( depth(i,j) == 1 ) THEN
1229	newbound = MAX(i - jpreci - nextra, pielb(iproc))
1230	#if defined TRIM_DEBUG
1231	IF ( lwp ) THEN
1232	WRITE (numout,*) 'Sub-domain ',iproc,': Low-i loop: row ',i &
1233	,' is land. New bound is ',newbound
1234	ENDIF
1235	#endif
1236	EXIT lowi
1237	ENDIF
1238	ENDDO
1239	ENDDO lowi
1240
1241	! Update with the new boundary and monitor.
1242
1243	IF ( newbound.NE.pielb(iproc) ) THEN
1244	#if defined TRIM_DEBUG
1245	IF ( lwp ) THEN
1246	IF ( newbound-pielb(iproc).GT.1 ) THEN
1247	WRITE(numout,'(a,i5,3(a,i6))') ' Process ',iproc-1 &
1248	,' trimmed ',newbound-pielb(iproc) &
1249	,' cols: ',pielb(iproc),' to ',newbound-1
1250	ELSE
1251	WRITE(numout,'(a,i5,3(a,i6))') ' Process ',iproc-1 &
1252	,' trimmed ',newbound-pielb(iproc) &
1253	,' col : ',pielb(iproc)
1254	ENDIF
1255	ENDIF
1256	#endif
1257	pielb(iproc) = newbound
1258	isTrimmed(4,iproc) = .TRUE.
1259	ENDIF
1260
1261	! Look at the high-i columns (Eastern edge of domain).
1262
1263	newbound = pieub(iproc)
1264	highi: DO i=pieub(iproc),pielb(iproc),-1
1265
1266	DO j=MAX(1,pjelb(iproc)-jprecj), MIN(jpjglo,pjeub(iproc)+jprecj)
1267
1268	IF ( depth(i,j) == 1 ) THEN
1269	! We've found a wet point in this column so this is as far
1270	! as we can trim.
1271	newbound = MIN(i + jpreci + nextra, pieub(iproc))
1272	#if defined TRIM_DEBUG
1273	IF ( lwp ) THEN
1274	WRITE (numout,"('Sub-domain ',I3,': High-i loop: col ',I3, &
1275	& ' contains water. New bound is ',I3)") iproc,i,newbound
1276	ENDIF
1277	#endif
1278	EXIT highi
1279
1280	ENDIF
1281	ENDDO
1282	ENDDO highi
1283
1284	! Update with the new boundary and monitor.
1285
1286	IF ( newbound.NE.pieub(iproc) ) THEN
1287	#if defined TRIM_DEBUG
1288	IF ( lwp ) THEN
1289	IF ( (pieub(iproc)-newbound) .GT.1 ) THEN
1290	WRITE (numout,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1291	,' trimmed ',pieub(iproc)-newbound &
1292	,' cols: ',newbound+1,' to ',pieub(iproc)
1293	ELSE
1294	WRITE (numout,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1295	,' trimmed ',pieub(iproc)-newbound &
1296	,' col : ',pieub(iproc)
1297	ENDIF
1298	ENDIF
1299	#endif
1300	pieub(iproc) = newbound
1301	isTrimmed(2,iproc) = .TRUE.
1302	ENDIF
1303
1304	! Look at the low-j rows (Southern edge of domain)
1305
1306	newbound = pjelb(iproc)
1307	lowj: DO j=pjelb(iproc),pjeub(iproc),1
1308
1309	! Scan this row for wet points
1310
1311	DO i=MAX(1,pielb(iproc)-jpreci),MIN(jpiglo,pieub(iproc)+jpreci)
1312	IF ( depth(i,j) == 1) THEN
1313	newbound = MAX(j - jpreci - nextra, pjelb(iproc))
1314	#if defined TRIM_DEBUG
1315	IF ( lwp ) THEN
1316	WRITE (numout,*) 'Sub-domain ',iproc,': Low-j loop: column ',j &
1317	,' is land. New bound is ',newbound
1318	ENDIF
1319	#endif
1320	EXIT lowj
1321	ENDIF
1322	ENDDO
1323
1324	ENDDO lowj
1325
1326
1327	! Update with the new boundary and monitor.
1328
1329	IF ( newbound.NE.pjelb(iproc) ) THEN
1330	#if defined TRIM_DEBUG
1331	IF ( lwp ) THEN
1332	IF ( (newbound-pjelb(iproc)) .GT.1 ) THEN
1333	WRITE (numout,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1334	,' trimmed ',newbound-pjelb(iproc) &
1335	,' rows: ',pjelb(iproc),' to ',newbound-1
1336	ELSE
1337	WRITE (numout,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1338	,' trimmed ',newbound-pjelb(iproc) &
1339	,' row : ',pjelb(iproc)
1340	ENDIF
1341	ENDIF
1342	#endif
1343	pjelb(iproc) = newbound
1344	isTrimmed(3,iproc) = .TRUE.
1345	ENDIF
1346
1347	! Look at the high-j rows (Northern edge of domain)
1348
1349	newbound = pjeub(iproc)
1350	highj: DO j=pjeub(iproc),pjelb(iproc),-1
1351
1352	! Scan this row for wet points
1353
1354	DO i=MAX(1,pielb(iproc)-jpreci),MIN(jpiglo,pieub(iproc)+jpreci)
1355	IF ( depth(i,j) == 1 ) THEN
1356	newbound = MIN(j + jpreci + nextra, pjeub(iproc))
1357	#if defined TRIM_DEBUG
1358	IF ( lwp ) then
1359	WRITE (numout,*) 'Sub-domain ',iproc,': High-j loop: column ',j &
1360	,' is land. New bound is ',newbound
1361	ENDIF
1362	#endif
1363	EXIT highj
1364	ENDIF
1365	ENDDO
1366	ENDDO highj
1367
1368	! Update with the new boundary and monitor.
1369
1370	IF ( newbound.NE.pjeub(iproc) ) THEN
1371	#if defined TRIM_DEBUG
1372	IF ( lwp ) THEN
1373	IF ( pjeub(iproc)-newbound.GT.1 ) THEN
1374	WRITE (numout,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1375	,' trimmed ',pjeub(iproc)-newbound &
1376	,' rows: ',newbound+1,' to ',pjeub(iproc)
1377	ELSE
1378	WRITE (numout,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1379	,' trimmed ',pjeub(iproc)-newbound &
1380	,' row : ',pjeub(iproc)
1381	ENDIF
1382	ENDIF
1383	#endif
1384	pjeub(iproc) = newbound
1385	isTrimmed(1,iproc) = .TRUE.
1386	ENDIF
1387
1388	! Reset the size of the sub-domain.
1389
1390	piesub(iproc) = pieub(iproc)-pielb(iproc)+1
1391	pjesub(iproc) = pjeub(iproc)-pjelb(iproc)+1
1392
1393	! endif active_subdomain
1394
1395	ENDDO subdomain_loop
1396
1397	END SUBROUTINE part_trim
1398
1399
1400	SUBROUTINE finish_partition(fromFile)
1401	USE mapcomm_mod, ONLY: ielb,ieub,pielb,pjelb,pieub,pjeub, &
1402	iesub,jesub,jeub,ilbext,iubext,jubext, &
1403	jlbext,pnactive,piesub,pjesub,jelb,pilbext, &
1404	piubext, pjlbext, pjubext, &
1405	trimmed, nidx,eidx,sidx,widx
1406	IMPLICIT NONE
1407	LOGICAL, INTENT(in), OPTIONAL :: fromFile
1408	! Locals
1409	INTEGER :: iproc, ierr
1410	LOGICAL :: lFromFile
1411
1412	! Check to see whether we're dealing with a partion that has been
1413	! read from file. If we are then there are some things we don't
1414	! calculate here.
1415	lFromFile = .FALSE.
1416	IF( PRESENT(fromFile) ) lFromFile = fromFile
1417
1418	IF(.NOT. lFromFile)THEN
1419	! Set the external boundary flags before boundaries are
1420	! altered by the trimming process and it becomes more difficult
1421	! to recognize which were the external boundaries.
1422
1423	DO iproc=1, nprocp, 1
1424	pilbext(iproc) = pielb(iproc).EQ.1
1425	piubext(iproc) = pieub(iproc).EQ.jpiglo
1426	pjlbext(iproc) = pjelb(iproc).EQ.1
1427	pjubext(iproc) = pjeub(iproc).EQ.jpjglo
1428	ENDDO
1429
1430	! Trim off redundant rows and columns containing all land.
1431	IF(.NOT. ALLOCATED(trimmed) )THEN
1432	ALLOCATE(trimmed(4,nprocp), Stat=ierr)
1433	IF(ierr /= 0)THEN
1434	CALL ctl_stop('STOP', &
1435	'Failed to allocate memory for domain trimming')
1436	END IF
1437	END IF
1438
1439	#if defined key_mpp_mpi
1440	IF ( nn_pttrim ) THEN
1441	nextra = 3
1442	CALL part_trim ( imask, trimmed, ierr )
1443	ELSE
1444	! Need non-zero nextra because otherwise hit trouble with fields
1445	! not being read from disk over land regions
1446	nextra = 2
1447	!nextra = 0 ! Don't need to back-off on message trimming
1448	! if we're not trimming the domains
1449	trimmed(1:4,1:nprocp) = .FALSE.
1450	ENDIF
1451	#else
1452	trimmed(1:4,1:nprocp) = .FALSE.
1453	#endif
1454	END IF ! not read from file
1455
1456	! Lower boundary (long.) of sub-domain, GLOBAL coords
1457	! before correction for global halos
1458	nimpp = pielb(narea)
1459
1460	! Is the domain on an external LONGITUDE boundary?
1461	nbondi = 0
1462	ilbext = pilbext(narea)
1463	IF(ilbext)THEN
1464	nbondi = -1
1465	END IF
1466
1467	IF( (.NOT. ilbext) .OR. (ilbext .AND. trimmed(widx,narea)) )THEN
1468	! It isn't, which means we must shift its lower boundary by
1469	! -jpreci to allow for the overlap of this domain with its
1470	! westerly neighbour.
1471	nimpp = nimpp - jpreci
1472	END IF
1473
1474	iubext = piubext(narea)
1475	IF(iubext)THEN
1476	nbondi = 1
1477	IF(ilbext)nbondi = 2 ! Both East and West boundaries are at
1478	! edges of global domain
1479	END IF
1480
1481	! Set local values for limits in global coords of the sub-domain
1482	! owned by this process.
1483	ielb = pielb (narea)
1484	ieub = pieub (narea)
1485	iesub = piesub(narea)
1486
1487	jpi = iesub + 2*jpreci ! jpi is the same for all domains - this is
1488	! what original decomposition did
1489	nlci = jpi
1490
1491	! If the domain is at the edge of the model domain and a cyclic
1492	! East-West b.c. is in effect then it already incorporates one
1493	! extra halo column (because of the way the model domain itself is
1494	! set-up). If we've trimmed-off dry rows and columns then, even if
1495	! a domain is on the model boundary, it may still need a halo so
1496	! we add one.
1497	IF( nbondi == -1 .AND. (.NOT. trimmed(widx,narea)) )THEN
1498	! Western boundary
1499	! First column of global domain is actually a halo but NEMO
1500	! still sets nldi to 1.
1501	nldi = 1 ! Lower bnd of int. region of sub-domain, LOCAL
1502	nlei = nldi + iesub - 1
1503	nlci = nlei + jpreci
1504	jpi = nlci
1505
1506	ELSE IF( nbondi == 1 .AND. (.NOT. trimmed(eidx,narea)) )THEN
1507	! Eastern boundary
1508	nldi = 1 + jpreci
1509	! Last column of global domain is actually a halo
1510	nlei = nldi + iesub - 1
1511	nlci = nlei
1512	jpi = nlei
1513
1514	ELSE IF( nbondi == 2)THEN
1515	! Both boundaries are on the edges of the global domain
1516	IF(trimmed(widx,narea))THEN
1517	nldi = 1 + jpreci
1518	ELSE
1519	nldi = 1
1520	END IF
1521	nlei = nldi + iesub - 1
1522
1523	IF(trimmed(eidx,narea))THEN
1524	nlci = nlei + jpreci
1525	ELSE
1526	nlci = nlei
1527	END IF
1528	jpi = nlci
1529
1530	ELSE
1531	! We have no external boundaries to worry about
1532	nldi = 1 + jpreci
1533	nlei = nldi + iesub - 1 !
1534	END IF
1535
1536	jpim1 = jpi - 1
1537
1538	! Lower ext. boundary (lat.) of sub-domain, global coords
1539	njmpp= pjelb (narea)
1540
1541	! Is the domain on an external LATITUDE boundary?
1542	nbondj = 0
1543	jlbext = pjlbext(narea)
1544	IF(jlbext)THEN
1545	nbondj = -1
1546	ENDIF
1547
1548	IF((.NOT. jlbext) .OR. (jlbext .AND. trimmed(sidx,narea)) )THEN
1549	! It isn't, which means we must shift its lower boundary by
1550	! -jprecj to allow for the overlap of this domain with its
1551	! southerly neighbour.
1552	njmpp = njmpp - jprecj
1553	END IF
1554	! ARPDBG - should we allow for trimming of northern edge of
1555	! sub-domains here?
1556	jubext = pjubext(narea)
1557	IF(jubext)THEN
1558	nbondj = 1
1559	IF(jlbext)nbondj = 2
1560	END IF
1561
1562	jelb = pjelb (narea) ! Lower bound of internal domain
1563	jeub = pjeub (narea) ! Upper bound of internal domain
1564	jesub = pjesub(narea) ! Extent of internal domain
1565
1566	jpj = jesub + 2*jprecj ! jpj is the same for all domains - this is
1567	! what original decomposition did
1568	nlcj = jpj
1569
1570	! Unlike the East-West boundaries, the global domain does not include
1571	! halo (rows) at the Northern and Southern edges. In fact, there is no
1572	! cyclic boundary condition in the North-South direction so there are no
1573	! halos at all at the edges of the global domain.
1574	IF( nbondj == -1 .AND. (.NOT. trimmed(sidx,narea)) )THEN
1575	! Southern edge
1576	nldj = 1 ! Start of internal region (local coords)
1577	nlej = nldj + jesub - 1 ! Upper bound of int. region of sub-domain, local
1578	nlcj = nlej + jprecj
1579	jpj = nlcj
1580	ELSE IF( nbondj == 1 .AND. (.NOT. trimmed(nidx,narea)) )THEN
1581	! Northern edge
1582	nldj = 1+jprecj ! Start of internal region (local coords)
1583	nlej = nldj + jesub - 1
1584	nlcj = nlej
1585	jpj = nlcj
1586	jpj = jpj + 1 ! Add one extra row for zero BC along N edge as
1587	! happens in orig. code when MOD(jpjglo,2)!=0
1588	! Many loops go up to j=jpjm1 so unless jpj>nlej
1589	! they won't cover the whole domain.
1590	ELSE IF( nbondj == 2)THEN
1591	! Both boundaries are on the edges of the global domain
1592
1593	IF( trimmed(sidx,narea) )THEN
1594	nldj = 1+jprecj
1595	ELSE
1596	nldj = 1
1597	END IF
1598	nlej = nldj + jesub - 1
1599
1600	IF( trimmed(nidx,narea) )THEN
1601	nlcj = nlej + jprecj
1602	jpj = nlcj
1603	ELSE
1604	nlcj = nlej
1605	jpj = nlcj
1606	END IF
1607	jpj = jpj + 1 ! Add one extra row for zero BC along N edge as
1608	! happens in orig. code when MOD(jpjglo,2)!=0
1609	ELSE
1610	! We have no external boundaries to worry about
1611	nldj = 1+jprecj ! Lower bound of int. region of sub-domain, local
1612	nlej = nldj + jesub - 1
1613	END IF
1614
1615	jpjm1 = jpj-1
1616	jpij = jpi*jpj
1617
1618
1619	END SUBROUTINE finish_partition
1620
1621
1622	SUBROUTINE mpp_ini_north
1623	!!----------------------------------------------------------------------
1624	!! * routine mpp_ini_north *
1625	!!
1626	!! ** Purpose : Initialize special communicator for north folding
1627	!! condition together with global variables needed in the mpp folding
1628	!!
1629	!! ** Method : - Look for northern processors
1630	!! - Put their number in nrank_north
1631	!! - Create groups for the world processors and the north
1632	!! processors
1633	!! - Create a communicator for northern processors
1634	!!
1635	!! ** output
1636	!! njmppmax = njmpp for northern procs
1637	!! ndim_rank_north = number of processors in the northern line
1638	!! nrank_north (ndim_rank_north) = number of the northern procs.
1639	!! ngrp_world = group ID for the world processors
1640	!! ngrp_north = group ID for the northern processors
1641	!! ncomm_north = communicator for the northern procs.
1642	!! north_root = number (in the world) of proc 0 in the northern comm.
1643	!! nwidthmax = width of widest northern domain
1644	!!
1645	!! History :
1646	!! ! 03-09 (J.M. Molines, MPI only )
1647	!! ! 08-09 (A.R. Porter - for new decomposition)
1648	!!----------------------------------------------------------------------
1649	USE par_oce, ONLY: jperio, jpni
1650	USE exchmod, ONLY: nrank_north, north_root, ndim_rank_north, &
1651	ncomm_north, ngrp_world, ngrp_north, &
1652	do_nfold, num_nfold_rows, nfold_npts
1653	USE dom_oce, ONLY: narea
1654	IMPLICIT none
1655	#ifdef key_mpp_shmem
1656	CALL ctl_stop('STOP', ' mpp_ini_north not available in SHMEM' )
1657	# elif key_mpp_mpi
1658	INTEGER :: ierr
1659	INTEGER :: jproc
1660	INTEGER :: ii,ji
1661	!!----------------------------------------------------------------------
1662
1663	! Look for how many procs on the northern boundary
1664	!
1665	ndim_rank_north = 0
1666	nwidthmax = 0
1667	do_nfold = .FALSE.
1668
1669	IF (.NOT. (jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1) ) THEN
1670	! No northern boundary to worry about
1671	RETURN
1672	END IF
1673
1674	DO jproc=1,mppsize,1
1675	IF ( pjubext(jproc) ) THEN
1676
1677	! If trimming of dry land from sub-domains is enabled
1678	! then check that this PE does actually have data to
1679	! contribute to the N-fold. If trimming is not enabled
1680	! then this condition will always be true for northern
1681	! PEs.
1682	IF( pjeub(jproc) > (jpjglo - num_nfold_rows) )THEN
1683
1684	ndim_rank_north = ndim_rank_north + 1
1685
1686	! and for the width of the widest northern domain...
1687	nwidthmax = MAX(nwidthmax, piesub(jproc))
1688	ENDIF
1689
1690	END IF
1691	END DO
1692	nwidthmax = nwidthmax + 2*jpreci ! Allow for halos
1693
1694	! Allocate the right size to nrank_north
1695	!
1696	ALLOCATE(nrank_north(ndim_rank_north), nfold_npts(ndim_rank_north), &
1697	Stat=ierr)
1698	IF( ierr /= 0 )THEN
1699	CALL ctl_stop('STOP','mpp_ini_north: failed to allocate arrays')
1700	END IF
1701
1702	#if defined PARTIT_DEBUG
1703	IF(lwp)THEN
1704	WRITE(,) 'mpp_ini_north: no. of northern PEs = ',ndim_rank_north
1705	WRITE(,) 'mpp_ini_north: nwidthmax = ',nwidthmax
1706	END IF
1707	#endif
1708	! Fill the nrank_north array with proc. number of northern procs.
1709	! Note : ranks start at 0 in MPI
1710	!
1711	ii=0
1712	DO ji = 1, mppsize, 1
1713	IF ( pjubext(ji) .AND. &
1714	(pjeub(ji) > (jpjglo - num_nfold_rows)) ) THEN
1715	ii=ii+1
1716	nrank_north(ii)=ji-1
1717
1718	! Flag that this PE does do North-fold (with trimming, checking
1719	! npolj is no longer sufficient)
1720	IF(ji == narea) do_nfold = .TRUE.
1721
1722	#if defined NO_NFOLD_GATHER
1723	! How many data points will this PE have to send for N-fold?
1724
1725	! No. of valid rows for n-fold = num_nfold_rows - <no. trimmed rows>
1726	! = num_nfold_rows - jpjglo + pjeub(ji)
1727
1728	! ARPDBG - could trim land-only rows/cols from this...
1729	nfold_npts(ii) = MAX(num_nfold_rows - jpjglo + pjeub(ji), 0) * &
1730	( nleit(ji) - nldit(ji) + 1 )
1731	#endif
1732	END IF
1733	END DO
1734	! create the world group
1735	!
1736	CALL MPI_COMM_GROUP(mpi_comm_opa,ngrp_world,ierr)
1737	!
1738	! Create the North group from the world group
1739	CALL MPI_GROUP_INCL(ngrp_world,ndim_rank_north,nrank_north, &
1740	ngrp_north,ierr)
1741
1742	! Create the North communicator , ie the pool of procs in the north group
1743	!
1744	CALL MPI_COMM_CREATE(mpi_comm_opa,ngrp_north,ncomm_north,ierr)
1745
1746
1747	! find proc number in the world of proc 0 in the north
1748	CALL MPI_GROUP_TRANSLATE_RANKS(ngrp_north,1,0,ngrp_world,north_root,ierr)
1749
1750	#endif
1751
1752	END SUBROUTINE mpp_ini_north
1753
1754
1755	SUBROUTINE eval_partition( nx, ny, mask, score )
1756
1757	! Compute the cost function for the current partition
1758	!
1759	! Assume that the time taken for a run is proportional
1760	! to the maximum over processors of:
1761	! w_processing * cost_processing
1762	! + w_communications * cost_communications
1763	! Assume further that cost_processing goes as
1764	! (number of wet points) + f_proc * (number of dry points)
1765	! (with f_proc << 1)
1766	! and that cost_communications goes as
1767	! (cost of intra-node communications) +
1768	! f_comm * (cost of inter-node communications)
1769	! (with f_comm << 1)
1770	!
1771	! However, because of the possiblity of network contention,
1772	! other factors may also matter, especially:
1773	! total over sub-domains of halo points with off-node neighbours
1774	! max over nodes of total off-node halo points and message counts
1775	!
1776	! With this in mind, we construct the ansatz
1777	! maximum over processors of {
1778	! w_1 * (number of wet points)
1779	! + w_2 * (number of dry points)
1780	! + w_3 * (halo points with off-node neighbours)
1781	! + w_4 * (halo points with on-node neighbours)
1782	! + ...
1783	! }
1784	#if defined key_mpp_mpi
1785	USE lib_mpp, ONLY: mppsize
1786	#endif
1787	USE mapcomm_mod, ONLY: iprocmap, land
1788	IMPLICIT NONE
1789	! Arguments
1790	INTEGER, INTENT(in) :: nx, ny
1791	INTEGER, INTENT(in) :: mask(nx,ny)
1792	REAL(wp), DIMENSION(pv_num_scores), INTENT(out) :: score
1793	! Local variables
1794	INTEGER :: iproc, inode, i, j
1795
1796	REAL(wp) :: score_wet, score_dry
1797	REAL(wp) :: score_pli, score_plx
1798	REAL(wp) :: score_pci, score_pcx
1799	REAL(wp) :: score_nli, score_nlx
1800	REAL(wp) :: score_nci, score_ncx
1801	REAL(wp) :: score_tli, score_tlx
1802	REAL(wp) :: score_tci, score_tcx
1803
1804	REAL(wp) :: score_too_narrow
1805
1806	REAL(wp) :: proc_overall_score
1807	REAL(wp) :: proc_comm_score
1808	REAL(wp) :: node_comm_score
1809
1810	REAL(wp), PARAMETER :: weight_wet = 1.00D0
1811	REAL(wp), PARAMETER :: weight_dry = 0.9D0
1812
1813	REAL(wp), PARAMETER :: weight_pli = 0.01D0
1814	REAL(wp), PARAMETER :: weight_plx = 0.20D0
1815	REAL(wp), PARAMETER :: weight_pci = 0.50D0
1816	REAL(wp), PARAMETER :: weight_pcx = 10.00D0
1817
1818	REAL(wp), PARAMETER :: weight_nli = 0.00D0
1819	REAL(wp), PARAMETER :: weight_nlx = 0.20D0
1820	REAL(wp), PARAMETER :: weight_nci = 0.00D0
1821	REAL(wp), PARAMETER :: weight_ncx = 10.00D0
1822
1823	REAL(wp), PARAMETER :: weight_tli = 0.00D0
1824	REAL(wp), PARAMETER :: weight_tlx = 0.01D0
1825	REAL(wp), PARAMETER :: weight_tci = 0.00D0
1826	REAL(wp), PARAMETER :: weight_tcx = 0.50D0
1827
1828	INTEGER :: peer, last_peer
1829
1830	! Which node is each process on?
1831	! Assumes that first nn_cpnode ranks are assigned to node 0,
1832	! next nn_cpnode ranks are assigned to node 1, etc
1833	INTEGER, ALLOCATABLE :: node(:)
1834
1835	#if defined key_mpp_mpi
1836
1837	ALLOCATE(node(nprocp))
1838
1839	DO iproc=1, nprocp
1840	node(iproc) = (iproc-1)/nn_cpnode
1841	END DO
1842
1843	! Calculate maximum per processor score
1844
1845	score(:) = 0.0_wp
1846
1847	! Total (over all processors) off node comms
1848	score_tli = 0.0_wp
1849	score_tlx = 0.0_wp
1850	score_tci = 0.0_wp
1851	score_tcx = 0.0_wp
1852
1853	! Set this to pv_awful if any sub-domain is too narrow.
1854	score_too_narrow = 0.0_wp
1855
1856	! loop over nodes in job, counting from 0
1857	node_loop: DO inode=0, (nprocp-1)/nn_cpnode
1858
1859	score_nli = 0.0_wp
1860	score_nlx = 0.0_wp
1861	score_nci = 0.0_wp
1862	score_ncx = 0.0_wp
1863
1864	! loop over processes in the node
1865	proc_loop: DO iproc=1+inode*nn_cpnode, &
1866	MIN(nprocp,(inode+1)*nn_cpnode)
1867
1868	score_wet = DBLE(pnactive(iproc))
1869	score_dry = DBLE(piesub(iproc)*pjesub(iproc)-score_wet)
1870
1871	score_pli = 0.0_wp
1872	score_plx = 0.0_wp
1873	score_pci = 0.0_wp
1874	score_pcx = 0.0_wp
1875
1876	! Things sometimes go wrong when a sub-domain has very
1877	! narrow partitions (2 grid points or less).
1878	! Prevent such problematic partitions from being selected.
1879	IF ( ((pieub(iproc)-pielb(iproc)) .LE. 2) &
1880	.OR. ((pjeub(iproc)-pjelb(iproc)) .LE. 2) ) THEN
1881	score_too_narrow = pv_awful
1882	END IF
1883
1884	IF (.NOT. pjlbext(iproc)) THEN
1885	j=pjelb(iproc)
1886	IF (j .GT. 1) THEN
1887	last_peer = -1
1888	DO i=pielb(iproc),pieub(iproc)
1889	IF ( (mask(i,j) .NE. land) &
1890	.AND. (mask(i,j-1) .NE. land)) THEN
1891	peer=iprocmap(i,j-1)
1892	! Total points involved in messages
1893	IF (node(peer) .EQ. inode) THEN
1894	score_pli = score_pli+1.0_wp
1895	ELSE
1896	score_plx = score_plx+1.0_wp
1897	END IF
1898	! Total number of messages
1899	IF (peer .NE. last_peer) THEN
1900	last_peer = peer
1901	IF (node(peer) .EQ. inode) THEN
1902	score_pci = score_pci+1.0_wp
1903	ELSE
1904	score_pcx = score_pcx+1.0_wp
1905	END IF
1906	END IF
1907	END IF
1908	END DO
1909	END IF
1910	END IF
1911
1912	IF (.NOT. pjubext(iproc)) THEN
1913	j=pjeub(iproc)
1914	IF (j .LT. ny) THEN
1915	last_peer = -1
1916	DO i=pielb(iproc),pieub(iproc)
1917	IF ( (mask(i,j) .NE. land) &
1918	.AND. (mask(i,j+1) .NE. land)) THEN
1919	peer=iprocmap(i,j+1)
1920	! Total points involved in messages
1921	IF (node(peer) .EQ. inode) THEN
1922	score_pli = score_pli+1.0_wp
1923	ELSE
1924	score_plx = score_plx+1.0_wp
1925	END IF
1926	! Total number of messages
1927	IF (peer .NE. last_peer) THEN
1928	last_peer = peer
1929	IF (node(peer) .EQ. inode) THEN
1930	score_pci = score_pci+1.0_wp
1931	ELSE
1932	score_pcx = score_pcx+1.0_wp
1933	END IF
1934	END IF
1935	END IF
1936	END DO
1937	END IF
1938	END IF
1939
1940	IF (.NOT. pilbext(iproc)) THEN
1941	i=pielb(iproc)
1942	IF (i .GT. 1) THEN
1943	last_peer = -1
1944	DO j=pjelb(iproc),pjeub(iproc)
1945	IF ( (mask(i,j) .NE. land) &
1946	.AND. (mask(i-1,j) .NE. land)) THEN
1947	peer=iprocmap(i-1,j)
1948	! Total points involved in messages
1949	IF (node(peer) .EQ. inode) THEN
1950	score_pli = score_pli+1.0_wp
1951	ELSE
1952	score_plx = score_plx+1.0_wp
1953	END IF
1954	! Total number of messages
1955	IF (peer .NE. last_peer) THEN
1956	last_peer = peer
1957	IF (node(peer) .EQ. inode) THEN
1958	score_pci = score_pci+1.0_wp
1959	ELSE
1960	score_pcx = score_pcx+1.0_wp
1961	END IF
1962	END IF
1963	END IF
1964	END DO
1965	END IF
1966	END IF
1967
1968	IF (.NOT. piubext(iproc)) THEN
1969	i=pieub(iproc)
1970	IF (i .LT. nx) THEN
1971	last_peer = -1
1972	DO j=pjelb(iproc),pjeub(iproc)
1973	IF ( (mask(i,j) .NE. land) &
1974	.AND. (mask(i+1,j) .NE. land)) THEN
1975	peer=iprocmap(i+1,j)
1976	! Total points involved in messages
1977	IF (node(peer) .EQ. inode) THEN
1978	score_pli = score_pli+1.0_wp
1979	ELSE
1980	score_plx = score_plx+1.0_wp
1981	END IF
1982	! Total number of messages
1983	IF (peer .NE. last_peer) THEN
1984	last_peer = peer
1985	IF (node(peer) .EQ. inode) THEN
1986	score_pci = score_pci+1.0_wp
1987	ELSE
1988	score_pcx = score_pcx+1.0_wp
1989	END IF
1990	END IF
1991	END IF
1992	END DO
1993	END IF
1994	END IF
1995
1996	score_nli = score_nli + score_pli
1997	score_nlx = score_nlx + score_plx
1998	score_nci = score_nci + score_pci
1999	score_ncx = score_ncx + score_pcx
2000
2001	proc_overall_score = weight_wet*score_wet &
2002	+ weight_dry*score_dry &
2003	+ weight_pli*score_pli &
2004	+ weight_plx*score_plx &
2005	+ weight_pci*score_pci &
2006	+ weight_pcx*score_pcx
2007
2008	proc_comm_score = weight_pli*score_pli &
2009	+ weight_plx*score_plx &
2010	+ weight_pci*score_pci &
2011	+ weight_pcx*score_pcx
2012
2013	IF (score(pv_index_overall) < proc_overall_score) &
2014	score(pv_index_overall) = proc_overall_score
2015	IF (score(pv_index_pcomms) < proc_comm_score) &
2016	score(pv_index_pcomms) = proc_comm_score
2017	IF (score(pv_index_wet) < score_wet) &
2018	score(pv_index_wet) = score_wet
2019	IF (score(pv_index_dry) < score_dry) &
2020	score(pv_index_dry) = score_dry
2021	IF (score(pv_index_pli) < score_pli) &
2022	score(pv_index_pli) = score_pli
2023	IF (score(pv_index_plx) < score_plx) &
2024	score(pv_index_plx) = score_plx
2025	IF (score(pv_index_pci) < score_pci) &
2026	score(pv_index_pci) = score_pci
2027	IF (score(pv_index_pcx) < score_pcx) &
2028	score(pv_index_pcx) = score_pcx
2029
2030	END DO proc_loop
2031
2032	score_tli = score_tli + score_nli
2033	score_tlx = score_tlx + score_nlx
2034	score_tci = score_tci + score_nci
2035	score_tcx = score_tcx + score_ncx
2036
2037	node_comm_score = weight_nli*score_nli &
2038	+ weight_nlx*score_nlx &
2039	+ weight_nci*score_nci &
2040	+ weight_ncx*score_ncx
2041
2042	IF (score(pv_index_ncomms) < node_comm_score) &
2043	score(pv_index_ncomms) = node_comm_score
2044	IF (score(pv_index_nli) < score_nli) &
2045	score(pv_index_nli) = score_nli
2046	IF (score(pv_index_nlx) < score_nlx) &
2047	score(pv_index_nlx) = score_nlx
2048	IF (score(pv_index_nci) < score_nci) &
2049	score(pv_index_nci) = score_nci
2050	IF (score(pv_index_ncx) < score_ncx) &
2051	score(pv_index_ncx) = score_ncx
2052
2053	END DO node_loop
2054
2055	score(pv_index_tcomms) = weight_tli*score_tli &
2056	+ weight_tlx*score_tlx &
2057	+ weight_tci*score_tci &
2058	+ weight_tcx*score_tcx
2059
2060	score(pv_index_tli) = score_tli
2061	score(pv_index_tlx) = score_tlx
2062	score(pv_index_tci) = score_tci
2063	score(pv_index_tcx) = score_tcx
2064
2065	score(pv_index_overall) = score(pv_index_overall) &
2066	+ score(pv_index_ncomms) &
2067	+ score(pv_index_tcomms) &
2068	+ score_too_narrow
2069
2070	DEALLOCATE(node)
2071
2072	#endif
2073
2074	END SUBROUTINE eval_partition
2075
2076
2077	SUBROUTINE calc_perms( nfactors, factors, &
2078	ndf, df, multiplicity, &
2079	nperms )
2080	IMPLICIT NONE
2081
2082	! Subroutine arguments
2083	!
2084	! Number of factors (including repetitions)
2085	INTEGER, INTENT(in) :: nfactors
2086
2087	! Factors (in descending order)
2088	INTEGER, INTENT(in) :: factors(nfactors)
2089
2090	! Number of distinct factors
2091	INTEGER, INTENT(out) :: ndf
2092
2093	! Distinct factors (in ascending order)
2094	INTEGER :: df(nfactors)
2095
2096	! Multiplicity of each distinct factor
2097	INTEGER :: multiplicity(nfactors)
2098
2099	! Number of distinct permutations
2100	INTEGER, INTENT(out) :: nperms
2101
2102	! Local variables
2103
2104	INTEGER :: product
2105	INTEGER :: i, j
2106
2107	product = factors(nfactors)
2108	ndf = 1
2109	df(:) = 0
2110	df(1) = factors(nfactors)
2111	multiplicity(:) = 0
2112	multiplicity(1) = 1
2113
2114	DO i=nfactors-1,1,-1
2115	product = product*factors(i)
2116	IF (factors(i) .EQ. df(ndf)) THEN
2117	multiplicity(ndf) = multiplicity(ndf)+1
2118	ELSE
2119	ndf = ndf+1
2120	df(ndf) = factors(i)
2121	multiplicity(ndf) = 1
2122	END IF
2123	END DO
2124	! write (,) 'number: ', product
2125
2126	! A careful code would try to avoid overflow here
2127	nperms = 1
2128	DO i=1, nfactors
2129	nperms = nperms*i
2130	END DO
2131	DO i=1, ndf
2132	DO j=1, multiplicity(i)
2133	nperms = nperms/j
2134	END DO
2135	END DO
2136
2137	! At this point, nperms is the number of distinct permutations
2138	! of the factors provided. But each of these permutations can
2139	! be split between x and y in (nfactors+1) ways, e.g.
2140	! x=(2,2,3), y=()
2141	! x=(2,3), y=(2)
2142	! x=(3), y=(2,2)
2143	! x=(), y=(2,2,3)
2144
2145	nperms = nperms*(nfactors+1)
2146	IF (lwp) THEN
2147	WRITE (numout,*) 'distinct factorisations: ', nperms
2148	END IF
2149
2150	END SUBROUTINE calc_perms
2151
2152
2153
2154	SUBROUTINE get_perm_factors( iperm, nf, ndf, df, m, &
2155	fx, nfx, fy, nfy, &
2156	prodx, prody, printit )
2157	USE dom_oce, ONLY: narea
2158	IMPLICIT NONE
2159	!!------------------------------------------------------------------
2160	! Our goal is to visit a particular permutation,
2161	! number perm ( 0 <= perm <= nperms-1 )
2162	! of nf things, of ndf distinct values (actually the prime
2163	! factors of number of MPI tasks), each of which can be repeated
2164	! with multiplicity m_i
2165	! assert nf = sum_{i=1..ndf} m(i)
2166	!
2167	! We don't care about lexicographic ordering, but we do
2168	! need to ensure that we don't visit any permutation twice,
2169	! in a loop over 0..nperms-1.
2170	! Textbook methods typically assume that all the things being
2171	! permuted are distinct.
2172	!
2173	! We use what I call a nested variable radix method.
2174	!
2175	! Stephen Pickles, STFC
2176	! Taken from POLCOMS code by Andrew Porter.
2177	!!------------------------------------------------------------------
2178	! Arguments
2179	INTEGER, INTENT(in) :: iperm, nf, ndf
2180	INTEGER, INTENT(in), DIMENSION(ndf) :: df, m
2181	INTEGER, INTENT(out), DIMENSION(nf) :: fx, fy
2182	INTEGER, INTENT(out) :: nfx, nfy
2183	INTEGER, INTENT(out) :: prodx, prody
2184	LOGICAL, INTENT(in) :: printit
2185	!
2186	! Local variables
2187	!
2188	INTEGER :: perm, split
2189	INTEGER, DIMENSION(nf) :: bin, a
2190	INTEGER, DIMENSION(ndf) :: ways
2191	INTEGER, DIMENSION(0:ndf) :: root, representation
2192	INTEGER :: i, j, k, v, p, q, r
2193	INTEGER :: unfilled, pm, rm
2194	INTEGER :: myinst
2195	LOGICAL, PARAMETER :: debug=.FALSE.
2196	!!------------------------------------------------------------------
2197
2198	! MPI rank of this process
2199	myinst = narea - 1
2200
2201	perm = iperm / (nf+1)
2202	! Where to split between x and y
2203	split = MOD(iperm, (nf+1))
2204
2205	! interpret ways(i) is the number of ways of distributing
2206	! m(i) copies of the i'th distinct factor into the remaining
2207	! bins
2208	k = nf
2209	DO i=1,ndf
2210	ways(i) = k
2211	DO j=2,m(i)
2212	ways(i) = ways(i)*(k-j+1)/j
2213	END DO
2214	k = k-m(i)
2215	END DO
2216
2217	! compute (outer) radices
2218	! root is the variable radix basis corresponding to ways
2219	! root(ndf) is always 1
2220	root(ndf) = 1
2221	root(0:ndf-1) = ways(1:ndf)
2222	DO i=ndf-1,0,-1
2223	root(i)=root(i)*root(i+1)
2224	END DO
2225
2226	bin(:) = 0
2227	unfilled = nf
2228
2229	r = perm
2230	! Next line is valid as long as perm < nperms
2231	representation(0) = 0
2232	DO i=1, ndf
2233	p = r/root(i)
2234	r = r - p*root(i)
2235	representation(i) = p
2236
2237	! At this point, we are considering distinct ways to
2238	! distribute m(i) copies of the i'th distinct factor into
2239	! the unfilled remaining bins. We want to select the p'th one.
2240
2241	DO j=1,unfilled
2242	a(j) = j
2243	END DO
2244	q = 0
2245	find_p: DO
2246	IF (q .GE. p) EXIT find_p
2247
2248	k=m(i)
2249	k_loop: DO
2250	IF (a(k) .EQ. (unfilled - m(i) + k)) THEN
2251	k=k-1
2252	ELSE
2253	EXIT k_loop
2254	END IF
2255	END DO k_loop
2256	a(k) = a(k) + 1
2257	DO v=k+1,m(i)
2258	a(v) = a(k) + v - k
2259	END DO
2260	q = q+1
2261	END DO find_p
2262
2263	IF (debug) THEN
2264	WRITE (*,'(A3)',advance='no') 'a=('
2265	DO j=1, m(i)-1
2266	WRITE (*,'(I3,A1)',advance='no') a(j), ','
2267	END DO
2268	WRITE (*,'(I3,A1)') a(m(i)), ')'
2269	END IF
2270
2271	DO j=1, m(i)
2272	pm=a(j)-j+1
2273
2274	! put this factor in the pm'th empty bin
2275	v = 1
2276	find_bin: DO k=1, nf
2277	IF (bin(k) .EQ. 0) THEN
2278	IF (v .EQ. pm) THEN
2279	bin(k) = df(i)
2280	unfilled = unfilled-1
2281	EXIT find_bin
2282	ELSE
2283	v=v+1
2284	END IF
2285	END IF
2286	END DO find_bin
2287
2288	END DO
2289
2290	END DO
2291
2292	! We have identified the (perm)th distinct permutation,
2293	! but we still need to split the factors between x and y.
2294	fx(:) = 0
2295	prodx = 1
2296	DO i=1,split
2297	fx(i)=bin(i)
2298	prodx=prodx*fx(i)
2299	END DO
2300	nfx=split
2301
2302	fy(:) = 0
2303	prody = 1
2304	j=1
2305	DO i=split+1,nf
2306	fy(j)=bin(i)
2307	prody=prody*fy(j)
2308	j=j+1
2309	END DO
2310	nfy=nf-nfx
2311
2312	IF (printit) THEN
2313
2314	WRITE (6,'(A14,I6,A1,I6,A2)',advance='no') &
2315	'factorisation[', myinst, ']', iperm, ' ('
2316	DO k=1,ndf-1
2317	WRITE (6,'(I4,A1)',advance='no') representation(k), ','
2318	END DO
2319	WRITE (6,'(I4,A1)',advance='no') representation(ndf), ')'
2320
2321	CALL print_factors(6,nfx,fx,nfy,fy)
2322
2323	END IF
2324
2325	END SUBROUTINE get_perm_factors
2326
2327
2328	SUBROUTINE print_factors(lu,nfx,fx,nfy,fy)
2329	IMPLICIT NONE
2330	INTEGER, INTENT(in) :: lu
2331	INTEGER, INTENT(in) :: nfx, nfy
2332	INTEGER, INTENT(in) :: fx(nfx), fy(nfy)
2333	INTEGER :: k
2334
2335	IF (nfx .GT. 0) THEN
2336	WRITE (lu,'(A1)',advance='no') ' '
2337	DO k=1,nfx-1
2338	WRITE (lu,'(I4,A1)',advance='no') fx(k), ','
2339	END DO
2340	WRITE (lu,'(I4)',advance='no') fx(nfx)
2341	END IF
2342	WRITE (lu,'(A1)',advance='no') ':'
2343	IF (nfy .GT. 0) THEN
2344	DO k=1,nfy-1
2345	WRITE (lu,'(I4,A1)',advance='no') fy(k), ','
2346	END DO
2347	WRITE (lu,'(I4)',advance='no') fy(nfy)
2348	END IF
2349	WRITE (lu,'(A1)') ' '
2350
2351	END SUBROUTINE print_factors
2352
2353
2354	CHARACTER(len=16) FUNCTION num_to_string(number)
2355	INTEGER, INTENT(in) :: number
2356
2357	CHARACTER*16 :: outs
2358
2359	WRITE (outs,'(I15)') number
2360	num_to_string = ADJUSTL(outs)
2361
2362	END FUNCTION num_to_string
2363
2364
2365	SUBROUTINE factor_string(fstr,nfx,fx,nfy,fy)
2366	IMPLICIT NONE
2367	CHARACTER*256, INTENT(out) :: fstr
2368	INTEGER, INTENT(in) :: nfx, nfy
2369	INTEGER, INTENT(in) :: fx(nfx), fy(nfy)
2370	!!----------------------------------------------------------------------
2371	!!----------------------------------------------------------------------
2372	INTEGER :: k
2373
2374	fstr = ' '
2375	IF (nfx .GT. 0) THEN
2376	DO k=1,nfx-1
2377	fstr = TRIM(fstr)//TRIM(num_to_string(fx(k)))//'x'
2378	END DO
2379	fstr = TRIM(fstr)//TRIM(num_to_string(fx(nfx)))
2380	END IF
2381	fstr = TRIM(fstr)//'-'
2382	IF (nfy .GT. 0) THEN
2383	DO k=1,nfy-1
2384	fstr = TRIM(fstr)//TRIM(num_to_string(fy(k)))//'x'
2385	END DO
2386	fstr = TRIM(fstr)//TRIM(num_to_string(fy(nfy)))
2387	END IF
2388	END SUBROUTINE factor_string
2389
2390
2391	SUBROUTINE write_partition_map(depth)
2392	IMPLICIT NONE
2393	INTEGER, DIMENSION(:,:), INTENT(in) :: depth
2394	!!----------------------------------------------------------------------
2395	!! Writes an ASCII and PPM format map of the domain decomposition
2396	!! to separate files.
2397	!!----------------------------------------------------------------------
2398	! Locals
2399	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: map
2400	INTEGER :: nx, ny
2401	INTEGER :: iproc, i, j, icol
2402	INTEGER :: lumapout
2403	CHARACTER(LEN=6) :: mode
2404	CHARACTER(LEN=40) :: mapout
2405	INTEGER, PARAMETER :: ncol=15
2406	INTEGER, DIMENSION(3,-2:ncol-1) :: rgbcol
2407	!!----------------------------------------------------------------------
2408
2409	nx = jpiglo
2410	ny = jpjglo
2411
2412	! Generate an integer map of the partitioning.
2413
2414	ALLOCATE (map(nx,ny))
2415	map = -1
2416	DO iproc=1,nprocp
2417	DO j=pjelb(iproc),pjeub(iproc)
2418	DO i=pielb(iproc),pieub(iproc)
2419	IF ( depth(i,j) == 0 ) THEN
2420
2421	! Identify the land using -2 which is set to black
2422	! in the colour map below.
2423	map(i,j) = -2
2424	ELSE
2425
2426	! Identify to which process the point is assigned.
2427	map(i,j) = iproc-1
2428	ENDIF
2429	ENDDO
2430	ENDDO
2431	ENDDO
2432
2433	! Write the map to a file for plotting.
2434
2435	IF ( lwp ) THEN
2436
2437	! ASCII format map file.
2438
2439	lumapout = 9
2440	mode = 'simple'
2441	IF ( nprocp.LT.10 ) THEN
2442	WRITE (mapout,'(''Map_'',a6,''_'',i1,''.dat'')') mode,nprocp
2443	ELSEIF ( nprocp.LT.100 ) THEN
2444	WRITE (mapout,'(''Map_'',a6,''_'',i2,''.dat'')') mode,nprocp
2445	ELSEIF ( nprocp.LT.1000 ) THEN
2446	WRITE (mapout,'(''Map_'',a6,''_'',i3,''.dat'')') mode,nprocp
2447	ELSE
2448	WRITE (mapout,'(''Map_'',a6,''_'',i4,''.dat'')') mode,nprocp
2449	ENDIF
2450	OPEN (lumapout,file=mapout)
2451	WRITE (lumapout,*) nx,ny
2452	DO j=1,ny
2453	! write (lumapout,'(15i5)') (map(i,j),i=1,ny)
2454	DO i=1,nx,1
2455	WRITE (lumapout,'(3i5)') i ,j, map(i,j)
2456	END DO
2457	ENDDO
2458	CLOSE (lumapout)
2459
2460	! PPM format map file.
2461
2462	lumapout = 10
2463	mode = 'partrk'
2464
2465	IF ( nprocp.LT.10 ) THEN
2466	WRITE (mapout,'(''Map_'',a6,''_'',i1,''.ppm'')') mode,nprocp
2467	ELSEIF ( nprocp.LT.100 ) THEN
2468	WRITE (mapout,'(''Map_'',a6,''_'',i2,''.ppm'')') mode,nprocp
2469	ELSEIF ( nprocp.LT.1000 ) THEN
2470	WRITE (mapout,'(''Map_'',a6,''_'',i3,''.ppm'')') mode,nprocp
2471	ELSE
2472	WRITE (mapout,'(''Map_'',a6,''_'',i4,''.ppm'')') mode,nprocp
2473	ENDIF
2474	OPEN (lumapout,file=mapout)
2475
2476	! PPM magic number.
2477	! Comment line
2478	! width and height of image (same as that of the domain)
2479	! Maximum colour value.
2480
2481	WRITE (lumapout,'(a)') 'P3'
2482	WRITE (lumapout,'(a)') '# '//mapout
2483	WRITE (lumapout,'(2i6)') nx,ny
2484	WRITE (lumapout,'(i6)') 255
2485
2486	! Define RGB colours. 0 is grey for the land. 1-16 for the sub-domains.
2487	! When there are more than 16 sub-domains the colours wrap around.
2488
2489	rgbcol(:,-2) = (/ 0, 0, 0 /)
2490	rgbcol(:,-1) = (/ 170, 170, 170 /)
2491	rgbcol(:, 0) = (/ 0, 0, 255 /) ! dark blue
2492	rgbcol(:, 1) = (/ 0, 85, 255 /) ! blue
2493	rgbcol(:, 2) = (/ 0, 170, 255 /) ! pale blue
2494	rgbcol(:, 3) = (/ 0, 255, 255 /) ! cyan
2495	rgbcol(:, 4) = (/ 0, 170, 0 /) ! dark green
2496	rgbcol(:, 5) = (/ 0, 255, 0 /) ! green
2497	rgbcol(:, 6) = (/ 0, 255, 170 /) ! blue-green
2498	rgbcol(:, 7) = (/ 128, 255, 0 /) ! yellow-green
2499	rgbcol(:, 8) = (/ 128, 170, 0 /) ! khaki
2500	rgbcol(:, 9) = (/ 255, 255, 0 /) ! yellow
2501	rgbcol(:,10) = (/ 255, 85, 0 /) ! orange
2502	rgbcol(:,11) = (/ 255, 0, 85 /) ! pink-ish
2503	rgbcol(:,12) = (/ 128, 0, 255 /) ! violet
2504	rgbcol(:,13) = (/ 255, 0, 255 /) ! magenta
2505	rgbcol(:,14) = (/ 170, 0, 128 /) ! purple
2506	!ma rgbcol(:,15) = (/ 255, 0, 85 /) ! red
2507
2508	! Write out the RGB pixels, one per point in the domain.
2509
2510	DO j=ny,1,-1
2511	DO i=1,nx
2512	IF ( map(i,j).LT.0 ) THEN
2513	icol = map(i,j)
2514	ELSE
2515	icol = MOD(map(i,j),ncol)
2516	ENDIF
2517	WRITE (lumapout,'(3i4)') &
2518	rgbcol(1,icol),rgbcol(2,icol),rgbcol(3,icol)
2519	ENDDO
2520	ENDDO
2521	CLOSE (lumapout)
2522	ENDIF ! (lwp)
2523
2524	DEALLOCATE (map)
2525
2526	END SUBROUTINE write_partition_map
2527
2528
2529	SUBROUTINE smooth_global_bathy(inbathy, imask)
2530	USE dom_oce
2531	USE domzgr, ONLY: rn_sbot_min, rn_sbot_max, rn_theta, rn_thetb, &
2532	rn_rmax, ln_s_sigma, rn_bb, rn_hc, fssig1, &
2533	namzgr_sco
2534	USE in_out_manager, ONLY: numnam
2535	IMPLICIT none
2536	!!----------------------------------------------------------------------
2537	!! Routine smooth_global_bathy
2538	!! Replicates the smoothing done on the decomposed domain in zgr_sco()
2539	!! in domzgr.F90. However, here the domain is NOT decomposed and
2540	!! every processor performs an identical smoothing on the whole model
2541	!! bathymetry. This is to ensure that the domain decomposition
2542	!! is done using a mask that is the same as that which is eventually
2543	!! computed after zgr_sco() has been called. (The smoothing process
2544	!! below can (erroneously) change whether grid points are wet or dry.)
2545	!!----------------------------------------------------------------------
2546	REAL(wp), INTENT(inout), DIMENSION(:,:) :: inbathy ! The bathymetry to
2547	! be smoothed
2548	INTEGER, INTENT(inout), DIMENSION(:,:) :: imask ! Mask holding index of
2549	! bottom level
2550	! Locals
2551	INTEGER :: ji, jj, jk, jl, ierr
2552	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
2553	INTEGER :: x_size, y_size
2554	REAL(wp) :: zrmax, zri, zrj, zcoeft
2555	REAL(wp), PARAMETER :: TOL_ZERO = 1.0E-20_wp ! Any value less than
2556	! this assumed zero
2557	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zenv, ztmp, zmsk, zbot, &
2558	zscosrf, zhbatt
2559	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgsigt3, zgdept
2560	!
2561	!!----------------------------------------------------------------------
2562
2563	! Do this because we've not decomposed the domain yet and therefore
2564	! jpi,jpj,nlc{i,j} etc. are not set.
2565	x_size = SIZE(inbathy, 1)
2566	y_size = SIZE(inbathy, 2)
2567
2568	ALLOCATE(zenv(x_size,y_size), ztmp(x_size,y_size), zmsk(x_size,y_size), &
2569	zbot(x_size,y_size), zgdept(x_size,y_size,jpkdta), zhbatt(x_size, y_size), &
2570	zscosrf(x_size,y_size), zgsigt3(x_size,y_size,jpkdta), Stat=ierr)
2571	IF( ierr /= 0 ) THEN
2572	CALL ctl_stop('smooth_global_bathy: ERROR - failed to allocate workspace arrays')
2573	RETURN
2574	ENDIF
2575
2576	REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching
2577	! parameters
2578	READ ( numnam, namzgr_sco )
2579
2580	zscosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
2581	zbot(:,:) = inbathy(:,:) ! ocean bottom depth
2582	! ! set maximum ocean depth
2583	inbathy(:,:) = MIN( rn_sbot_max, inbathy(:,:) )
2584
2585	WHERE( inbathy(:,:) > TOL_ZERO ) inbathy(:,:) = MAX( rn_sbot_min, inbathy(:,:) )
2586
2587	! use r-value to create hybrid coordinates
2588	zenv(:,:) = MAX( inbathy(:,:), rn_sbot_min )
2589	!
2590	! Smooth the bathymetry (if required)
2591	!
2592	jl = 0
2593	zrmax = 1._wp
2594	! ! ================ !
2595	DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop !
2596	! ! ================ !
2597	jl = jl + 1
2598	zrmax = 0._wp
2599	zmsk(:,:) = 0._wp
2600
2601	DO jj = 1, y_size
2602	DO ji = 1, x_size
2603	iip1 = MIN( ji+1, x_size ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
2604	ijp1 = MIN( jj+1, y_size ) ! force zrj = 0 on last row (jj=nclj+1 to jpj)
2605	zri = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
2606	zrj = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
2607	zrmax = MAX( zrmax, zri, zrj )
2608
2609	IF( zri > rn_rmax ) zmsk(ji ,jj ) = 1._wp
2610	IF( zri > rn_rmax ) zmsk(iip1,jj ) = 1._wp
2611	IF( zrj > rn_rmax ) zmsk(ji ,jj ) = 1._wp
2612	IF( zrj > rn_rmax ) zmsk(ji ,ijp1) = 1._wp
2613	END DO
2614	END DO
2615
2616	!
2617	IF(lwp)WRITE(numout,"('smooth_global_bathy : iter=',I5,' rmax=',F8.4,' nb of pt= ',I8)") &
2618	jl, zrmax, INT( SUM(zmsk(:,:) ) )
2619	!
2620
2621	! Copy current surface before next smoothing iteration
2622	ztmp(:,:) = zenv(:,:)
2623
2624	DO jj = 1, y_size
2625	DO ji = 1, x_size
2626	iip1 = MIN( ji+1, x_size ) ! last line (ji=nlci)
2627	ijp1 = MIN( jj+1, y_size ) ! last raw (jj=nlcj)
2628	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
2629	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
2630	IF( zmsk(ji,jj) == 1._wp ) THEN
2631	ztmp(ji,jj) = ( &
2632	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
2633	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) &
2634	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
2635	& ) / ( &
2636	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
2637	& + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) &
2638	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
2639	& )
2640	ENDIF
2641	END DO
2642	END DO
2643	!
2644	DO jj = 1,y_size
2645	DO ji = 1,x_size
2646	IF( zmsk(ji,jj) >= 1._wp-TOL_ZERO ) zenv(ji,jj) = MAX( ztmp(ji,jj), inbathy(ji,jj) )
2647	END DO
2648	END DO
2649	!
2650	! ! ================ !
2651	END DO ! End loop !
2652	! ! ================ !
2653	!
2654	! ! envelop bathymetry saved in zhbatt
2655	zhbatt(:,:) = zenv(:,:)
2656	! gphit calculated in nemo_init->dom_init->dom_hgr and dom_hgr requires that
2657	! partitioning already done. Could repeat its calculation here but since AMM doesn't
2658	! require it we leave it out for the moment ARPDBG
2659	CALL ctl_warn( ' ARPDBG - NOT checking whether s-coordinates are tapered in vicinity of the Equator' )
2660	!!$ IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
2661	!!$ CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
2662	!!$ DO jj = 1, jpj
2663	!!$ DO ji = 1, jpi
2664	!!$ ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 )
2665	!!$ hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
2666	!!$ END DO
2667	!!$ END DO
2668	!!$ ENDIF
2669
2670	! Subtract off rn_sbot_min so can check for land using zenv = LAND (0)
2671	inbathy(:,:) = zenv(:,:) - rn_sbot_min
2672
2673
2674	! ! =======================
2675	! ! s-ccordinate fields (gdep., e3.)
2676	! ! =======================
2677	!
2678	! non-dimensional "sigma" for model level depth at w- and t-levels
2679
2680	IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths
2681	! ! below rn_hc, with uniform sigma in shallower waters
2682	DO ji = 1, x_size
2683	DO jj = 1, y_size
2684
2685	IF( zhbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
2686	DO jk = 1, jpk
2687	zgsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
2688	END DO
2689	ELSE ! shallow water, uniform sigma
2690	DO jk = 1, jpk
2691	zgsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
2692	END DO
2693	ENDIF
2694	!
2695	DO jk = 1, jpk
2696	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
2697	zgdept (ji,jj,jk) = zscosrf(ji,jj) + (zhbatt(ji,jj)-rn_hc)zgsigt3(ji,jj,jk)+rn_hczcoeft
2698	END DO
2699	!
2700	END DO ! for all jj's
2701	END DO ! for all ji's
2702	ELSE
2703	CALL ctl_stop('STOP', &
2704	'partition_mod::smooth_global_bathy() only supports ln_s_sigma = .TRUE. currently!')
2705	END IF
2706
2707	! HYBRID scheme
2708	DO jj = 1, y_size
2709	DO ji = 1, x_size
2710	DO jk = 1, jpkm1
2711	IF( zbot(ji,jj) >= zgdept(ji,jj,jk) ) imask(ji,jj) = MAX( 2, jk )
2712	IF( zbot(ji,jj) == 0._wp ) imask(ji,jj) = 0
2713	END DO
2714	END DO
2715	END DO
2716
2717	! Dump to file for debugging ARPDBG
2718	IF(lwp)THEN
2719	OPEN(UNIT=1098, FILE='smoothed_bathy.dat', STATUS='REPLACE', &
2720	ACTION='WRITE', IOSTAT=jj)
2721	IF(jj == 0)THEN
2722	DO jj = 1, y_size
2723	DO ji = 1, x_size
2724	WRITE (1098,"(I4,1x,I4,3(E14.4,1x),I4)") ji, jj, &
2725	inbathy(ji,jj), zbot(ji,jj), &
2726	inbathy(ji,jj)-zbot(ji,jj), imask(ji,jj)
2727	END DO
2728	WRITE (1098,*)
2729	END DO
2730	CLOSE(1098)
2731	END IF
2732	END IF
2733
2734	END SUBROUTINE smooth_global_bathy
2735
2736
2737	SUBROUTINE global_bot_level(imask)
2738	USE par_oce, ONLY: jperio
2739	IMPLICIT none
2740	!!----------------------------------------------------------------------
2741	!! Compute the deepest level for any of the u,v,w or T grids. (Code
2742	!! taken from zgr_bot_level() and intermediate arrays for U and V
2743	!! removed.)
2744	!!----------------------------------------------------------------------
2745	INTEGER, DIMENSION(:,:), INTENT(inout) :: imask
2746	! Locals
2747	INTEGER :: ji, jj
2748	INTEGER :: x_size, y_size
2749
2750	! Do this because we've not decomposed the domain yet and therefore
2751	! jpi,jpj,nlc{i,j} etc. are not set.
2752	x_size = SIZE(imask, 1)
2753	y_size = SIZE(imask, 2)
2754
2755	imask(:,:) = MAX( imask(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
2756
2757	!
2758	! Compute and store the deepest bottom k-index of any grid-type at
2759	! each grid point.
2760	! For use in removing data below ocean floor from halo exchanges.
2761	DO jj = 1, y_size-1
2762	DO ji = 1, x_size-1
2763	imask(ji,jj) = MAX(imask(ji,jj)+1, & ! W (= T-level + 1)
2764	MAX(1, MIN(imask(ji+1,jj), imask(ji,jj)) ), & ! U
2765	MAX(1, MIN(imask(ji ,jj+1), imask(ji,jj))) ) ! V
2766	END DO
2767	imask(x_size,jj) = imask(x_size-1,jj)
2768	END DO
2769
2770	! Check on jperio because we've not set cyclic_bc in mapcomms yet
2771	IF(jperio == 0)THEN
2772	! In zgr_bat_ctl mbathy is zeroed along E/W boundaries if there are no
2773	! cyclic boundary conditions.
2774	imask(1,:) =0
2775	imask(x_size,:)=0
2776	ELSE IF(jperio == 1 .OR. jperio == 4 .OR. jperio == 6)THEN
2777	! Impose global E-W cyclic boundary conditions on the array holding the
2778	! deepest level
2779	imask(1,:) = imask(x_size - 1, :)
2780	imask(x_size,:) = imask(2,:)
2781	END IF
2782
2783	! Dump to file for debugging ARPDBG
2784	IF(lwp)THEN
2785	OPEN(UNIT=1098, FILE='bathy_bottom.dat', STATUS='REPLACE', &
2786	ACTION='WRITE', IOSTAT=jj)
2787	IF(jj == 0)THEN
2788	DO jj = 1, y_size
2789	DO ji = 1, x_size
2790	WRITE (1098,"(I4,1x,I4,1x,I4)") ji, jj, imask(ji,jj)
2791	END DO
2792	WRITE (1098,*)
2793	END DO
2794	CLOSE(1098)
2795	END IF
2796	END IF
2797
2798	END SUBROUTINE global_bot_level
2799
2800
2801	SUBROUTINE read_partition(ierr)
2802	USE par_oce, ONLY: jpni, jpnj, jpnij
2803	USE mapcomm_mod, ONLY: eidx, widx, sidx, nidx, trimmed, &
2804	pilbext, piubext, pjlbext, pjubext
2805	IMPLICIT none
2806	INTEGER, INTENT(out) :: ierr
2807	! Locals
2808	INTEGER, PARAMETER :: funit = 1099
2809	INTEGER :: idom, ndom
2810	CHARACTER(len=200) :: linein
2811	!======================================================================
2812
2813	ierr = 0
2814
2815	OPEN(UNIT=funit, file='partition.dat', status='OLD', &
2816	ACTION='READ', IOSTAT=ierr)
2817	IF(ierr /= 0)THEN
2818	CALL ctl_warn('read_partition: failed to read partitioning from file - will calculate it instead.')
2819	RETURN
2820	END IF
2821
2822	! Number of procs in x and y
2823	CALL read_next_line(funit, linein, ierr)
2824	READ(linein,FMT=*) jpni, jpnj
2825
2826	! Store their product
2827	jpnij = jpni*jpnj
2828
2829	! Check that the implied number of PEs matches that
2830	! in our MPI world
2831	ndom = jpni*jpnj
2832	IF(ndom /= mppsize)THEN
2833	CALL ctl_stop('STOP', &
2834	'read_partition: no. of PEs specified in partition.dat does not match no. of PEs in use by this job.')
2835	END IF
2836
2837	! Read the description of each sub-domain
2838	domains: DO idom = 1, ndom, 1
2839
2840	! Coordinates of bottom-left (SW) corner of domain
2841	CALL read_next_line(funit, linein, ierr)
2842	READ(linein,FMT=*) pielb(idom), pjelb(idom)
2843	! Top-right (NE) corner
2844	CALL read_next_line(funit, linein, ierr)
2845	READ(linein,FMT=*) pieub(idom), pjeub(idom)
2846	! Whether this domain has external boundaries and has been trimmed
2847	CALL read_next_line(funit, linein, ierr)
2848	READ(linein,FMT=*) pilbext(idom), trimmed(widx,idom)
2849	CALL read_next_line(funit, linein, ierr)
2850	READ(linein,FMT=*) piubext(idom), trimmed(eidx,idom)
2851	CALL read_next_line(funit, linein, ierr)
2852	READ(linein,FMT=*) pjlbext(idom), trimmed(sidx,idom)
2853	CALL read_next_line(funit, linein, ierr)
2854	READ(linein,FMT=*) pjubext(idom), trimmed(nidx,idom)
2855
2856	piesub(idom) = pieub(idom) - pielb(idom) + 1
2857	pjesub(idom) = pjeub(idom) - pjelb(idom) + 1
2858
2859	END DO domains
2860
2861	! All done - close the file
2862	CLOSE(UNIT=funit)
2863
2864	CALL finish_partition(fromFile=.TRUE.)
2865
2866	END SUBROUTINE read_partition
2867
2868	!===================================================================
2869
2870	SUBROUTINE write_partition
2871	USE par_oce, ONLY: jpni, jpnj
2872	USE mapcomm_mod, ONLY: eidx, widx, sidx, nidx, trimmed, &
2873	pjubext, pjlbext, piubext, pilbext, &
2874	pielb, pieub, pjelb, pjeub
2875	IMPLICIT none
2876	INTEGER, PARAMETER :: funit = 1099
2877	INTEGER :: ierr
2878	INTEGER :: idom
2879
2880	! Only PE 0 (narea==1) writes this file
2881	IF(narea /= 1) RETURN
2882
2883	OPEN(UNIT=funit, file='partition.dat.new', status='REPLACE', &
2884	ACTION='WRITE', IOSTAT=ierr)
2885	IF(ierr /= 0)THEN
2886	CALL ctl_warn('write_partition: failed to write partition description to file.')
2887	RETURN
2888	END IF
2889	WRITE(UNIT=funit,FMT="('# jpni jpnj')")
2890	WRITE(UNIT=funit,FMT="(I5,1x,I5)") jpni, jpnj
2891
2892	DO idom = 1, mppsize, 1
2893	WRITE(UNIT=funit,FMT="('# Domain: ',I5)") idom
2894	IF(idom==1)WRITE(UNIT=funit,FMT="('# Lower bounds: x y')")
2895	WRITE(UNIT=funit,FMT="(I5,1x,I5)") pielb(idom), pjelb(idom)
2896	IF(idom==1)WRITE(UNIT=funit,FMT="('# Upper bounds: x y')")
2897	WRITE(UNIT=funit,FMT="(I5,1x,I5)") pieub(idom), pjeub(idom)
2898	IF(idom==1)WRITE(UNIT=funit,FMT="('# x: Lower bound external, trimmed')")
2899	WRITE(UNIT=funit,FMT="(L5,1x,L5)") pilbext(idom), trimmed(widx,idom)
2900	IF(idom==1)WRITE(UNIT=funit,FMT="('# x: Upper bound external, trimmed')")
2901	WRITE(UNIT=funit,FMT="(L5,1x,L5)") piubext(idom), trimmed(eidx,idom)
2902	IF(idom==1)WRITE(UNIT=funit,FMT="('# y: Lower bound external, trimmed')")
2903	WRITE(UNIT=funit,FMT="(L5,1x,L5)") pjlbext(idom), trimmed(sidx,idom)
2904	IF(idom==1)WRITE(UNIT=funit,FMT="('# y: Upper bound external, trimmed')")
2905	WRITE(UNIT=funit,FMT="(L5,1x,L5)") pjubext(idom), trimmed(nidx,idom)
2906	END DO
2907
2908	CLOSE(UNIT=funit)
2909
2910	END SUBROUTINE write_partition
2911
2912	SUBROUTINE read_next_line(funit, linein, ierr)
2913	IMPLICIT none
2914	!!------------------------------------------------------------------
2915	INTEGER, INTENT( in) :: funit ! Unit no. to read
2916	CHARACTER(len=200), INTENT(out) :: linein ! String containing next
2917	! non-comment line in file
2918	INTEGER, INTENT(out) :: ierr ! Error flag (0==OK)
2919	!!------------------------------------------------------------------
2920
2921	ierr = 0
2922
2923	READ(UNIT=funit,FMT="(200A)") linein
2924
2925	! Comment lines begin with '#'. Skip those plus any blank
2926	! lines...
2927	DO WHILE( INDEX( TRIM(ADJUSTL(linein)),'#') /= 0 .OR. &
2928	LEN_TRIM(linein) == 0 )
2929	READ(UNIT=funit,FMT="(200A)") linein
2930	END DO
2931
2932	WRITE(,)'returning linein >>'//linein//'<<'
2933
2934	END SUBROUTINE read_next_line
2935
2936	END MODULE partition_mod

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

Download in other formats: