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

Last change on this file since 10397 was 10397, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 5 introduce mpp_delay_sum, 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, ALLOCATABLE, DIMENSION(:)     ::   iin, ii_nono, ii_noea          ! 1D workspace
152      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ijn, ii_noso, ii_nowe          !  -     -
153      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ilci, ibondi, ipproc   ! 2D workspace
154      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ilcj, ibondj, ipolj    !  -     -
155      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilei, ildi, iono, ioea         !  -     -
156      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilej, ildj, ioso, iowe         !  -     -
157      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   llisoce                        !  -     -
158      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           &
159           &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
160           &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             & 
161           &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
162           &             cn_ice, nn_ice_dta,                                     &
163           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
164           &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
165      !!----------------------------------------------------------------------
166
167      ! do we need to take into account bdy_msk?
168      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY
169      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
170903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp )
171      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY
172      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
173904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp )
174      !
175      IF(               ln_read_cfg ) CALL iom_open( cn_domcfg,    numbot )
176      IF( ln_bdy .AND. ln_mask_file ) CALL iom_open( cn_mask_file, numbdy )
177      !
178      !  1. Dimension arrays for subdomains
179      ! -----------------------------------
180      !
181      ! If dimensions of processor grid weren't specified in the namelist file
182      ! then we calculate them here now that we have our communicator size
183      IF( jpni < 1 .OR. jpnj < 1 )   CALL mpp_init_bestpartition( mppsize, jpni, jpnj )
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 chossen domain decomposition ', jpni, jpnj
192!!$            WRITE(ctmp2,*) '   exceeds the maximum number of ocean subdomains = ', inijmin
193!!$            WRITE(ctmp3,*) '   we suppressed ', jpni*jpnj - mppsize, ' land subdomains '
194!!$            WRITE(ctmp4,*) '   BUT we had to keep ', mppsize - inijmin, ' land subdomains that are useless...'
195!!$            CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 )
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), INTENT(  out) ::   kimppt, kjmppt
658      INTEGER, DIMENSION(knbi,knbj), 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      jpimax = ( nx_global+2-2*nn_hls + (knbi-1) ) / knbi + 2*nn_hls    ! first  dim.
667      jpjmax = ( 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      !
673      !  1. Dimension arrays for subdomains
674      ! -----------------------------------
675      !  Computation of local domain sizes klci() klcj()
676      !  These dimensions depend on global sizes knbi,knbj and jpiglo,jpjglo
677      !  The subdomains are squares lesser than or equal to the global
678      !  dimensions divided by the number of processors minus the overlap array.
679      !
680      ireci = 2 * nn_hls
681      irecj = 2 * nn_hls
682      iresti = 1 + MOD( jpiglo - ireci -1 , knbi )
683      irestj = 1 + MOD( jpjglo - irecj -1 , knbj )
684      !
685      !  Need to use kimax and kjmax here since jpi and jpj not yet defined
686#if defined key_nemocice_decomp
687      ! Change padding to be consistent with CICE
688      klci(1:knbi-1      ,:) = kimax
689      klci(knbi          ,:) = jpiglo - (knbi - 1) * (kimax - nreci)
690      klcj(:,      1:knbj-1) = kjmax
691      klcj(:,          knbj) = jpjglo - (knbj - 1) * (kjmax - nrecj)
692#else
693      klci(1:iresti      ,:) = kimax
694      klci(iresti+1:knbi ,:) = kimax-1
695      klcj(:,      1:irestj) = kjmax
696      klcj(:, irestj+1:knbj) = kjmax-1
697#endif
698
699      !  2. Index arrays for subdomains
700      ! -------------------------------
701      kimppt(:,:) = 1
702      kjmppt(:,:) = 1
703      !
704      IF( knbi > 1 ) THEN
705         DO jj = 1, knbj
706            DO ji = 2, knbi
707               kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - ireci
708            END DO
709         END DO
710      ENDIF
711      !
712      IF( knbj > 1 )THEN
713         DO jj = 2, knbj
714            DO ji = 1, knbi
715               kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - irecj
716            END DO
717         END DO
718      ENDIF
719     
720   END SUBROUTINE mpp_basic_decomposition
721
722
723   SUBROUTINE mpp_init_bestpartition( knbij, knbi, knbj, ldlist )
724      !!----------------------------------------------------------------------
725      !!                 ***  ROUTINE mpp_init_bestpartition  ***
726      !!
727      !! ** Purpose :
728      !!
729      !! ** Method  :
730      !!----------------------------------------------------------------------
731      INTEGER,           INTENT(in   ) ::   knbij         ! total number if subdomains               (knbi*knbj)
732      INTEGER, OPTIONAL, INTENT(  out) ::   knbi, knbj    ! number if subdomains along i and j (knbi and knbj)
733      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldlist        ! .true.: print the list the best domain decompositions (with land)
734      !
735      INTEGER :: ji, jj, ii, iitarget
736      INTEGER :: iszitst, iszjtst
737      INTEGER :: isziref, iszjref
738      INTEGER :: inbij, iszij
739      INTEGER :: inbimax, inbjmax, inbijmax
740      INTEGER :: isz0, isz1
741      INTEGER, DIMENSION(  :), ALLOCATABLE :: indexok
742      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi0, inbj0, inbij0   ! number of subdomains along i,j
743      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi0, iszj0, iszij0   ! max size of the subdomains along i,j
744      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi1, inbj1, inbij1   ! number of subdomains along i,j
745      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi1, iszj1, iszij1   ! max size of the subdomains along i,j
746      LOGICAL :: llist
747      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d                 ! max size of the subdomains along i,j
748      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisoce              !  -     -
749      REAL(wp)::   zpropland
750      !!----------------------------------------------------------------------
751      !
752      llist = .FALSE.
753      IF( PRESENT(ldlist) ) llist = ldlist
754
755      CALL mpp_init_landprop( zpropland )                      ! get the proportion of land point over the gloal domain
756      inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) )    ! define the largest possible value for jpni*jpnj
757      !
758      IF( llist ) THEN   ;   inbijmax = inbij*2
759      ELSE               ;   inbijmax = inbij
760      ENDIF
761      !
762      ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax))
763      !
764      inbimax = 0
765      inbjmax = 0
766      isziref = jpiglo*jpjglo+1
767      iszjref = jpiglo*jpjglo+1
768      !
769      ! get the list of knbi that gives a smaller jpimax than knbi-1
770      ! get the list of knbj that gives a smaller jpjmax than knbj-1
771      DO ji = 1, inbijmax     
772#if defined key_nemocice_decomp
773         iszitst = ( nx_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
774#else
775         iszitst = ( jpiglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls
776#endif
777         IF( iszitst < isziref ) THEN
778            isziref = iszitst
779            inbimax = inbimax + 1
780            inbi0(inbimax) = ji
781            iszi0(inbimax) = isziref
782         ENDIF
783#if defined key_nemocice_decomp
784         iszjtst = ( ny_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
785#else
786         iszjtst = ( jpjglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls
787#endif
788         IF( iszjtst < iszjref ) THEN
789            iszjref = iszjtst
790            inbjmax = inbjmax + 1
791            inbj0(inbjmax) = ji
792            iszj0(inbjmax) = iszjref
793         ENDIF
794      END DO
795
796      ! combine these 2 lists to get all possible knbi*knbj <  inbijmax
797      ALLOCATE( llmsk2d(inbimax,inbjmax) )
798      DO jj = 1, inbjmax
799         DO ji = 1, inbimax
800            IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN   ;   llmsk2d(ji,jj) = .TRUE.
801            ELSE                                            ;   llmsk2d(ji,jj) = .FALSE.
802            ENDIF
803         END DO
804      END DO
805      isz1 = COUNT(llmsk2d)
806      ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) )
807      ii = 0
808      DO jj = 1, inbjmax
809         DO ji = 1, inbimax
810            IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN
811               ii = ii + 1
812               inbi1(ii) = inbi0(ji)
813               inbj1(ii) = inbj0(jj)
814               iszi1(ii) = iszi0(ji)
815               iszj1(ii) = iszj0(jj)
816            END IF
817         END DO
818      END DO
819      DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
820      DEALLOCATE( llmsk2d )
821
822      ALLOCATE( inbij1(isz1), iszij1(isz1) )
823      inbij1(:) = inbi1(:) * inbj1(:)
824      iszij1(:) = iszi1(:) * iszj1(:)
825
826      ! if therr is no land and no print
827      IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN
828         ! get the smaller partition which gives the smallest subdomain size
829         ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1)
830         knbi = inbi1(ii)
831         knbj = inbj1(ii)
832         DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 )
833         RETURN
834      ENDIF
835
836      ! extract only the partitions which reduce the subdomain size in comparison with smaller partitions
837      ALLOCATE( indexok(isz1) )                                 ! to store indices of the best partitions
838      isz0 = 0                                                  ! number of best partitions     
839      inbij = 1                                                 ! start with the min value of inbij1 => 1
840      iszij = jpiglo*jpjglo+1                                   ! default: larger than global domain
841      DO WHILE( inbij <= inbijmax )                             ! if we did not reach the max of inbij1
842         ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1)   ! warning: send back the first occurence if multiple results
843         IF ( iszij1(ii) < iszij ) THEN
844            isz0 = isz0 + 1
845            indexok(isz0) = ii
846            iszij = iszij1(ii)
847         ENDIF
848         inbij = MINVAL(inbij1, mask = inbij1 > inbij)   ! warning: return largest integer value if mask = .false. everywhere
849      END DO
850      DEALLOCATE( inbij1, iszij1 )
851
852      ! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size)
853      ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) )
854      DO ji = 1, isz0
855         ii = indexok(ji)
856         inbi0(ji) = inbi1(ii)
857         inbj0(ji) = inbj1(ii)
858         iszi0(ji) = iszi1(ii)
859         iszj0(ji) = iszj1(ii)
860      END DO
861      DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 )
862
863      IF( llist ) THEN  ! we print about 21 best partitions
864         IF(lwp) THEN
865            WRITE(numout,*)
866            WRITE(numout,         *) '                  For your information:'
867            WRITE(numout,'(a,i5,a)') '  list of the best partitions around ',   knbij, ' mpi processes'
868            WRITE(numout,         *) '  --------------------------------------', '-----', '--------------'
869            WRITE(numout,*)
870         END IF
871         iitarget = MINLOC( inbi0(:)*inbj0(:), mask = inbi0(:)*inbj0(:) >= knbij, dim = 1 )
872         DO ji = MAX(1,iitarget-10), MIN(isz0,iitarget+10)
873            ALLOCATE( llisoce(inbi0(ji), inbj0(ji)) )
874            CALL mpp_init_isoce( inbi0(ji), inbj0(ji), llisoce ) ! Warning: must be call by all cores (call mpp_sum)
875            inbij = COUNT(llisoce)
876            DEALLOCATE( llisoce )
877            IF(lwp) WRITE(numout,'(a, i5, a, i5, a, i4, a, i4, a, i9, a, i5, a, i5, a)')    &
878               &     'nb_cores ' , inbij,' oce + ', inbi0(ji)*inbj0(ji) - inbij             &
879               &                                , ' land ( ', inbi0(ji),' x ', inbj0(ji),   &
880               & ' ), nb_points ', iszi0(ji)*iszj0(ji),' ( ', iszi0(ji),' x ', iszj0(ji),' )'
881         END DO
882         DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
883         RETURN
884      ENDIF
885     
886      DEALLOCATE( iszi0, iszj0 )
887      inbij = inbijmax + 1        ! default: larger than possible
888      ii = isz0+1                 ! start from the end of the list (smaller subdomains)
889      DO WHILE( inbij > knbij )   ! while the number of ocean subdomains exceed the number of procs
890         ii = ii -1 
891         ALLOCATE( llisoce(inbi0(ii), inbj0(ii)) )
892         CALL mpp_init_isoce( inbi0(ii), inbj0(ii), llisoce )            ! must be done by all core
893         inbij = COUNT(llisoce)
894         DEALLOCATE( llisoce )
895      END DO
896      knbi = inbi0(ii)
897      knbj = inbj0(ii)
898      DEALLOCATE( inbi0, inbj0 )
899      !
900   END SUBROUTINE mpp_init_bestpartition
901   
902   
903   SUBROUTINE mpp_init_landprop( propland )
904      !!----------------------------------------------------------------------
905      !!                  ***  ROUTINE mpp_init_landprop  ***
906      !!
907      !! ** Purpose : the the proportion of land points in the surface land-sea mask
908      !!
909      !! ** Method  : read iproc strips (of length jpiglo) of the land-sea mask
910      !!----------------------------------------------------------------------
911      REAL(wp), INTENT(  out) :: propland    ! proportion of land points (between 0 and 1)
912      !
913      INTEGER, DIMENSION(jpni*jpnj) ::   kusedom_1d
914      INTEGER :: inboce 
915      INTEGER :: iproc, idiv, ijsz
916      INTEGER :: ijstr, ijend, ijcnt
917      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce
918      !!----------------------------------------------------------------------
919      ! do nothing if there is no land-sea mask
920      IF( numbot == -1 .and. numbdy == -1 ) THEN
921         propland = 0.
922         RETURN
923      ENDIF
924
925      ! number of processes reading the bathymetry file
926      iproc = MINVAL( (/mppsize, jpjglo/2, 100/) )  ! read a least 2 lines, no more that 100 processes reading at the same time
927     
928      ! we want to read iproc strips of the land-sea mask. -> pick up iproc processes among mppsize processes
929      IF( iproc == 1 ) THEN   ;   idiv = mppsize
930      ELSE                    ;   idiv = ( mppsize - 1 ) / ( iproc - 1 )
931      ENDIF
932      ijsz = jpjglo / iproc
933      IF( narea <= MOD(jpjglo,iproc) ) ijsz = ijsz + 1
934     
935      IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= iproc ) ) THEN
936         !
937         ijstr = (narea-1)*(jpjglo/iproc) + MIN(narea-1, MOD(jpjglo,iproc)) + 1
938         ijend = ijstr + ijsz - 1
939         ijcnt = ijend - ijstr + 1
940         !
941         ALLOCATE( lloce(jpiglo, ijcnt) )   ! allocate the strip
942         CALL mpp_init_readbot_strip( ijstr, ijcnt, lloce )
943         inboce = COUNT(lloce)
944         DEALLOCATE(lloce)
945         !
946      ELSE
947         inboce = 0
948      ENDIF
949      CALL mpp_sum( 'mppini', inboce )
950      !
951      propland = REAL( jpiglo*jpjglo - inboce, wp ) / REAL( jpiglo*jpjglo, wp ) 
952      !
953   END SUBROUTINE mpp_init_landprop
954   
955   
956   SUBROUTINE mpp_init_isoce( knbi, knbj, ldisoce )
957      !!----------------------------------------------------------------------
958      !!                  ***  ROUTINE mpp_init_nboce  ***
959      !!
960      !! ** Purpose : check for a mpi domain decomposition knbi x knbj which
961      !!              subdomains contain at least 1 ocean point
962      !!
963      !! ** Method  : read knbj strips (of length jpiglo) of the land-sea mask
964      !!----------------------------------------------------------------------
965      INTEGER,                       INTENT(in   ) ::   knbi, knbj
966      LOGICAL, DIMENSION(knbi,knbj), INTENT(  out) ::   ldisoce   !
967      !
968      INTEGER, DIMENSION(knbi,knbj) ::   inboce
969      INTEGER, DIMENSION(knbi*knbj) ::   inboce_1d
970      INTEGER :: idiv, i2read, inj
971      INTEGER :: iimax, ijmax
972      INTEGER :: ji,jj
973      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce
974      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ilci
975      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ilcj
976      !!----------------------------------------------------------------------
977      ! do nothing if there is no land-sea mask
978      IF( numbot == -1 .AND. numbdy == -1 ) THEN
979         ldisoce(:,:) = .TRUE.
980         RETURN
981      ENDIF
982
983      ! we want to read knbj strips of the land-sea mask. -> pick up knbj processes among mppsize processes
984      IF( knbj == 1 ) THEN   ;   idiv = mppsize
985      ELSE                   ;   idiv = ( mppsize - 1 ) / ( knbj - 1 )
986      ENDIF
987      inboce(:,:) = 0
988      IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= knbj ) ) THEN
989         !
990         ALLOCATE( iimppt(knbi,knbj), ijmppt(knbi,knbj), ilci(knbi,knbj), ilcj(knbi,knbj) )
991         CALL mpp_basic_decomposition( knbi, knbj, iimax, ijmax, iimppt, ijmppt, ilci, ilcj )
992         !
993         i2read = knbj / mppsize    ! strip number to be read by this process
994         IF( ( narea - 1 ) / idiv < MOD(knbj,mppsize) ) i2read = i2read + 1
995         DO jj = 1, i2read
996            ! strip number to be read (from 1 to knbj)
997            inj = ( narea - 1 ) * ( knbj / mppsize ) + MIN( MOD(knbj,mppsize), ( narea - 1 ) / idiv ) + jj
998            ALLOCATE( lloce(jpiglo, ilcj(1,inj)) )                              ! allocate the strip
999            CALL mpp_init_readbot_strip( ijmppt(1,inj), ilcj(1,inj), lloce )   ! read the strip
1000            DO  ji = 1, knbi
1001               inboce(ji,inj) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ilci(ji,1)-1,:) )
1002            END DO
1003            DEALLOCATE(lloce)
1004         END DO
1005         !
1006         DEALLOCATE(iimppt, ijmppt, ilci, ilcj)
1007         !
1008      ENDIF
1009     
1010      inboce_1d = RESHAPE(inboce, (/ knbi*knbj /))
1011      CALL mpp_sum( 'mppini', inboce_1d )
1012      inboce = RESHAPE(inboce_1d, (/knbi, knbj/))
1013      ldisoce = inboce /= 0
1014      !
1015   END SUBROUTINE mpp_init_isoce
1016   
1017   
1018   SUBROUTINE mpp_init_readbot_strip( kjstr, kjcnt, ldoce )
1019      !!----------------------------------------------------------------------
1020      !!                  ***  ROUTINE mpp_init_readbot_strip  ***
1021      !!
1022      !! ** Purpose : Read relevant bathymetric information in order to
1023      !!              provide a land/sea mask used for the elimination
1024      !!              of land domains, in an mpp computation.
1025      !!
1026      !! ** Method  : read stipe of size (jpiglo,...)
1027      !!----------------------------------------------------------------------
1028      INTEGER                         , INTENT(in   ) :: kjstr       !
1029      INTEGER                         , INTENT(in   ) :: kjcnt       !
1030      LOGICAL, DIMENSION(jpiglo,kjcnt), INTENT(  out) :: ldoce       !
1031      !
1032      INTEGER                           ::   inumsave                     ! local logical unit
1033      REAL(wp), DIMENSION(jpiglo,kjcnt) ::   zbot, zbdy 
1034      !!----------------------------------------------------------------------
1035      !
1036      inumsave = numout   ;   numout = numnul   !   redirect all print to /dev/null
1037      !
1038      IF( numbot /= -1 ) THEN
1039         CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
1040      ELSE
1041         zbot(:,:) = 1.   ! put a non-null value
1042      ENDIF
1043
1044       IF( numbdy /= -1 ) THEN   ! Adjust  with bdy_msk if it exists   
1045         CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
1046         zbot(:,:) = zbot(:,:) * zbdy(:,:)
1047      ENDIF
1048      !
1049      ldoce = zbot > 0.
1050      numout = inumsave
1051      !
1052   END SUBROUTINE mpp_init_readbot_strip
1053
1054
1055   SUBROUTINE mpp_init_ioipsl
1056      !!----------------------------------------------------------------------
1057      !!                  ***  ROUTINE mpp_init_ioipsl  ***
1058      !!
1059      !! ** Purpose :   
1060      !!
1061      !! ** Method  :   
1062      !!
1063      !! History :
1064      !!   9.0  !  04-03  (G. Madec )  MPP-IOIPSL
1065      !!   " "  !  08-12  (A. Coward)  addition in case of jpni*jpnj < jpnij
1066      !!----------------------------------------------------------------------
1067      INTEGER, DIMENSION(2) ::   iglo, iloc, iabsf, iabsl, ihals, ihale, idid
1068      !!----------------------------------------------------------------------
1069
1070      ! The domain is split only horizontally along i- or/and j- direction
1071      ! So we need at the most only 1D arrays with 2 elements.
1072      ! Set idompar values equivalent to the jpdom_local_noextra definition
1073      ! used in IOM. This works even if jpnij .ne. jpni*jpnj.
1074      iglo(1) = jpiglo
1075      iglo(2) = jpjglo
1076      iloc(1) = nlci
1077      iloc(2) = nlcj
1078      iabsf(1) = nimppt(narea)
1079      iabsf(2) = njmppt(narea)
1080      iabsl(:) = iabsf(:) + iloc(:) - 1
1081      ihals(1) = nldi - 1
1082      ihals(2) = nldj - 1
1083      ihale(1) = nlci - nlei
1084      ihale(2) = nlcj - nlej
1085      idid(1) = 1
1086      idid(2) = 2
1087
1088      IF(lwp) THEN
1089          WRITE(numout,*)
1090          WRITE(numout,*) 'mpp_init_ioipsl :   iloc  = ', iloc (1), iloc (2)
1091          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf(1), iabsf(2)
1092          WRITE(numout,*) '                    ihals = ', ihals(1), ihals(2)
1093          WRITE(numout,*) '                    ihale = ', ihale(1), ihale(2)
1094      ENDIF
1095      !
1096      CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
1097      !
1098   END SUBROUTINE mpp_init_ioipsl 
1099
1100
1101   SUBROUTINE mpp_init_nfdcom
1102      !!----------------------------------------------------------------------
1103      !!                     ***  ROUTINE  mpp_init_nfdcom  ***
1104      !! ** Purpose :   Setup for north fold exchanges with explicit
1105      !!                point-to-point messaging
1106      !!
1107      !! ** Method :   Initialization of the northern neighbours lists.
1108      !!----------------------------------------------------------------------
1109      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
1110      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
1111      !!----------------------------------------------------------------------
1112      INTEGER  ::   sxM, dxM, sxT, dxT, jn
1113      INTEGER  ::   njmppmax
1114      !!----------------------------------------------------------------------
1115      !
1116      njmppmax = MAXVAL( njmppt )
1117      !
1118      !initializes the north-fold communication variables
1119      isendto(:) = 0
1120      nsndto     = 0
1121      !
1122      IF ( njmpp == njmppmax ) THEN      ! if I am a process in the north
1123         !
1124         !sxM is the first point (in the global domain) needed to compute the north-fold for the current process
1125         sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
1126         !dxM is the last point (in the global domain) needed to compute the north-fold for the current process
1127         dxM = jpiglo - nimppt(narea) + 2
1128         !
1129         ! loop over the other north-fold processes to find the processes
1130         ! managing the points belonging to the sxT-dxT range
1131         !
1132         DO jn = 1, jpni
1133            !
1134            sxT = nfiimpp(jn, jpnj)                            ! sxT = 1st  point (in the global domain) of the jn process
1135            dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1    ! dxT = last point (in the global domain) of the jn process
1136            !
1137            IF    ( sxT < sxM  .AND.  sxM < dxT ) THEN
1138               nsndto          = nsndto + 1
1139               isendto(nsndto) = jn
1140            ELSEIF( sxM <= sxT  .AND.  dxM >= dxT ) THEN
1141               nsndto          = nsndto + 1
1142               isendto(nsndto) = jn
1143            ELSEIF( dxM <  dxT  .AND.  sxT <  dxM ) THEN
1144               nsndto          = nsndto + 1
1145               isendto(nsndto) = jn
1146            ENDIF
1147            !
1148         END DO
1149         nfsloop = 1
1150         nfeloop = nlci
1151         DO jn = 2,jpni-1
1152            IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
1153               IF( nfipproc(jn-1,jpnj) == -1 )   nfsloop = nldi
1154               IF( nfipproc(jn+1,jpnj) == -1 )   nfeloop = nlei
1155            ENDIF
1156         END DO
1157         !
1158      ENDIF
1159      l_north_nogather = .TRUE.
1160      !
1161   END SUBROUTINE mpp_init_nfdcom
1162
1163
1164#endif
1165
1166   !!======================================================================
1167END MODULE mppini
Note: See TracBrowser for help on using the repository browser.