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

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

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

Download in other formats: