source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/partition_mod.F90 @ 4477

Last change on this file since 4477 was 4477, checked in by trackstand2, 7 years ago

Add check to catch case where run out of space for domains while partitioning

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