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.
mppini.F90 in NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/LBC – NEMO

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/LBC/mppini.F90 @ 10409

Last change on this file since 10409 was 10409, checked in by smasson, 5 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 6 add print when using non-optimal domain decomposition, see #2133

  • Property svn:keywords set to Id
File size: 53.9 KB
Line 
1MODULE mppini
2   !!======================================================================
3   !!                       ***  MODULE mppini   ***
4   !! Ocean initialization : distributed memory computing initialization
5   !!======================================================================
6   !! History :  6.0  !  1994-11  (M. Guyon)  Original code
7   !!  OPA       7.0  !  1995-04  (J. Escobar, M. Imbard)
8   !!            8.0  !  1998-05  (M. Imbard, J. Escobar, L. Colombet )  SHMEM and MPI versions
9   !!  NEMO      1.0  !  2004-01  (G. Madec, J.M Molines)  F90 : free form , north fold jpni > 1
10   !!            3.4  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE) add mpp_init_nfdcom
11   !!            3.   ! 2013-06  (I. Epicoco, S. Mocavero, CMCC) mpp_init_nfdcom: setup avoiding MPI communication
12   !!            4.0  !  2016-06  (G. Madec)  use domain configuration file instead of bathymetry file
13   !!            4.0  !  2017-06  (J.M. Molines, T. Lovato) merge of mppini and mppini_2
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!  mpp_init          : Lay out the global domain over processors with/without land processor elimination
18   !!  mpp_init_mask     : Read global bathymetric information to facilitate land suppression
19   !!  mpp_init_ioipsl   : IOIPSL initialization in mpp
20   !!  mpp_init_partition: Calculate MPP domain decomposition
21   !!  factorise         : Calculate the factors of the no. of MPI processes
22   !!  mpp_init_nfdcom   : Setup for north fold exchanges with explicit point-to-point messaging
23   !!----------------------------------------------------------------------
24   USE dom_oce        ! ocean space and time domain
25   USE bdy_oce        ! open BounDarY 
26   !
27   USE lbcnfd  , ONLY : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges
28   USE lib_mpp        ! distribued memory computing library
29   USE iom            ! nemo I/O library
30   USE ioipsl         ! I/O IPSL library
31   USE in_out_manager ! I/O Manager
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC mpp_init       ! called by opa.F90
37
38   INTEGER :: numbot = -1  ! 'bottom_level' local logical unit
39   INTEGER :: numbdy = -1  ! 'bdy_msk'      local logical unit
40   
41   !!----------------------------------------------------------------------
42   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
43   !! $Id$
44   !! Software governed by the CeCILL license (see ./LICENSE)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48#if ! defined key_mpp_mpi
49   !!----------------------------------------------------------------------
50   !!   Default option :                            shared memory computing
51   !!----------------------------------------------------------------------
52
53   SUBROUTINE mpp_init
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE mpp_init  ***
56      !!
57      !! ** Purpose :   Lay out the global domain over processors.
58      !!
59      !! ** Method  :   Shared memory computing, set the local processor
60      !!              variables to the value of the global domain
61      !!----------------------------------------------------------------------
62      !
63      jpimax = jpiglo
64      jpjmax = jpjglo
65      jpi    = jpiglo
66      jpj    = jpjglo
67      jpk    = jpkglo
68      jpim1  = jpi-1                                            ! inner domain indices
69      jpjm1  = jpj-1                                            !   "           "
70      jpkm1  = MAX( 1, jpk-1 )                                  !   "           "
71      jpij   = jpi*jpj
72      jpni   = 1
73      jpnj   = 1
74      jpnij  = jpni*jpnj
75      nimpp  = 1           !
76      njmpp  = 1
77      nlci   = jpi
78      nlcj   = jpj
79      nldi   = 1
80      nldj   = 1
81      nlei   = jpi
82      nlej   = jpj
83      nbondi = 2
84      nbondj = 2
85      nidom  = FLIO_DOM_NONE
86      npolj = jperio
87      l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
88      l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
89      !
90      IF(lwp) THEN
91         WRITE(numout,*)
92         WRITE(numout,*) 'mpp_init : NO massively parallel processing'
93         WRITE(numout,*) '~~~~~~~~ '
94         WRITE(numout,*) '   l_Iperio = ', l_Iperio, '    l_Jperio = ', l_Jperio 
95         WRITE(numout,*) '     npolj  = ',   npolj , '      njmpp  = ', njmpp
96      ENDIF
97      !
98      IF(  jpni /= 1 .OR. jpnj /= 1 .OR. jpnij /= 1 )                                     &
99         CALL ctl_stop( 'mpp_init: equality  jpni = jpnj = jpnij = 1 is not satisfied',   &
100            &           'the domain is lay out for distributed memory computing!' )
101         !
102   END SUBROUTINE mpp_init
103
104#else
105   !!----------------------------------------------------------------------
106   !!   'key_mpp_mpi'                     MPI massively parallel processing
107   !!----------------------------------------------------------------------
108
109
110   SUBROUTINE mpp_init
111      !!----------------------------------------------------------------------
112      !!                  ***  ROUTINE mpp_init  ***
113      !!                   
114      !! ** Purpose :   Lay out the global domain over processors.
115      !!      If land processors are to be eliminated, this program requires the
116      !!      presence of the domain configuration file. Land processors elimination
117      !!      is performed if jpni x jpnj /= jpnij. In this case, using the MPP_PREP
118      !!      preprocessing tool, help for defining the best cutting out.
119      !!
120      !! ** Method  :   Global domain is distributed in smaller local domains.
121      !!      Periodic condition is a function of the local domain position
122      !!      (global boundary or neighbouring domain) and of the global
123      !!      periodic
124      !!      Type :         jperio global periodic condition
125      !!
126      !! ** Action : - set domain parameters
127      !!                    nimpp     : longitudinal index
128      !!                    njmpp     : latitudinal  index
129      !!                    narea     : number for local area
130      !!                    nlci      : first dimension
131      !!                    nlcj      : second dimension
132      !!                    nbondi    : mark for "east-west local boundary"
133      !!                    nbondj    : mark for "north-south local boundary"
134      !!                    nproc     : number for local processor
135      !!                    noea      : number for local neighboring processor
136      !!                    nowe      : number for local neighboring processor
137      !!                    noso      : number for local neighboring processor
138      !!                    nono      : number for local neighboring processor
139      !!----------------------------------------------------------------------
140      INTEGER ::   ji, jj, jn, jproc, jarea   ! dummy loop indices
141      INTEGER ::   inijmin
142      INTEGER ::   i2add
143      INTEGER ::   inum                       ! local logical unit
144      INTEGER ::   idir, ifreq, icont         ! local integers
145      INTEGER ::   ii, il1, ili, imil         !   -       -
146      INTEGER ::   ij, il2, ilj, ijm1         !   -       -
147      INTEGER ::   iino, ijno, iiso, ijso     !   -       -
148      INTEGER ::   iiea, ijea, iiwe, ijwe     !   -       -
149      INTEGER ::   iarea0                     !   -       -
150      INTEGER ::   ierr, ios                  !
151      INTEGER ::   inbi, inbj, iimax,  ijmax  !
152      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   iin, ii_nono, ii_noea          ! 1D workspace
153      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ijn, ii_noso, ii_nowe          !  -     -
154      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ilci, ibondi, ipproc   ! 2D workspace
155      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ilcj, ibondj, ipolj    !  -     -
156      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilei, ildi, iono, ioea         !  -     -
157      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilej, ildj, ioso, iowe         !  -     -
158      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   llisoce                        !  -     -
159      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           &
160           &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
161           &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             & 
162           &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
163           &             cn_ice, nn_ice_dta,                                     &
164           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
165           &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
166      !!----------------------------------------------------------------------
167
168      ! do we need to take into account bdy_msk?
169      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY
170      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
171903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp )
172      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY
173      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
174904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp )
175      !
176      IF(               ln_read_cfg ) CALL iom_open( cn_domcfg,    numbot )
177      IF( ln_bdy .AND. ln_mask_file ) CALL iom_open( cn_mask_file, numbdy )
178      !
179      !  1. Dimension arrays for subdomains
180      ! -----------------------------------
181      !
182      ! If dimensions of processor grid weren't specified in the namelist file
183      ! then we calculate them here now that we have our communicator size
184      IF( jpni < 1 .OR. jpnj < 1 ) THEN
185         CALL mpp_init_bestpartition( mppsize, jpni, jpnj )
186      ELSE
187         CALL mpp_init_bestpartition( mppsize, inbi, inbj )
188         CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax )
189         CALL mpp_basic_decomposition( inbi, inbj,  iimax,  ijmax )
190         IF( iimax*ijmax < jpimax*jpjmax ) THEN
191            WRITE(ctmp1,*) '   The chosen domain decomposition ', jpni, jpnj
192            WRITE(ctmp2,*) '   has larger MPI subdomains ', jpimax, jpjmax, jpimax*jpjmax
193            WRITE(ctmp3,*) '   than the following domain decompostion ', inbi, inbj
194            WRITE(ctmp4,*) '   which MPI subdomains size is ', iimax, ijmax, iimax*ijmax
195            CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4, '    --- YOU ARE WASTING CPU... ---' )
196         ENDIF
197      ENDIF
198     
199      ! look for land mpi subdomains...
200      ALLOCATE( llisoce(jpni,jpnj) )
201      CALL mpp_init_isoce( jpni, jpnj, llisoce )
202      inijmin = COUNT( llisoce )   ! number of oce subdomains
203
204      IF( mppsize < inijmin ) THEN
205         WRITE(ctmp1,*) '   With this specified domain decomposition: jpni =', jpni, ' jpnj =', jpnj
206         WRITE(ctmp2,*) '   we can eliminate only ', jpni*jpnj - inijmin, ' land mpi subdomains therefore '
207         WRITE(ctmp3,*) '   the number of ocean mpi subdomains (', inijmin,') exceed the number of MPI processes:', mppsize
208         WRITE(ctmp4,*) '   ==>>> There is the list of best domain decompositions you should use: '
209         CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 )
210         CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core
211         CALL ctl_stop( 'STOP' )
212      ENDIF
213
214      IF( mppsize > jpni*jpnj ) THEN
215         WRITE(ctmp1,*) '   The number of mpi processes: ', mppsize
216         WRITE(ctmp2,*) '   exceeds the maximum number of subdomains (ocean+land) = ', jpni*jpnj
217         WRITE(ctmp3,*) '   defined by the following domain decomposition: jpni =', jpni, ' jpnj =', jpnj
218         WRITE(ctmp4,*) '   ==>>> There is the list of best domain decompositions you should use: '
219         CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 )
220         CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core
221         CALL ctl_stop( 'STOP' )
222      ENDIF
223
224      jpnij = mppsize   ! force jpnij definition <-- remove as much land subdomains as needed to reach this condition
225      IF( mppsize > inijmin ) THEN
226         WRITE(ctmp1,*) '   The number of mpi processes: ', mppsize
227         WRITE(ctmp2,*) '   exceeds the maximum number of ocean subdomains = ', inijmin
228         WRITE(ctmp3,*) '   we suppressed ', jpni*jpnj - mppsize, ' land subdomains '
229         WRITE(ctmp4,*) '   BUT we had to keep ', mppsize - inijmin, ' land subdomains that are useless...'
230         CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 )
231      ELSE   ! mppsize = inijmin
232         IF(lwp) THEN
233            WRITE(numout,*) 'mpp_init: You use an optimal domain decomposition'
234            WRITE(numout,*) '~~~~~~~~ '
235            WRITE(numout,*) '   Number of mpi processes: ', mppsize
236            WRITE(numout,*) '   Number of ocean subdomains = ', inijmin
237            WRITE(numout,*) '   Number of suppressed land subdomains = ', jpni*jpnj - inijmin
238         ENDIF
239      ENDIF
240
241      IF( numbot /= -1 )   CALL iom_close( numbot )
242      IF( numbdy /= -1 )   CALL iom_close( numbdy )
243   
244      ALLOCATE(  nfiimpp(jpni,jpnj), nfipproc(jpni,jpnj), nfilcit(jpni,jpnj) ,    &
245         &       nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) ,    &
246         &       njmppt(jpnij) , ibonjt(jpnij) , nldit(jpnij) , nldjt(jpnij) ,    &
247         &                                       nleit(jpnij) , nlejt(jpnij) ,    &
248         &       iin(jpnij), ii_nono(jpnij), ii_noea(jpnij),   &
249         &       ijn(jpnij), ii_noso(jpnij), ii_nowe(jpnij),   &
250         &       iimppt(jpni,jpnj), ilci(jpni,jpnj), ibondi(jpni,jpnj), ipproc(jpni,jpnj),   &
251         &       ijmppt(jpni,jpnj), ilcj(jpni,jpnj), ibondj(jpni,jpnj), ipolj(jpni,jpnj),   &
252         &       ilei(jpni,jpnj), ildi(jpni,jpnj), iono(jpni,jpnj), ioea(jpni,jpnj),   &
253         &       ilej(jpni,jpnj), ildj(jpni,jpnj), ioso(jpni,jpnj), iowe(jpni,jpnj),   &
254         &       STAT=ierr )
255      CALL mpp_sum( 'mppini', ierr )
256      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' )
257     
258#if defined key_agrif
259      IF( .NOT. Agrif_Root() ) THEN       ! AGRIF children: specific setting (cf. agrif_user.F90)
260         IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells )   &
261            CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells' )
262         IF( jpjglo /= nbcellsy + 2 + 2*nbghostcells )   &
263            CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpjglo == nbcellsy + 2 + 2*nbghostcells' )
264         IF( ln_use_jattr )   CALL ctl_stop( 'STOP', 'mpp_init:Agrif children requires ln_use_jattr = .false. ' )
265      ENDIF
266#endif
267      !
268      !  2. Index arrays for subdomains
269      ! -----------------------------------
270      !
271      nreci = 2 * nn_hls
272      nrecj = 2 * nn_hls
273      CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ilci, ilcj )
274      nfiimpp(:,:) = iimppt(:,:)
275      nfilcit(:,:) = ilci(:,:)
276      !
277      IF(lwp) THEN
278         WRITE(numout,*)
279         WRITE(numout,*) 'MPI Message Passing MPI - domain lay out over processors'
280         WRITE(numout,*)
281         WRITE(numout,*) '   defines mpp subdomains'
282         WRITE(numout,*) '      jpni = ', jpni 
283         WRITE(numout,*) '      jpnj = ', jpnj
284         WRITE(numout,*)
285         WRITE(numout,*) '      sum ilci(i,1) = ', sum(ilci(:,1)), ' jpiglo = ', jpiglo
286         WRITE(numout,*) '      sum ilcj(1,j) = ', sum(ilcj(1,:)), ' jpjglo = ', jpjglo
287      ENDIF
288     
289      ! 3. Subdomain description in the Regular Case
290      ! --------------------------------------------
291      ! specific cases where there is no communication -> must do the periodicity by itself
292      ! Warning: because of potential land-area suppression, do not use nbond[ij] == 2 
293      l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
294      l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
295     
296      DO jarea = 1, jpni*jpnj
297         !
298         iarea0 = jarea - 1
299         ii = 1 + MOD(iarea0,jpni)
300         ij = 1 +     iarea0/jpni
301         ili = ilci(ii,ij)
302         ilj = ilcj(ii,ij)
303         ibondi(ii,ij) = 0                         ! default: has e-w neighbours
304         IF( ii   ==    1 )   ibondi(ii,ij) = -1   ! first column, has only e neighbour
305         IF( ii   == jpni )   ibondi(ii,ij) =  1   ! last column,  has only w neighbour
306         IF( jpni ==    1 )   ibondi(ii,ij) =  2   ! has no e-w neighbour
307         ibondj(ii,ij) = 0                         ! default: has n-s neighbours
308         IF( ij   ==    1 )   ibondj(ii,ij) = -1   ! first row, has only n neighbour
309         IF( ij   == jpnj )   ibondj(ii,ij) =  1   ! last row,  has only s neighbour
310         IF( jpnj ==    1 )   ibondj(ii,ij) =  2   ! has no n-s neighbour
311
312         ! Subdomain neighbors (get their zone number): default definition
313         ioso(ii,ij) = iarea0 - jpni
314         iowe(ii,ij) = iarea0 - 1
315         ioea(ii,ij) = iarea0 + 1
316         iono(ii,ij) = iarea0 + jpni
317         ildi(ii,ij) =  1  + nn_hls
318         ilei(ii,ij) = ili - nn_hls
319         ildj(ii,ij) =  1  + nn_hls
320         ilej(ii,ij) = ilj - nn_hls
321
322         ! East-West periodicity: change ibondi, ioea, iowe
323         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7 ) THEN
324            IF( jpni  /= 1 )   ibondi(ii,ij) = 0                        ! redefine: all have e-w neighbours
325            IF( ii ==    1 )   iowe(ii,ij) = iarea0 +        (jpni-1)   ! redefine: first column, address of w neighbour
326            IF( ii == jpni )   ioea(ii,ij) = iarea0 -        (jpni-1)   ! redefine: last column,  address of e neighbour
327         ENDIF
328
329         ! Simple North-South periodicity: change ibondj, ioso, iono
330         IF( jperio == 2 .OR. jperio == 7 ) THEN
331            IF( jpnj  /= 1 )   ibondj(ii,ij) = 0                        ! redefine: all have n-s neighbours
332            IF( ij ==    1 )   ioso(ii,ij) = iarea0 + jpni * (jpnj-1)   ! redefine: first row, address of s neighbour
333            IF( ij == jpnj )   iono(ii,ij) = iarea0 - jpni * (jpnj-1)   ! redefine: last row,  address of n neighbour
334         ENDIF
335
336         ! North fold: define ipolj, change iono. Warning: we do not change ibondj...
337         ipolj(ii,ij) = 0
338         IF( jperio == 3 .OR. jperio == 4 ) THEN
339            ijm1 = jpni*(jpnj-1)
340            imil = ijm1+(jpni+1)/2
341            IF( jarea > ijm1 ) ipolj(ii,ij) = 3
342            IF( MOD(jpni,2) == 1 .AND. jarea == imil ) ipolj(ii,ij) = 4
343            IF( ipolj(ii,ij) == 3 ) iono(ii,ij) = jpni*jpnj-jarea+ijm1   ! MPI rank of northern neighbour
344         ENDIF
345         IF( jperio == 5 .OR. jperio == 6 ) THEN
346            ijm1 = jpni*(jpnj-1)
347            imil = ijm1+(jpni+1)/2
348            IF( jarea > ijm1) ipolj(ii,ij) = 5
349            IF( MOD(jpni,2) == 1 .AND. jarea == imil ) ipolj(ii,ij) = 6
350            IF( ipolj(ii,ij) == 5) iono(ii,ij) = jpni*jpnj-jarea+ijm1    ! MPI rank of northern neighbour
351         ENDIF
352         !
353      END DO
354
355      ! 4. deal with land subdomains
356      ! ----------------------------
357      !
358      ! specify which subdomains are oce subdomains; other are land subdomains
359      ipproc(:,:) = -1
360      icont = -1
361      DO jarea = 1, jpni*jpnj
362         iarea0 = jarea - 1
363         ii = 1 + MOD(iarea0,jpni)
364         ij = 1 +     iarea0/jpni
365         IF( llisoce(ii,ij) ) THEN
366            icont = icont + 1
367            ipproc(ii,ij) = icont
368            iin(icont+1) = ii
369            ijn(icont+1) = ij
370         ENDIF
371      END DO
372      ! if needed add some land subdomains to reach jpnij active subdomains
373      i2add = jpnij - inijmin
374      DO jarea = 1, jpni*jpnj
375         iarea0 = jarea - 1
376         ii = 1 + MOD(iarea0,jpni)
377         ij = 1 +     iarea0/jpni
378         IF( .NOT. llisoce(ii,ij) .AND. i2add > 0 ) THEN
379            icont = icont + 1
380            ipproc(ii,ij) = icont
381            iin(icont+1) = ii
382            ijn(icont+1) = ij
383            i2add = i2add - 1
384         ENDIF
385      END DO
386      nfipproc(:,:) = ipproc(:,:)
387
388      ! neighbour treatment: change ibondi, ibondj if next to a land zone
389      DO jarea = 1, jpni*jpnj
390         ii = 1 + MOD( jarea-1  , jpni )
391         ij = 1 +     (jarea-1) / jpni
392         ! land-only area with an active n neigbour
393         IF ( ipproc(ii,ij) == -1 .AND. 0 <= iono(ii,ij) .AND. iono(ii,ij) <= jpni*jpnj-1 ) THEN
394            iino = 1 + MOD( iono(ii,ij) , jpni )                    ! ii index of this n neigbour
395            ijno = 1 +      iono(ii,ij) / jpni                      ! ij index of this n neigbour
396            ! In case of north fold exchange: I am the n neigbour of my n neigbour!! (#1057)
397            ! --> for northern neighbours of northern row processors (in case of north-fold)
398            !     need to reverse the LOGICAL direction of communication
399            idir = 1                                           ! we are indeed the s neigbour of this n neigbour
400            IF( ij == jpnj .AND. ijno == jpnj )   idir = -1    ! both are on the last row, we are in fact the n neigbour
401            IF( ibondj(iino,ijno) == idir     )   ibondj(iino,ijno) =   2     ! this n neigbour had only a s/n neigbour -> no more
402            IF( ibondj(iino,ijno) == 0        )   ibondj(iino,ijno) = -idir   ! this n neigbour had both, n-s neighbours -> keep 1
403         ENDIF
404         ! land-only area with an active s neigbour
405         IF( ipproc(ii,ij) == -1 .AND. 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= jpni*jpnj-1 ) THEN
406            iiso = 1 + MOD( ioso(ii,ij) , jpni )                    ! ii index of this s neigbour
407            ijso = 1 +      ioso(ii,ij) / jpni                      ! ij index of this s neigbour
408            IF( ibondj(iiso,ijso) == -1 )   ibondj(iiso,ijso) = 2   ! this s neigbour had only a n neigbour    -> no more neigbour
409            IF( ibondj(iiso,ijso) ==  0 )   ibondj(iiso,ijso) = 1   ! this s neigbour had both, n-s neighbours -> keep s neigbour
410         ENDIF
411         ! land-only area with an active e neigbour
412         IF( ipproc(ii,ij) == -1 .AND. 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= jpni*jpnj-1 ) THEN
413            iiea = 1 + MOD( ioea(ii,ij) , jpni )                    ! ii index of this e neigbour
414            ijea = 1 +      ioea(ii,ij) / jpni                      ! ij index of this e neigbour
415            IF( ibondi(iiea,ijea) == 1 )   ibondi(iiea,ijea) =  2   ! this e neigbour had only a w neigbour    -> no more neigbour
416            IF( ibondi(iiea,ijea) == 0 )   ibondi(iiea,ijea) = -1   ! this e neigbour had both, e-w neighbours -> keep e neigbour
417         ENDIF
418         ! land-only area with an active w neigbour
419         IF( ipproc(ii,ij) == -1 .AND. 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= jpni*jpnj-1) THEN
420            iiwe = 1 + MOD( iowe(ii,ij) , jpni )                    ! ii index of this w neigbour
421            ijwe = 1 +      iowe(ii,ij) / jpni                      ! ij index of this w neigbour
422            IF( ibondi(iiwe,ijwe) == -1 )   ibondi(iiwe,ijwe) = 2   ! this w neigbour had only a e neigbour    -> no more neigbour
423            IF( ibondi(iiwe,ijwe) ==  0 )   ibondi(iiwe,ijwe) = 1   ! this w neigbour had both, e-w neighbours -> keep w neigbour
424         ENDIF
425      END DO
426
427      ! Update il[de][ij] according to modified ibond[ij]
428      ! ----------------------
429      DO jproc = 1, jpnij
430         ii = iin(jproc)
431         ij = ijn(jproc)
432         IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) =  1
433         IF( ibondi(ii,ij) ==  1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij)
434         IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) =  1
435         IF( ibondj(ii,ij) ==  1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij)
436      END DO
437     
438      ! 5. Subdomain print
439      ! ------------------
440      IF(lwp) THEN
441         ifreq = 4
442         il1 = 1
443         DO jn = 1, (jpni-1)/ifreq+1
444            il2 = MIN(jpni,il1+ifreq-1)
445            WRITE(numout,*)
446            WRITE(numout,9400) ('***',ji=il1,il2-1)
447            DO jj = jpnj, 1, -1
448               WRITE(numout,9403) ('   ',ji=il1,il2-1)
449               WRITE(numout,9402) jj, (ilci(ji,jj),ilcj(ji,jj),ji=il1,il2)
450               WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2)
451               WRITE(numout,9403) ('   ',ji=il1,il2-1)
452               WRITE(numout,9400) ('***',ji=il1,il2-1)
453            END DO
454            WRITE(numout,9401) (ji,ji=il1,il2)
455            il1 = il1+ifreq
456         END DO
457 9400    FORMAT('           ***'   ,20('*************',a3)    )
458 9403    FORMAT('           *     ',20('         *   ',a3)    )
459 9401    FORMAT('              '   ,20('   ',i3,'          ') )
460 9402    FORMAT('       ',i3,' *  ',20(i3,'  x',i3,'   *   ') )
461 9404    FORMAT('           *  '   ,20('      ',i3,'   *   ') )
462      ENDIF
463         
464      ! just to save nono etc for all proc
465      ! warning ii*ij (zone) /= nproc (processors)!
466      ! ioso = zone number, ii_noso = proc number
467      ii_noso(:) = -1
468      ii_nono(:) = -1
469      ii_noea(:) = -1
470      ii_nowe(:) = -1 
471      DO jproc = 1, jpnij
472         ii = iin(jproc)
473         ij = ijn(jproc)
474         IF( 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= (jpni*jpnj-1) ) THEN
475            iiso = 1 + MOD( ioso(ii,ij) , jpni )
476            ijso = 1 +      ioso(ii,ij) / jpni
477            ii_noso(jproc) = ipproc(iiso,ijso)
478         ENDIF
479         IF( 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= (jpni*jpnj-1) ) THEN
480          iiwe = 1 + MOD( iowe(ii,ij) , jpni )
481          ijwe = 1 +      iowe(ii,ij) / jpni
482          ii_nowe(jproc) = ipproc(iiwe,ijwe)
483         ENDIF
484         IF( 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= (jpni*jpnj-1) ) THEN
485            iiea = 1 + MOD( ioea(ii,ij) , jpni )
486            ijea = 1 +      ioea(ii,ij) / jpni
487            ii_noea(jproc)= ipproc(iiea,ijea)
488         ENDIF
489         IF( 0 <= iono(ii,ij) .AND. iono(ii,ij) <= (jpni*jpnj-1) ) THEN
490            iino = 1 + MOD( iono(ii,ij) , jpni )
491            ijno = 1 +      iono(ii,ij) / jpni
492            ii_nono(jproc)= ipproc(iino,ijno)
493         ENDIF
494      END DO
495   
496      ! 6. Change processor name
497      ! ------------------------
498      ii = iin(narea)
499      ij = ijn(narea)
500      !
501      ! set default neighbours
502      noso = ii_noso(narea)
503      nowe = ii_nowe(narea)
504      noea = ii_noea(narea)
505      nono = ii_nono(narea)
506      nlci = ilci(ii,ij) 
507      nldi = ildi(ii,ij)
508      nlei = ilei(ii,ij)
509      nlcj = ilcj(ii,ij) 
510      nldj = ildj(ii,ij)
511      nlej = ilej(ii,ij)
512      nbondi = ibondi(ii,ij)
513      nbondj = ibondj(ii,ij)
514      nimpp = iimppt(ii,ij) 
515      njmpp = ijmppt(ii,ij)
516      jpi = nlci
517      jpj = nlcj
518      jpk = jpkglo                                             ! third dim
519#if defined key_agrif
520      ! simple trick to use same vertical grid as parent but different number of levels:
521      ! Save maximum number of levels in jpkglo, then define all vertical grids with this number.
522      ! Suppress once vertical online interpolation is ok
523!!$      IF(.NOT.Agrif_Root())   jpkglo = Agrif_Parent( jpkglo )
524#endif
525      jpim1 = jpi-1                                            ! inner domain indices
526      jpjm1 = jpj-1                                            !   "           "
527      jpkm1 = MAX( 1, jpk-1 )                                  !   "           "
528      jpij  = jpi*jpj                                          !  jpi x j
529      DO jproc = 1, jpnij
530         ii = iin(jproc)
531         ij = ijn(jproc)
532         nlcit(jproc) = ilci(ii,ij)
533         nldit(jproc) = ildi(ii,ij)
534         nleit(jproc) = ilei(ii,ij)
535         nlcjt(jproc) = ilcj(ii,ij)
536         nldjt(jproc) = ildj(ii,ij)
537         nlejt(jproc) = ilej(ii,ij)
538         ibonit(jproc) = ibondi(ii,ij)
539         ibonjt(jproc) = ibondj(ii,ij)
540         nimppt(jproc) = iimppt(ii,ij) 
541         njmppt(jproc) = ijmppt(ii,ij) 
542      END DO
543
544      ! Save processor layout in ascii file
545      IF (lwp) THEN
546         CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
547         WRITE(inum,'(a)') '   jpnij   jpimax  jpjmax    jpk  jpiglo  jpjglo'//&
548   &           ' ( local:    narea     jpi     jpj )'
549         WRITE(inum,'(6i8,a,3i8,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,&
550   &           ' ( local: ',narea,jpi,jpj,' )'
551         WRITE(inum,'(a)') 'nproc nlci nlcj nldi nldj nlei nlej nimp njmp nono noso nowe noea nbondi nbondj '
552
553         DO jproc = 1, jpnij
554            WRITE(inum,'(13i5,2i7)')   jproc-1, nlcit  (jproc), nlcjt  (jproc),   &
555               &                                nldit  (jproc), nldjt  (jproc),   &
556               &                                nleit  (jproc), nlejt  (jproc),   &
557               &                                nimppt (jproc), njmppt (jproc),   & 
558               &                                ii_nono(jproc), ii_noso(jproc),   &
559               &                                ii_nowe(jproc), ii_noea(jproc),   &
560               &                                ibonit (jproc), ibonjt (jproc) 
561         END DO
562      END IF
563
564      !                          ! north fold parameter
565      ! Defined npolj, either 0, 3 , 4 , 5 , 6
566      ! In this case the important thing is that npolj /= 0
567      ! Because if we go through these line it is because jpni >1 and thus
568      ! we must use lbcnorthmpp, which tests only npolj =0 or npolj /= 0
569      npolj = 0
570      ij = ijn(narea)
571      IF( jperio == 3 .OR. jperio == 4 ) THEN
572         IF( ij == jpnj )   npolj = 3
573      ENDIF
574      IF( jperio == 5 .OR. jperio == 6 ) THEN
575         IF( ij == jpnj )   npolj = 5
576      ENDIF
577      !
578      nproc = narea-1
579      IF(lwp) THEN
580         WRITE(numout,*)
581         WRITE(numout,*) '   resulting internal parameters : '
582         WRITE(numout,*) '      nproc  = ', nproc
583         WRITE(numout,*) '      nowe   = ', nowe  , '   noea  =  ', noea
584         WRITE(numout,*) '      nono   = ', nono  , '   noso  =  ', noso
585         WRITE(numout,*) '      nbondi = ', nbondi
586         WRITE(numout,*) '      nbondj = ', nbondj
587         WRITE(numout,*) '      npolj  = ', npolj
588         WRITE(numout,*) '    l_Iperio = ', l_Iperio
589         WRITE(numout,*) '    l_Jperio = ', l_Jperio
590         WRITE(numout,*) '      nlci   = ', nlci
591         WRITE(numout,*) '      nlcj   = ', nlcj
592         WRITE(numout,*) '      nimpp  = ', nimpp
593         WRITE(numout,*) '      njmpp  = ', njmpp
594         WRITE(numout,*) '      nreci  = ', nreci 
595         WRITE(numout,*) '      nrecj  = ', nrecj 
596         WRITE(numout,*) '      nn_hls = ', nn_hls 
597      ENDIF
598
599      !                          ! Prepare mpp north fold
600      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
601         CALL mpp_ini_north
602         IF (lwp) THEN
603            WRITE(numout,*)
604            WRITE(numout,*) '   ==>>>   North fold boundary prepared for jpni >1'
605            ! additional prints in layout.dat
606            WRITE(inum,*)
607            WRITE(inum,*)
608            WRITE(inum,*) 'number of subdomains located along the north fold : ', ndim_rank_north
609            WRITE(inum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
610            DO jproc = 1, ndim_rank_north, 5
611               WRITE(inum,*) nrank_north( jproc:MINVAL( (/jproc+4,ndim_rank_north/) ) )
612            END DO
613         ENDIF
614      ENDIF
615      !
616      CALL mpp_init_ioipsl       ! Prepare NetCDF output file (if necessary)
617      !
618      IF( ln_nnogather ) THEN
619         CALL mpp_init_nfdcom     ! northfold neighbour lists
620         IF (lwp) THEN
621            WRITE(inum,*)
622            WRITE(inum,*)
623            WRITE(inum,*) 'north fold exchanges with explicit point-to-point messaging :'
624            WRITE(inum,*) 'nfsloop : ', nfsloop
625            WRITE(inum,*) 'nfeloop : ', nfeloop
626            WRITE(inum,*) 'nsndto : ', nsndto
627            WRITE(inum,*) 'isendto : ', isendto
628         ENDIF
629      ENDIF
630      !
631      IF (lwp) CLOSE(inum)   
632      !
633      DEALLOCATE(iin, ijn, ii_nono, ii_noea, ii_noso, ii_nowe,    &
634         &       iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj,   &
635         &       ilci, ilcj, ilei, ilej, ildi, ildj,              &
636         &       iono, ioea, ioso, iowe, llisoce)
637      !
638    END SUBROUTINE mpp_init
639
640
641    SUBROUTINE mpp_basic_decomposition( knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj)
642      !!----------------------------------------------------------------------
643      !!                  ***  ROUTINE mpp_basic_decomposition  ***
644      !!                   
645      !! ** Purpose :   Lay out the global domain over processors.
646      !!
647      !! ** Method  :   Global domain is distributed in smaller local domains.
648      !!
649      !! ** Action : - set for all knbi*knbj domains:
650      !!                    kimppt     : longitudinal index
651      !!                    kjmppt     : latitudinal  index
652      !!                    klci       : first dimension
653      !!                    klcj       : second dimension
654      !!----------------------------------------------------------------------
655      INTEGER,                                 INTENT(in   ) ::   knbi, knbj
656      INTEGER,                                 INTENT(  out) ::   kimax, kjmax
657      INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT(  out) ::   kimppt, kjmppt
658      INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT(  out) ::   klci, klcj
659      !
660      INTEGER ::   ji, jj
661      INTEGER ::   iresti, irestj
662      INTEGER ::   ireci, irecj
663      !!----------------------------------------------------------------------
664      !
665#if defined key_nemocice_decomp
666      kimax = ( nx_global+2-2*nn_hls + (knbi-1) ) / knbi + 2*nn_hls    ! first  dim.
667      kjmax = ( ny_global+2-2*nn_hls + (knbj-1) ) / knbj + 2*nn_hls    ! second dim.
668#else
669      kimax = ( jpiglo - 2*nn_hls + (knbi-1) ) / knbi + 2*nn_hls    ! first  dim.
670      kjmax = ( jpjglo - 2*nn_hls + (knbj-1) ) / knbj + 2*nn_hls    ! second dim.
671#endif
672      IF( .NOT. PRESENT(kimppt) ) RETURN
673      !
674      !  1. Dimension arrays for subdomains
675      ! -----------------------------------
676      !  Computation of local domain sizes klci() klcj()
677      !  These dimensions depend on global sizes knbi,knbj and jpiglo,jpjglo
678      !  The subdomains are squares lesser than or equal to the global
679      !  dimensions divided by the number of processors minus the overlap array.
680      !
681      ireci = 2 * nn_hls
682      irecj = 2 * nn_hls
683      iresti = 1 + MOD( jpiglo - ireci -1 , knbi )
684      irestj = 1 + MOD( jpjglo - irecj -1 , knbj )
685      !
686      !  Need to use kimax and kjmax here since jpi and jpj not yet defined
687#if defined key_nemocice_decomp
688      ! Change padding to be consistent with CICE
689      klci(1:knbi-1      ,:) = kimax
690      klci(knbi          ,:) = jpiglo - (knbi - 1) * (kimax - nreci)
691      klcj(:,      1:knbj-1) = kjmax
692      klcj(:,          knbj) = jpjglo - (knbj - 1) * (kjmax - nrecj)
693#else
694      klci(1:iresti      ,:) = kimax
695      klci(iresti+1:knbi ,:) = kimax-1
696      klcj(:,      1:irestj) = kjmax
697      klcj(:, irestj+1:knbj) = kjmax-1
698#endif
699
700      !  2. Index arrays for subdomains
701      ! -------------------------------
702      kimppt(:,:) = 1
703      kjmppt(:,:) = 1
704      !
705      IF( knbi > 1 ) THEN
706         DO jj = 1, knbj
707            DO ji = 2, knbi
708               kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - ireci
709            END DO
710         END DO
711      ENDIF
712      !
713      IF( knbj > 1 )THEN
714         DO jj = 2, knbj
715            DO ji = 1, knbi
716               kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - irecj
717            END DO
718         END DO
719      ENDIF
720     
721   END SUBROUTINE mpp_basic_decomposition
722
723
724   SUBROUTINE mpp_init_bestpartition( knbij, knbi, knbj, ldlist )
725      !!----------------------------------------------------------------------
726      !!                 ***  ROUTINE mpp_init_bestpartition  ***
727      !!
728      !! ** Purpose :
729      !!
730      !! ** Method  :
731      !!----------------------------------------------------------------------
732      INTEGER,           INTENT(in   ) ::   knbij         ! total number if subdomains               (knbi*knbj)
733      INTEGER, OPTIONAL, INTENT(  out) ::   knbi, knbj    ! number if subdomains along i and j (knbi and knbj)
734      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldlist        ! .true.: print the list the best domain decompositions (with land)
735      !
736      INTEGER :: ji, jj, ii, iitarget
737      INTEGER :: iszitst, iszjtst
738      INTEGER :: isziref, iszjref
739      INTEGER :: inbij, iszij
740      INTEGER :: inbimax, inbjmax, inbijmax
741      INTEGER :: isz0, isz1
742      INTEGER, DIMENSION(  :), ALLOCATABLE :: indexok
743      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi0, inbj0, inbij0   ! number of subdomains along i,j
744      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi0, iszj0, iszij0   ! max size of the subdomains along i,j
745      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi1, inbj1, inbij1   ! number of subdomains along i,j
746      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi1, iszj1, iszij1   ! max size of the subdomains along i,j
747      LOGICAL :: llist
748      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d                 ! max size of the subdomains along i,j
749      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisoce              !  -     -
750      REAL(wp)::   zpropland
751      !!----------------------------------------------------------------------
752      !
753      llist = .FALSE.
754      IF( PRESENT(ldlist) ) llist = ldlist
755
756      CALL mpp_init_landprop( zpropland )                      ! get the proportion of land point over the gloal domain
757      inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) )    ! define the largest possible value for jpni*jpnj
758      !
759      IF( llist ) THEN   ;   inbijmax = inbij*2
760      ELSE               ;   inbijmax = inbij
761      ENDIF
762      !
763      ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax))
764      !
765      inbimax = 0
766      inbjmax = 0
767      isziref = jpiglo*jpjglo+1
768      iszjref = jpiglo*jpjglo+1
769      !
770      ! get the list of knbi that gives a smaller jpimax than knbi-1
771      ! get the list of knbj that gives a smaller jpjmax than knbj-1
772      DO ji = 1, inbijmax     
773#if defined key_nemocice_decomp
774         iszitst = ( nx_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
775#else
776         iszitst = ( jpiglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls
777#endif
778         IF( iszitst < isziref ) THEN
779            isziref = iszitst
780            inbimax = inbimax + 1
781            inbi0(inbimax) = ji
782            iszi0(inbimax) = isziref
783         ENDIF
784#if defined key_nemocice_decomp
785         iszjtst = ( ny_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
786#else
787         iszjtst = ( jpjglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls
788#endif
789         IF( iszjtst < iszjref ) THEN
790            iszjref = iszjtst
791            inbjmax = inbjmax + 1
792            inbj0(inbjmax) = ji
793            iszj0(inbjmax) = iszjref
794         ENDIF
795      END DO
796
797      ! combine these 2 lists to get all possible knbi*knbj <  inbijmax
798      ALLOCATE( llmsk2d(inbimax,inbjmax) )
799      DO jj = 1, inbjmax
800         DO ji = 1, inbimax
801            IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN   ;   llmsk2d(ji,jj) = .TRUE.
802            ELSE                                            ;   llmsk2d(ji,jj) = .FALSE.
803            ENDIF
804         END DO
805      END DO
806      isz1 = COUNT(llmsk2d)
807      ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) )
808      ii = 0
809      DO jj = 1, inbjmax
810         DO ji = 1, inbimax
811            IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN
812               ii = ii + 1
813               inbi1(ii) = inbi0(ji)
814               inbj1(ii) = inbj0(jj)
815               iszi1(ii) = iszi0(ji)
816               iszj1(ii) = iszj0(jj)
817            END IF
818         END DO
819      END DO
820      DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
821      DEALLOCATE( llmsk2d )
822
823      ALLOCATE( inbij1(isz1), iszij1(isz1) )
824      inbij1(:) = inbi1(:) * inbj1(:)
825      iszij1(:) = iszi1(:) * iszj1(:)
826
827      ! if therr is no land and no print
828      IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN
829         ! get the smaller partition which gives the smallest subdomain size
830         ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1)
831         knbi = inbi1(ii)
832         knbj = inbj1(ii)
833         DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 )
834         RETURN
835      ENDIF
836
837      ! extract only the partitions which reduce the subdomain size in comparison with smaller partitions
838      ALLOCATE( indexok(isz1) )                                 ! to store indices of the best partitions
839      isz0 = 0                                                  ! number of best partitions     
840      inbij = 1                                                 ! start with the min value of inbij1 => 1
841      iszij = jpiglo*jpjglo+1                                   ! default: larger than global domain
842      DO WHILE( inbij <= inbijmax )                             ! if we did not reach the max of inbij1
843         ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1)   ! warning: send back the first occurence if multiple results
844         IF ( iszij1(ii) < iszij ) THEN
845            isz0 = isz0 + 1
846            indexok(isz0) = ii
847            iszij = iszij1(ii)
848         ENDIF
849         inbij = MINVAL(inbij1, mask = inbij1 > inbij)   ! warning: return largest integer value if mask = .false. everywhere
850      END DO
851      DEALLOCATE( inbij1, iszij1 )
852
853      ! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size)
854      ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) )
855      DO ji = 1, isz0
856         ii = indexok(ji)
857         inbi0(ji) = inbi1(ii)
858         inbj0(ji) = inbj1(ii)
859         iszi0(ji) = iszi1(ii)
860         iszj0(ji) = iszj1(ii)
861      END DO
862      DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 )
863
864      IF( llist ) THEN  ! we print about 21 best partitions
865         IF(lwp) THEN
866            WRITE(numout,*)
867            WRITE(numout,         *) '                  For your information:'
868            WRITE(numout,'(a,i5,a)') '  list of the best partitions around ',   knbij, ' mpi processes'
869            WRITE(numout,         *) '  --------------------------------------', '-----', '--------------'
870            WRITE(numout,*)
871         END IF
872         iitarget = MINLOC( inbi0(:)*inbj0(:), mask = inbi0(:)*inbj0(:) >= knbij, dim = 1 )
873         DO ji = MAX(1,iitarget-10), MIN(isz0,iitarget+10)
874            ALLOCATE( llisoce(inbi0(ji), inbj0(ji)) )
875            CALL mpp_init_isoce( inbi0(ji), inbj0(ji), llisoce ) ! Warning: must be call by all cores (call mpp_sum)
876            inbij = COUNT(llisoce)
877            DEALLOCATE( llisoce )
878            IF(lwp) WRITE(numout,'(a, i5, a, i5, a, i4, a, i4, a, i9, a, i5, a, i5, a)')    &
879               &     'nb_cores ' , inbij,' oce + ', inbi0(ji)*inbj0(ji) - inbij             &
880               &                                , ' land ( ', inbi0(ji),' x ', inbj0(ji),   &
881               & ' ), nb_points ', iszi0(ji)*iszj0(ji),' ( ', iszi0(ji),' x ', iszj0(ji),' )'
882         END DO
883         DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
884         RETURN
885      ENDIF
886     
887      DEALLOCATE( iszi0, iszj0 )
888      inbij = inbijmax + 1        ! default: larger than possible
889      ii = isz0+1                 ! start from the end of the list (smaller subdomains)
890      DO WHILE( inbij > knbij )   ! while the number of ocean subdomains exceed the number of procs
891         ii = ii -1 
892         ALLOCATE( llisoce(inbi0(ii), inbj0(ii)) )
893         CALL mpp_init_isoce( inbi0(ii), inbj0(ii), llisoce )            ! must be done by all core
894         inbij = COUNT(llisoce)
895         DEALLOCATE( llisoce )
896      END DO
897      knbi = inbi0(ii)
898      knbj = inbj0(ii)
899      DEALLOCATE( inbi0, inbj0 )
900      !
901   END SUBROUTINE mpp_init_bestpartition
902   
903   
904   SUBROUTINE mpp_init_landprop( propland )
905      !!----------------------------------------------------------------------
906      !!                  ***  ROUTINE mpp_init_landprop  ***
907      !!
908      !! ** Purpose : the the proportion of land points in the surface land-sea mask
909      !!
910      !! ** Method  : read iproc strips (of length jpiglo) of the land-sea mask
911      !!----------------------------------------------------------------------
912      REAL(wp), INTENT(  out) :: propland    ! proportion of land points (between 0 and 1)
913      !
914      INTEGER, DIMENSION(jpni*jpnj) ::   kusedom_1d
915      INTEGER :: inboce 
916      INTEGER :: iproc, idiv, ijsz
917      INTEGER :: ijstr, ijend, ijcnt
918      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce
919      !!----------------------------------------------------------------------
920      ! do nothing if there is no land-sea mask
921      IF( numbot == -1 .and. numbdy == -1 ) THEN
922         propland = 0.
923         RETURN
924      ENDIF
925
926      ! number of processes reading the bathymetry file
927      iproc = MINVAL( (/mppsize, jpjglo/2, 100/) )  ! read a least 2 lines, no more that 100 processes reading at the same time
928     
929      ! we want to read iproc strips of the land-sea mask. -> pick up iproc processes among mppsize processes
930      IF( iproc == 1 ) THEN   ;   idiv = mppsize
931      ELSE                    ;   idiv = ( mppsize - 1 ) / ( iproc - 1 )
932      ENDIF
933      ijsz = jpjglo / iproc
934      IF( narea <= MOD(jpjglo,iproc) ) ijsz = ijsz + 1
935     
936      IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= iproc ) ) THEN
937         !
938         ijstr = (narea-1)*(jpjglo/iproc) + MIN(narea-1, MOD(jpjglo,iproc)) + 1
939         ijend = ijstr + ijsz - 1
940         ijcnt = ijend - ijstr + 1
941         !
942         ALLOCATE( lloce(jpiglo, ijcnt) )   ! allocate the strip
943         CALL mpp_init_readbot_strip( ijstr, ijcnt, lloce )
944         inboce = COUNT(lloce)
945         DEALLOCATE(lloce)
946         !
947      ELSE
948         inboce = 0
949      ENDIF
950      CALL mpp_sum( 'mppini', inboce )
951      !
952      propland = REAL( jpiglo*jpjglo - inboce, wp ) / REAL( jpiglo*jpjglo, wp ) 
953      !
954   END SUBROUTINE mpp_init_landprop
955   
956   
957   SUBROUTINE mpp_init_isoce( knbi, knbj, ldisoce )
958      !!----------------------------------------------------------------------
959      !!                  ***  ROUTINE mpp_init_nboce  ***
960      !!
961      !! ** Purpose : check for a mpi domain decomposition knbi x knbj which
962      !!              subdomains contain at least 1 ocean point
963      !!
964      !! ** Method  : read knbj strips (of length jpiglo) of the land-sea mask
965      !!----------------------------------------------------------------------
966      INTEGER,                       INTENT(in   ) ::   knbi, knbj
967      LOGICAL, DIMENSION(knbi,knbj), INTENT(  out) ::   ldisoce   !
968      !
969      INTEGER, DIMENSION(knbi,knbj) ::   inboce
970      INTEGER, DIMENSION(knbi*knbj) ::   inboce_1d
971      INTEGER :: idiv, i2read, inj
972      INTEGER :: iimax, ijmax
973      INTEGER :: ji,jj
974      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce
975      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ilci
976      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ilcj
977      !!----------------------------------------------------------------------
978      ! do nothing if there is no land-sea mask
979      IF( numbot == -1 .AND. numbdy == -1 ) THEN
980         ldisoce(:,:) = .TRUE.
981         RETURN
982      ENDIF
983
984      ! we want to read knbj strips of the land-sea mask. -> pick up knbj processes among mppsize processes
985      IF( knbj == 1 ) THEN   ;   idiv = mppsize
986      ELSE                   ;   idiv = ( mppsize - 1 ) / ( knbj - 1 )
987      ENDIF
988      inboce(:,:) = 0
989      IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= knbj ) ) THEN
990         !
991         ALLOCATE( iimppt(knbi,knbj), ijmppt(knbi,knbj), ilci(knbi,knbj), ilcj(knbi,knbj) )
992         CALL mpp_basic_decomposition( knbi, knbj, iimax, ijmax, iimppt, ijmppt, ilci, ilcj )
993         !
994         i2read = knbj / mppsize    ! strip number to be read by this process
995         IF( ( narea - 1 ) / idiv < MOD(knbj,mppsize) ) i2read = i2read + 1
996         DO jj = 1, i2read
997            ! strip number to be read (from 1 to knbj)
998            inj = ( narea - 1 ) * ( knbj / mppsize ) + MIN( MOD(knbj,mppsize), ( narea - 1 ) / idiv ) + jj
999            ALLOCATE( lloce(jpiglo, ilcj(1,inj)) )                              ! allocate the strip
1000            CALL mpp_init_readbot_strip( ijmppt(1,inj), ilcj(1,inj), lloce )   ! read the strip
1001            DO  ji = 1, knbi
1002               inboce(ji,inj) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ilci(ji,1)-1,:) )
1003            END DO
1004            DEALLOCATE(lloce)
1005         END DO
1006         !
1007         DEALLOCATE(iimppt, ijmppt, ilci, ilcj)
1008         !
1009      ENDIF
1010     
1011      inboce_1d = RESHAPE(inboce, (/ knbi*knbj /))
1012      CALL mpp_sum( 'mppini', inboce_1d )
1013      inboce = RESHAPE(inboce_1d, (/knbi, knbj/))
1014      ldisoce = inboce /= 0
1015      !
1016   END SUBROUTINE mpp_init_isoce
1017   
1018   
1019   SUBROUTINE mpp_init_readbot_strip( kjstr, kjcnt, ldoce )
1020      !!----------------------------------------------------------------------
1021      !!                  ***  ROUTINE mpp_init_readbot_strip  ***
1022      !!
1023      !! ** Purpose : Read relevant bathymetric information in order to
1024      !!              provide a land/sea mask used for the elimination
1025      !!              of land domains, in an mpp computation.
1026      !!
1027      !! ** Method  : read stipe of size (jpiglo,...)
1028      !!----------------------------------------------------------------------
1029      INTEGER                         , INTENT(in   ) :: kjstr       !
1030      INTEGER                         , INTENT(in   ) :: kjcnt       !
1031      LOGICAL, DIMENSION(jpiglo,kjcnt), INTENT(  out) :: ldoce       !
1032      !
1033      INTEGER                           ::   inumsave                     ! local logical unit
1034      REAL(wp), DIMENSION(jpiglo,kjcnt) ::   zbot, zbdy 
1035      !!----------------------------------------------------------------------
1036      !
1037      inumsave = numout   ;   numout = numnul   !   redirect all print to /dev/null
1038      !
1039      IF( numbot /= -1 ) THEN
1040         CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
1041      ELSE
1042         zbot(:,:) = 1.   ! put a non-null value
1043      ENDIF
1044
1045       IF( numbdy /= -1 ) THEN   ! Adjust  with bdy_msk if it exists   
1046         CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
1047         zbot(:,:) = zbot(:,:) * zbdy(:,:)
1048      ENDIF
1049      !
1050      ldoce = zbot > 0.
1051      numout = inumsave
1052      !
1053   END SUBROUTINE mpp_init_readbot_strip
1054
1055
1056   SUBROUTINE mpp_init_ioipsl
1057      !!----------------------------------------------------------------------
1058      !!                  ***  ROUTINE mpp_init_ioipsl  ***
1059      !!
1060      !! ** Purpose :   
1061      !!
1062      !! ** Method  :   
1063      !!
1064      !! History :
1065      !!   9.0  !  04-03  (G. Madec )  MPP-IOIPSL
1066      !!   " "  !  08-12  (A. Coward)  addition in case of jpni*jpnj < jpnij
1067      !!----------------------------------------------------------------------
1068      INTEGER, DIMENSION(2) ::   iglo, iloc, iabsf, iabsl, ihals, ihale, idid
1069      !!----------------------------------------------------------------------
1070
1071      ! The domain is split only horizontally along i- or/and j- direction
1072      ! So we need at the most only 1D arrays with 2 elements.
1073      ! Set idompar values equivalent to the jpdom_local_noextra definition
1074      ! used in IOM. This works even if jpnij .ne. jpni*jpnj.
1075      iglo(1) = jpiglo
1076      iglo(2) = jpjglo
1077      iloc(1) = nlci
1078      iloc(2) = nlcj
1079      iabsf(1) = nimppt(narea)
1080      iabsf(2) = njmppt(narea)
1081      iabsl(:) = iabsf(:) + iloc(:) - 1
1082      ihals(1) = nldi - 1
1083      ihals(2) = nldj - 1
1084      ihale(1) = nlci - nlei
1085      ihale(2) = nlcj - nlej
1086      idid(1) = 1
1087      idid(2) = 2
1088
1089      IF(lwp) THEN
1090          WRITE(numout,*)
1091          WRITE(numout,*) 'mpp_init_ioipsl :   iloc  = ', iloc (1), iloc (2)
1092          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf(1), iabsf(2)
1093          WRITE(numout,*) '                    ihals = ', ihals(1), ihals(2)
1094          WRITE(numout,*) '                    ihale = ', ihale(1), ihale(2)
1095      ENDIF
1096      !
1097      CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
1098      !
1099   END SUBROUTINE mpp_init_ioipsl 
1100
1101
1102   SUBROUTINE mpp_init_nfdcom
1103      !!----------------------------------------------------------------------
1104      !!                     ***  ROUTINE  mpp_init_nfdcom  ***
1105      !! ** Purpose :   Setup for north fold exchanges with explicit
1106      !!                point-to-point messaging
1107      !!
1108      !! ** Method :   Initialization of the northern neighbours lists.
1109      !!----------------------------------------------------------------------
1110      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
1111      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
1112      !!----------------------------------------------------------------------
1113      INTEGER  ::   sxM, dxM, sxT, dxT, jn
1114      INTEGER  ::   njmppmax
1115      !!----------------------------------------------------------------------
1116      !
1117      njmppmax = MAXVAL( njmppt )
1118      !
1119      !initializes the north-fold communication variables
1120      isendto(:) = 0
1121      nsndto     = 0
1122      !
1123      IF ( njmpp == njmppmax ) THEN      ! if I am a process in the north
1124         !
1125         !sxM is the first point (in the global domain) needed to compute the north-fold for the current process
1126         sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
1127         !dxM is the last point (in the global domain) needed to compute the north-fold for the current process
1128         dxM = jpiglo - nimppt(narea) + 2
1129         !
1130         ! loop over the other north-fold processes to find the processes
1131         ! managing the points belonging to the sxT-dxT range
1132         !
1133         DO jn = 1, jpni
1134            !
1135            sxT = nfiimpp(jn, jpnj)                            ! sxT = 1st  point (in the global domain) of the jn process
1136            dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1    ! dxT = last point (in the global domain) of the jn process
1137            !
1138            IF    ( sxT < sxM  .AND.  sxM < dxT ) THEN
1139               nsndto          = nsndto + 1
1140               isendto(nsndto) = jn
1141            ELSEIF( sxM <= sxT  .AND.  dxM >= dxT ) THEN
1142               nsndto          = nsndto + 1
1143               isendto(nsndto) = jn
1144            ELSEIF( dxM <  dxT  .AND.  sxT <  dxM ) THEN
1145               nsndto          = nsndto + 1
1146               isendto(nsndto) = jn
1147            ENDIF
1148            !
1149         END DO
1150         nfsloop = 1
1151         nfeloop = nlci
1152         DO jn = 2,jpni-1
1153            IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
1154               IF( nfipproc(jn-1,jpnj) == -1 )   nfsloop = nldi
1155               IF( nfipproc(jn+1,jpnj) == -1 )   nfeloop = nlei
1156            ENDIF
1157         END DO
1158         !
1159      ENDIF
1160      l_north_nogather = .TRUE.
1161      !
1162   END SUBROUTINE mpp_init_nfdcom
1163
1164
1165#endif
1166
1167   !!======================================================================
1168END MODULE mppini
Note: See TracBrowser for help on using the repository browser.