source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/LBC/mppini.F90 @ 12586

Last change on this file since 12586 was 12586, checked in by francesca, 9 months ago

Add extra-halo support (jperio 3,4) - ticket #2366

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