source: NEMO/releases/release-4.0/src/OCE/LBC/mppini.F90 @ 10616

Last change on this file since 10616 was 10616, checked in by smasson, 19 months ago

v4: improve information/error messages in mppini

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