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

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

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

Download in other formats: