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

Last change on this file since 3837 was 3837, checked in by trackstand2, 11 years ago
Merge of finiss
File size: 104.6 KB

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

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

Download in other formats: