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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
mppini.F90 in NEMO/trunk/src/OCE/LBC – NEMO

source: NEMO/trunk/src/OCE/LBC/mppini.F90

Last change on this file was 15302, checked in by smasson, 3 years ago

trunk: minor improvment in pt2pt (east-west comm only send/recv inner values)

  • Property svn:keywords set to Id
File size: 71.6 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 init_nfdcom
11   !!            3.   !  2013-06  (I. Epicoco, S. Mocavero, CMCC)  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   !!      init_ioipsl: IOIPSL initialization in mpp
19   !!      init_nfdcom: Setup for north fold exchanges with explicit point-to-point messaging
20   !!      init_doloop: set the starting/ending indices of DO-loop used in do_loop_substitute
21   !!----------------------------------------------------------------------
22   USE dom_oce        ! ocean space and time domain
23   USE bdy_oce        ! open BounDarY
24   !
25   USE lbcnfd         ! Setup of north fold exchanges
26   USE lib_mpp        ! distribued memory computing library
27   USE iom            ! nemo I/O library
28   USE ioipsl         ! I/O IPSL library
29   USE in_out_manager ! I/O Manager
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   mpp_init       ! called by nemogcm.F90
35   PUBLIC   mpp_getnum     ! called by prtctl
36   PUBLIC   mpp_basesplit  ! called by prtctl
37   PUBLIC   mpp_is_ocean   ! called by prtctl
38
39   INTEGER ::   numbot = -1   ! 'bottom_level' local logical unit
40   INTEGER ::   numbdy = -1   ! 'bdy_msk'      local logical unit
41
42   !!----------------------------------------------------------------------
43   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
44   !! $Id$
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49#if defined key_mpi_off
50   !!----------------------------------------------------------------------
51   !!   Default option :                            shared memory computing
52   !!----------------------------------------------------------------------
53
54   SUBROUTINE mpp_init
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE mpp_init  ***
57      !!
58      !! ** Purpose :   Lay out the global domain over processors.
59      !!
60      !! ** Method  :   Shared memory computing, set the local processor
61      !!              variables to the value of the global domain
62      !!----------------------------------------------------------------------
63      !
64      nn_comm = 1
65      nn_hls  = 1
66      jpiglo  = Ni0glo + 2 * nn_hls
67      jpjglo  = Nj0glo + 2 * nn_hls
68      jpimax  = jpiglo
69      jpjmax  = jpjglo
70      jpi     = jpiglo
71      jpj     = jpjglo
72      jpk     = MAX( 2, jpkglo )
73      jpij    = jpi*jpj
74      jpni    = 1
75      jpnj    = 1
76      jpnij   = jpni*jpnj
77      nimpp   = 1
78      njmpp   = 1
79      nidom   = FLIO_DOM_NONE
80      !
81      mpiSnei(:,:) = -1
82      mpiRnei(:,:) = -1
83      l_SelfPerio(1:2) = l_Iperio                  !  west,  east periodicity by itself
84      l_SelfPerio(3:4) = l_Jperio                  ! south, north periodicity by itself
85      l_SelfPerio(5:8) = l_Iperio .AND. l_Jperio   ! corners bi-periodicity by itself
86      l_IdoNFold = l_NFold                         ! is this process doing North fold?
87      !
88      CALL init_doloop                       ! set start/end indices or do-loop depending on the halo width value (nn_hls)
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,*) '     njmpp  = ', njmpp
96      ENDIF
97      !
98#if defined key_agrif
99      call agrif_nemo_init()
100#endif
101   END SUBROUTINE mpp_init
102
103#else
104   !!----------------------------------------------------------------------
105   !!                   MPI massively parallel processing
106   !!----------------------------------------------------------------------
107
108
109   SUBROUTINE mpp_init
110      !!----------------------------------------------------------------------
111      !!                  ***  ROUTINE mpp_init  ***
112      !!
113      !! ** Purpose :   Lay out the global domain over processors.
114      !!      If land processors are to be eliminated, this program requires the
115      !!      presence of the domain configuration file. Land processors elimination
116      !!      is performed if jpni x jpnj /= jpnij. In this case, using the MPP_PREP
117      !!      preprocessing tool, help for defining the best cutting out.
118      !!
119      !! ** Method  :   Global domain is distributed in smaller local domains.
120      !!      Periodic condition is a function of the local domain position
121      !!      (global boundary or neighbouring domain) and of the global periodic
122      !!
123      !! ** Action : - set domain parameters
124      !!                    nimpp     : longitudinal index
125      !!                    njmpp     : latitudinal  index
126      !!                    narea     : number for local area
127      !!                    mpinei    : number of neighboring domains (starting at 0, -1 if no neighbourg)
128      !!----------------------------------------------------------------------
129      INTEGER ::   ji, jj, jn, jp, jh
130      INTEGER ::   ii, ij, ii2, ij2
131      INTEGER ::   inijmin   ! number of oce subdomains
132      INTEGER ::   inum, inum0
133      INTEGER ::   ifreq, il1, imil, il2, ijm1
134      INTEGER ::   ierr, ios
135      INTEGER ::   inbi, inbj, iimax, ijmax, icnt1, icnt2
136      INTEGER, DIMENSION(16*n_hlsmax) :: ichanged
137      INTEGER, ALLOCATABLE, DIMENSION(:    ) ::   iin, ijn
138      INTEGER, ALLOCATABLE, DIMENSION(:,:  ) ::   iimppt, ijpi, ipproc
139      INTEGER, ALLOCATABLE, DIMENSION(:,:  ) ::   ijmppt, ijpj
140      INTEGER, ALLOCATABLE, DIMENSION(:,:  ) ::   impi
141      INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::   inei
142      LOGICAL ::   llbest, llauto
143      LOGICAL ::   llwrtlay
144      LOGICAL ::   llmpi_Iperio, llmpi_Jperio, llmpiNFold
145      LOGICAL ::   ln_listonly
146      LOGICAL, ALLOCATABLE, DIMENSION(:,:  ) ::   llisOce  ! is not land-domain only?
147      LOGICAL, ALLOCATABLE, DIMENSION(:,:,:) ::   llnei    ! are neighbourgs existing?
148      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           &
149           &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
150           &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &
151           &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
152           &             cn_ice, nn_ice_dta,                                     &
153           &             ln_vol, nn_volctl, nn_rimwidth
154      NAMELIST/nammpp/ jpni, jpnj, nn_hls, ln_nnogather, ln_listonly, nn_comm
155      !!----------------------------------------------------------------------
156      !
157      llwrtlay = lwm .OR. sn_cfctl%l_layout
158      !
159      !  0. read namelists parameters
160      ! -----------------------------------
161      !
162      READ  ( numnam_ref, nammpp, IOSTAT = ios, ERR = 901 )
163901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nammpp in reference namelist' )
164      READ  ( numnam_cfg, nammpp, IOSTAT = ios, ERR = 902 )
165902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nammpp in configuration namelist' )
166      !
167      nn_hls = MAX(1, nn_hls)   ! nn_hls must be > 0
168      IF(lwp) THEN
169            WRITE(numout,*) '   Namelist nammpp'
170         IF( jpni < 1 .OR. jpnj < 1 ) THEN
171            WRITE(numout,*) '      jpni and jpnj will be calculated automatically'
172         ELSE
173            WRITE(numout,*) '      processor grid extent in i                            jpni = ', jpni
174            WRITE(numout,*) '      processor grid extent in j                            jpnj = ', jpnj
175         ENDIF
176            WRITE(numout,*) '      avoid use of mpi_allgather at the north fold  ln_nnogather = ', ln_nnogather
177            WRITE(numout,*) '      halo width (applies to both rows and columns)       nn_hls = ', nn_hls
178      ENDIF
179      !
180      IF(lwm)   WRITE( numond, nammpp )
181      !
182      jpiglo = Ni0glo + 2 * nn_hls
183      jpjglo = Nj0glo + 2 * nn_hls
184      !
185      ! do we need to take into account bdy_msk?
186      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
187903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)' )
188      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
189904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)' )
190      !
191      IF(               ln_read_cfg ) CALL iom_open( cn_domcfg,    numbot )
192      IF( ln_bdy .AND. ln_mask_file ) CALL iom_open( cn_mask_file, numbdy )
193      !
194      IF( ln_listonly )   CALL bestpartition( MAX(mppsize,jpni*jpnj), ldlist = .TRUE. )   ! must be done by all core
195      !
196      !  1. Dimension arrays for subdomains
197      ! -----------------------------------
198      !
199      ! If dimensions of MPI processes grid weren't specified in the namelist file
200      ! then we calculate them here now that we have our communicator size
201      IF(lwp) THEN
202         WRITE(numout,*)
203         WRITE(numout,*) 'mpp_init:'
204         WRITE(numout,*) '~~~~~~~~ '
205      ENDIF
206      IF( jpni < 1 .OR. jpnj < 1 ) THEN
207         CALL bestpartition( mppsize, jpni, jpnj )           ! best mpi decomposition for mppsize mpi processes
208         llauto = .TRUE.
209         llbest = .TRUE.
210      ELSE
211         llauto = .FALSE.
212         CALL bestpartition( mppsize, inbi, inbj, icnt2 )    ! best mpi decomposition for mppsize mpi processes
213         ! largest subdomain size for mpi decoposition jpni*jpnj given in the namelist
214         CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, jpni, jpnj, jpimax, jpjmax )
215         ! largest subdomain size for mpi decoposition inbi*inbj given by bestpartition
216         CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, inbi, inbj,  iimax,  ijmax )
217         icnt1 = jpni*jpnj - mppsize   ! number of land subdomains that should be removed to use mppsize mpi processes
218         IF(lwp) THEN
219            WRITE(numout,9000) '   The chosen domain decomposition ', jpni, ' x ', jpnj, ' with ', icnt1, ' land subdomains'
220            WRITE(numout,9002) '      - uses a total of ',  mppsize,' mpi process'
221            WRITE(numout,9000) '      - has mpi subdomains with a maximum size of (jpi = ', jpimax, ', jpj = ', jpjmax,   &
222               &                                                                ', jpi*jpj = ', jpimax*jpjmax, ')'
223            WRITE(numout,9000) '   The best domain decompostion ', inbi, ' x ', inbj, ' with ', icnt2, ' land subdomains'
224            WRITE(numout,9002) '      - uses a total of ',  inbi*inbj-icnt2,' mpi process'
225            WRITE(numout,9000) '      - has mpi subdomains with a maximum size of (jpi = ',  iimax, ', jpj = ',  ijmax,   &
226               &                                                             ', jpi*jpj = ',  iimax* ijmax, ')'
227         ENDIF
228         IF( iimax*ijmax < jpimax*jpjmax ) THEN   ! chosen subdomain size is larger that the best subdomain size
229            llbest = .FALSE.
230            IF ( inbi*inbj-icnt2 < mppsize ) THEN
231               WRITE(ctmp1,*) '   ==> You could therefore have smaller mpi subdomains with less mpi processes'
232            ELSE
233               WRITE(ctmp1,*) '   ==> You could therefore have smaller mpi subdomains with the same number of mpi processes'
234            ENDIF
235            CALL ctl_warn( ' ', ctmp1, ' ', '    ---   YOU ARE WASTING CPU...   ---', ' ' )
236         ELSE IF ( iimax*ijmax == jpimax*jpjmax .AND. (inbi*inbj-icnt2) <  mppsize) THEN
237            llbest = .FALSE.
238            WRITE(ctmp1,*) '   ==> You could therefore have the same mpi subdomains size with less mpi processes'
239            CALL ctl_warn( ' ', ctmp1, ' ', '    ---   YOU ARE WASTING CPU...   ---', ' ' )
240         ELSE
241            llbest = .TRUE.
242         ENDIF
243      ENDIF
244
245      ! look for land mpi subdomains...
246      ALLOCATE( llisOce(jpni,jpnj) )
247      CALL mpp_is_ocean( llisOce )
248      inijmin = COUNT( llisOce )   ! number of oce subdomains
249
250      IF( mppsize < inijmin ) THEN   ! too many oce subdomains: can happen only if jpni and jpnj are prescribed...
251         WRITE(ctmp1,9001) '   With this specified domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj
252         WRITE(ctmp2,9002) '   we can eliminate only ', jpni*jpnj - inijmin, ' land mpi subdomains therefore '
253         WRITE(ctmp3,9001) '   the number of ocean mpi subdomains (', inijmin,') exceed the number of MPI processes:', mppsize
254         WRITE(ctmp4,*) '   ==>>> There is the list of best domain decompositions you should use: '
255         CALL ctl_stop( ctmp1, ctmp2, ctmp3, ' ', ctmp4, ' ' )
256         CALL bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core
257      ENDIF
258
259      IF( mppsize > jpni*jpnj ) THEN   ! not enough mpi subdomains for the total number of mpi processes
260         IF(lwp) THEN
261            WRITE(numout,9003) '   The number of mpi processes: ', mppsize
262            WRITE(numout,9003) '   exceeds the maximum number of subdomains (ocean+land) = ', jpni*jpnj
263            WRITE(numout,9001) '   defined by the following domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj
264            WRITE(numout,   *) '   You should: '
265           IF( llauto ) THEN
266               WRITE(numout,*) '     - either prescribe your domain decomposition with the namelist variables'
267               WRITE(numout,*) '       jpni and jpnj to match the number of mpi process you want to use, '
268               WRITE(numout,*) '       even IF it not the best choice...'
269               WRITE(numout,*) '     - or keep the automatic and optimal domain decomposition by picking up one'
270               WRITE(numout,*) '       of the number of mpi process proposed in the list bellow'
271            ELSE
272               WRITE(numout,*) '     - either properly prescribe your domain decomposition with jpni and jpnj'
273               WRITE(numout,*) '       in order to be consistent with the number of mpi process you want to use'
274               WRITE(numout,*) '       even IF it not the best choice...'
275               WRITE(numout,*) '     - or use the automatic and optimal domain decomposition and pick up one of'
276               WRITE(numout,*) '       the domain decomposition proposed in the list bellow'
277            ENDIF
278            WRITE(numout,*)
279         ENDIF
280         CALL bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core
281      ENDIF
282
283      jpnij = mppsize   ! force jpnij definition <-- remove as much land subdomains as needed to reach this condition
284      IF( mppsize > inijmin ) THEN
285         WRITE(ctmp1,9003) '   The number of mpi processes: ', mppsize
286         WRITE(ctmp2,9003) '   exceeds the maximum number of ocean subdomains = ', inijmin
287         WRITE(ctmp3,9002) '   we suppressed ', jpni*jpnj - mppsize, ' land subdomains '
288         WRITE(ctmp4,9002) '   BUT we had to keep ', mppsize - inijmin, ' land subdomains that are useless...'
289         CALL ctl_warn( ctmp1, ctmp2, ctmp3, ctmp4, ' ', '    --- YOU ARE WASTING CPU... ---', ' ' )
290      ELSE   ! mppsize = inijmin
291         IF(lwp) THEN
292            IF(llbest) WRITE(numout,*) '   ==> you use the best mpi decomposition'
293            WRITE(numout,*)
294            WRITE(numout,9003) '   Number of mpi processes: ', mppsize
295            WRITE(numout,9003) '   Number of ocean subdomains = ', inijmin
296            WRITE(numout,9003) '   Number of suppressed land subdomains = ', jpni*jpnj - inijmin
297            WRITE(numout,*)
298         ENDIF
299      ENDIF
3009000  FORMAT (a, i4, a, i4, a, i7, a)
3019001  FORMAT (a, i4, a, i4)
3029002  FORMAT (a, i4, a)
3039003  FORMAT (a, i5)
304
305      ALLOCATE( nfimpp(jpni), nfproc(jpni), nfjpi(jpni),   &
306         &      iin(jpnij), ijn(jpnij),   &
307         &      iimppt(jpni,jpnj), ijmppt(jpni,jpnj), ijpi(jpni,jpnj), ijpj(jpni,jpnj), ipproc(jpni,jpnj),   &
308         &      inei(8,jpni,jpnj), llnei(8,jpni,jpnj),   &
309         &      impi(8,jpnij),   &
310         &      STAT=ierr )
311      CALL mpp_sum( 'mppini', ierr )
312      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' )
313
314#if defined key_agrif
315         CALL agrif_nemo_init()
316#endif
317      !
318      !  2. Index arrays for subdomains
319      ! -----------------------------------
320      !
321      CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ijpi, ijpj )
322      CALL mpp_getnum( llisOce, ipproc, iin, ijn )
323      !
324      ii = iin(narea)
325      ij = ijn(narea)
326      jpi   = ijpi(ii,ij)
327      jpj   = ijpj(ii,ij)
328      jpk   = MAX( 2, jpkglo )
329      jpij  = jpi*jpj
330      nimpp = iimppt(ii,ij)
331      njmpp = ijmppt(ii,ij)
332      !
333      CALL init_doloop    ! set start/end indices of do-loop, depending on the halo width value (nn_hls)
334      CALL init_locglo    ! define now functions needed to convert indices from/to global to/from local domains
335      !
336      IF(lwp) THEN
337         WRITE(numout,*)
338         WRITE(numout,*) 'MPI Message Passing MPI - domain lay out over processors'
339         WRITE(numout,*)
340         WRITE(numout,*) '   defines mpp subdomains'
341         WRITE(numout,*) '      jpni = ', jpni
342         WRITE(numout,*) '      jpnj = ', jpnj
343         WRITE(numout,*) '     jpnij = ', jpnij
344         WRITE(numout,*) '     nimpp = ', nimpp
345         WRITE(numout,*) '     njmpp = ', njmpp
346         WRITE(numout,*)
347         WRITE(numout,*) '      sum ijpi(i,1) = ', sum(ijpi(:,1)), ' jpiglo = ', jpiglo
348         WRITE(numout,*) '      sum ijpj(1,j) = ', SUM(ijpj(1,:)), ' jpjglo = ', jpjglo
349         
350         ! Subdomain grid print
351         ifreq = 4
352         il1 = 1
353         DO jn = 1, (jpni-1)/ifreq+1
354            il2 = MIN(jpni,il1+ifreq-1)
355            WRITE(numout,*)
356            WRITE(numout,9400) ('***',ji=il1,il2-1)
357            DO jj = jpnj, 1, -1
358               WRITE(numout,9403) ('   ',ji=il1,il2-1)
359               WRITE(numout,9402) jj, (ijpi(ji,jj),ijpj(ji,jj),ji=il1,il2)
360               WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2)
361               WRITE(numout,9403) ('   ',ji=il1,il2-1)
362               WRITE(numout,9400) ('***',ji=il1,il2-1)
363            END DO
364            WRITE(numout,9401) (ji,ji=il1,il2)
365            il1 = il1+ifreq
366         END DO
367 9400    FORMAT('           ***'   ,20('*************',a3)    )
368 9403    FORMAT('           *     ',20('         *   ',a3)    )
369 9401    FORMAT('              '   ,20('   ',i3,'          ') )
370 9402    FORMAT('       ',i3,' *  ',20(i3,'  x',i3,'   *   ') )
371 9404    FORMAT('           *  '   ,20('     ' ,i4,'   *   ') )
372      ENDIF
373      !
374      ! Store informations for the north pole folding communications
375      nfproc(:) = ipproc(:,jpnj)
376      nfimpp(:) = iimppt(:,jpnj)
377      nfjpi (:) =   ijpi(:,jpnj)
378      !
379      ! 3. Define Western, Eastern, Southern and Northern neighbors + corners in the subdomain grid reference
380      ! ------------------------------------------------------------------------------------------------------
381      !
382      ! note that North fold is has specific treatment for its MPI communications.
383      ! This must not be treated as a "usual" communication with a northern neighbor.
384      !    -> North fold processes have no Northern neighbor in the definition done bellow
385      !
386      llmpi_Iperio = jpni > 1 .AND. l_Iperio                         ! do i-periodicity with an MPI communication?
387      llmpi_Jperio = jpnj > 1 .AND. l_Jperio                         ! do j-periodicity with an MPI communication?
388      !
389      l_SelfPerio(1:2) = l_Iperio .AND. jpni == 1                    !  west,  east periodicity by itself
390      l_SelfPerio(3:4) = l_Jperio .AND. jpnj == 1                    ! south, north periodicity by itself
391      l_SelfPerio(5:8) = l_SelfPerio(jpwe) .AND. l_SelfPerio(jpso)   ! corners bi-periodicity by itself
392      !
393      ! define neighbors mapping (1/2): default definition: ignore if neighbours are land-only subdomains or not
394      DO jj = 1, jpnj
395         DO ji = 1, jpni
396            !
397            IF ( llisOce(ji,jj) ) THEN                     ! this subdomain has some ocean: it has neighbours
398               !
399               inum0 = ji - 1 + ( jj - 1 ) * jpni             ! index in the subdomains grid. start at 0
400               !
401               ! Is there a neighbor?
402               llnei(jpwe,ji,jj) = ji >   1  .OR. llmpi_Iperio           ! West  nei exists if not the first column or llmpi_Iperio
403               llnei(jpea,ji,jj) = ji < jpni .OR. llmpi_Iperio           ! East  nei exists if not the last  column or llmpi_Iperio
404               llnei(jpso,ji,jj) = jj >   1  .OR. llmpi_Jperio           ! South nei exists if not the first line   or llmpi_Jperio
405               llnei(jpno,ji,jj) = jj < jpnj .OR. llmpi_Jperio           ! North nei exists if not the last  line   or llmpi_Jperio
406               llnei(jpsw,ji,jj) = llnei(jpwe,ji,jj) .AND. llnei(jpso,ji,jj)   ! So-We nei exists if both South and West nei exist
407               llnei(jpse,ji,jj) = llnei(jpea,ji,jj) .AND. llnei(jpso,ji,jj)   ! So-Ea nei exists if both South and East nei exist
408               llnei(jpnw,ji,jj) = llnei(jpwe,ji,jj) .AND. llnei(jpno,ji,jj)   ! No-We nei exists if both North and West nei exist
409               llnei(jpne,ji,jj) = llnei(jpea,ji,jj) .AND. llnei(jpno,ji,jj)   ! No-Ea nei exists if both North and East nei exist
410               !
411               ! Which index (starting at 0) have neighbors in the subdomains grid?
412               IF( llnei(jpwe,ji,jj) )   inei(jpwe,ji,jj) =            inum0 -    1 + jpni        * COUNT( (/ ji ==    1 /) )
413               IF( llnei(jpea,ji,jj) )   inei(jpea,ji,jj) =            inum0 +    1 - jpni        * COUNT( (/ ji == jpni /) )
414               IF( llnei(jpso,ji,jj) )   inei(jpso,ji,jj) =            inum0 - jpni + jpni * jpnj * COUNT( (/ jj ==    1 /) )
415               IF( llnei(jpno,ji,jj) )   inei(jpno,ji,jj) =            inum0 + jpni - jpni * jpnj * COUNT( (/ jj == jpnj /) )
416               IF( llnei(jpsw,ji,jj) )   inei(jpsw,ji,jj) = inei(jpso,ji,jj) -    1 + jpni        * COUNT( (/ ji ==    1 /) )
417               IF( llnei(jpse,ji,jj) )   inei(jpse,ji,jj) = inei(jpso,ji,jj) +    1 - jpni        * COUNT( (/ ji == jpni /) )
418               IF( llnei(jpnw,ji,jj) )   inei(jpnw,ji,jj) = inei(jpno,ji,jj) -    1 + jpni        * COUNT( (/ ji ==    1 /) )
419               IF( llnei(jpne,ji,jj) )   inei(jpne,ji,jj) = inei(jpno,ji,jj) +    1 - jpni        * COUNT( (/ ji == jpni /) )
420               !
421            ELSE                                           ! land-only domain has no neighbour
422               llnei(:,ji,jj) = .FALSE.
423            ENDIF
424            !
425         END DO
426      END DO
427      !
428      ! define neighbors mapping (2/2): check if neighbours are not land-only subdomains
429      DO jj = 1, jpnj
430         DO ji = 1, jpni
431            DO jn = 1, 8
432               IF( llnei(jn,ji,jj) ) THEN   ! if a neighbour is existing -> this should not be a land-only domain
433                  ii = 1 + MOD( inei(jn,ji,jj) , jpni )
434                  ij = 1 +      inei(jn,ji,jj) / jpni
435                  llnei(jn,ji,jj) = llisOce( ii, ij )
436               ENDIF
437            END DO
438         END DO
439      END DO
440      !
441      ! update index of the neighbours in the subdomains grid
442      WHERE( .NOT. llnei )   inei = -1
443      !
444      ! Save processor layout in ascii file
445      IF (llwrtlay) THEN
446         CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
447         WRITE(inum,'(a)') '  jpnij jpimax jpjmax    jpk jpiglo jpjglo ( local:   narea    jpi    jpj )'
448         WRITE(inum,'(6i7,a,3i7,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,' ( local: ',narea,jpi,jpj,' )'
449         WRITE(inum,*) 
450         WRITE(inum,       *) '------------------------------------'
451         WRITE(inum,'(a,i2)') ' Mapping of the default neighnourgs '
452         WRITE(inum,       *) '------------------------------------'
453         WRITE(inum,*) 
454         WRITE(inum,'(a)') '  rank    ii    ij   jpi   jpj nimpp njmpp mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine'
455         DO jp = 1, jpnij
456            ii = iin(jp)
457            ij = ijn(jp)
458            WRITE(inum,'(15i6)')  jp-1, ii, ij, ijpi(ii,ij), ijpj(ii,ij), iimppt(ii,ij), ijmppt(ii,ij), inei(:,ii,ij)
459         END DO
460      ENDIF
461
462      !
463      ! 4. Define Western, Eastern, Southern and Northern neighbors + corners for each mpi process
464      ! ------------------------------------------------------------------------------------------
465      !
466      ! rewrite information from "subdomain grid" to mpi process list
467      ! Warning, for example:
468      !    position of the northern neighbor in the "subdomain grid"
469      !    position of the northern neighbor in the "mpi process list"
470     
471      ! default definition: no neighbors
472      impi(:,:) = -1   ! (starting at 0, -1 if no neighbourg)
473     
474      DO jp = 1, jpnij
475         ii = iin(jp)
476         ij = ijn(jp)
477         DO jn = 1, 8
478            IF( llnei(jn,ii,ij) ) THEN   ! must be tested as some land-domain can be kept to fit mppsize
479               ii2 = 1 + MOD( inei(jn,ii,ij) , jpni )
480               ij2 = 1 +      inei(jn,ii,ij) / jpni
481               impi(jn,jp) = ipproc( ii2, ij2 )
482            ENDIF
483         END DO
484      END DO
485
486      !
487      ! 4. keep information for the local process
488      ! -----------------------------------------
489      !
490      ! set default neighbours
491      mpinei(:) = impi(:,narea)
492      DO jh = 1, n_hlsmax
493         mpiSnei(jh,:) = impi(:,narea)   ! default definition
494         mpiRnei(jh,:) = impi(:,narea)
495      END DO
496      !
497      IF(lwp) THEN
498         WRITE(numout,*)
499         WRITE(numout,*) '   resulting internal parameters : '
500         WRITE(numout,*) '      narea = ', narea
501         WRITE(numout,*) '      mpi nei  west = ', mpinei(jpwe)  , '   mpi nei  east = ', mpinei(jpea)
502         WRITE(numout,*) '      mpi nei south = ', mpinei(jpso)  , '   mpi nei north = ', mpinei(jpno)
503         WRITE(numout,*) '      mpi nei so-we = ', mpinei(jpsw)  , '   mpi nei so-ea = ', mpinei(jpse)
504         WRITE(numout,*) '      mpi nei no-we = ', mpinei(jpnw)  , '   mpi nei no-ea = ', mpinei(jpne)
505      ENDIF
506      !
507      CALL mpp_ini_nc(nn_hls)    ! Initialize communicator for neighbourhood collective communications
508      DO jh = 1, n_hlsmax
509         mpi_nc_com4(jh) = mpi_nc_com4(nn_hls)   ! default definition
510         mpi_nc_com8(jh) = mpi_nc_com8(nn_hls)
511      END DO
512      !                          ! Exclude exchanges which contain only land points
513      !
514      IF( jpnij > 1 ) CALL init_excl_landpt
515      !
516      !                          ! Prepare mpp north fold
517      !
518      llmpiNFold =          jpni  > 1 .AND. l_NFold           ! is the North fold done with an MPI communication?
519      l_IdoNFold = ijn(narea) == jpnj .AND. l_NFold           ! is this process doing North fold?
520      !
521      IF( llmpiNFold )   CALL init_nfdcom( llwrtlay, inum )   ! init northfold communication, must be done after init_excl_landpt
522      !
523      !                          !  Save processor layout changes in ascii file
524      !
525      DO jh = 1, n_hlsmax    ! different halo size
526         DO ji = 1, 8
527            ichanged(16*(jh-1)  +ji) = COUNT( mpinei(ji:ji) /= mpiSnei(jh,ji:ji) )
528            ichanged(16*(jh-1)+8+ji) = COUNT( mpinei(ji:ji) /= mpiRnei(jh,ji:ji) )
529         END DO
530      END DO
531      CALL mpp_sum( "mpp_init", ichanged )   ! must be called by all processes
532      IF (llwrtlay) THEN
533         WRITE(inum,*) 
534         WRITE(inum,       *) '----------------------------------------------------------------------'
535         WRITE(inum,'(a,i2)') ' Mapping of the neighnourgs once excluding comm with only land points '
536         WRITE(inum,       *) '----------------------------------------------------------------------'
537         DO jh = 1, n_hlsmax    ! different halo size
538            WRITE(inum,*) 
539            WRITE(inum,'(a,i2)') 'halo size: ', jh
540            WRITE(inum,       *) '---------'
541            WRITE(inum,'(a)') '  rank    ii    ij mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine'
542            WRITE(inum,   '(11i6,a)')  narea-1, iin(narea), ijn(narea),   mpinei(:), ' <- Org'
543            WRITE(inum,'(18x,8i6,a,i1,a)')   mpiSnei(jh,:), ' <- Send ', COUNT( mpinei(:) /= mpiSnei(jh,:) ), ' modif'
544            WRITE(inum,'(18x,8i6,a,i1,a)')   mpiRnei(jh,:), ' <- Recv ', COUNT( mpinei(:) /= mpiRnei(jh,:) ), ' modif'
545            WRITE(inum,*) ' total changes among all mpi tasks:'
546            WRITE(inum,*) '       mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine'
547            WRITE(inum,'(a,8i6)') ' Send: ', ichanged(jh*16-15:jh*16-8)
548            WRITE(inum,'(a,8i6)') ' Recv: ', ichanged(jh*16 -7:jh*16  )
549         END DO
550      ENDIF
551      !
552      CALL init_ioipsl           ! Prepare NetCDF output file (if necessary)
553      !
554      IF (llwrtlay) CLOSE(inum)
555      !
556      DEALLOCATE(iin, ijn, iimppt, ijmppt, ijpi, ijpj, ipproc, inei, llnei, impi, llisOce)
557      !
558    END SUBROUTINE mpp_init
559
560#endif
561
562    SUBROUTINE mpp_basesplit( kiglo, kjglo, khls, knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj)
563      !!----------------------------------------------------------------------
564      !!                  ***  ROUTINE mpp_basesplit  ***
565      !!
566      !! ** Purpose :   Lay out the global domain over processors.
567      !!
568      !! ** Method  :   Global domain is distributed in smaller local domains.
569      !!
570      !! ** Action : - set for all knbi*knbj domains:
571      !!                    kimppt     : longitudinal index
572      !!                    kjmppt     : latitudinal  index
573      !!                    klci       : first dimension
574      !!                    klcj       : second dimension
575      !!----------------------------------------------------------------------
576      INTEGER,                                 INTENT(in   ) ::   kiglo, kjglo
577      INTEGER,                                 INTENT(in   ) ::   khls
578      INTEGER,                                 INTENT(in   ) ::   knbi, knbj
579      INTEGER,                                 INTENT(  out) ::   kimax, kjmax
580      INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT(  out) ::   kimppt, kjmppt
581      INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT(  out) ::   klci, klcj
582      !
583      INTEGER ::   ji, jj
584      INTEGER ::   i2hls
585      INTEGER ::   iresti, irestj, irm, ijpjmin
586      !!----------------------------------------------------------------------
587      i2hls = 2*khls
588      !
589#if defined key_nemocice_decomp
590      kimax = ( nx_global+2-i2hls + (knbi-1) ) / knbi + i2hls    ! first  dim.
591      kjmax = ( ny_global+2-i2hls + (knbj-1) ) / knbj + i2hls    ! second dim.
592#else
593      kimax = ( kiglo - i2hls + (knbi-1) ) / knbi + i2hls    ! first  dim.
594      kjmax = ( kjglo - i2hls + (knbj-1) ) / knbj + i2hls    ! second dim.
595#endif
596      IF( .NOT. PRESENT(kimppt) ) RETURN
597      !
598      !  1. Dimension arrays for subdomains
599      ! -----------------------------------
600      !  Computation of local domain sizes klci() klcj()
601      !  These dimensions depend on global sizes knbi,knbj and kiglo,kjglo
602      !  The subdomains are squares lesser than or equal to the global
603      !  dimensions divided by the number of processors minus the overlap array.
604      !
605      iresti = 1 + MOD( kiglo - i2hls - 1 , knbi )
606      irestj = 1 + MOD( kjglo - i2hls - 1 , knbj )
607      !
608      !  Need to use kimax and kjmax here since jpi and jpj not yet defined
609#if defined key_nemocice_decomp
610      ! Change padding to be consistent with CICE
611      klci(1:knbi-1,:       ) = kimax
612      klci(  knbi  ,:       ) = kiglo - (knbi - 1) * (kimax - i2hls)
613      klcj(:       ,1:knbj-1) = kjmax
614      klcj(:       ,  knbj  ) = kjglo - (knbj - 1) * (kjmax - i2hls)
615#else
616      klci(1:iresti      ,:) = kimax
617      klci(iresti+1:knbi ,:) = kimax-1
618      IF( MINVAL(klci) < 3*khls ) THEN
619         WRITE(ctmp1,*) '   mpp_basesplit: minimum value of jpi must be >= ', 3*khls
620         WRITE(ctmp2,*) '   We have ', MINVAL(klci)
621         CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
622      ENDIF
623      IF( l_NFold ) THEN
624         ! minimize the size of the last row to compensate for the north pole folding coast
625         IF( c_NFtype == 'T' )   ijpjmin = 2+3*khls   ! V and F folding must be outside of southern halos
626         IF( c_NFtype == 'F' )   ijpjmin = 1+3*khls   ! V and F folding must be outside of southern halos
627         irm = knbj - irestj                          ! total number of lines to be removed
628         klcj(:,knbj) = MAX( ijpjmin, kjmax-irm )     ! we must have jpj >= ijpjmin in the last row
629         irm = irm - ( kjmax - klcj(1,knbj) )         ! remaining number of lines to remove
630         irestj = knbj - 1 - irm
631         klcj(:, irestj+1:knbj-1) = kjmax-1
632      ELSE
633         klcj(:, irestj+1:knbj  ) = kjmax-1
634      ENDIF
635      klcj(:,1:irestj) = kjmax
636      IF( MINVAL(klcj) < 3*khls ) THEN
637         WRITE(ctmp1,*) '   mpp_basesplit: minimum value of jpj must be >= ', 3*khls
638         WRITE(ctmp2,*) '   We have ', MINVAL(klcj)
639         CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
640      ENDIF
641#endif
642
643      !  2. Index arrays for subdomains
644      ! -------------------------------
645      kimppt(:,:) = 1
646      kjmppt(:,:) = 1
647      !
648      IF( knbi > 1 ) THEN
649         DO jj = 1, knbj
650            DO ji = 2, knbi
651               kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - i2hls
652            END DO
653         END DO
654      ENDIF
655      !
656      IF( knbj > 1 )THEN
657         DO jj = 2, knbj
658            DO ji = 1, knbi
659               kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - i2hls
660            END DO
661         END DO
662      ENDIF
663
664   END SUBROUTINE mpp_basesplit
665
666
667   SUBROUTINE bestpartition( knbij, knbi, knbj, knbcnt, ldlist )
668      !!----------------------------------------------------------------------
669      !!                 ***  ROUTINE bestpartition  ***
670      !!
671      !! ** Purpose :
672      !!
673      !! ** Method  :
674      !!----------------------------------------------------------------------
675      INTEGER,           INTENT(in   ) ::   knbij         ! total number of subdomains (knbi*knbj)
676      INTEGER, OPTIONAL, INTENT(  out) ::   knbi, knbj    ! number if subdomains along i and j (knbi and knbj)
677      INTEGER, OPTIONAL, INTENT(  out) ::   knbcnt        ! number of land subdomains
678      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldlist        ! .true.: print the list the best domain decompositions (with land)
679      !
680      INTEGER :: ji, jj, ii, iitarget
681      INTEGER :: iszitst, iszjtst
682      INTEGER :: isziref, iszjref
683      INTEGER :: iszimin, iszjmin
684      INTEGER :: inbij, iszij
685      INTEGER :: inbimax, inbjmax, inbijmax, inbijold
686      INTEGER :: isz0, isz1
687      INTEGER, DIMENSION(  :), ALLOCATABLE :: indexok
688      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi0, inbj0, inbij0   ! number of subdomains along i,j
689      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi0, iszj0, iszij0   ! max size of the subdomains along i,j
690      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi1, inbj1, inbij1   ! number of subdomains along i,j
691      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi1, iszj1, iszij1   ! max size of the subdomains along i,j
692      LOGICAL :: llist
693      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d                 ! max size of the subdomains along i,j
694      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisOce              !  -     -
695      REAL(wp)::   zpropland
696      !!----------------------------------------------------------------------
697      !
698      llist = .FALSE.
699      IF( PRESENT(ldlist) ) llist = ldlist
700
701      CALL mpp_init_landprop( zpropland )                      ! get the proportion of land point over the gloal domain
702      inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) )    ! define the largest possible value for jpni*jpnj
703      !
704      IF( llist ) THEN   ;   inbijmax = inbij*2
705      ELSE               ;   inbijmax = inbij
706      ENDIF
707      !
708      ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax))
709      !
710      inbimax = 0
711      inbjmax = 0
712      isziref = jpiglo*jpjglo+1   ! define a value that is larger than the largest possible
713      iszjref = jpiglo*jpjglo+1
714      !
715      ! WARNING, see also init_excl_landpt: minimum subdomain size defined here according to nn_hls (and not n_hlsmax)
716      ! --> If, one day, we want to use local halos largers than nn_hls, we must replace nn_hls by n_hlsmax
717      !
718      iszimin = 3*nn_hls          ! minimum size of the MPI subdomain so halos are always adressing neighbor inner domain
719      iszjmin = 3*nn_hls
720      IF( c_NFtype == 'T' )   iszjmin = MAX(iszjmin, 2+3*nn_hls)   ! V and F folding must be outside of southern halos
721      IF( c_NFtype == 'F' )   iszjmin = MAX(iszjmin, 1+3*nn_hls)   ! V and F folding must be outside of southern halos
722      !
723      ! get the list of knbi that gives a smaller jpimax than knbi-1
724      ! get the list of knbj that gives a smaller jpjmax than knbj-1
725      DO ji = 1, inbijmax
726#if defined key_nemocice_decomp
727         iszitst = ( nx_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
728#else
729         iszitst = ( Ni0glo + (ji-1) ) / ji + 2*nn_hls   ! max subdomain i-size
730#endif
731         IF( iszitst < isziref .AND. iszitst >= iszimin ) THEN
732            isziref = iszitst
733            inbimax = inbimax + 1
734            inbi0(inbimax) = ji
735            iszi0(inbimax) = isziref
736         ENDIF
737#if defined key_nemocice_decomp
738         iszjtst = ( ny_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
739#else
740         iszjtst = ( Nj0glo + (ji-1) ) / ji + 2*nn_hls   ! max subdomain j-size
741#endif
742         IF( iszjtst < iszjref .AND. iszjtst >= iszjmin ) THEN
743            iszjref = iszjtst
744            inbjmax = inbjmax + 1
745            inbj0(inbjmax) = ji
746            iszj0(inbjmax) = iszjref
747         ENDIF
748      END DO
749      IF( inbimax == 0 ) THEN
750         WRITE(ctmp1,'(a,i2,a,i2)') '   mpp_ini bestpartition: Ni0glo (', Ni0glo, ') is too small to be used with nn_hls = ', nn_hls
751         CALL ctl_stop( 'STOP', ctmp1 )
752      ENDIF
753      IF( inbjmax == 0 ) THEN
754         WRITE(ctmp1,'(a,i2,a,i2)') '   mpp_ini bestpartition: Nj0glo (', Nj0glo, ') is too small to be used with nn_hls = ', nn_hls
755         CALL ctl_stop( 'STOP', ctmp1 )
756      ENDIF
757
758      ! combine these 2 lists to get all possible knbi*knbj <  inbijmax
759      ALLOCATE( llmsk2d(inbimax,inbjmax) )
760      DO jj = 1, inbjmax
761         DO ji = 1, inbimax
762            IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN   ;   llmsk2d(ji,jj) = .TRUE.
763            ELSE                                            ;   llmsk2d(ji,jj) = .FALSE.
764            ENDIF
765         END DO
766      END DO
767      isz1 = COUNT(llmsk2d)
768      ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) )
769      ii = 0
770      DO jj = 1, inbjmax
771         DO ji = 1, inbimax
772            IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN
773               ii = ii + 1
774               inbi1(ii) = inbi0(ji)
775               inbj1(ii) = inbj0(jj)
776               iszi1(ii) = iszi0(ji)
777               iszj1(ii) = iszj0(jj)
778            ENDIF
779         END DO
780      END DO
781      DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
782      DEALLOCATE( llmsk2d )
783
784      ALLOCATE( inbij1(isz1), iszij1(isz1) )
785      inbij1(:) = inbi1(:) * inbj1(:)
786      iszij1(:) = iszi1(:) * iszj1(:)
787
788      ! if there is no land and no print
789      IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN
790         ! get the smaller partition which gives the smallest subdomain size
791         ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1)
792         knbi = inbi1(ii)
793         knbj = inbj1(ii)
794         IF(PRESENT(knbcnt))   knbcnt = 0
795         DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 )
796         RETURN
797      ENDIF
798
799      ! extract only the partitions which reduce the subdomain size in comparison with smaller partitions
800      ALLOCATE( indexok(isz1) )                                 ! to store indices of the best partitions
801      isz0 = 0                                                  ! number of best partitions
802      inbij = 1                                                 ! start with the min value of inbij1 => 1
803      iszij = jpiglo*jpjglo+1                                   ! default: larger than global domain
804      DO WHILE( inbij <= inbijmax )                             ! if we did not reach the max of inbij1
805         ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1)   ! warning: send back the first occurence if multiple results
806         IF ( iszij1(ii) < iszij ) THEN
807            ii = MINLOC( iszi1+iszj1, mask = iszij1 == iszij1(ii) .AND. inbij1 == inbij, dim = 1)  ! select the smaller perimeter if multiple min
808            isz0 = isz0 + 1
809            indexok(isz0) = ii
810            iszij = iszij1(ii)
811         ENDIF
812         inbij = MINVAL(inbij1, mask = inbij1 > inbij)   ! warning: return largest integer value if mask = .false. everywhere
813      END DO
814      DEALLOCATE( inbij1, iszij1 )
815
816      ! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size)
817      ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) )
818      DO ji = 1, isz0
819         ii = indexok(ji)
820         inbi0(ji) = inbi1(ii)
821         inbj0(ji) = inbj1(ii)
822         iszi0(ji) = iszi1(ii)
823         iszj0(ji) = iszj1(ii)
824      END DO
825      DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 )
826
827      IF( llist ) THEN
828         IF(lwp) THEN
829            WRITE(numout,*)
830            WRITE(numout,*) '                  For your information:'
831            WRITE(numout,*) '  list of the best partitions including land supression'
832            WRITE(numout,*) '  -----------------------------------------------------'
833            WRITE(numout,*)
834         ENDIF
835         ji = isz0   ! initialization with the largest value
836         ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) )
837         CALL mpp_is_ocean( llisOce )   ! Warning: must be call by all cores (call mpp_sum)
838         inbijold = COUNT(llisOce)
839         DEALLOCATE( llisOce )
840         DO ji =isz0-1,1,-1
841            ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) )
842            CALL mpp_is_ocean( llisOce )   ! Warning: must be call by all cores (call mpp_sum)
843            inbij = COUNT(llisOce)
844            DEALLOCATE( llisOce )
845            IF(lwp .AND. inbij < inbijold) THEN
846               WRITE(numout,'(a, i6, a, i6, a, f4.1, a, i9, a, i6, a, i6, a)')                                 &
847                  &   'nb_cores oce: ', inbij, ', land domains excluded: ', inbi0(ji)*inbj0(ji) - inbij,       &
848                  &   ' (', REAL(inbi0(ji)*inbj0(ji) - inbij,wp) / REAL(inbi0(ji)*inbj0(ji),wp) *100.,         &
849                  &   '%), largest oce domain: ', iszi0(ji)*iszj0(ji), ' ( ', iszi0(ji),' x ', iszj0(ji), ' )'
850               inbijold = inbij
851            ENDIF
852         END DO
853         DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
854         IF(lwp) THEN
855            WRITE(numout,*)
856            WRITE(numout,*)  '  -----------------------------------------------------------'
857         ENDIF
858         CALL mppsync
859         CALL mppstop( ld_abort = .TRUE. )
860      ENDIF
861
862      DEALLOCATE( iszi0, iszj0 )
863      inbij = inbijmax + 1        ! default: larger than possible
864      ii = isz0+1                 ! start from the end of the list (smaller subdomains)
865      DO WHILE( inbij > knbij )   ! while the number of ocean subdomains exceed the number of procs
866         ii = ii -1
867         ALLOCATE( llisOce(inbi0(ii), inbj0(ii)) )
868         CALL mpp_is_ocean( llisOce )            ! must be done by all core
869         inbij = COUNT(llisOce)
870         DEALLOCATE( llisOce )
871      END DO
872      knbi = inbi0(ii)
873      knbj = inbj0(ii)
874      IF(PRESENT(knbcnt))   knbcnt = knbi * knbj - inbij
875      DEALLOCATE( inbi0, inbj0 )
876      !
877   END SUBROUTINE bestpartition
878
879
880   SUBROUTINE mpp_init_landprop( propland )
881      !!----------------------------------------------------------------------
882      !!                  ***  ROUTINE mpp_init_landprop  ***
883      !!
884      !! ** Purpose : the the proportion of land points in the surface land-sea mask
885      !!
886      !! ** Method  : read iproc strips (of length Ni0glo) of the land-sea mask
887      !!----------------------------------------------------------------------
888      REAL(wp), INTENT(  out) :: propland    ! proportion of land points in the global domain (between 0 and 1)
889      !
890      INTEGER, DIMENSION(jpni*jpnj) ::   kusedom_1d
891      INTEGER :: inboce, iarea
892      INTEGER :: iproc, idiv, ijsz
893      INTEGER :: ijstr
894      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce
895      !!----------------------------------------------------------------------
896      ! do nothing if there is no land-sea mask
897      IF( numbot == -1 .and. numbdy == -1 ) THEN
898         propland = 0.
899         RETURN
900      ENDIF
901
902      ! number of processes reading the bathymetry file
903      iproc = MINVAL( (/mppsize, Nj0glo/2, 100/) )  ! read a least 2 lines, no more that 100 processes reading at the same time
904
905      ! we want to read iproc strips of the land-sea mask. -> pick up iproc processes every idiv processes starting at 1
906      IF( iproc == 1 ) THEN   ;   idiv = mppsize
907      ELSE                    ;   idiv = ( mppsize - 1 ) / ( iproc - 1 )
908      ENDIF
909
910      iarea = (narea-1)/idiv   ! involed process number (starting counting at 0)
911      IF( MOD( narea-1, idiv ) == 0 .AND. iarea < iproc ) THEN   ! beware idiv can be = to 1
912         !
913         ijsz = Nj0glo / iproc                                               ! width of the stripe to read
914         IF( iarea < MOD(Nj0glo,iproc) ) ijsz = ijsz + 1
915         ijstr = iarea*(Nj0glo/iproc) + MIN(iarea, MOD(Nj0glo,iproc)) + 1    ! starting j position of the reading
916         !
917         ALLOCATE( lloce(Ni0glo, ijsz) )                                     ! allocate the strip
918         CALL read_mask( 1, ijstr, Ni0glo, ijsz, lloce )
919         inboce = COUNT(lloce)                                               ! number of ocean point in the stripe
920         DEALLOCATE(lloce)
921         !
922      ELSE
923         inboce = 0
924      ENDIF
925      CALL mpp_sum( 'mppini', inboce )   ! total number of ocean points over the global domain
926      !
927      propland = REAL( Ni0glo*Nj0glo - inboce, wp ) / REAL( Ni0glo*Nj0glo, wp )
928      !
929   END SUBROUTINE mpp_init_landprop
930
931
932   SUBROUTINE mpp_is_ocean( ldIsOce )
933      !!----------------------------------------------------------------------
934      !!                  ***  ROUTINE mpp_is_ocean  ***
935      !!
936      !! ** Purpose : Check for a mpi domain decomposition inbi x inbj which
937      !!              subdomains, including 1 halo (even if nn_hls>1), contain
938      !!              at least 1 ocean point.
939      !!              We must indeed ensure that each subdomain that is a neighbour
940      !!              of a land subdomain, has only land points on its boundary
941      !!              (inside the inner subdomain) with the land subdomain.
942      !!              This is needed to get the proper bondary conditions on
943      !!              a subdomain with a closed boundary.
944      !!
945      !! ** Method  : read inbj strips (of length Ni0glo) of the land-sea mask
946      !!----------------------------------------------------------------------
947      LOGICAL, DIMENSION(:,:), INTENT(  out) ::   ldIsOce        ! .true. if a sub domain constains 1 ocean point
948      !
949      INTEGER :: idiv, iimax, ijmax, iarea
950      INTEGER :: inbi, inbj, inx, iny, inry, isty
951      INTEGER :: ji, jn
952      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   inboce           ! number oce oce pint in each mpi subdomain
953      INTEGER, ALLOCATABLE, DIMENSION(:  ) ::   inboce_1d
954      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ijpi
955      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ijpj
956      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce            ! lloce(i,j) = .true. if the point (i,j) is ocean
957      !!----------------------------------------------------------------------
958      ! do nothing if there is no land-sea mask
959      IF( numbot == -1 .AND. numbdy == -1 ) THEN
960         ldIsOce(:,:) = .TRUE.
961         RETURN
962      ENDIF
963      !
964      inbi = SIZE( ldIsOce, dim = 1 )
965      inbj = SIZE( ldIsOce, dim = 2 )
966      !
967      ! we want to read inbj strips of the land-sea mask. -> pick up inbj processes every idiv processes starting at 1
968      IF           ( inbj == 1 ) THEN   ;   idiv = mppsize
969      ELSE IF ( mppsize < inbj ) THEN   ;   idiv = 1
970      ELSE                              ;   idiv = ( mppsize - 1 ) / ( inbj - 1 )
971      ENDIF
972      !
973      ALLOCATE( inboce(inbi,inbj), inboce_1d(inbi*inbj) )
974      inboce(:,:) = 0          ! default no ocean point found
975      !
976      DO jn = 0, (inbj-1)/mppsize   ! if mppsize < inbj : more strips than mpi processes (because of potential land domains)
977         !
978         iarea = (narea-1)/idiv + jn * mppsize + 1                     ! involed process number (starting counting at 1)
979         IF( MOD( narea-1, idiv ) == 0 .AND. iarea <= inbj ) THEN      ! beware idiv can be = to 1
980            !
981            ALLOCATE( iimppt(inbi,inbj), ijmppt(inbi,inbj), ijpi(inbi,inbj), ijpj(inbi,inbj) )
982            CALL mpp_basesplit( Ni0glo, Nj0glo, 0, inbi, inbj, iimax, ijmax, iimppt, ijmppt, ijpi, ijpj )
983            !
984            inx = Ni0glo + 2   ;   iny = ijpj(1,iarea) + 2             ! strip size + 1 halo on each direction (even if nn_hls>1)
985            ALLOCATE( lloce(inx, iny) )                                ! allocate the strip
986            inry = iny - COUNT( (/ iarea == 1, iarea == inbj /) )      ! number of point to read in y-direction
987            isty = 1 + COUNT( (/ iarea == 1 /) )                       ! read from the first or the second line?
988            CALL read_mask( 1, ijmppt(1,iarea) - 2 + isty, Ni0glo, inry, lloce(2:inx-1, isty:inry+isty-1) )   ! read the strip
989            !
990            IF( iarea == 1    ) THEN                                   ! the first line was not read
991               IF( l_Jperio ) THEN                                     !   north-south periodocity
992                  CALL read_mask( 1, Nj0glo, Ni0glo, 1, lloce(2:inx-1, 1) )   !   read the last line -> first line of lloce
993               ELSE
994                  lloce(2:inx-1,  1) = .FALSE.                         !   closed boundary
995               ENDIF
996            ENDIF
997            IF( iarea == inbj ) THEN                                   ! the last line was not read
998               IF( l_Jperio ) THEN                                     !   north-south periodocity
999                  CALL read_mask( 1, 1, Ni0glo, 1, lloce(2:inx-1,iny) )   !      read the first line -> last line of lloce
1000               ELSEIF( c_NFtype == 'T' ) THEN                          !   north-pole folding T-pivot, T-point
1001                  lloce(2,iny) = lloce(2,iny-2)                        !      here we have 1 halo (even if nn_hls>1)
1002                  DO ji = 3,inx-1
1003                     lloce(ji,iny  ) = lloce(inx-ji+2,iny-2)           !      ok, we have at least 3 lines
1004                  END DO
1005                  DO ji = inx/2+2,inx-1
1006                     lloce(ji,iny-1) = lloce(inx-ji+2,iny-1)
1007                  END DO
1008               ELSEIF( c_NFtype == 'F' ) THEN                          !   north-pole folding F-pivot, T-point, 1 halo
1009                  lloce(inx/2+1,iny-1) = lloce(inx/2,iny-1)            !      here we have 1 halo (even if nn_hls>1)
1010                  lloce(inx  -1,iny-1) = lloce(2    ,iny-1)
1011                  DO ji = 2,inx-1
1012                     lloce(ji,iny) = lloce(inx-ji+1,iny-1)
1013                  END DO
1014               ELSE                                                    !   closed boundary
1015                  lloce(2:inx-1,iny) = .FALSE.
1016               ENDIF
1017            ENDIF
1018            !                                                          ! first and last column were not read
1019            IF( l_Iperio ) THEN
1020               lloce(1,:) = lloce(inx-1,:)   ;   lloce(inx,:) = lloce(2,:)   ! east-west periodocity
1021            ELSE
1022               lloce(1,:) = .FALSE.          ;   lloce(inx,:) = .FALSE.      ! closed boundary
1023            ENDIF
1024            !
1025            DO  ji = 1, inbi
1026               inboce(ji,iarea) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ijpi(ji,1)+1,:) )   ! lloce as 2 points more than Ni0glo
1027            END DO
1028            !
1029            DEALLOCATE(lloce)
1030            DEALLOCATE(iimppt, ijmppt, ijpi, ijpj)
1031            !
1032         ENDIF
1033      END DO
1034
1035      inboce_1d = RESHAPE(inboce, (/ inbi*inbj /))
1036      CALL mpp_sum( 'mppini', inboce_1d )
1037      inboce = RESHAPE(inboce_1d, (/inbi, inbj/))
1038      ldIsOce(:,:) = inboce(:,:) /= 0
1039      DEALLOCATE(inboce, inboce_1d)
1040      !
1041   END SUBROUTINE mpp_is_ocean
1042
1043
1044   SUBROUTINE read_mask( kistr, kjstr, kicnt, kjcnt, ldoce )
1045      !!----------------------------------------------------------------------
1046      !!                  ***  ROUTINE read_mask  ***
1047      !!
1048      !! ** Purpose : Read relevant bathymetric information in order to
1049      !!              provide a land/sea mask used for the elimination
1050      !!              of land domains, in an mpp computation.
1051      !!
1052      !! ** Method  : read stipe of size (Ni0glo,...)
1053      !!----------------------------------------------------------------------
1054      INTEGER                        , INTENT(in   ) ::   kistr, kjstr   ! starting i and j position of the reading
1055      INTEGER                        , INTENT(in   ) ::   kicnt, kjcnt   ! number of points to read in i and j directions
1056      LOGICAL, DIMENSION(kicnt,kjcnt), INTENT(  out) ::   ldoce          ! ldoce(i,j) = .true. if the point (i,j) is ocean
1057      !
1058      INTEGER                          ::   inumsave                     ! local logical unit
1059      REAL(wp), DIMENSION(kicnt,kjcnt) ::   zbot, zbdy
1060      !!----------------------------------------------------------------------
1061      !
1062      inumsave = numout   ;   numout = numnul   !   redirect all print to /dev/null
1063      !
1064      IF( numbot /= -1 ) THEN
1065         CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) )
1066      ELSE
1067         zbot(:,:) = 1._wp                      ! put a non-null value
1068      ENDIF
1069      !
1070      IF( numbdy /= -1 ) THEN                   ! Adjust with bdy_msk if it exists
1071         CALL iom_get ( numbdy, jpdom_unknown,     'bdy_msk', zbdy, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) )
1072         zbot(:,:) = zbot(:,:) * zbdy(:,:)
1073      ENDIF
1074      !
1075      ldoce(:,:) = NINT(zbot(:,:)) > 0
1076      numout = inumsave
1077      !
1078   END SUBROUTINE read_mask
1079
1080
1081   SUBROUTINE mpp_getnum( ldIsOce, kproc, kipos, kjpos )
1082      !!----------------------------------------------------------------------
1083      !!                  ***  ROUTINE mpp_getnum  ***
1084      !!
1085      !! ** Purpose : give a number to each MPI subdomains (starting at 0)
1086      !!
1087      !! ** Method  : start from bottom left. First skip land subdomain, and finally use them if needed
1088      !!----------------------------------------------------------------------
1089      LOGICAL, DIMENSION(:,:), INTENT(in   ) ::   ldIsOce     ! F if land process
1090      INTEGER, DIMENSION(:,:), INTENT(  out) ::   kproc       ! subdomain number (-1 if not existing, starting at 0)
1091      INTEGER, DIMENSION(  :), INTENT(  out) ::   kipos       ! i-position of the subdomain (from 1 to jpni)
1092      INTEGER, DIMENSION(  :), INTENT(  out) ::   kjpos       ! j-position of the subdomain (from 1 to jpnj)
1093      !
1094      INTEGER :: ii, ij, jarea, iarea0
1095      INTEGER :: icont, i2add , ini, inj, inij
1096      !!----------------------------------------------------------------------
1097      !
1098      ini = SIZE(ldIsOce, dim = 1)
1099      inj = SIZE(ldIsOce, dim = 2)
1100      inij = SIZE(kipos)
1101      !
1102      ! specify which subdomains are oce subdomains; other are land subdomains
1103      kproc(:,:) = -1
1104      icont = -1
1105      DO jarea = 1, ini*inj
1106         iarea0 = jarea - 1
1107         ii = 1 + MOD(iarea0,ini)
1108         ij = 1 +     iarea0/ini
1109         IF( ldIsOce(ii,ij) ) THEN
1110            icont = icont + 1
1111            kproc(ii,ij) = icont
1112            kipos(icont+1) = ii
1113            kjpos(icont+1) = ij
1114         ENDIF
1115      END DO
1116      ! if needed add some land subdomains to reach inij active subdomains
1117      i2add = inij - COUNT( ldIsOce )
1118      DO jarea = 1, ini*inj
1119         iarea0 = jarea - 1
1120         ii = 1 + MOD(iarea0,ini)
1121         ij = 1 +     iarea0/ini
1122         IF( .NOT. ldIsOce(ii,ij) .AND. i2add > 0 ) THEN
1123            icont = icont + 1
1124            kproc(ii,ij) = icont
1125            kipos(icont+1) = ii
1126            kjpos(icont+1) = ij
1127            i2add = i2add - 1
1128         ENDIF
1129      END DO
1130      !
1131   END SUBROUTINE mpp_getnum
1132
1133
1134   SUBROUTINE init_excl_landpt
1135      !!----------------------------------------------------------------------
1136      !!                  ***  ROUTINE   ***
1137      !!
1138      !! ** Purpose : exclude exchanges which contain only land points
1139      !!
1140      !! ** Method  : if a send or receive buffer constains only land point we
1141      !!              flag off the corresponding communication
1142      !!              Warning: this selection depend on the halo size -> loop on halo size
1143      !!
1144      !!----------------------------------------------------------------------
1145      INTEGER ::   inumsave
1146      INTEGER ::   jh
1147      INTEGER ::   ipi, ipj
1148      INTEGER ::   iiwe, iiea, iist, iisz 
1149      INTEGER ::   ijso, ijno, ijst, ijsz 
1150      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zmsk
1151      LOGICAL , DIMENSION(Ni_0,Nj_0,1)      ::   lloce
1152      !!----------------------------------------------------------------------
1153      !
1154      ! read the land-sea mask on the inner domain
1155      CALL read_mask( nimpp, njmpp, Ni_0, Nj_0, lloce(:,:,1) )
1156      !
1157      ! Here we look only at communications excluding the NP folding.
1158      !   --> we switch off lbcnfd at this stage (init_nfdcom called after init_excl_landpt)...
1159      l_IdoNFold = .FALSE.
1160      !
1161      ! WARNING, see also bestpartition: minimum subdomain size defined in bestpartition according to nn_hls.
1162      ! If, one day, we want to use local halos largers than nn_hls, we must replace nn_hls by n_hlsmax in bestpartition
1163      !
1164      DO jh = 1, MIN(nn_hls, n_hlsmax)   ! different halo size
1165         !
1166         ipi = Ni_0 + 2*jh   ! local domain size
1167         ipj = Nj_0 + 2*jh
1168         !
1169         ALLOCATE( zmsk(ipi,ipj) )
1170         zmsk(jh+1:jh+Ni_0,jh+1:jh+Nj_0) = REAL(COUNT(lloce, dim = 3), wp)   ! define inner domain -> need REAL to use lbclnk
1171         CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh)                 ! fill halos
1172         ! Beware, coastal F points can be used in the code -> we may need communications for these points F points even if tmask = 0
1173         ! -> the mask we must use here is equal to 1 as soon as one of the 4 neighbours is oce (sum of the mask, not multiplication)
1174         zmsk(jh+1:jh+Ni_0,jh+1:jh+Nj_0) = zmsk(jh+1:jh+Ni_0,jh+1  :jh+Nj_0  ) + zmsk(jh+1+1:jh+Ni_0+1,jh+1  :jh+Nj_0  )   &
1175            &                            + zmsk(jh+1:jh+Ni_0,jh+1+1:jh+Nj_0+1) + zmsk(jh+1+1:jh+Ni_0+1,jh+1+1:jh+Nj_0+1)
1176         CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh)                 ! fill halos again!
1177         !       
1178         iiwe = jh   ;   iiea = Ni_0   ! bottom-left corner - 1 of the sent data
1179         ijso = jh   ;   ijno = Nj_0
1180         IF( nn_comm == 1 ) THEN
1181            iist =  0   ;   iisz = ipi
1182            ijst = jh   ;   ijsz = Nj_0
1183         ELSE
1184            iist = jh   ;   iisz = Ni_0
1185            ijst = jh   ;   ijsz = Nj_0
1186         ENDIF
1187IF( nn_comm == 1 ) THEN       ! SM: NOT WORKING FOR NEIGHBOURHOOD COLLECTIVE COMMUNICATIONS, I DON'T KNOW WHY...
1188         ! do not send if we send only land points
1189         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiSnei(jh,jpwe) = -1
1190         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiSnei(jh,jpea) = -1
1191         IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh  ) )) == 0 )   mpiSnei(jh,jpso) = -1
1192         IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh  ) )) == 0 )   mpiSnei(jh,jpno) = -1
1193         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijso+1:ijso+jh  ) )) == 0 )   mpiSnei(jh,jpsw) = -1
1194         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijso+1:ijso+jh  ) )) == 0 )   mpiSnei(jh,jpse) = -1
1195         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijno+1:ijno+jh  ) )) == 0 )   mpiSnei(jh,jpnw) = -1
1196         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijno+1:ijno+jh  ) )) == 0 )   mpiSnei(jh,jpne) = -1
1197         !
1198         iiwe = iiwe-jh   ;   iiea = iiea+jh   ! bottom-left corner - 1 of the received data
1199         ijso = ijso-jh   ;   ijno = ijno+jh
1200         ! do not send if we send only land points
1201         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiRnei(jh,jpwe) = -1
1202         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiRnei(jh,jpea) = -1
1203         IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh  ) )) == 0 )   mpiRnei(jh,jpso) = -1
1204         IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh  ) )) == 0 )   mpiRnei(jh,jpno) = -1
1205         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijso+1:ijso+jh  ) )) == 0 )   mpiRnei(jh,jpsw) = -1
1206         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijso+1:ijso+jh  ) )) == 0 )   mpiRnei(jh,jpse) = -1
1207         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijno+1:ijno+jh  ) )) == 0 )   mpiRnei(jh,jpnw) = -1
1208         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijno+1:ijno+jh  ) )) == 0 )   mpiRnei(jh,jpne) = -1
1209ENDIF
1210         !
1211         ! Specific (and rare) problem in corner treatment because we do 1st West-East comm, next South-North comm
1212         IF( nn_comm == 1 ) THEN
1213            IF( mpiSnei(jh,jpwe) > -1 )   mpiSnei(jh, (/jpsw,jpnw/) ) = -1   ! SW and NW corners already sent through West nei
1214            IF( mpiSnei(jh,jpea) > -1 )   mpiSnei(jh, (/jpse,jpne/) ) = -1   ! SE and NE corners already sent through East nei
1215            IF( mpiRnei(jh,jpso) > -1 )   mpiRnei(jh, (/jpsw,jpse/) ) = -1   ! SW and SE corners will be received through South nei
1216            IF( mpiRnei(jh,jpno) > -1 )   mpiRnei(jh, (/jpnw,jpne/) ) = -1   ! NW and NE corners will be received through North nei
1217        ENDIF
1218         !
1219         DEALLOCATE( zmsk )
1220         !
1221         CALL mpp_ini_nc(jh)    ! Initialize/Update communicator for neighbourhood collective communications
1222         !
1223      END DO
1224
1225   END SUBROUTINE init_excl_landpt
1226
1227
1228   SUBROUTINE init_ioipsl
1229      !!----------------------------------------------------------------------
1230      !!                  ***  ROUTINE init_ioipsl  ***
1231      !!
1232      !! ** Purpose :
1233      !!
1234      !! ** Method  :
1235      !!
1236      !! History :
1237      !!   9.0  !  04-03  (G. Madec )  MPP-IOIPSL
1238      !!   " "  !  08-12  (A. Coward)  addition in case of jpni*jpnj < jpnij
1239      !!----------------------------------------------------------------------
1240      INTEGER, DIMENSION(2) ::   iglo, iloc, iabsf, iabsl, ihals, ihale, idid
1241      !!----------------------------------------------------------------------
1242
1243      ! The domain is split only horizontally along i- or/and j- direction
1244      ! So we need at the most only 1D arrays with 2 elements.
1245      ! Set idompar values equivalent to the jpdom_local_noextra definition
1246      ! used in IOM. This works even if jpnij .ne. jpni*jpnj.
1247      iglo( :) = (/ Ni0glo, Nj0glo /)
1248      iloc( :) = (/ Ni_0  , Nj_0   /)
1249      iabsf(:) = (/ Nis0  , Njs0   /) + (/ nimpp, njmpp /) - 1 - nn_hls   ! corresponds to mig0(Nis0) but mig0 is not yet defined!
1250      iabsl(:) = iabsf(:) + iloc(:) - 1
1251      ihals(:) = (/ 0     , 0      /)
1252      ihale(:) = (/ 0     , 0      /)
1253      idid( :) = (/ 1     , 2      /)
1254
1255      IF(lwp) THEN
1256          WRITE(numout,*)
1257          WRITE(numout,*) 'mpp init_ioipsl :   iloc  = ', iloc
1258          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf
1259          WRITE(numout,*) '                    ihals = ', ihals
1260          WRITE(numout,*) '                    ihale = ', ihale
1261      ENDIF
1262      !
1263      CALL flio_dom_set ( jpnij, narea-1, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
1264      !
1265   END SUBROUTINE init_ioipsl
1266
1267
1268   SUBROUTINE init_nfdcom( ldwrtlay, knum )
1269      !!----------------------------------------------------------------------
1270      !!                     ***  ROUTINE  init_nfdcom  ***
1271      !! ** Purpose :   Setup for north fold exchanges with explicit
1272      !!                point-to-point messaging
1273      !!
1274      !! ** Method :   Initialization of the northern neighbours lists.
1275      !!----------------------------------------------------------------------
1276      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
1277      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
1278      !!    3.0  ! 2021-09 complete rewrite using informations from gather north fold
1279      !!----------------------------------------------------------------------
1280      LOGICAL, INTENT(in   ) ::   ldwrtlay   ! true if additional prints in layout.dat
1281      INTEGER, INTENT(in   ) ::   knum       ! layout.dat unit
1282      !
1283      REAL(wp), DIMENSION(jpi,jpj,2,4) ::   zinfo
1284      INTEGER , DIMENSION(10) ::   irknei ! too many elements but safe...
1285      INTEGER                 ::   ji, jj, jg, jn   ! dummy loop indices
1286      LOGICAL                 ::   lnew
1287      !!----------------------------------------------------------------------
1288      !
1289      IF (lwp) THEN
1290         WRITE(numout,*)
1291         WRITE(numout,*) '   ==>>>   North fold boundary prepared for jpni >1'
1292      ENDIF
1293      !
1294      CALL mpp_ini_northgather   ! we need to init the nfd with gathering in all cases as it is used to define the no-gather case
1295      !
1296      IF(ldwrtlay) THEN      ! additional prints in layout.dat
1297         WRITE(knum,*)
1298         WRITE(knum,*)
1299         WRITE(knum,*) 'Number of subdomains located along the north fold : ', ndim_rank_north
1300         WRITE(knum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
1301         DO jn = 1, ndim_rank_north, 5
1302            WRITE(knum,*) nrank_north( jn:MINVAL( (/jn+4,ndim_rank_north/) ) )
1303         END DO
1304      ENDIF
1305     
1306      nfd_nbnei = 0   ! defaul def (useless?)
1307      IF( ln_nnogather ) THEN
1308         !
1309         ! Use the "gather nfd" to know how to do the nfd: for ji point, which process send data from which of its ji-index?
1310         ! Note that nfd is perfectly symetric: I receive data from X <=> I send data to X  (-> no deadlock)
1311         !
1312         zinfo(:,:,:,:) = HUGE(0._wp)   ! default def to make sure we don't use the halos
1313         DO jg = 1, 4   ! grid type: T, U, V, F
1314            DO jj = nn_hls+1, jpj-nn_hls                ! inner domain (warning do_loop_substitute not yet defined)
1315               DO ji = nn_hls+1, jpi-nn_hls             ! inner domain (warning do_loop_substitute not yet defined)
1316                  zinfo(ji,jj,1,jg) = REAL(narea, wp)   ! mpi_rank + 1 (as default lbc_lnk fill is 0
1317                  zinfo(ji,jj,2,jg) = REAL(ji, wp)      ! ji of this proc
1318               END DO
1319            END DO
1320         END DO
1321         !
1322         ln_nnogather = .FALSE.   ! force "classical" North pole folding to fill all halos -> should be no more HUGE values...
1323         CALL lbc_lnk( 'mppini', zinfo(:,:,:,1), 'T', 1._wp )   ! Do 4 calls instead of 1 to save memory as the nogather version
1324         CALL lbc_lnk( 'mppini', zinfo(:,:,:,2), 'U', 1._wp )   ! creates buffer arrays with jpiglo as the first dimension
1325         CALL lbc_lnk( 'mppini', zinfo(:,:,:,3), 'V', 1._wp )   !
1326         CALL lbc_lnk( 'mppini', zinfo(:,:,:,4), 'F', 1._wp )   !
1327         ln_nnogather = .TRUE.
1328         
1329         IF( l_IdoNFold ) THEN   ! only the procs involed in the NFD must take care of this
1330           
1331            ALLOCATE( nfd_rksnd(jpi,4), nfd_jisnd(jpi,4) )          ! neighbour rand and remote ji-index for each grid (T, U, V, F)
1332            nfd_rksnd(:,:) = NINT( zinfo(:, jpj, 1, :) ) - 1        ! neighbour MPI rank
1333            nfd_jisnd(:,:) = NINT( zinfo(:, jpj, 2, :) ) - nn_hls   ! neighbour ji index (shifted as we don't send the halos)
1334            WHERE( nfd_rksnd == -1 )   nfd_jisnd = 1                ! use ji=1 if no neighbour, see mpp_nfd_generic.h90
1335           
1336            nfd_nbnei = 1                ! Number of neighbour sending data for the nfd. We have at least 1 neighbour!
1337            irknei(1) = nfd_rksnd(1,1)   ! which is the 1st one (I can be neighbour of myself, exclude land-proc are also ok)
1338            DO jg = 1, 4
1339               DO ji = 1, jpi     ! we must be able to fill the full line including halos
1340                  lnew = .TRUE.   ! new neighbour?
1341                  DO jn = 1, nfd_nbnei
1342                     IF( irknei(jn) == nfd_rksnd(ji,jg) )   lnew = .FALSE.   ! already found
1343                  END DO
1344                  IF( lnew ) THEN
1345                     jn = nfd_nbnei + 1
1346                     nfd_nbnei = jn
1347                     irknei(jn) = nfd_rksnd(ji,jg)
1348                  ENDIF
1349               END DO
1350            END DO
1351           
1352            ALLOCATE( nfd_rknei(nfd_nbnei) )
1353            nfd_rknei(:) = irknei(1:nfd_nbnei)
1354            ! re-number nfd_rksnd according to the indexes of nfd_rknei
1355            DO jn = 1, nfd_nbnei
1356               WHERE( nfd_rksnd == nfd_rknei(jn) )   nfd_rksnd = jn
1357            END DO
1358           
1359            IF( ldwrtlay ) THEN
1360               WRITE(knum,*)
1361               WRITE(knum,*) 'north fold exchanges with explicit point-to-point messaging :'
1362               WRITE(knum,*) '   number of neighbours for the NF: nfd_nbnei : ', nfd_nbnei
1363               IF( nfd_nbnei > 0 )   WRITE(knum,*) '   neighbours MPI ranks                       : ', nfd_rknei
1364            ENDIF
1365           
1366         ENDIF   ! l_IdoNFold
1367         !
1368      ENDIF   ! ln_nnogather
1369      !
1370   END SUBROUTINE init_nfdcom
1371
1372
1373   SUBROUTINE init_doloop
1374      !!----------------------------------------------------------------------
1375      !!                  ***  ROUTINE init_doloop  ***
1376      !!
1377      !! ** Purpose :   set the starting/ending indices of DO-loop
1378      !!              These indices are used in do_loop_substitute.h90
1379      !!----------------------------------------------------------------------
1380      !
1381      Nis0 =   1+nn_hls
1382      Njs0 =   1+nn_hls
1383      Nie0 = jpi-nn_hls
1384      Nje0 = jpj-nn_hls
1385      !
1386      Ni_0 = Nie0 - Nis0 + 1
1387      Nj_0 = Nje0 - Njs0 + 1
1388      !
1389      jpkm1 = jpk-1                             !   "           "
1390      !
1391   END SUBROUTINE init_doloop
1392
1393   
1394   SUBROUTINE init_locglo
1395      !!----------------------------------------------------------------------
1396      !!                     ***  ROUTINE init_locglo  ***
1397      !!
1398      !! ** Purpose :   initialization of global domain <--> local domain indices
1399      !!
1400      !! ** Method  :
1401      !!
1402      !! ** Action  : - mig , mjg : local  domain indices ==> global domain, including halos, indices
1403      !!              - mig0, mjg0: local  domain indices ==> global domain, excluding halos, indices
1404      !!              - mi0 , mi1 : global domain indices ==> local  domain indices
1405      !!              - mj0 , mj1   (if global point not in the local domain ==> mi0>mi1 and/or mj0>mj1)
1406      !!----------------------------------------------------------------------
1407      INTEGER ::   ji, jj   ! dummy loop argument
1408      !!----------------------------------------------------------------------
1409      !
1410      ALLOCATE( mig(jpi), mjg(jpj), mig0(jpi), mjg0(jpj) )
1411      ALLOCATE( mi0(jpiglo), mi1(jpiglo), mj0(jpjglo), mj1(jpjglo) )
1412      !
1413      DO ji = 1, jpi                 ! local domain indices ==> global domain indices, including halos
1414        mig(ji) = ji + nimpp - 1
1415      END DO
1416      DO jj = 1, jpj
1417        mjg(jj) = jj + njmpp - 1
1418      END DO
1419      !                              ! local domain indices ==> global domain indices, excluding halos
1420      !
1421      mig0(:) = mig(:) - nn_hls
1422      mjg0(:) = mjg(:) - nn_hls
1423      !                              ! global domain, including halos, indices ==> local domain indices
1424      !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the
1425      !                                   ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.
1426      DO ji = 1, jpiglo
1427        mi0(ji) = MAX( 1 , MIN( ji - nimpp + 1, jpi+1 ) )
1428        mi1(ji) = MAX( 0 , MIN( ji - nimpp + 1, jpi   ) )
1429      END DO
1430      DO jj = 1, jpjglo
1431        mj0(jj) = MAX( 1 , MIN( jj - njmpp + 1, jpj+1 ) )
1432        mj1(jj) = MAX( 0 , MIN( jj - njmpp + 1, jpj   ) )
1433      END DO
1434      !
1435   END SUBROUTINE init_locglo
1436   
1437   !!======================================================================
1438END MODULE mppini
Note: See TracBrowser for help on using the repository browser.