New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
partition_mod.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC – NEMO

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

Last change on this file since 4485 was 4485, checked in by trackstand2, 10 years ago

Undo early setting of mbkmax because not needed for halo swaps anymore

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