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/TOOLS/RK_PARTITION/src – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/partition_mod.f90 @ 4396

Last change on this file since 4396 was 4396, checked in by trackstand2, 10 years ago
WIP on reading grid dims from netcdf
File size: 82.4 KB

Line
1	MODULE partition_mod
2	USE par_kind
3	IMPLICIT none
4
5	PRIVATE
6
7	LOGICAL, PARAMETER :: PARTIT_DEBUG = .TRUE. ! Ctrl of debug messages from
8	! partitioning code
9	LOGICAL, PARAMETER :: TRIM_DEBUG = .FALSE. ! Ctrl of debug messages from
10	! trimming code.
11
12	INTEGER, PUBLIC, SAVE :: jpni, jpnj, jpnij
13	INTEGER, PUBLIC, SAVE :: jpiglo, jpjglo
14
15	! Value representing land in the mask used for partitioning and message
16	! trimming
17	INTEGER, PUBLIC, PARAMETER :: LAND = 0
18
19	INTEGER, PARAMETER :: MAX_FACTORS = 20 ! Max no. of factors factor() can return
20	! Most of the following parameters should be input vars - TODO ARPDEBUG
21	INTEGER, PARAMETER :: jpreci = 1, jprecj = 1
22	INTEGER, PARAMETER :: nextra = 2 ! Extent to back off from coastlines when
23	! trimming - causes sea-ice problems if 0.
24	LOGICAL, PARAMETER :: trim_land = .TRUE. ! Whether to trim land from sub domains
25	LOGICAL, PUBLIC, SAVE :: cyclic_bc ! Whether we have cyclic bound. cond.'s
26
27	INTEGER, SAVE :: nprocp ! No. of PEs to partition over
28
29	! Information held by the ensemble leader about all processes.
30	! pielb Lower (west) longitude bound index.
31	! pieub Upper (east) longitude bound index.
32	! piesub Number of longitude gridpoints.
33	! pilbext True if the lower longitude boundary is external.
34	! piubext True if the upper longitude boundary is external.
35	! pjelb Lower (south) latitude bound index.
36	! pjeub Upper (north) latitude bound index.
37	! pjesub Number of latitude gridpoints.
38	! pjlbext True if the lower latitude boundary is external.
39	! pjubext True if the upper latitude boundary is external.
40
41	INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: pielb, pieub, piesub, pjelb, pjeub, pjesub, pnactive
42	LOGICAL, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: pilbext, piubext, pjlbext, pjubext
43
44	! Process dependent partitioning information.
45
46	! ielb Lower (west) longitude bound index.
47	! ieub Upper (east) longitude bound index.
48	! iesub Number of longitude gridpoints.
49	! ilbext True if the lower longitude boundary is external.
50	! iubext True if the upper longitude boundary is external.
51	! jelb Lower (south) latitude bound index.
52	! jeub Upper (north) latitude bound index.
53	! jesub Number of latitude gridpoints.
54	! jlbext True if the lower latitude boundary is external.
55	! jubext True if the upper latitude boundary is external.
56	! trimmed True if boundary on given proc trimmed.
57
58	INTEGER, SAVE :: ielb, ieub, iesub
59	INTEGER, SAVE :: jelb, jeub, jesub
60	LOGICAL, SAVE :: ilbext, iubext, jlbext, jubext
61
62	! Stores whether a domain's boundaries have been trimmed as
63	! trimmed(boundary, PE) where boundary is one of {n,e,s,w}idx
64	! for the Northern, Eastern, Southern or Western boundary, respectively.
65	LOGICAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: trimmed
66	INTEGER, PARAMETER :: nidx = 1, eidx= 2, sidx = 3, widx = 4
67
68	! Parameters for the cost function used when evaluating different
69	! possible domain partitions
70
71	! Mnemonics:
72	! p = process (i.e. core)
73	! n = node
74	! l = length (number of points)
75	! c = count (number of messages)
76	! i = internal (intra-node, or on-node)
77	! x = external (inter-node, or off-node)
78
79	INTEGER, PARAMETER :: pv_index_overall = 1
80	INTEGER, PARAMETER :: pv_index_wet = 2
81	INTEGER, PARAMETER :: pv_index_dry = 3
82	INTEGER, PARAMETER :: pv_index_pli = 4
83	INTEGER, PARAMETER :: pv_index_plx = 5
84	INTEGER, PARAMETER :: pv_index_pci = 6
85	INTEGER, PARAMETER :: pv_index_pcx = 7
86	INTEGER, PARAMETER :: pv_index_nli = 8
87	INTEGER, PARAMETER :: pv_index_nlx = 9
88	INTEGER, PARAMETER :: pv_index_nci = 10
89	INTEGER, PARAMETER :: pv_index_ncx = 11
90	INTEGER, PARAMETER :: pv_index_tli = 12
91	INTEGER, PARAMETER :: pv_index_tlx = 13
92	INTEGER, PARAMETER :: pv_index_tci = 14
93	INTEGER, PARAMETER :: pv_index_tcx = 15
94
95	INTEGER, PARAMETER :: pv_index_pcomms = 16
96	INTEGER, PARAMETER :: pv_index_ncomms = 17
97	INTEGER, PARAMETER :: pv_index_tcomms = 18
98
99	INTEGER, PARAMETER :: pv_num_scores = 18
100
101	REAL(wp),PARAMETER :: pv_awful = 1.0e20
102
103	PUBLIC allocate_arrays, nemo_partition, write_partition, partition_mca_rk
104
105	CONTAINS
106
107	SUBROUTINE allocate_arrays(num_pes, ierr)
108	!!----------------------------------------------------------------------
109	! Allocate partitioning arrays.
110	!!----------------------------------------------------------------------
111	IMPLICIT none
112	INTEGER, INTENT(In) :: num_pes
113	INTEGER, INTENT(out) :: ierr
114	!!----------------------------------------------------------------------
115
116	IF ( .NOT. ALLOCATED(pielb) ) THEN
117	ALLOCATE (pielb(num_pes), pieub(num_pes), piesub(num_pes), &
118	pilbext(num_pes), piubext(num_pes), &
119	pjelb(num_pes), pjeub(num_pes), pjesub(num_pes), &
120	pjlbext(num_pes), pjubext(num_pes), pnactive(num_pes), &
121	trimmed(4,num_pes), Stat = ierr)
122	IF(ierr /= 0)THEN
123	WRITE(,)('ERROR: failed to allocate partitioning arrays')
124	RETURN
125	END IF
126	ENDIF
127
128	END SUBROUTINE allocate_arrays
129
130
131	SUBROUTINE nemo_partition( num_pes )
132	!!----------------------------------------------------------------------
133	!! * ROUTINE nemo_partition *
134	!!
135	!! ** Purpose : Work out a sensible factorisation of the number of
136	!! processors for the x and y dimensions.
137	!! ** Method :
138	!!----------------------------------------------------------------------
139	IMPLICIT None
140	INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have
141	!
142	INTEGER :: ii, jj, jn, ierr
143	INTEGER :: ifact1, ifact2 ! factors of num_pes, ifact1 <= ifact2
144	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilcit, ilcjt, iimppt, ijmppt
145	INTEGER :: nreci, nrecj, iresti, irestj
146	INTEGER :: jpi, jpj, nbondi, nbondj
147	INTEGER :: nldi, nldj, nlei, nlej, nlci, nlcj
148	INTEGER :: zidom, zjdom
149	!!----------------------------------------------------------------------
150
151	! Copy no. of PEs into module variable
152	nprocp = num_pes
153
154	! Factorise the number of processors into ifact1*ifact2, such that
155	! ifact1 and ifact2 are as nearly equal as possible.
156
157	CALL sqfact( num_pes, ifact1, ifact2 )
158
159	! Make sure that the smaller dimension of the processor grid
160	! is given the smaller dimension of the global domain
161	IF( jpiglo <= jpjglo) THEN
162	jpni = ifact1
163	jpnj = ifact2
164	ELSE
165	jpni = ifact2
166	jpnj = ifact1
167	ENDIF
168
169	! Calc jp{i,j} as done in nemo_init
170	jpi = ( jpiglo-2jpreci + (jpni-1) ) / jpni + 2jpreci ! first dim.
171	jpj = ( jpjglo-2jprecj + (jpnj-1) ) / jpnj + 2jprecj ! second dim.
172
173	! This should never happen
174	IF( (jpni*jpnj) /= num_pes) THEN
175	WRITE (, ) 'ERROR: internal error - factorisation of number of PEs failed'
176	RETURN
177	ENDIF
178
179	! This should only happen if num_pes is prime
180	IF( ifact1 <= 1 ) THEN
181	WRITE (, ) 'ERROR: factorisation of number of PEs failed'
182	RETURN
183	ENDIF
184	!
185	jpnij = jpni*jpnj
186
187	WRITE(*,FMT="(/'Generating REGULAR partition on ',I4,'x',I4,' process grid'/)") &
188	jpni, jpnj
189	!
190	! 1. Dimension arrays for subdomains
191	! -----------------------------------
192	! Computation of local domain sizes ilcit() ilcjt()
193	! These dimensions depend on global sizes jpni,jpnj and jpiglo,jpjglo
194	! The subdomains are squares lesser than or equal to the global
195	! dimensions divided by the number of processors minus the overlap
196	! array (cf. par_oce.F90).
197	ALLOCATE(ilcit(jpni,jpnj), ilcjt(jpni,jpnj), &
198	iimppt(jpni,jpnj), ijmppt(jpni,jpnj), Stat=ierr)
199
200	nreci = 2 * jpreci
201	nrecj = 2 * jprecj
202	iresti = MOD( jpiglo - nreci , jpni )
203	irestj = MOD( jpjglo - nrecj , jpnj )
204
205	IF( iresti == 0 ) iresti = jpni
206	DO jj = 1, jpnj
207	DO ii = 1, iresti
208	ilcit(ii,jj) = jpi
209	END DO
210	DO ii = iresti+1, jpni
211	ilcit(ii,jj) = jpi -1
212	END DO
213	END DO
214
215	IF( irestj == 0 ) irestj = jpnj
216	DO ii = 1, jpni
217	DO jj = 1, irestj
218	ilcjt(ii,jj) = jpj
219	END DO
220	DO jj = irestj+1, jpnj
221	ilcjt(ii,jj) = jpj -1
222	END DO
223	END DO
224
225	zidom = nreci
226	DO ii = 1, jpni
227	zidom = zidom + ilcit(ii,1) - nreci
228	END DO
229	!WRITE(,)
230	!WRITE(,)' sum ilcit(i,1) = ', zidom, ' jpiglo = ', jpiglo
231
232	zjdom = nrecj
233	DO jj = 1, jpnj
234	zjdom = zjdom + ilcjt(1,jj) - nrecj
235	END DO
236	!WRITE(,)' sum ilcit(1,j) = ', zjdom, ' jpjglo = ', jpjglo
237	!WRITE(,)
238
239
240	! 2. Index arrays for subdomains
241	! -------------------------------
242
243	iimppt(:,:) = 1
244	ijmppt(:,:) = 1
245
246	IF( jpni > 1 ) THEN
247	DO jj = 1, jpnj
248	DO ii = 2, jpni
249	iimppt(ii,jj) = iimppt(ii-1,jj) + ilcit(ii-1,jj) - nreci
250	END DO
251	END DO
252	ENDIF
253
254	IF( jpnj > 1 ) THEN
255	DO jj = 2, jpnj
256	DO ii = 1, jpni
257	ijmppt(ii,jj) = ijmppt(ii,jj-1)+ilcjt(ii,jj-1)-nrecj
258	END DO
259	END DO
260	ENDIF
261
262	! 3. Subdomain description
263	! ------------------------
264
265	DO jn = 1, jpnij
266	ii = 1 + MOD( jn-1, jpni )
267	jj = 1 + (jn-1) / jpni
268	nlci = ilcit (ii,jj)
269	nlcj = ilcjt (ii,jj)
270
271	pielb(jn) = iimppt(ii,jj)
272	pjelb(jn) = ijmppt(ii,jj)
273
274	nbondj = -1 ! general case
275	IF( jn > jpni ) nbondj = 0 ! first row of processor
276	IF( jn > (jpnj-1)*jpni ) nbondj = 1 ! last row of processor
277	IF( jpnj == 1 ) nbondj = 2 ! one processor only in j-direction
278	nbondi = 0 !
279	IF( MOD( jn, jpni ) == 1 ) nbondi = -1 !
280	IF( MOD( jn, jpni ) == 0 ) nbondi = 1 !
281	IF( jpni == 1 ) nbondi = 2 ! one processor only in i-direction
282	nldi = 1 + jpreci
283	nlei = nlci - jpreci
284	pilbext(jn) = .FALSE.
285	piubext(jn) = .FALSE.
286	IF( nbondi == -1 .OR. nbondi == 2 )THEN
287	nldi = 1
288	pilbext(jn) = .TRUE.
289	END IF
290	IF( nbondi == 1 .OR. nbondi == 2 )THEN
291	nlei = nlci
292	piubext(jn) = .TRUE.
293	END IF
294
295	nldj = 1 + jprecj
296	nlej = nlcj - jprecj
297	pjlbext(jn) = .FALSE.
298	pjubext(jn) = .FALSE.
299	IF( nbondj == -1 .OR. nbondj == 2 )THEN
300	nldj = 1
301	pjlbext(jn) = .TRUE.
302	END IF
303	IF( nbondj == 1 .OR. nbondj == 2 )THEN
304	nlej = nlcj
305	pjubext(jn) = .TRUE.
306	END IF
307
308	! nldj corresponds to domain lower bound and nlej the upper bound
309	pieub(jn) = pielb(jn) + (nlei - nldi + 1) - 1
310	pjeub(jn) = pjelb(jn) + (nlej - nldj + 1) - 1
311
312	! No trimming in this partitioning
313	trimmed(:,jn) = .FALSE.
314	END DO
315
316	END SUBROUTINE nemo_partition
317
318	!==========================================================================
319
320	SUBROUTINE sqfact ( kn, kna, knb )
321	!!----------------------------------------------------------------------
322	!! * ROUTINE sqfact *
323	!!
324	!! ** Purpose : return factors (kna, knb) of kn, such that
325	!! (1) kna*knb=kn
326	!! (2) kna and knb are as near equal as possible
327	!! (3) kna < knb
328	!! ** Method : Search backwards from the square root of kn,
329	!! until we find an integer that cleanly divides kn
330	!! ** Preconditions : kn must be positive
331	!!----------------------------------------------------------------------
332	INTEGER, INTENT(in ) :: kn
333	INTEGER, INTENT( out) :: kna, knb
334
335	! Search backwards from the square root of n.
336
337	fact_loop: DO kna=INT(SQRT(REAL(kn))),1,-1
338	IF ( kn/kna*kna == kn ) THEN
339	EXIT fact_loop
340	ENDIF
341	END DO fact_loop
342
343	IF( kna < 1 ) kna = 1
344
345	! kna divides kn cleanly. Work out the other factor.
346	knb = kn/kna
347
348	END SUBROUTINE sqfact
349
350	!==========================================================================
351
352	SUBROUTINE partition_mca_rk(mask, nproc, ncpn, ierr)
353	IMPLICIT NONE
354	!!------------------------------------------------------------------
355	!! Multi-Core Aware recursive partitioning of the domain. As for
356	!! partition_rk but choose the partion
357	!! so as to minimize off-node MPI communication
358	!!
359	!! Original by Stephen Pickles for POLCOMS code.
360	!! Implementation in NEMO by Andrew Porter, 26/01/2012
361	!!------------------------------------------------------------------
362	! Subroutine arguments.
363	INTEGER, INTENT(in) :: mask(:,:)
364	INTEGER, INTENT(in) :: nproc
365	INTEGER, INTENT(in) :: ncpn ! No. of cores per node on target machine
366	INTEGER, INTENT(out) :: ierr
367	! Local variables
368	INTEGER :: ii
369	INTEGER, DIMENSION(MAX_FACTORS) :: fx, fy, factors
370	INTEGER, DIMENSION(MAX_FACTORS) :: df, multiplicity
371	INTEGER :: nfx, nfy, nfactors, ndf, nperms
372	!INTEGER :: check_nprocx, check_nprocy, check_nprocp
373	INTEGER :: iperm
374	CHARACTER(LEN=256) :: fstr
375	INTEGER :: nprocx, nprocy
376
377	! A high score is bad
378	REAL(wp), DIMENSION(pv_num_scores) :: score
379	REAL(wp) :: best_score
380	INTEGER :: best_perm
381	REAL(wp), DIMENSION(2,pv_num_scores) :: best, wrst
382	!!------------------------------------------------------------------
383
384	! Copy no. of PEs into module variable
385	nprocp = nproc
386
387	! Factorise the total number of MPI processes that we have
388	CALL factor (factors,MAX_FACTORS,nfactors,nproc,ierr)
389
390	IF ( ierr.NE.0 ) THEN
391	WRITE (*,FMT="(/'partition_mca_rk: factorisation failed ',I3)") ierr
392	RETURN
393	ELSE
394	WRITE (*,FMT="(/'partition_mca_rk: factors are: ',(I3,1x))") factors(:nfactors)
395	WRITE (,)
396	ENDIF
397
398	CALL calc_perms( nfactors, factors, &
399	ndf, df, multiplicity, &
400	nperms )
401
402	DO ii=1,pv_num_scores
403	best(1,ii) = pv_awful
404	best(2,ii) = -1.0_wp
405	END DO
406	DO ii=1,pv_num_scores
407	wrst(1,ii) = 0.0_wp
408	wrst(2,ii) = -1.0_wp
409	END DO
410
411	WRITE(*,"('% Partn',2X,10(A4,2X),4(A9,1X),A7)") &
412	'Wet', 'Dry', &
413	'pli', 'plx', 'pci', 'pcx', &
414	'nlx', 'ncx', 'tlx', 'tcx', &
415	'P comms', 'N comms', 'T comms', 'Overall', &
416	'Factors'
417
418	! perm_loop: DO iperm=myinst, nperms-1, nprocp
419	! Now serial so we have to do them all ourselves - could use
420	! OpenMP?
421	perm_loop: DO iperm=1, nperms-1, 1
422
423	CALL get_perm_factors( iperm, nfactors, ndf, df, multiplicity, &
424	fx, nfx, fy, nfy, &
425	nprocx, nprocy)
426
427	CALL partition_rk_core(mask, jpiglo, jpjglo, &
428	MAX_FACTORS, fx, nfx, fy, nfy, nprocp, ierr)
429
430	IF (ierr .NE. 0) CYCLE perm_loop
431	CALL finish_partition(mask)
432
433	! Compute the cost function for this partition
434	!
435	CALL eval_partition(ncpn, jpiglo, jpjglo, mask, score)
436	CALL factor_string(fstr,nfx,fx,nfy,fy)
437
438	WRITE (6,'(''%'',I6,1X,10(I5,1X),3(F9.2,1X),E12.4,1x,(A))') &
439	iperm, &
440	INT(score(pv_index_wet)), &
441	INT(score(pv_index_dry)), &
442	INT(score(pv_index_pli)), &
443	INT(score(pv_index_plx)), &
444	INT(score(pv_index_pci)), &
445	INT(score(pv_index_pcx)), &
446	INT(score(pv_index_nlx)), &
447	INT(score(pv_index_ncx)), &
448	INT(score(pv_index_tlx)), &
449	INT(score(pv_index_tcx)), &
450	score(pv_index_pcomms), &
451	score(pv_index_ncomms), &
452	score(pv_index_tcomms), &
453	score(pv_index_overall), &
454	TRIM(fstr)
455
456	DO ii=1,pv_num_scores
457	IF (score(ii) .LT. best(1,ii)) THEN
458	best(1,ii) = score(ii)
459	best(2,ii) = iperm
460	END IF
461	IF (score(ii) .GT. wrst(1,ii)) THEN
462	wrst(1,ii) = score(ii)
463	wrst(2,ii) = iperm
464	END IF
465	END DO
466
467	END DO perm_loop
468
469	! Now choose the "best" partition
470
471	best_score = best(1,pv_index_overall)
472	best_perm = INT( best(2,pv_index_overall) )
473
474	WRITE (*,'(/A32,A20,A20)') &
475	' ',' best score / perm ',' worst score / perm'
476	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') 'overall: ', &
477	best(1,pv_index_overall), best(2,pv_index_overall), &
478	wrst(1,pv_index_overall), wrst(2,pv_index_overall)
479	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') 'wet points: ', &
480	best(1,pv_index_wet), best(2,pv_index_wet), &
481	wrst(1,pv_index_wet), wrst(2,pv_index_wet)
482	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') 'dry points: ', &
483	best(1,pv_index_dry), best(2,pv_index_dry), &
484	wrst(1,pv_index_dry), wrst(2,pv_index_dry)
485	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
486	'proc max on-node wet points: ', &
487	best(1,pv_index_pli), best(2,pv_index_pli), &
488	wrst(1,pv_index_pli), wrst(2,pv_index_pli)
489	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
490	'proc max off-node wet points: ', &
491	best(1,pv_index_plx), best(2,pv_index_plx), &
492	wrst(1,pv_index_plx), wrst(2,pv_index_plx)
493	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
494	'proc max on-node messages: ', &
495	best(1,pv_index_pci), best(2,pv_index_pci), &
496	wrst(1,pv_index_pci), wrst(2,pv_index_pci)
497	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
498	'proc max off-node messages: ', &
499	best(1,pv_index_pcx), best(2,pv_index_pcx), &
500	wrst(1,pv_index_pcx), wrst(2,pv_index_pcx)
501	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
502	'node max off-node wet points: ', &
503	best(1,pv_index_nlx), best(2,pv_index_nlx), &
504	wrst(1,pv_index_nlx), wrst(2,pv_index_nlx)
505	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
506	'node max off-node messages: ', &
507	best(1,pv_index_ncx), best(2,pv_index_ncx), &
508	wrst(1,pv_index_ncx), wrst(2,pv_index_ncx)
509	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
510	'total off-node wet points: ', &
511	best(1,pv_index_tlx), best(2,pv_index_tlx), &
512	wrst(1,pv_index_tlx), wrst(2,pv_index_tlx)
513	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
514	'per core communications cost: ', &
515	best(1,pv_index_pcomms), best(2,pv_index_pcomms), &
516	wrst(1,pv_index_pcomms), wrst(2,pv_index_pcomms)
517	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
518	'per node communications cost: ', &
519	best(1,pv_index_ncomms), best(2,pv_index_ncomms), &
520	wrst(1,pv_index_ncomms), wrst(2,pv_index_ncomms)
521	WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') &
522	'network communications cost: ', &
523	best(1,pv_index_tcomms), best(2,pv_index_tcomms), &
524	wrst(1,pv_index_tcomms), wrst(2,pv_index_tcomms)
525
526	WRITE (*,"(/'partition_mca_rk: overall best perm is ',I6,', score=',F12.3/)") &
527	best_perm, best_score
528
529	! Use the same partition on all processes
530
531	! If a particular factorisation has been forced by
532	! the nn_{x,y}factors fields in the nammpp section of the namelist
533	! then use that one instead
534
535	!!$ IF ((nxfactors + nyfactors) > 0) THEN
536	!!$ check_nprocx = 1
537	!!$ check_nprocy = 1
538	!!$ DO ii=1,nxfactors
539	!!$ check_nprocx = check_nprocx*xfactors(ii)
540	!!$ END DO
541	!!$ DO ii=1,nyfactors
542	!!$ check_nprocy = check_nprocy*yfactors(ii)
543	!!$ END DO
544	!!$ check_nprocp = check_nprocx*check_nprocy
545	!!$ IF (check_nprocp .EQ. nprocp) THEN
546	!!$ nprocx = check_nprocx
547	!!$ nprocy = check_nprocy
548	!!$ nfx = nxfactors
549	!!$ nfy = nyfactors
550	!!$ fx(1:nfx) = xfactors(1:nfx)
551	!!$ fy(1:nfy) = yfactors(1:nfy)
552	!!$ ELSE
553	!!$ CALL ctl_stop('STOP', 'part_mca_rk: '// &
554	!!$ 'requested factorisation does not match nprocp')
555	!!$ END IF
556	!!$ ELSE
557	! Use the best factorisation found above
558	IF (best_perm < 0.0_wp) THEN
559	WRITE (,) 'partition_mca_rk: no feasible partition found'
560	ierr = 1
561	RETURN
562	END IF
563	CALL get_perm_factors(best_perm, nfactors, ndf, df, multiplicity, &
564	fx, nfx, fy, nfy, &
565	nprocx, nprocy)
566	!!$ END IF
567
568	! Set corresponding NEMO variables for PE grid, even though it is now
569	! rather irregular
570	jpni = nprocx
571	jpnj = nprocy
572
573	WRITE (,)
574	WRITE (*,'(A39)',advance='no') &
575	'partition_mca_rk: partitioning with factors '
576	CALL print_factors(nfx, fx, nfy, fy)
577
578	CALL partition_rk_core ( mask, jpiglo, jpjglo, &
579	MAX_FACTORS, &
580	fx, nfx, fy, nfy, &
581	nproc, ierr )
582
583	IF (ierr .NE. 0) THEN
584	WRITE (,) 'partition_mca_rk: failed to apply selection partition'
585	RETURN
586	END IF
587	CALL finish_partition(mask)
588
589	CALL eval_partition(ncpn, jpiglo, jpjglo, mask, score)
590	CALL factor_string(fstr,nfx,fx,nfy,fy)
591	WRITE(*,'(10(A7,1X),4(A9,1X),A7)') &
592	'Wet', 'Dry', &
593	'pli', 'plx', 'pci', 'pcx', &
594	'nlx', 'ncx', 'tlx', 'tcx', &
595	'P comms', 'N comms', 'T comms', 'Overall', &
596	'Factors'
597	WRITE (*,'(10(F7.0,1X),4(F9.2,1X),A50)') &
598	score(pv_index_wet), &
599	score(pv_index_dry), &
600	score(pv_index_pli), &
601	score(pv_index_plx), &
602	score(pv_index_pci), &
603	score(pv_index_pcx), &
604	score(pv_index_nlx), &
605	score(pv_index_ncx), &
606	score(pv_index_tlx), &
607	score(pv_index_tcx), &
608	score(pv_index_pcomms), &
609	score(pv_index_ncomms), &
610	score(pv_index_tcomms), &
611	score(pv_index_overall), &
612	fstr
613
614	END SUBROUTINE partition_mca_rk
615
616
617	SUBROUTINE partition_rk_core( mask, nx, ny, maxfax, fx, nfx, fy, nfy, &
618	nproc, ierr )
619	IMPLICIT NONE
620	!!------------------------------------------------------------------
621	!!------------------------------------------------------------------
622	INTEGER, INTENT(in) :: nx, ny
623	INTEGER, INTENT(in) :: mask(nx,ny)
624	INTEGER, INTENT(in) :: maxfax, nfx, nfy
625	INTEGER, INTENT(in) :: fx(maxfax), fy(maxfax)
626	INTEGER, INTENT(in) :: nproc
627	INTEGER, INTENT(out) :: ierr
628	! Local variables
629	INTEGER :: istart, jstart, istop, jstop
630	INTEGER :: f,gnactive
631	INTEGER :: i,ifax,nfax,ifin,ifx,ify,ilb,iproc
632	INTEGER :: ist,isub,isub_old,isub_new,iub
633	INTEGER :: j,jfin,jlb,jst,jub,line
634	INTEGER :: ngone,nsub_old,nsub_new,ntarget,ntry
635	LOGICAL :: partx
636
637	! Zero the error flag
638	ierr = 0
639
640	! Count the significant (non-zero) factors.
641	nfax = nfx+nfy
642	IF( PARTIT_DEBUG )THEN
643	WRITE (*,"('jpni = ',I3,1x,'nfx = ',I3,/'Factorized into: ',(I3,1x))") &
644	jpni, nfx, fx(:nfx)
645	WRITE (*,"('jpnj = ',I3,1x,'nfy = ',I3,/'Factorized into: ',(I3,1x))") &
646	jpnj, nfy, fy(:nfy)
647	WRITE (*,"('Partitioning a total of ',I3,' times')") nfax
648	END IF
649
650	! Define the full domain as the one and only sub-domain.
651	istart = 1
652	istop = nx
653	jstart = 1
654	jstop = ny
655
656	nsub_old = 1
657	pielb(nsub_old) = istart
658	pjelb(nsub_old) = jstart
659	pieub(nsub_old) = istop
660	pjeub(nsub_old) = jstop
661	! Count the number of active points.
662	gnactive = 0
663	DO j=jstart,jstop,1
664	DO i=istart,istop,1
665	IF ( mask(i,j) == 1 ) gnactive = gnactive+1
666	ENDDO
667	ENDDO
668
669	IF( PARTIT_DEBUG )THEN
670	WRITE (,) 'Total wet points ',gnactive
671	END IF
672
673	pnactive(nsub_old) = gnactive
674
675	! Partition for each factor.
676	DO ifax=1,nfax
677	IF ( ifax.EQ.1 ) THEN
678	! Start by partitioning the dimension with more factors.
679	partx = nfx.GE.nfy
680	ifx = 0
681	ify = 0
682	ELSE
683	! Try to toggle the partitioning direction.
684	partx = .NOT.partx
685	IF ( partx ) THEN
686	! If we have run out of factors in x, switch to y.
687	partx = .NOT. ifx+1.GT.nfx
688	ELSE
689	! If we have run out of factors in y, switch to x.
690	partx = ify+1.GT.nfy
691	ENDIF
692	ENDIF
693
694	IF( PARTIT_DEBUG )THEN
695	WRITE (,)
696	WRITE (,) '###########################################'
697	WRITE (,)
698	WRITE (,) 'Partition ',ifax,'partx = ',partx
699	END IF
700
701	IF ( partx ) THEN
702	! ============================================================
703	! Partition in x.
704	! ============================================================
705	ifx = ifx+1
706	f = fx(ifx)
707	nsub_new = nsub_old*f
708
709	IF( PARTIT_DEBUG )THEN
710	WRITE (,) 'Partitioning in x from ',nsub_old &
711	,' to ',nsub_new
712	END IF
713
714	DO isub_old=1, nproc, nproc/nsub_old
715	! Check that there are sufficient points to factorise.
716	IF ( pieub(isub_old)-pielb(isub_old)+1.LT.f ) THEN
717	WRITE (,) &
718	'partition_rk: insufficient points to partition'
719	ierr = 1
720	RETURN
721	ENDIF
722
723	! Set the target number of points in the new sub-domains.
724	ntarget = NINT(REAL(pnactive(isub_old))/REAL(f))
725	ist = pielb(isub_old)
726	iub = pieub(isub_old)
727	jlb = pjelb(isub_old)
728	jub = pjeub(isub_old)
729
730	IF( PARTIT_DEBUG )THEN
731	WRITE (*,"(/'=======================================')")
732	WRITE (*,"(/'Partitioning sub-domain ',I3,': (',I3,':',I3,')(',I3,':',I3,')')") &
733	isub_old,ist,iub,jlb,jub
734	WRITE (*,"('Old partition has ',I8,' points')") pnactive(isub_old)
735	WRITE (*,"('Target is ',I8,' points')") ntarget
736	END IF
737	! Create the new sub-domains.
738	ngone = 0
739	DO isub=1,f,1
740	isub_new = isub_old + (isub-1)*nproc/nsub_new
741
742	IF( PARTIT_DEBUG )THEN
743	WRITE (,)
744	WRITE (,) 'Creating new sub-domain ',isub_new
745	END IF
746	! The last new sub-domain takes the remaining data.
747	IF ( isub.EQ.f ) THEN
748	ifin = iub
749	ELSE
750	! Scan lines in x counting active grid points until we
751	! exceed the target.
752	ntry = 0
753	ifin = -1
754	scanrows: DO i=ist,iub,1
755	! Count up the contribution from the next line.
756	line = 0
757	DO j=jlb,jub,1
758	IF ( mask(i,j) == 1 ) line = line+1
759	ENDDO
760
761	! Check against the target.
762	IF ( ntry+line.GT.ntarget ) THEN
763	! Is it nearer the target to include this line or not ?
764	IF ( ntry+line-ntarget.GT.ntarget-ntry ) THEN
765	! Nearer start than end. Finish at previous line.
766	ifin = i-1
767	ELSE
768	! Nearer end than start. Include this line.
769	ifin = i
770	ntry = ntry + line
771	ENDIF
772
773	IF( PARTIT_DEBUG )THEN
774	WRITE (,) 'Reached target at ',ifin &
775	,' ntry = ',ntry
776	END IF
777	EXIT scanrows
778	ENDIF
779	! Add in the line count to the running total.
780	ntry = ntry + line
781	ENDDO scanrows
782	IF ( ifin.LT.0 ) ifin = iub
783	ENDIF
784
785	! Set up the parameters of the new sub-domain.
786	pnactive(isub_new) = 0
787	DO j=jlb,jub
788	DO i=ist,ifin
789	IF ( mask(i,j) == 1 ) THEN
790	pnactive(isub_new) = pnactive(isub_new)+1
791	ENDIF
792	ENDDO
793	ENDDO
794	pielb(isub_new) = ist
795	pieub(isub_new) = ifin
796	pjelb(isub_new) = jlb
797	pjeub(isub_new) = jub
798
799	IF( PARTIT_DEBUG )THEN
800	WRITE (,) 'New subdomain is ',ist,ifin,jlb,jub &
801	,' with ',pnactive(isub_new),' points'
802	END IF
803	! Reset the target based on the actual number of points
804	! remaining, split between the partitions remaining.
805	ngone = ngone+ntry
806	! Start the next search at the next point.
807	ist = ifin+1
808	ENDDO
809	ENDDO
810
811	ELSE
812
813	! ============================================================
814	! Partition in y.
815	! ============================================================
816	ify = ify+1
817	f = fy(ify)
818	nsub_new = nsub_old*f
819
820	IF( PARTIT_DEBUG )THEN
821	WRITE (,) 'Partitioning in y from ',nsub_old &
822	,' to ',nsub_new
823	END IF
824
825	DO isub_old=1, nproc, nproc/nsub_old
826	! Check that there are sufficient points to factorise.
827	IF ( pjeub(isub_old)-pjelb(isub_old)+1.LT.f ) THEN
828	WRITE (,) &
829	'partition_rk: insufficient points to partition'
830	ierr = 1
831	RETURN
832	ENDIF
833	! Set the target number of points in the new sub-domains.
834	ntarget = NINT(REAL(pnactive(isub_old))/REAL(f))
835	ilb = pielb(isub_old)
836	iub = pieub(isub_old)
837	jst = pjelb(isub_old)
838	jub = pjeub(isub_old)
839
840	IF( PARTIT_DEBUG )THEN
841	WRITE (,)
842	WRITE (,) '======================================='
843	WRITE (,)
844	WRITE (,) 'Partitioning sub-domain ',isub_old,' : ' &
845	,ilb,iub,jst,jub
846	WRITE (,) 'Old partition has ',pnactive(isub_old) &
847	,' points'
848	WRITE (,) 'Target is ',ntarget
849	END IF
850	! Create the new sub-domains.
851	ngone = 0
852	DO isub=1,f
853	isub_new = isub_old + (isub-1)*nproc/nsub_new
854
855	IF( PARTIT_DEBUG )THEN
856	WRITE (,)
857	WRITE (,) 'Creating new sub-domain ',isub_new
858	END IF
859
860	! The last new sub-domain takes the remaining data.
861	IF ( isub == f ) THEN
862	jfin = jub
863	ELSE
864	! Scan lines in y counting active grid points until we
865	! exceed the target.
866	ntry = 0
867	jfin = -1
868	scancols: DO j=jst,jub
869	! Count up the contribution from the next line.
870	line = 0
871	DO i=ilb,iub
872	IF ( mask(i,j) == 1 ) line = line+1
873	ENDDO
874
875	! Check against the target.
876	IF ( ntry+line > ntarget ) THEN
877	! Is it nearer the target to include this line or not ?
878	IF ( ntry+line-ntarget > ntarget-ntry ) THEN
879	! Nearer start than end. Finish at previous line.
880	jfin = j-1
881	ELSE
882	! Nearer end than start. Include this line.
883	jfin = j
884	ntry = ntry + line
885	ENDIF
886
887	IF( PARTIT_DEBUG )THEN
888	WRITE (,) 'Reached target at ',jfin &
889	,' ntry = ',ntry
890	END IF
891
892	EXIT scancols
893	ENDIF
894	! Add in the line count to the running total.
895	ntry = ntry + line
896	ENDDO scancols
897	IF ( jfin.LT.0 ) jfin = jub
898	ENDIF
899	! Set up the parameters of the new sub-domain.
900	pnactive(isub_new) = 0
901	DO j=jst,jfin
902	DO i=ilb,iub
903	IF ( mask(i,j) == 1 ) THEN
904	pnactive(isub_new) = pnactive(isub_new)+1
905	ENDIF
906	ENDDO
907	ENDDO
908	pielb(isub_new) = ilb
909	pieub(isub_new) = iub
910	pjelb(isub_new) = jst
911	pjeub(isub_new) = jfin
912
913	IF( PARTIT_DEBUG )THEN
914	WRITE (,) 'New subdomain is ',ilb,iub,jst,jfin &
915	,' with ',pnactive(isub_new),' points'
916	END IF
917
918	! Reset the target based on the actual number of points
919	! remaining, split between the partitions remaining.
920	ngone = ngone+ntry
921
922	IF( PARTIT_DEBUG ) WRITE (,) 'Target reset to ',ntarget
923
924	! Start the next search at the next point.
925	jst = jfin+1
926	ENDDO
927	ENDDO
928	ENDIF
929	! Update the number of sub-domains ready for the next loop.
930	nsub_old = nsub_new
931	ENDDO
932
933	IF( PARTIT_DEBUG )THEN
934	WRITE (,)
935	WRITE (,) '========================================='
936	WRITE (,)
937	END IF
938
939	! Set the size of each sub-domain.
940	DO iproc=1, nproc, 1
941	piesub(iproc) = pieub(iproc)-pielb(iproc)+1
942	pjesub(iproc) = pjeub(iproc)-pjelb(iproc)+1
943
944	IF( PARTIT_DEBUG )THEN
945	WRITE (,) 'Domain ',iproc,' has ',pnactive(iproc), ' ocean points'
946	WRITE(*,"(I3,': [',I3,':',I3,'][',I3,':',I3,']')") &
947	iproc, pielb(iproc), pieub(iproc), pjelb(iproc), pjeub(iproc)
948	END IF
949	ENDDO
950
951	END SUBROUTINE partition_rk_core
952
953
954	SUBROUTINE factor ( ifax, maxfax, nfax, n, ierr )
955
956	!!------------------------------------------------------------------
957	! Subroutine to return the prime factors of n.
958	! nfax factors are returned in array ifax which is of maximum
959	! dimension maxfax.
960	!!------------------------------------------------------------------
961	IMPLICIT NONE
962	! Subroutine arguments
963	INTEGER ierr, n, nfax, maxfax
964	INTEGER ifax(maxfax)
965	! Local variables.
966	INTEGER i, ifac, l, nu
967	INTEGER lfax(20)
968	! lfax contains the set of allowed factors.
969	DATA (lfax(i),i=1,20) / 71,67,59,53,47,43,41,37,31,29 &
970	,23,19,17,13,11, 7, 5, 3, 2, 1 /
971	! Clear the error flag.
972	ierr = 0
973	! Find the factors of n.
974	IF ( n.EQ.1 ) THEN
975	nfax = 0
976	GOTO 20
977	ENDIF
978	! nu holds the unfactorised part of the number.
979	! nfax holds the number of factors found.
980	! l points to the allowed factor list.
981	! ifac holds the current factor.
982	nu = n
983	nfax = 0
984	l = 1
985	ifac = lfax(l)
986	! Label 10 is the start of the factor search loop.
987	10 CONTINUE
988	! Test whether the factor will divide.
989	IF ( MOD(nu,ifac).EQ.0 ) THEN
990	! Add the factor to the list.
991	nfax = nfax+1
992	IF ( nfax.GT.maxfax ) THEN
993	ierr = 6
994	WRITE (,) &
995	'FACTOR: insufficient space in factor array ',nfax
996	RETURN
997	END IF
998	ifax(nfax) = ifac
999	! Divide out the factor from the remaining number.
1000	! If unity remains, the number has been
1001	! successfully factorized.
1002	nu = nu/ifac
1003	IF ( nu.EQ.1 ) GOTO 20
1004	! Loop to try the factor again.
1005	GOTO 10
1006	ENDIF
1007	! Move on to the next factor in the list.
1008	l = l+1
1009	ifac = lfax(l)
1010	! Unless the end of the factor list has been reached, loop.
1011	IF ( ifac.GT.1 ) GOTO 10
1012	! There is a factor in n which is not in the lfac list.
1013	! Add the remaining number onto the end of the list.
1014	nfax = nfax+1
1015	IF ( nfax.GT.maxfax ) THEN
1016	ierr = 6
1017	WRITE (,) 'FACTOR: insufficient space in factor array ',nfax
1018	RETURN
1019	ENDIF
1020	ifax(nfax) = nu
1021	! Label 20 is the exit point from the factor search loop.
1022	20 CONTINUE
1023
1024	END SUBROUTINE factor
1025
1026
1027	SUBROUTINE part_trim ( depth, isTrimmed, ierr )
1028	!!------------------------------------------------------------------
1029	!!
1030	!! Examines all the sub-domains and trims off boundary rows and
1031	!! columns which are all land.
1032	!!
1033	!! depth real input Depth map.
1034	!! isTrimmed logical output Whether N,E,S,W edge
1035	!! of domain has been
1036	!! trimmed
1037	!! ierr int output Error flag.
1038	!!
1039	!! Mike Ashworth, CLRC Daresbury Laboratory, July 1999
1040	!!------------------------------------------------------------------
1041	IMPLICIT NONE
1042
1043	! Subroutine arguments.
1044	INTEGER, INTENT(in) :: depth(jpiglo,jpjglo)
1045	LOGICAL, DIMENSION(:,:), INTENT(out) :: isTrimmed
1046	INTEGER, INTENT(out) :: ierr
1047	! Local variables.
1048	INTEGER :: i, iproc, j, newbound
1049
1050	! Clear the error code.
1051	ierr = 0
1052
1053	! Clear all flags
1054	! N E S W
1055	! 1 2 3 4
1056	isTrimmed(1:4,1:nprocp) = .FALSE.
1057
1058	! Examine each sub-domain in turn.
1059
1060	subdomain_loop: DO iproc=1,nprocp
1061
1062	! Do not trim if there are no active points at all.
1063	! Otherwise we will trim away the whole sub-domain and we
1064	! will be in big trouble.
1065
1066	! Look at the low-i columns (Western edge of domain)
1067
1068	newbound = pielb(iproc)
1069	lowi: DO i=pielb(iproc),pieub(iproc)
1070
1071	! Scan this column for wet points
1072
1073	DO j=MAX(1,pjelb(iproc)-jprecj),MIN(jpjglo,pjeub(iproc)+jprecj)
1074
1075	IF ( depth(i,j) == 1 ) THEN
1076	newbound = MAX(i - jpreci - nextra, pielb(iproc))
1077
1078	IF( PARTIT_DEBUG )THEN
1079	WRITE (,) 'Sub-domain ',iproc,': Low-i loop: row ',i &
1080	,' is land. New bound is ',newbound
1081	END IF
1082
1083	EXIT lowi
1084	ENDIF
1085	ENDDO
1086	ENDDO lowi
1087
1088	! Update with the new boundary and monitor.
1089
1090	IF ( newbound.NE.pielb(iproc) ) THEN
1091
1092	IF( TRIM_DEBUG )THEN
1093	IF ( newbound-pielb(iproc).GT.1 ) THEN
1094	WRITE(*,'(a,i5,3(a,i6))') ' Process ',iproc-1 &
1095	,' trimmed ',newbound-pielb(iproc) &
1096	,' cols: ',pielb(iproc),' to ',newbound-1
1097	ELSE
1098	WRITE(*,'(a,i5,3(a,i6))') ' Process ',iproc-1 &
1099	,' trimmed ',newbound-pielb(iproc) &
1100	,' col : ',pielb(iproc)
1101	ENDIF
1102	END IF
1103
1104	pielb(iproc) = newbound
1105	isTrimmed(4,iproc) = .TRUE.
1106	ENDIF
1107
1108	! Look at the high-i columns (Eastern edge of domain).
1109
1110	newbound = pieub(iproc)
1111	highi: DO i=pieub(iproc),pielb(iproc),-1
1112
1113	DO j=MAX(1,pjelb(iproc)-jprecj), MIN(jpjglo,pjeub(iproc)+jprecj)
1114
1115	IF ( depth(i,j) == 1 ) THEN
1116	! We've found a wet point in this column so this is as far
1117	! as we can trim.
1118	newbound = MIN(i + jpreci + nextra, pieub(iproc))
1119
1120	IF( PARTIT_DEBUG )THEN
1121	WRITE (*,"('Sub-domain ',I3,': High-i loop: col ',I3, &
1122	& ' contains water. New bound is ',I3)") iproc,i,newbound
1123	END IF
1124
1125	EXIT highi
1126
1127	ENDIF
1128	ENDDO
1129	ENDDO highi
1130
1131	! Update with the new boundary and monitor.
1132
1133	IF ( newbound.NE.pieub(iproc) ) THEN
1134
1135	IF( TRIM_DEBUG )THEN
1136	IF ( (pieub(iproc)-newbound) .GT.1 ) THEN
1137	WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1138	,' trimmed ',pieub(iproc)-newbound &
1139	,' cols: ',newbound+1,' to ',pieub(iproc)
1140	ELSE
1141	WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1142	,' trimmed ',pieub(iproc)-newbound &
1143	,' col : ',pieub(iproc)
1144	ENDIF
1145	END IF
1146
1147	pieub(iproc) = newbound
1148	isTrimmed(2,iproc) = .TRUE.
1149	ENDIF
1150
1151	! Look at the low-j rows (Southern edge of domain)
1152
1153	newbound = pjelb(iproc)
1154	lowj: DO j=pjelb(iproc),pjeub(iproc),1
1155
1156	! Scan this row for wet points
1157
1158	DO i=MAX(1,pielb(iproc)-jpreci),MIN(jpiglo,pieub(iproc)+jpreci)
1159	IF ( depth(i,j) == 1) THEN
1160	newbound = MAX(j - jpreci - nextra, pjelb(iproc))
1161
1162	IF( PARTIT_DEBUG )THEN
1163	WRITE (,) 'Sub-domain ',iproc,': Low-j loop: column ',j &
1164	,' is land. New bound is ',newbound
1165	END IF
1166
1167	EXIT lowj
1168	ENDIF
1169	ENDDO
1170
1171	ENDDO lowj
1172
1173
1174	! Update with the new boundary and monitor.
1175
1176	IF ( newbound.NE.pjelb(iproc) ) THEN
1177
1178	IF( TRIM_DEBUG )THEN
1179	IF ( (newbound-pjelb(iproc)) .GT.1 ) THEN
1180	WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1181	,' trimmed ',newbound-pjelb(iproc) &
1182	,' rows: ',pjelb(iproc),' to ',newbound-1
1183	ELSE
1184	WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1185	,' trimmed ',newbound-pjelb(iproc) &
1186	,' row : ',pjelb(iproc)
1187	ENDIF
1188	END IF
1189
1190	pjelb(iproc) = newbound
1191	isTrimmed(3,iproc) = .TRUE.
1192	ENDIF
1193
1194	! Look at the high-j rows (Northern edge of domain)
1195
1196	newbound = pjeub(iproc)
1197	highj: DO j=pjeub(iproc),pjelb(iproc),-1
1198
1199	! Scan this row for wet points
1200
1201	DO i=MAX(1,pielb(iproc)-jpreci),MIN(jpiglo,pieub(iproc)+jpreci)
1202	IF ( depth(i,j) == 1 ) THEN
1203	newbound = MIN(j + jpreci + nextra, pjeub(iproc))
1204
1205	IF( TRIM_DEBUG )THEN
1206	WRITE(,) 'Sub-domain ',iproc,': High-j loop: column ',j &
1207	,' is land. New bound is ',newbound
1208	END IF
1209
1210	EXIT highj
1211	ENDIF
1212	ENDDO
1213	ENDDO highj
1214
1215	! Update with the new boundary and monitor.
1216
1217	IF ( newbound.NE.pjeub(iproc) ) THEN
1218
1219	IF( TRIM_DEBUG )THEN
1220	IF ( pjeub(iproc)-newbound.GT.1 ) THEN
1221	WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1222	,' trimmed ',pjeub(iproc)-newbound &
1223	,' rows: ',newbound+1,' to ',pjeub(iproc)
1224	ELSE
1225	WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
1226	,' trimmed ',pjeub(iproc)-newbound &
1227	,' row : ',pjeub(iproc)
1228	ENDIF
1229	ENDIF
1230	pjeub(iproc) = newbound
1231	isTrimmed(1,iproc) = .TRUE.
1232	ENDIF
1233
1234	! Reset the size of the sub-domain.
1235
1236	piesub(iproc) = pieub(iproc)-pielb(iproc)+1
1237	pjesub(iproc) = pjeub(iproc)-pjelb(iproc)+1
1238
1239	! endif active_subdomain
1240
1241	ENDDO subdomain_loop
1242
1243	END SUBROUTINE part_trim
1244
1245
1246	SUBROUTINE finish_partition(imask)
1247	IMPLICIT NONE
1248	INTEGER, DIMENSION(:,:), INTENT(in) :: imask
1249	! Locals
1250	INTEGER :: iproc, ierr
1251
1252	! Set the external boundary flags before boundaries are
1253	! altered by the trimming process and it becomes more difficult
1254	! to recognize which were the external boundaries.
1255
1256	DO iproc=1, nprocp, 1
1257	pilbext(iproc) = pielb(iproc).EQ.1
1258	piubext(iproc) = pieub(iproc).EQ.jpiglo
1259	pjlbext(iproc) = pjelb(iproc).EQ.1
1260	pjubext(iproc) = pjeub(iproc).EQ.jpjglo
1261	ENDDO
1262
1263	! Trim off redundant rows and columns containing all land.
1264
1265	IF ( trim_land ) THEN
1266	CALL part_trim ( imask, trimmed, ierr )
1267	ELSE
1268	trimmed(1:4,1:nprocp) = .FALSE.
1269	ENDIF
1270
1271	!!$
1272	!!$ ! Lower boundary (long.) of sub-domain, GLOBAL coords
1273	!!$ ! before correction for global halos
1274	!!$ nimpp = pielb(narea)
1275	!!$
1276	!!$ ! Is the domain on an external LONGITUDE boundary?
1277	!!$ nbondi = 0
1278	!!$ ilbext = pilbext(narea)
1279	!!$ IF(ilbext)THEN
1280	!!$ nbondi = -1
1281	!!$ END IF
1282	!!$
1283	!!$ IF( (.NOT. ilbext) .OR. (ilbext .AND. trimmed(widx,narea)) )THEN
1284	!!$ ! It isn't, which means we must shift its lower boundary by
1285	!!$ ! -jpreci to allow for the overlap of this domain with its
1286	!!$ ! westerly neighbour.
1287	!!$ nimpp = nimpp - jpreci
1288	!!$ END IF
1289	!!$
1290	!!$ iubext = piubext(narea)
1291	!!$ IF(iubext)THEN
1292	!!$ nbondi = 1
1293	!!$ IF(ilbext)nbondi = 2 ! Both East and West boundaries are at
1294	!!$ ! edges of global domain
1295	!!$ END IF
1296	!!$
1297	!!$ ! Set local values for limits in global coords of the sub-domain
1298	!!$ ! owned by this process.
1299	!!$ ielb = pielb (narea)
1300	!!$ ieub = pieub (narea)
1301	!!$ iesub = piesub(narea)
1302	!!$
1303	!!$ jpi = iesub + 2*jpreci ! jpi is the same for all domains - this is
1304	!!$ ! what original decomposition did
1305	!!$ nlci = jpi
1306	!!$
1307	!!$ ! If the domain is at the edge of the model domain and a cyclic
1308	!!$ ! East-West b.c. is in effect then it already incorporates one
1309	!!$ ! extra halo column (because of the way the model domain itself is
1310	!!$ ! set-up). If we've trimmed-off dry rows and columns then, even if
1311	!!$ ! a domain is on the model boundary, it may still need a halo so
1312	!!$ ! we add one.
1313	!!$ IF( nbondi == -1 .AND. (.NOT. trimmed(widx,narea)) )THEN
1314	!!$ ! Western boundary
1315	!!$ ! First column of global domain is actually a halo but NEMO
1316	!!$ ! still sets nldi to 1.
1317	!!$ nldi = 1 ! Lower bnd of int. region of sub-domain, LOCAL
1318	!!$ nlei = nldi + iesub - 1
1319	!!$ nlci = nlei + jpreci
1320	!!$ jpi = nlci
1321	!!$
1322	!!$ ELSE IF( nbondi == 1 .AND. (.NOT. trimmed(eidx,narea)) )THEN
1323	!!$ ! Eastern boundary
1324	!!$ nldi = 1 + jpreci
1325	!!$ ! Last column of global domain is actually a halo
1326	!!$ nlei = nldi + iesub - 1
1327	!!$ nlci = nlei
1328	!!$ jpi = nlei
1329	!!$
1330	!!$ ELSE IF( nbondi == 2)THEN
1331	!!$ ! Both boundaries are on the edges of the global domain
1332	!!$ IF(trimmed(widx,narea))THEN
1333	!!$ nldi = 1 + jpreci
1334	!!$ ELSE
1335	!!$ nldi = 1
1336	!!$ END IF
1337	!!$ nlei = nldi + iesub - 1
1338	!!$
1339	!!$ IF(trimmed(eidx,narea))THEN
1340	!!$ nlci = nlei + jpreci
1341	!!$ ELSE
1342	!!$ nlci = nlei
1343	!!$ END IF
1344	!!$ jpi = nlci
1345	!!$
1346	!!$ ELSE
1347	!!$ ! We have no external boundaries to worry about
1348	!!$ nldi = 1 + jpreci
1349	!!$ nlei = nldi + iesub - 1 !
1350	!!$ END IF
1351	!!$
1352	!!$ jpim1 = jpi - 1
1353	!!$
1354	!!$ ! Lower ext. boundary (lat.) of sub-domain, global coords
1355	!!$ njmpp= pjelb (narea)
1356	!!$
1357	!!$ ! Is the domain on an external LATITUDE boundary?
1358	!!$ nbondj = 0
1359	!!$ jlbext = pjlbext(narea)
1360	!!$ IF(jlbext)THEN
1361	!!$ nbondj = -1
1362	!!$ ENDIF
1363	!!$
1364	!!$ IF((.NOT. jlbext) .OR. (jlbext .AND. trimmed(sidx,narea)) )THEN
1365	!!$ ! It isn't, which means we must shift its lower boundary by
1366	!!$ ! -jprecj to allow for the overlap of this domain with its
1367	!!$ ! southerly neighbour.
1368	!!$ njmpp = njmpp - jprecj
1369	!!$ END IF
1370	!!$ jubext = pjubext(narea)
1371	!!$ IF(jubext)THEN
1372	!!$ nbondj = 1
1373	!!$ IF(jlbext)nbondj = 2
1374	!!$ END IF
1375	!!$
1376	!!$ jelb = pjelb (narea) ! Lower bound of internal domain
1377	!!$ jeub = pjeub (narea) ! Upper bound of internal domain
1378	!!$ jesub = pjesub(narea) ! Extent of internal domain
1379	!!$
1380	!!$ jpj = jesub + 2*jprecj ! jpj is the same for all domains - this is
1381	!!$ ! what original decomposition did
1382	!!$ nlcj = jpj
1383	!!$
1384	!!$ ! Unlike the East-West boundaries, the global domain does not include
1385	!!$ ! halo (rows) at the Northern and Southern edges. In fact, there is no
1386	!!$ ! cyclic boundary condition in the North-South direction so there are no
1387	!!$ ! halos at all at the edges of the global domain.
1388	!!$ IF( nbondj == -1 .AND. (.NOT. trimmed(sidx,narea)) )THEN
1389	!!$ ! Southern edge
1390	!!$ nldj = 1 ! Start of internal region (local coords)
1391	!!$ nlej = nldj + jesub - 1 ! Upper bound of int. region of sub-domain, local
1392	!!$ nlcj = nlej + jprecj
1393	!!$ jpj = nlcj
1394	!!$ ELSE IF( nbondj == 1 .AND. (.NOT. trimmed(nidx,narea)) )THEN
1395	!!$ ! Northern edge
1396	!!$ nldj = 1+jprecj ! Start of internal region (local coords)
1397	!!$ nlej = nldj + jesub - 1
1398	!!$ nlcj = nlej
1399	!!$ jpj = nlcj
1400	!!$ jpj = jpj + 1 ! Add one extra row for zero BC along N edge as
1401	!!$ ! happens in orig. code when MOD(jpjglo,2)!=0
1402	!!$ ! Many loops go up to j=jpjm1 so unless jpj>nlej
1403	!!$ ! they won't cover the whole domain.
1404	!!$ ELSE IF( nbondj == 2)THEN
1405	!!$ ! Both boundaries are on the edges of the global domain
1406	!!$
1407	!!$ IF( trimmed(sidx,narea) )THEN
1408	!!$ nldj = 1+jprecj
1409	!!$ ELSE
1410	!!$ nldj = 1
1411	!!$ END IF
1412	!!$ nlej = nldj + jesub - 1
1413	!!$
1414	!!$ IF( trimmed(nidx,narea) )THEN
1415	!!$ nlcj = nlej + jprecj
1416	!!$ jpj = nlcj
1417	!!$ ELSE
1418	!!$ nlcj = nlej
1419	!!$ jpj = nlcj
1420	!!$ END IF
1421	!!$ jpj = jpj + 1 ! Add one extra row for zero BC along N edge as
1422	!!$ ! happens in orig. code when MOD(jpjglo,2)!=0
1423	!!$ ELSE
1424	!!$ ! We have no external boundaries to worry about
1425	!!$ nldj = 1+jprecj ! Lower bound of int. region of sub-domain, local
1426	!!$ nlej = nldj + jesub - 1
1427	!!$ END IF
1428	!!$
1429	!!$ jpjm1 = jpj-1
1430	!!$ jpij = jpi*jpj
1431
1432	END SUBROUTINE finish_partition
1433
1434
1435	SUBROUTINE eval_partition( cores_per_node, nx, ny, mask, score )
1436
1437	! Compute the cost function for the current partition
1438	!
1439	! Assume that the time taken for a run is proportional
1440	! to the maximum over processors of:
1441	! w_processing * cost_processing
1442	! + w_communications * cost_communications
1443	! Assume further that cost_processing goes as
1444	! (number of wet points) + f_proc * (number of dry points)
1445	! (with f_proc << 1)
1446	! and that cost_communications goes as
1447	! (cost of intra-node communications) +
1448	! f_comm * (cost of inter-node communications)
1449	! (with f_comm << 1)
1450	!
1451	! However, because of the possiblity of network contention,
1452	! other factors may also matter, especially:
1453	! total over sub-domains of halo points with off-node neighbours
1454	! max over nodes of total off-node halo points and message counts
1455	!
1456	! With this in mind, we construct the ansatz
1457	! maximum over processors of {
1458	! w_1 * (number of wet points)
1459	! + w_2 * (number of dry points)
1460	! + w_3 * (halo points with off-node neighbours)
1461	! + w_4 * (halo points with on-node neighbours)
1462	! + ...
1463	! }
1464	IMPLICIT NONE
1465	! Arguments
1466	INTEGER, INTENT(in) :: cores_per_node
1467	INTEGER, INTENT(in) :: nx, ny
1468	INTEGER, INTENT(in) :: mask(nx,ny)
1469	REAL(wp), DIMENSION(pv_num_scores), INTENT(out) :: score
1470	! Local variables
1471	INTEGER :: iproc, inode, i, j
1472
1473	REAL(wp) :: score_wet, score_dry
1474	REAL(wp) :: score_pli, score_plx
1475	REAL(wp) :: score_pci, score_pcx
1476	REAL(wp) :: score_nli, score_nlx
1477	REAL(wp) :: score_nci, score_ncx
1478	REAL(wp) :: score_tli, score_tlx
1479	REAL(wp) :: score_tci, score_tcx
1480
1481	REAL(wp) :: score_too_narrow
1482
1483	REAL(wp) :: proc_overall_score
1484	REAL(wp) :: proc_comm_score
1485	REAL(wp) :: node_comm_score
1486
1487	REAL(wp), PARAMETER :: weight_wet = 1.00D0
1488	REAL(wp), PARAMETER :: weight_dry = 0.9D0
1489
1490	REAL(wp), PARAMETER :: weight_pli = 0.01D0
1491	REAL(wp), PARAMETER :: weight_plx = 0.20D0
1492	REAL(wp), PARAMETER :: weight_pci = 0.50D0
1493	REAL(wp), PARAMETER :: weight_pcx = 10.00D0
1494
1495	REAL(wp), PARAMETER :: weight_nli = 0.00D0
1496	REAL(wp), PARAMETER :: weight_nlx = 0.20D0
1497	REAL(wp), PARAMETER :: weight_nci = 0.00D0
1498	REAL(wp), PARAMETER :: weight_ncx = 10.00D0
1499
1500	REAL(wp), PARAMETER :: weight_tli = 0.00D0
1501	REAL(wp), PARAMETER :: weight_tlx = 0.01D0
1502	REAL(wp), PARAMETER :: weight_tci = 0.00D0
1503	REAL(wp), PARAMETER :: weight_tcx = 0.50D0
1504
1505	INTEGER :: peer, last_peer
1506
1507	! Which node is each process on?
1508	! Assumes that first cores_per_node ranks are assigned to node 0,
1509	! next cores_per_node ranks are assigned to node 1, etc
1510	INTEGER, ALLOCATABLE :: node(:)
1511
1512	ALLOCATE(node(nprocp))
1513
1514	DO iproc=1, nprocp
1515	node(iproc) = (iproc-1)/cores_per_node
1516	END DO
1517
1518	! Calculate maximum per processor score
1519
1520	score(:) = 0.0_wp
1521
1522	! Total (over all processors) off node comms
1523	score_tli = 0.0_wp
1524	score_tlx = 0.0_wp
1525	score_tci = 0.0_wp
1526	score_tcx = 0.0_wp
1527
1528	! Set this to pv_awful if any sub-domain is too narrow.
1529	score_too_narrow = 0.0_wp
1530
1531	! loop over nodes in job, counting from 0
1532	node_loop: DO inode=0, (nprocp-1)/cores_per_node
1533
1534	score_nli = 0.0_wp
1535	score_nlx = 0.0_wp
1536	score_nci = 0.0_wp
1537	score_ncx = 0.0_wp
1538
1539	! loop over processes in the node
1540	proc_loop: DO iproc=1+inode*cores_per_node, &
1541	MIN(nprocp,(inode+1)*cores_per_node)
1542
1543	score_wet = DBLE(pnactive(iproc))
1544	score_dry = DBLE(piesub(iproc)*pjesub(iproc)-score_wet)
1545
1546	score_pli = 0.0_wp
1547	score_plx = 0.0_wp
1548	score_pci = 0.0_wp
1549	score_pcx = 0.0_wp
1550
1551	! Things sometimes go wrong when a sub-domain has very
1552	! narrow partitions (2 grid points or less).
1553	! Prevent such problematic partitions from being selected.
1554	IF ( ((pieub(iproc)-pielb(iproc)) .LE. 2) &
1555	.OR. ((pjeub(iproc)-pjelb(iproc)) .LE. 2) ) THEN
1556	score_too_narrow = pv_awful
1557	END IF
1558
1559	IF (.NOT. pjlbext(iproc)) THEN
1560	j=pjelb(iproc)
1561	IF (j .GT. 1) THEN
1562	last_peer = -1
1563	DO i=pielb(iproc),pieub(iproc)
1564	IF ( (mask(i,j) .NE. land) &
1565	.AND. (mask(i,j-1) .NE. land)) THEN
1566	peer=iprocmap(i,j-1)
1567	! Total points involved in messages
1568	IF (node(peer) .EQ. inode) THEN
1569	score_pli = score_pli+1.0_wp
1570	ELSE
1571	score_plx = score_plx+1.0_wp
1572	END IF
1573	! Total number of messages
1574	IF (peer .NE. last_peer) THEN
1575	last_peer = peer
1576	IF (node(peer) .EQ. inode) THEN
1577	score_pci = score_pci+1.0_wp
1578	ELSE
1579	score_pcx = score_pcx+1.0_wp
1580	END IF
1581	END IF
1582	END IF
1583	END DO
1584	END IF
1585	END IF
1586
1587	IF (.NOT. pjubext(iproc)) THEN
1588	j=pjeub(iproc)
1589	IF (j .LT. ny) THEN
1590	last_peer = -1
1591	DO i=pielb(iproc),pieub(iproc)
1592	IF ( (mask(i,j) .NE. land) &
1593	.AND. (mask(i,j+1) .NE. land)) THEN
1594	peer=iprocmap(i,j+1)
1595	! Total points involved in messages
1596	IF (node(peer) .EQ. inode) THEN
1597	score_pli = score_pli+1.0_wp
1598	ELSE
1599	score_plx = score_plx+1.0_wp
1600	END IF
1601	! Total number of messages
1602	IF (peer .NE. last_peer) THEN
1603	last_peer = peer
1604	IF (node(peer) .EQ. inode) THEN
1605	score_pci = score_pci+1.0_wp
1606	ELSE
1607	score_pcx = score_pcx+1.0_wp
1608	END IF
1609	END IF
1610	END IF
1611	END DO
1612	END IF
1613	END IF
1614
1615	IF (.NOT. pilbext(iproc)) THEN
1616	i=pielb(iproc)
1617	IF (i .GT. 1) THEN
1618	last_peer = -1
1619	DO j=pjelb(iproc),pjeub(iproc)
1620	IF ( (mask(i,j) .NE. land) &
1621	.AND. (mask(i-1,j) .NE. land)) THEN
1622	peer=iprocmap(i-1,j)
1623	! Total points involved in messages
1624	IF (node(peer) .EQ. inode) THEN
1625	score_pli = score_pli+1.0_wp
1626	ELSE
1627	score_plx = score_plx+1.0_wp
1628	END IF
1629	! Total number of messages
1630	IF (peer .NE. last_peer) THEN
1631	last_peer = peer
1632	IF (node(peer) .EQ. inode) THEN
1633	score_pci = score_pci+1.0_wp
1634	ELSE
1635	score_pcx = score_pcx+1.0_wp
1636	END IF
1637	END IF
1638	END IF
1639	END DO
1640	END IF
1641	END IF
1642
1643	IF (.NOT. piubext(iproc)) THEN
1644	i=pieub(iproc)
1645	IF (i .LT. nx) THEN
1646	last_peer = -1
1647	DO j=pjelb(iproc),pjeub(iproc)
1648	IF ( (mask(i,j) .NE. land) &
1649	.AND. (mask(i+1,j) .NE. land)) THEN
1650	peer=iprocmap(i+1,j)
1651	! Total points involved in messages
1652	IF (node(peer) .EQ. inode) THEN
1653	score_pli = score_pli+1.0_wp
1654	ELSE
1655	score_plx = score_plx+1.0_wp
1656	END IF
1657	! Total number of messages
1658	IF (peer .NE. last_peer) THEN
1659	last_peer = peer
1660	IF (node(peer) .EQ. inode) THEN
1661	score_pci = score_pci+1.0_wp
1662	ELSE
1663	score_pcx = score_pcx+1.0_wp
1664	END IF
1665	END IF
1666	END IF
1667	END DO
1668	END IF
1669	END IF
1670
1671	score_nli = score_nli + score_pli
1672	score_nlx = score_nlx + score_plx
1673	score_nci = score_nci + score_pci
1674	score_ncx = score_ncx + score_pcx
1675
1676	proc_overall_score = weight_wet*score_wet &
1677	+ weight_dry*score_dry &
1678	+ weight_pli*score_pli &
1679	+ weight_plx*score_plx &
1680	+ weight_pci*score_pci &
1681	+ weight_pcx*score_pcx
1682
1683	proc_comm_score = weight_pli*score_pli &
1684	+ weight_plx*score_plx &
1685	+ weight_pci*score_pci &
1686	+ weight_pcx*score_pcx
1687
1688	IF (score(pv_index_overall) < proc_overall_score) &
1689	score(pv_index_overall) = proc_overall_score
1690	IF (score(pv_index_pcomms) < proc_comm_score) &
1691	score(pv_index_pcomms) = proc_comm_score
1692	IF (score(pv_index_wet) < score_wet) &
1693	score(pv_index_wet) = score_wet
1694	IF (score(pv_index_dry) < score_dry) &
1695	score(pv_index_dry) = score_dry
1696	IF (score(pv_index_pli) < score_pli) &
1697	score(pv_index_pli) = score_pli
1698	IF (score(pv_index_plx) < score_plx) &
1699	score(pv_index_plx) = score_plx
1700	IF (score(pv_index_pci) < score_pci) &
1701	score(pv_index_pci) = score_pci
1702	IF (score(pv_index_pcx) < score_pcx) &
1703	score(pv_index_pcx) = score_pcx
1704
1705	END DO proc_loop
1706
1707	score_tli = score_tli + score_nli
1708	score_tlx = score_tlx + score_nlx
1709	score_tci = score_tci + score_nci
1710	score_tcx = score_tcx + score_ncx
1711
1712	node_comm_score = weight_nli*score_nli &
1713	+ weight_nlx*score_nlx &
1714	+ weight_nci*score_nci &
1715	+ weight_ncx*score_ncx
1716
1717	IF (score(pv_index_ncomms) < node_comm_score) &
1718	score(pv_index_ncomms) = node_comm_score
1719	IF (score(pv_index_nli) < score_nli) &
1720	score(pv_index_nli) = score_nli
1721	IF (score(pv_index_nlx) < score_nlx) &
1722	score(pv_index_nlx) = score_nlx
1723	IF (score(pv_index_nci) < score_nci) &
1724	score(pv_index_nci) = score_nci
1725	IF (score(pv_index_ncx) < score_ncx) &
1726	score(pv_index_ncx) = score_ncx
1727
1728	END DO node_loop
1729
1730	score(pv_index_tcomms) = weight_tli*score_tli &
1731	+ weight_tlx*score_tlx &
1732	+ weight_tci*score_tci &
1733	+ weight_tcx*score_tcx
1734
1735	score(pv_index_tli) = score_tli
1736	score(pv_index_tlx) = score_tlx
1737	score(pv_index_tci) = score_tci
1738	score(pv_index_tcx) = score_tcx
1739
1740	score(pv_index_overall) = score(pv_index_overall) &
1741	+ score(pv_index_ncomms) &
1742	+ score(pv_index_tcomms) &
1743	+ score_too_narrow
1744
1745	DEALLOCATE(node)
1746
1747	END SUBROUTINE eval_partition
1748
1749
1750	SUBROUTINE calc_perms( nfactors, factors, &
1751	ndf, df, multiplicity, &
1752	nperms )
1753	IMPLICIT NONE
1754
1755	! Subroutine arguments
1756	! Number of factors (including repetitions)
1757	INTEGER, INTENT(in) :: nfactors
1758	! Factors (in descending order)
1759	INTEGER, INTENT(in) :: factors(nfactors)
1760	! Number of distinct factors
1761	INTEGER, INTENT(out) :: ndf
1762	! Distinct factors (in ascending order)
1763	INTEGER :: df(nfactors)
1764	! Multiplicity of each distinct factor
1765	INTEGER :: multiplicity(nfactors)
1766	! Number of distinct permutations
1767	INTEGER, INTENT(out) :: nperms
1768
1769	! Local variables
1770	INTEGER :: product
1771	INTEGER :: i, j
1772
1773	product = factors(nfactors)
1774	ndf = 1
1775	df(:) = 0
1776	df(1) = factors(nfactors)
1777	multiplicity(:) = 0
1778	multiplicity(1) = 1
1779
1780	DO i=nfactors-1,1,-1
1781	product = product*factors(i)
1782	IF (factors(i) .EQ. df(ndf)) THEN
1783	multiplicity(ndf) = multiplicity(ndf)+1
1784	ELSE
1785	ndf = ndf+1
1786	df(ndf) = factors(i)
1787	multiplicity(ndf) = 1
1788	END IF
1789	END DO
1790
1791	! A careful code would try to avoid overflow here
1792	nperms = 1
1793	DO i=1, nfactors
1794	nperms = nperms*i
1795	END DO
1796	DO i=1, ndf
1797	DO j=1, multiplicity(i)
1798	nperms = nperms/j
1799	END DO
1800	END DO
1801
1802	! At this point, nperms is the number of distinct permutations
1803	! of the factors provided. But each of these permutations can
1804	! be split between x and y in (nfactors+1) ways, e.g.
1805	! x=(2,2,3), y=()
1806	! x=(2,3), y=(2)
1807	! x=(3), y=(2,2)
1808	! x=(), y=(2,2,3)
1809
1810	nperms = nperms*(nfactors+1)
1811
1812	WRITE (*,FMT="('Distinct factorisations: ', I4/)") nperms
1813
1814	END SUBROUTINE calc_perms
1815
1816
1817
1818	SUBROUTINE get_perm_factors( iperm, nf, ndf, df, m, &
1819	fx, nfx, fy, nfy, &
1820	prodx, prody )
1821	IMPLICIT NONE
1822	!!------------------------------------------------------------------
1823	! Our goal is to visit a particular permutation,
1824	! number perm ( 0 <= perm <= nperms-1 )
1825	! of nf things, of ndf distinct values (actually the prime
1826	! factors of number of MPI tasks), each of which can be repeated
1827	! with multiplicity m_i
1828	! assert nf = sum_{i=1..ndf} m(i)
1829	!
1830	! We don't care about lexicographic ordering, but we do
1831	! need to ensure that we don't visit any permutation twice,
1832	! in a loop over 0..nperms-1.
1833	! Textbook methods typically assume that all the things being
1834	! permuted are distinct.
1835	!
1836	! We use what I call a nested variable radix method.
1837	!
1838	! Stephen Pickles, STFC
1839	! Taken from POLCOMS code by Andrew Porter.
1840	!!------------------------------------------------------------------
1841	! Arguments
1842	INTEGER, INTENT(in) :: iperm, nf, ndf
1843	INTEGER, INTENT(in), DIMENSION(ndf) :: df, m
1844	INTEGER, INTENT(out), DIMENSION(nf) :: fx, fy
1845	INTEGER, INTENT(out) :: nfx, nfy
1846	INTEGER, INTENT(out) :: prodx, prody
1847	!
1848	! Local variables
1849	!
1850	INTEGER :: perm, split
1851	INTEGER, DIMENSION(nf) :: bin, a
1852	INTEGER, DIMENSION(ndf) :: ways
1853	INTEGER, DIMENSION(0:ndf) :: root, representation
1854	INTEGER :: i, j, k, v, p, q, r
1855	INTEGER :: unfilled, pm
1856	INTEGER, SAVE :: myinst = 0
1857	LOGICAL, PARAMETER :: debug=.FALSE.
1858	!!------------------------------------------------------------------
1859
1860	! This is a serial version so instead of using MPI rank
1861	! use the count of how many times we've been called
1862	myinst = myinst + 1
1863
1864	perm = iperm / (nf+1)
1865	! Where to split between x and y
1866	split = MOD(iperm, (nf+1))
1867
1868	! interpret ways(i) is the number of ways of distributing
1869	! m(i) copies of the i'th distinct factor into the remaining
1870	! bins
1871	k = nf
1872	DO i=1,ndf
1873	ways(i) = k
1874	DO j=2,m(i)
1875	ways(i) = ways(i)*(k-j+1)/j
1876	END DO
1877	k = k-m(i)
1878	END DO
1879
1880	! compute (outer) radices
1881	! root is the variable radix basis corresponding to ways
1882	! root(ndf) is always 1
1883	root(ndf) = 1
1884	root(0:ndf-1) = ways(1:ndf)
1885	DO i=ndf-1,0,-1
1886	root(i)=root(i)*root(i+1)
1887	END DO
1888
1889	bin(:) = 0
1890	unfilled = nf
1891
1892	r = perm
1893	! Next line is valid as long as perm < nperms
1894	representation(0) = 0
1895	DO i=1, ndf
1896	p = r/root(i)
1897	r = r - p*root(i)
1898	representation(i) = p
1899
1900	! At this point, we are considering distinct ways to
1901	! distribute m(i) copies of the i'th distinct factor into
1902	! the unfilled remaining bins. We want to select the p'th one.
1903
1904	DO j=1,unfilled
1905	a(j) = j
1906	END DO
1907	q = 0
1908	find_p: DO
1909	IF (q .GE. p) EXIT find_p
1910
1911	k=m(i)
1912	k_loop: DO
1913	IF (a(k) .EQ. (unfilled - m(i) + k)) THEN
1914	k=k-1
1915	ELSE
1916	EXIT k_loop
1917	END IF
1918	END DO k_loop
1919	a(k) = a(k) + 1
1920	DO v=k+1,m(i)
1921	a(v) = a(k) + v - k
1922	END DO
1923	q = q+1
1924	END DO find_p
1925
1926	IF (debug) THEN
1927	WRITE (*,'(A3)',advance='no') 'a=('
1928	DO j=1, m(i)-1
1929	WRITE (*,'(I3,A1)',advance='no') a(j), ','
1930	END DO
1931	WRITE (*,'(I3,A1)') a(m(i)), ')'
1932	END IF
1933
1934	DO j=1, m(i)
1935	pm=a(j)-j+1
1936
1937	! put this factor in the pm'th empty bin
1938	v = 1
1939	find_bin: DO k=1, nf
1940	IF (bin(k) .EQ. 0) THEN
1941	IF (v .EQ. pm) THEN
1942	bin(k) = df(i)
1943	unfilled = unfilled-1
1944	EXIT find_bin
1945	ELSE
1946	v=v+1
1947	END IF
1948	END IF
1949	END DO find_bin
1950
1951	END DO
1952
1953	END DO
1954
1955	! We have identified the (perm)th distinct permutation,
1956	! but we still need to split the factors between x and y.
1957	fx(:) = 0
1958	prodx = 1
1959	DO i=1,split
1960	fx(i)=bin(i)
1961	prodx=prodx*fx(i)
1962	END DO
1963	nfx=split
1964
1965	fy(:) = 0
1966	prody = 1
1967	j=1
1968	DO i=split+1,nf
1969	fy(j)=bin(i)
1970	prody=prody*fy(j)
1971	j=j+1
1972	END DO
1973	nfy=nf-nfx
1974
1975	WRITE (6,'(A14,I6,A1,I6,A2)',advance='no') &
1976	'factorisation[', myinst, ']', iperm, ' ('
1977	DO k=1,ndf-1
1978	WRITE (6,'(I4,A1)',advance='no') representation(k), ','
1979	END DO
1980	WRITE (6,'(I4,A1)',advance='no') representation(ndf), ')'
1981
1982	CALL print_factors(nfx,fx,nfy,fy)
1983
1984	END SUBROUTINE get_perm_factors
1985
1986
1987	SUBROUTINE print_factors(nfx,fx,nfy,fy)
1988	IMPLICIT NONE
1989	INTEGER, INTENT(in) :: nfx, nfy
1990	INTEGER, INTENT(in) :: fx(nfx), fy(nfy)
1991	INTEGER :: k
1992
1993	IF (nfx .GT. 0) THEN
1994	WRITE (*,'(A1)',advance='no') ' '
1995	DO k=1,nfx-1
1996	WRITE (*,'(I4,A1)',advance='no') fx(k), ','
1997	END DO
1998	WRITE (*,'(I4)',advance='no') fx(nfx)
1999	END IF
2000	WRITE (*,'(A1)',advance='no') ':'
2001	IF (nfy .GT. 0) THEN
2002	DO k=1,nfy-1
2003	WRITE (*,'(I4,A1)',advance='no') fy(k), ','
2004	END DO
2005	WRITE (*,'(I4)',advance='no') fy(nfy)
2006	END IF
2007	WRITE (*,'(A1)') ' '
2008
2009	END SUBROUTINE print_factors
2010
2011
2012	CHARACTER(len=16) FUNCTION num_to_string(number)
2013	INTEGER, INTENT(in) :: number
2014
2015	CHARACTER*16 :: outs
2016
2017	WRITE (outs,'(I15)') number
2018	num_to_string = ADJUSTL(outs)
2019
2020	END FUNCTION num_to_string
2021
2022
2023	SUBROUTINE factor_string(fstr,nfx,fx,nfy,fy)
2024	IMPLICIT NONE
2025	CHARACTER*256, INTENT(out) :: fstr
2026	INTEGER, INTENT(in) :: nfx, nfy
2027	INTEGER, INTENT(in) :: fx(nfx), fy(nfy)
2028	! Locals
2029	INTEGER :: k
2030
2031	fstr = ' '
2032	IF (nfx .GT. 0) THEN
2033	DO k=1,nfx-1
2034	fstr = TRIM(fstr)//TRIM(num_to_string(fx(k)))//'x'
2035	END DO
2036	fstr = TRIM(fstr)//TRIM(num_to_string(fx(nfx)))
2037	END IF
2038	fstr = TRIM(fstr)//'-'
2039	IF (nfy .GT. 0) THEN
2040	DO k=1,nfy-1
2041	fstr = TRIM(fstr)//TRIM(num_to_string(fy(k)))//'x'
2042	END DO
2043	fstr = TRIM(fstr)//TRIM(num_to_string(fy(nfy)))
2044	END IF
2045	END SUBROUTINE factor_string
2046
2047	FUNCTION iprocmap(ia, ja)
2048	!!------------------------------------------------------------------
2049	! Returns the process number (1...N) of the process whose sub-domain
2050	! contains the point with global coordinates (i,j).
2051	! If no process owns the point, returns zero.
2052
2053	! i int input global x-coordinate
2054	! j int input global y-coordinate
2055
2056	! Mike Ashworth, CLRC Daresbury Laboratory, July 1999
2057	! Andrew Porter, STFC Daresbury Laboratory, May 2008
2058	!!------------------------------------------------------------------
2059	IMPLICIT NONE
2060
2061	! Function arguments.
2062	INTEGER :: iprocmap
2063	INTEGER, INTENT(in) :: ia, ja
2064	! Local variables.
2065	INTEGER :: iproc, i, j, iwidth
2066
2067	iprocmap = 0
2068
2069	! Make sure we don't change variable values in calling
2070	! routine...
2071	i = ia
2072	j = ja
2073	IF(cyclic_bc)THEN
2074	! Allow for fact that first and last columns in global domain
2075	! are actually halos
2076	iwidth = jpiglo - 2*jpreci
2077	IF(i >= jpiglo) i = ia - iwidth
2078	IF(i <= 1 ) i = ia + iwidth
2079	! No cyclic BCs in North/South direction
2080	!IF(j > jpjglo) j = ja - jpjglo
2081	!IF(j < 1 ) j = ja + jpjglo
2082	END IF
2083
2084	! Search the processes for the one owning point (i,j).
2085
2086	DO iproc=1,nprocp
2087	IF ( pielb(iproc).LE.i .AND. i.LE.pieub(iproc) .AND. &
2088	pjelb(iproc).LE.j .AND. j.LE.pjeub(iproc) ) THEN
2089	iprocmap = iproc
2090	EXIT
2091	ENDIF
2092	ENDDO
2093	END FUNCTION iprocmap
2094
2095	!===================================================================
2096
2097	SUBROUTINE write_partition(imask)
2098	IMPLICIT none
2099	INTEGER, DIMENSION(:,:), INTENT(in) :: imask
2100	! Locals
2101	INTEGER, PARAMETER :: funit = 1099
2102	INTEGER :: i, j, iproc, ierr
2103	INTEGER :: idom
2104	CHARACTER(LEN=4) :: intStr
2105	! Space between plot and (each) axis
2106	INTEGER, PARAMETER :: axis_offset = 4
2107	INTEGER, PARAMETER :: tic_len = 4
2108	INTEGER, PARAMETER :: tic_label_space = 2
2109	INTEGER, PARAMETER :: xtic_incr = 50
2110	INTEGER, PARAMETER :: ytic_incr = 50
2111
2112	! Write the file that NEMO can read
2113
2114	OPEN(UNIT=funit, file='partition.dat.new', status='REPLACE', &
2115	ACTION='WRITE', IOSTAT=ierr)
2116	IF(ierr /= 0)THEN
2117	WRITE(,) 'ERROR: write_partition: failed to write partition description to file.'
2118	RETURN
2119	END IF
2120
2121	WRITE(UNIT=funit,FMT="(I5,1x,I5)") jpni, jpnj
2122
2123	DO idom = 1,nprocp, 1
2124	WRITE(UNIT=funit,FMT="('Domain: ',I5)") idom
2125	WRITE(UNIT=funit,FMT="(I5,1x,I5)") pielb(idom), pjelb(idom)
2126	WRITE(UNIT=funit,FMT="(I5,1x,I5)") pieub(idom), pjeub(idom)
2127	WRITE(UNIT=funit,FMT="(L5,1x,L5)") pilbext(idom), trimmed(widx,idom)
2128	WRITE(UNIT=funit,FMT="(L5,1x,L5)") piubext(idom), trimmed(eidx,idom)
2129	WRITE(UNIT=funit,FMT="(L5,1x,L5)") pjlbext(idom), trimmed(sidx,idom)
2130	WRITE(UNIT=funit,FMT="(L5,1x,L5)") pjubext(idom), trimmed(nidx,idom)
2131	END DO
2132
2133	CLOSE(UNIT=funit)
2134
2135	! Write out domains in postscript
2136
2137	OPEN(UNIT=997, FILE="domain_decomp.ps", &
2138	STATUS='REPLACE', ACTION='WRITE', IOSTAT=ierr)
2139
2140	IF(ierr /= 0)THEN
2141	WRITE(,) 'WARNING: failed to open domain_decomp.ps for writing'
2142	RETURN
2143	END IF
2144
2145	! Header of ps file
2146	WRITE (997,FMT='("%!PS-Adobe-1.0")')
2147	WRITE (997,FMT='("/Helvetica findfont %load the font dictionary")')
2148	WRITE (997,FMT='("12 scalefont %scale to 12pt")')
2149	WRITE (997,FMT='("setfont %make this the current font")')
2150	WRITE (997,FMT='("/u { 2 mul } def %set up a scaling factor")')
2151
2152	! Put green circles at dry points
2153	WRITE (997,FMT='("% Filled green circles at dry points")')
2154	WRITE (997,FMT='("0.1 setlinewidth")') ! Thin green line
2155	WRITE (997,FMT='("0 1 0 setrgbcolor")')
2156	WRITE (997,FMT='("newpath")')
2157	DO j = 1,jpjglo,1
2158	DO i = 1,jpiglo,1
2159	IF(imask(i,j) /= 1)THEN
2160	! Radius of each circle is unity
2161	WRITE (997,FMT='(I3," u ",I3," u 1 0 360 arc closepath")') i, j
2162	! Use 'fill' instead of 'stroke' to get filled circles
2163	WRITE (997,FMT='("fill")')
2164	END IF
2165	END DO
2166	END DO
2167
2168	! Draw the outline of the global domain in red
2169	WRITE (997,FMT='("% Draw the outline of the global domain in red")')
2170	WRITE (997,FMT='("3.0 setlinewidth")') ! Fat line of 3mm for global dom.
2171	WRITE (997,FMT='("1 0 0 setrgbcolor")')
2172	WRITE (997,FMT='("newpath")')
2173	WRITE (997,FMT='("% Cursor initialization")')
2174	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') 1,1
2175	WRITE (997,FMT='("% Drawing the rectangle")')
2176	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo,1
2177	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo,jpjglo
2178	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') 1,jpjglo
2179	WRITE (997,FMT='("closepath stroke")')
2180
2181	! Axes
2182	WRITE (997,FMT='("0 0 0 setrgbcolor")')! back to black
2183	WRITE (997,FMT='("% Draw axes")')
2184	WRITE (997,FMT='("newpath")')
2185	WRITE (997,FMT='("% y axis")')
2186	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset,1
2187	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo+axis_offset,jpjglo
2188	WRITE (997,FMT='("closepath stroke")')
2189	! Tic marks
2190	! At y==1
2191	j = 1
2192	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset,j
2193	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo+axis_offset+tic_len,j
2194	WRITE (997,FMT='("closepath stroke")')
2195	! Tic label
2196	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset+tic_len+tic_label_space,j
2197	WRITE (997,FMT='("(1) dup stringwidth pop 2 div neg 0 rmoveto show")')
2198
2199	DO j = ytic_incr, jpjglo, ytic_incr
2200	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset,j
2201	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo+axis_offset+tic_len,j
2202	WRITE (997,FMT='("closepath stroke")')
2203	! Tic label
2204	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset+tic_len+tic_label_space,j
2205	WRITE (intStr, FMT='(I4)') j
2206	WRITE (997,FMT='("(",(A),") show")') TRIM(ADJUSTL(intStr))
2207	END DO
2208
2209	WRITE (997,FMT='("% x axis")')
2210	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') 1,jpjglo+axis_offset
2211	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo,jpjglo+axis_offset
2212	WRITE (997,FMT='("closepath stroke")')
2213	WRITE (997,FMT='("closepath stroke")')
2214	! Tic marks
2215	! For x==1
2216	i = 1
2217	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') i,jpjglo+axis_offset
2218	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') i,jpjglo+axis_offset+tic_len
2219	WRITE (997,FMT='("closepath stroke")')
2220	! Tic label
2221	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') i,jpjglo+axis_offset+tic_len+tic_label_space
2222	WRITE (997,FMT='("(1) dup stringwidth pop 2 div neg 0 rmoveto show")')
2223
2224	DO i = xtic_incr, jpiglo, xtic_incr
2225	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') i,jpjglo+axis_offset
2226	WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') i,jpjglo+axis_offset+tic_len
2227	WRITE (997,FMT='("closepath stroke")')
2228	! Tic label
2229	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') i,jpjglo+axis_offset+tic_len+tic_label_space
2230	WRITE (intStr, FMT='(I4)') i
2231	WRITE (997,FMT='("(",(A),") dup stringwidth pop 2 div neg 0 rmoveto show")') TRIM(ADJUSTL(intStr))
2232	END DO
2233
2234	! Now draw the outline of each individual domain, nprocp should have been
2235	! set at the very beginning of the recursive partioning process.
2236	WRITE (997,FMT='("0.7 setlinewidth")') ! Back to default width
2237	WRITE (997,FMT='("0 0 0 setrgbcolor")')! and colour
2238	DO iproc=1,nprocp
2239	WRITE (997,FMT='("newpath")')
2240	WRITE (997,FMT='("% Cursor initialization")')
2241	WRITE (997,FMT='(I3," u ",1x,I3," u moveto %BL Corner")') pielb(iproc),pjelb(iproc)
2242	WRITE (997,FMT='("% Drawing the rectangle")')
2243	WRITE (997,FMT='(I3," u ",1x,I3," u lineto %BR Corner")') pieub(iproc),pjelb(iproc)
2244	WRITE (997,FMT='(I3," u ",1x,I3," u lineto %TR Corner")') pieub(iproc),pjeub(iproc)
2245	WRITE (997,FMT='(I3," u ",1x,I3," u lineto %TL Corner")') pielb(iproc),pjeub(iproc)
2246	WRITE (997,FMT='("closepath stroke")')
2247	WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') &
2248	INT(0.5*(pieub(iproc)+pielb(iproc))), &
2249	INT(0.5*(pjeub(iproc)+pjelb(iproc)-1))
2250	! Write processor label, left justified
2251	WRITE (intStr, FMT='(I4)') iproc-1
2252	! Use postscipt operator 'stringwidth' to get the approximate width of the
2253	! string containing the PE id and then adjust position so it is centred
2254	! about point we've just moveto'd. This doesn't attempt to deal with
2255	! vertical centering.
2256	WRITE (997,FMT='("(",(A),") dup stringwidth pop 2 div neg 0 rmoveto show")') TRIM(ADJUSTL(intStr))
2257	END DO
2258	WRITE (997,FMT='("showpage")')
2259	WRITE (997,FMT='("%%EOF")')
2260	CLOSE(997)
2261
2262	END SUBROUTINE write_partition
2263
2264	END MODULE partition_mod

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

Download in other formats: