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 @ 15267

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

trunk: new nogather nolding, #2724

  • Property svn:keywords set to Id
File size: 70.5 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      iszimin = 3*nn_hls          ! minimum size of the MPI subdomain so halos are always adressing neighbor inner domain
716      iszjmin = 3*nn_hls
717      IF( c_NFtype == 'T' )   iszjmin = MAX(iszjmin, 2+3*nn_hls)   ! V and F folding must be outside of southern halos
718      IF( c_NFtype == 'F' )   iszjmin = MAX(iszjmin, 1+3*nn_hls)   ! V and F folding must be outside of southern halos
719      !
720      ! get the list of knbi that gives a smaller jpimax than knbi-1
721      ! get the list of knbj that gives a smaller jpjmax than knbj-1
722      DO ji = 1, inbijmax
723#if defined key_nemocice_decomp
724         iszitst = ( nx_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
725#else
726         iszitst = ( Ni0glo + (ji-1) ) / ji + 2*nn_hls   ! max subdomain i-size
727#endif
728         IF( iszitst < isziref .AND. iszitst >= iszimin ) THEN
729            isziref = iszitst
730            inbimax = inbimax + 1
731            inbi0(inbimax) = ji
732            iszi0(inbimax) = isziref
733         ENDIF
734#if defined key_nemocice_decomp
735         iszjtst = ( ny_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
736#else
737         iszjtst = ( Nj0glo + (ji-1) ) / ji + 2*nn_hls   ! max subdomain j-size
738#endif
739         IF( iszjtst < iszjref .AND. iszjtst >= iszjmin ) THEN
740            iszjref = iszjtst
741            inbjmax = inbjmax + 1
742            inbj0(inbjmax) = ji
743            iszj0(inbjmax) = iszjref
744         ENDIF
745      END DO
746      IF( inbimax == 0 ) THEN
747         WRITE(ctmp1,'(a,i2,a,i2)') '   mpp_ini bestpartition: Ni0glo (', Ni0glo, ') is too small to be used with nn_hls = ', nn_hls
748         CALL ctl_stop( 'STOP', ctmp1 )
749      ENDIF
750      IF( inbjmax == 0 ) THEN
751         WRITE(ctmp1,'(a,i2,a,i2)') '   mpp_ini bestpartition: Nj0glo (', Nj0glo, ') is too small to be used with nn_hls = ', nn_hls
752         CALL ctl_stop( 'STOP', ctmp1 )
753      ENDIF
754
755      ! combine these 2 lists to get all possible knbi*knbj <  inbijmax
756      ALLOCATE( llmsk2d(inbimax,inbjmax) )
757      DO jj = 1, inbjmax
758         DO ji = 1, inbimax
759            IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN   ;   llmsk2d(ji,jj) = .TRUE.
760            ELSE                                            ;   llmsk2d(ji,jj) = .FALSE.
761            ENDIF
762         END DO
763      END DO
764      isz1 = COUNT(llmsk2d)
765      ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) )
766      ii = 0
767      DO jj = 1, inbjmax
768         DO ji = 1, inbimax
769            IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN
770               ii = ii + 1
771               inbi1(ii) = inbi0(ji)
772               inbj1(ii) = inbj0(jj)
773               iszi1(ii) = iszi0(ji)
774               iszj1(ii) = iszj0(jj)
775            ENDIF
776         END DO
777      END DO
778      DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
779      DEALLOCATE( llmsk2d )
780
781      ALLOCATE( inbij1(isz1), iszij1(isz1) )
782      inbij1(:) = inbi1(:) * inbj1(:)
783      iszij1(:) = iszi1(:) * iszj1(:)
784
785      ! if there is no land and no print
786      IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN
787         ! get the smaller partition which gives the smallest subdomain size
788         ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1)
789         knbi = inbi1(ii)
790         knbj = inbj1(ii)
791         IF(PRESENT(knbcnt))   knbcnt = 0
792         DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 )
793         RETURN
794      ENDIF
795
796      ! extract only the partitions which reduce the subdomain size in comparison with smaller partitions
797      ALLOCATE( indexok(isz1) )                                 ! to store indices of the best partitions
798      isz0 = 0                                                  ! number of best partitions
799      inbij = 1                                                 ! start with the min value of inbij1 => 1
800      iszij = jpiglo*jpjglo+1                                   ! default: larger than global domain
801      DO WHILE( inbij <= inbijmax )                             ! if we did not reach the max of inbij1
802         ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1)   ! warning: send back the first occurence if multiple results
803         IF ( iszij1(ii) < iszij ) THEN
804            ii = MINLOC( iszi1+iszj1, mask = iszij1 == iszij1(ii) .AND. inbij1 == inbij, dim = 1)  ! select the smaller perimeter if multiple min
805            isz0 = isz0 + 1
806            indexok(isz0) = ii
807            iszij = iszij1(ii)
808         ENDIF
809         inbij = MINVAL(inbij1, mask = inbij1 > inbij)   ! warning: return largest integer value if mask = .false. everywhere
810      END DO
811      DEALLOCATE( inbij1, iszij1 )
812
813      ! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size)
814      ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) )
815      DO ji = 1, isz0
816         ii = indexok(ji)
817         inbi0(ji) = inbi1(ii)
818         inbj0(ji) = inbj1(ii)
819         iszi0(ji) = iszi1(ii)
820         iszj0(ji) = iszj1(ii)
821      END DO
822      DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 )
823
824      IF( llist ) THEN
825         IF(lwp) THEN
826            WRITE(numout,*)
827            WRITE(numout,*) '                  For your information:'
828            WRITE(numout,*) '  list of the best partitions including land supression'
829            WRITE(numout,*) '  -----------------------------------------------------'
830            WRITE(numout,*)
831         ENDIF
832         ji = isz0   ! initialization with the largest value
833         ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) )
834         CALL mpp_is_ocean( llisOce )   ! Warning: must be call by all cores (call mpp_sum)
835         inbijold = COUNT(llisOce)
836         DEALLOCATE( llisOce )
837         DO ji =isz0-1,1,-1
838            ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) )
839            CALL mpp_is_ocean( llisOce )   ! Warning: must be call by all cores (call mpp_sum)
840            inbij = COUNT(llisOce)
841            DEALLOCATE( llisOce )
842            IF(lwp .AND. inbij < inbijold) THEN
843               WRITE(numout,'(a, i6, a, i6, a, f4.1, a, i9, a, i6, a, i6, a)')                                 &
844                  &   'nb_cores oce: ', inbij, ', land domains excluded: ', inbi0(ji)*inbj0(ji) - inbij,       &
845                  &   ' (', REAL(inbi0(ji)*inbj0(ji) - inbij,wp) / REAL(inbi0(ji)*inbj0(ji),wp) *100.,         &
846                  &   '%), largest oce domain: ', iszi0(ji)*iszj0(ji), ' ( ', iszi0(ji),' x ', iszj0(ji), ' )'
847               inbijold = inbij
848            ENDIF
849         END DO
850         DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
851         IF(lwp) THEN
852            WRITE(numout,*)
853            WRITE(numout,*)  '  -----------------------------------------------------------'
854         ENDIF
855         CALL mppsync
856         CALL mppstop( ld_abort = .TRUE. )
857      ENDIF
858
859      DEALLOCATE( iszi0, iszj0 )
860      inbij = inbijmax + 1        ! default: larger than possible
861      ii = isz0+1                 ! start from the end of the list (smaller subdomains)
862      DO WHILE( inbij > knbij )   ! while the number of ocean subdomains exceed the number of procs
863         ii = ii -1
864         ALLOCATE( llisOce(inbi0(ii), inbj0(ii)) )
865         CALL mpp_is_ocean( llisOce )            ! must be done by all core
866         inbij = COUNT(llisOce)
867         DEALLOCATE( llisOce )
868      END DO
869      knbi = inbi0(ii)
870      knbj = inbj0(ii)
871      IF(PRESENT(knbcnt))   knbcnt = knbi * knbj - inbij
872      DEALLOCATE( inbi0, inbj0 )
873      !
874   END SUBROUTINE bestpartition
875
876
877   SUBROUTINE mpp_init_landprop( propland )
878      !!----------------------------------------------------------------------
879      !!                  ***  ROUTINE mpp_init_landprop  ***
880      !!
881      !! ** Purpose : the the proportion of land points in the surface land-sea mask
882      !!
883      !! ** Method  : read iproc strips (of length Ni0glo) of the land-sea mask
884      !!----------------------------------------------------------------------
885      REAL(wp), INTENT(  out) :: propland    ! proportion of land points in the global domain (between 0 and 1)
886      !
887      INTEGER, DIMENSION(jpni*jpnj) ::   kusedom_1d
888      INTEGER :: inboce, iarea
889      INTEGER :: iproc, idiv, ijsz
890      INTEGER :: ijstr
891      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce
892      !!----------------------------------------------------------------------
893      ! do nothing if there is no land-sea mask
894      IF( numbot == -1 .and. numbdy == -1 ) THEN
895         propland = 0.
896         RETURN
897      ENDIF
898
899      ! number of processes reading the bathymetry file
900      iproc = MINVAL( (/mppsize, Nj0glo/2, 100/) )  ! read a least 2 lines, no more that 100 processes reading at the same time
901
902      ! we want to read iproc strips of the land-sea mask. -> pick up iproc processes every idiv processes starting at 1
903      IF( iproc == 1 ) THEN   ;   idiv = mppsize
904      ELSE                    ;   idiv = ( mppsize - 1 ) / ( iproc - 1 )
905      ENDIF
906
907      iarea = (narea-1)/idiv   ! involed process number (starting counting at 0)
908      IF( MOD( narea-1, idiv ) == 0 .AND. iarea < iproc ) THEN   ! beware idiv can be = to 1
909         !
910         ijsz = Nj0glo / iproc                                               ! width of the stripe to read
911         IF( iarea < MOD(Nj0glo,iproc) ) ijsz = ijsz + 1
912         ijstr = iarea*(Nj0glo/iproc) + MIN(iarea, MOD(Nj0glo,iproc)) + 1    ! starting j position of the reading
913         !
914         ALLOCATE( lloce(Ni0glo, ijsz) )                                     ! allocate the strip
915         CALL read_mask( 1, ijstr, Ni0glo, ijsz, lloce )
916         inboce = COUNT(lloce)                                               ! number of ocean point in the stripe
917         DEALLOCATE(lloce)
918         !
919      ELSE
920         inboce = 0
921      ENDIF
922      CALL mpp_sum( 'mppini', inboce )   ! total number of ocean points over the global domain
923      !
924      propland = REAL( Ni0glo*Nj0glo - inboce, wp ) / REAL( Ni0glo*Nj0glo, wp )
925      !
926   END SUBROUTINE mpp_init_landprop
927
928
929   SUBROUTINE mpp_is_ocean( ldIsOce )
930      !!----------------------------------------------------------------------
931      !!                  ***  ROUTINE mpp_is_ocean  ***
932      !!
933      !! ** Purpose : Check for a mpi domain decomposition inbi x inbj which
934      !!              subdomains, including 1 halo (even if nn_hls>1), contain
935      !!              at least 1 ocean point.
936      !!              We must indeed ensure that each subdomain that is a neighbour
937      !!              of a land subdomain, has only land points on its boundary
938      !!              (inside the inner subdomain) with the land subdomain.
939      !!              This is needed to get the proper bondary conditions on
940      !!              a subdomain with a closed boundary.
941      !!
942      !! ** Method  : read inbj strips (of length Ni0glo) of the land-sea mask
943      !!----------------------------------------------------------------------
944      LOGICAL, DIMENSION(:,:), INTENT(  out) ::   ldIsOce        ! .true. if a sub domain constains 1 ocean point
945      !
946      INTEGER :: idiv, iimax, ijmax, iarea
947      INTEGER :: inbi, inbj, inx, iny, inry, isty
948      INTEGER :: ji, jn
949      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   inboce           ! number oce oce pint in each mpi subdomain
950      INTEGER, ALLOCATABLE, DIMENSION(:  ) ::   inboce_1d
951      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ijpi
952      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ijpj
953      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce            ! lloce(i,j) = .true. if the point (i,j) is ocean
954      !!----------------------------------------------------------------------
955      ! do nothing if there is no land-sea mask
956      IF( numbot == -1 .AND. numbdy == -1 ) THEN
957         ldIsOce(:,:) = .TRUE.
958         RETURN
959      ENDIF
960      !
961      inbi = SIZE( ldIsOce, dim = 1 )
962      inbj = SIZE( ldIsOce, dim = 2 )
963      !
964      ! we want to read inbj strips of the land-sea mask. -> pick up inbj processes every idiv processes starting at 1
965      IF           ( inbj == 1 ) THEN   ;   idiv = mppsize
966      ELSE IF ( mppsize < inbj ) THEN   ;   idiv = 1
967      ELSE                              ;   idiv = ( mppsize - 1 ) / ( inbj - 1 )
968      ENDIF
969      !
970      ALLOCATE( inboce(inbi,inbj), inboce_1d(inbi*inbj) )
971      inboce(:,:) = 0          ! default no ocean point found
972      !
973      DO jn = 0, (inbj-1)/mppsize   ! if mppsize < inbj : more strips than mpi processes (because of potential land domains)
974         !
975         iarea = (narea-1)/idiv + jn * mppsize + 1                     ! involed process number (starting counting at 1)
976         IF( MOD( narea-1, idiv ) == 0 .AND. iarea <= inbj ) THEN      ! beware idiv can be = to 1
977            !
978            ALLOCATE( iimppt(inbi,inbj), ijmppt(inbi,inbj), ijpi(inbi,inbj), ijpj(inbi,inbj) )
979            CALL mpp_basesplit( Ni0glo, Nj0glo, 0, inbi, inbj, iimax, ijmax, iimppt, ijmppt, ijpi, ijpj )
980            !
981            inx = Ni0glo + 2   ;   iny = ijpj(1,iarea) + 2             ! strip size + 1 halo on each direction (even if nn_hls>1)
982            ALLOCATE( lloce(inx, iny) )                                ! allocate the strip
983            inry = iny - COUNT( (/ iarea == 1, iarea == inbj /) )      ! number of point to read in y-direction
984            isty = 1 + COUNT( (/ iarea == 1 /) )                       ! read from the first or the second line?
985            CALL read_mask( 1, ijmppt(1,iarea) - 2 + isty, Ni0glo, inry, lloce(2:inx-1, isty:inry+isty-1) )   ! read the strip
986            !
987            IF( iarea == 1    ) THEN                                   ! the first line was not read
988               IF( l_Jperio ) THEN                                     !   north-south periodocity
989                  CALL read_mask( 1, Nj0glo, Ni0glo, 1, lloce(2:inx-1, 1) )   !   read the last line -> first line of lloce
990               ELSE
991                  lloce(2:inx-1,  1) = .FALSE.                         !   closed boundary
992               ENDIF
993            ENDIF
994            IF( iarea == inbj ) THEN                                   ! the last line was not read
995               IF( l_Jperio ) THEN                                     !   north-south periodocity
996                  CALL read_mask( 1, 1, Ni0glo, 1, lloce(2:inx-1,iny) )   !      read the first line -> last line of lloce
997               ELSEIF( c_NFtype == 'T' ) THEN                          !   north-pole folding T-pivot, T-point
998                  lloce(2,iny) = lloce(2,iny-2)                        !      here we have 1 halo (even if nn_hls>1)
999                  DO ji = 3,inx-1
1000                     lloce(ji,iny  ) = lloce(inx-ji+2,iny-2)           !      ok, we have at least 3 lines
1001                  END DO
1002                  DO ji = inx/2+2,inx-1
1003                     lloce(ji,iny-1) = lloce(inx-ji+2,iny-1)
1004                  END DO
1005               ELSEIF( c_NFtype == 'F' ) THEN                          !   north-pole folding F-pivot, T-point, 1 halo
1006                  lloce(inx/2+1,iny-1) = lloce(inx/2,iny-1)            !      here we have 1 halo (even if nn_hls>1)
1007                  lloce(inx  -1,iny-1) = lloce(2    ,iny-1)
1008                  DO ji = 2,inx-1
1009                     lloce(ji,iny) = lloce(inx-ji+1,iny-1)
1010                  END DO
1011               ELSE                                                    !   closed boundary
1012                  lloce(2:inx-1,iny) = .FALSE.
1013               ENDIF
1014            ENDIF
1015            !                                                          ! first and last column were not read
1016            IF( l_Iperio ) THEN
1017               lloce(1,:) = lloce(inx-1,:)   ;   lloce(inx,:) = lloce(2,:)   ! east-west periodocity
1018            ELSE
1019               lloce(1,:) = .FALSE.          ;   lloce(inx,:) = .FALSE.      ! closed boundary
1020            ENDIF
1021            !
1022            DO  ji = 1, inbi
1023               inboce(ji,iarea) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ijpi(ji,1)+1,:) )   ! lloce as 2 points more than Ni0glo
1024            END DO
1025            !
1026            DEALLOCATE(lloce)
1027            DEALLOCATE(iimppt, ijmppt, ijpi, ijpj)
1028            !
1029         ENDIF
1030      END DO
1031
1032      inboce_1d = RESHAPE(inboce, (/ inbi*inbj /))
1033      CALL mpp_sum( 'mppini', inboce_1d )
1034      inboce = RESHAPE(inboce_1d, (/inbi, inbj/))
1035      ldIsOce(:,:) = inboce(:,:) /= 0
1036      DEALLOCATE(inboce, inboce_1d)
1037      !
1038   END SUBROUTINE mpp_is_ocean
1039
1040
1041   SUBROUTINE read_mask( kistr, kjstr, kicnt, kjcnt, ldoce )
1042      !!----------------------------------------------------------------------
1043      !!                  ***  ROUTINE read_mask  ***
1044      !!
1045      !! ** Purpose : Read relevant bathymetric information in order to
1046      !!              provide a land/sea mask used for the elimination
1047      !!              of land domains, in an mpp computation.
1048      !!
1049      !! ** Method  : read stipe of size (Ni0glo,...)
1050      !!----------------------------------------------------------------------
1051      INTEGER                        , INTENT(in   ) ::   kistr, kjstr   ! starting i and j position of the reading
1052      INTEGER                        , INTENT(in   ) ::   kicnt, kjcnt   ! number of points to read in i and j directions
1053      LOGICAL, DIMENSION(kicnt,kjcnt), INTENT(  out) ::   ldoce          ! ldoce(i,j) = .true. if the point (i,j) is ocean
1054      !
1055      INTEGER                          ::   inumsave                     ! local logical unit
1056      REAL(wp), DIMENSION(kicnt,kjcnt) ::   zbot, zbdy
1057      !!----------------------------------------------------------------------
1058      !
1059      inumsave = numout   ;   numout = numnul   !   redirect all print to /dev/null
1060      !
1061      IF( numbot /= -1 ) THEN
1062         CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) )
1063      ELSE
1064         zbot(:,:) = 1._wp                      ! put a non-null value
1065      ENDIF
1066      !
1067      IF( numbdy /= -1 ) THEN                   ! Adjust with bdy_msk if it exists
1068         CALL iom_get ( numbdy, jpdom_unknown,     'bdy_msk', zbdy, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) )
1069         zbot(:,:) = zbot(:,:) * zbdy(:,:)
1070      ENDIF
1071      !
1072      ldoce(:,:) = NINT(zbot(:,:)) > 0
1073      numout = inumsave
1074      !
1075   END SUBROUTINE read_mask
1076
1077
1078   SUBROUTINE mpp_getnum( ldIsOce, kproc, kipos, kjpos )
1079      !!----------------------------------------------------------------------
1080      !!                  ***  ROUTINE mpp_getnum  ***
1081      !!
1082      !! ** Purpose : give a number to each MPI subdomains (starting at 0)
1083      !!
1084      !! ** Method  : start from bottom left. First skip land subdomain, and finally use them if needed
1085      !!----------------------------------------------------------------------
1086      LOGICAL, DIMENSION(:,:), INTENT(in   ) ::   ldIsOce     ! F if land process
1087      INTEGER, DIMENSION(:,:), INTENT(  out) ::   kproc       ! subdomain number (-1 if not existing, starting at 0)
1088      INTEGER, DIMENSION(  :), INTENT(  out) ::   kipos       ! i-position of the subdomain (from 1 to jpni)
1089      INTEGER, DIMENSION(  :), INTENT(  out) ::   kjpos       ! j-position of the subdomain (from 1 to jpnj)
1090      !
1091      INTEGER :: ii, ij, jarea, iarea0
1092      INTEGER :: icont, i2add , ini, inj, inij
1093      !!----------------------------------------------------------------------
1094      !
1095      ini = SIZE(ldIsOce, dim = 1)
1096      inj = SIZE(ldIsOce, dim = 2)
1097      inij = SIZE(kipos)
1098      !
1099      ! specify which subdomains are oce subdomains; other are land subdomains
1100      kproc(:,:) = -1
1101      icont = -1
1102      DO jarea = 1, ini*inj
1103         iarea0 = jarea - 1
1104         ii = 1 + MOD(iarea0,ini)
1105         ij = 1 +     iarea0/ini
1106         IF( ldIsOce(ii,ij) ) THEN
1107            icont = icont + 1
1108            kproc(ii,ij) = icont
1109            kipos(icont+1) = ii
1110            kjpos(icont+1) = ij
1111         ENDIF
1112      END DO
1113      ! if needed add some land subdomains to reach inij active subdomains
1114      i2add = inij - COUNT( ldIsOce )
1115      DO jarea = 1, ini*inj
1116         iarea0 = jarea - 1
1117         ii = 1 + MOD(iarea0,ini)
1118         ij = 1 +     iarea0/ini
1119         IF( .NOT. ldIsOce(ii,ij) .AND. i2add > 0 ) THEN
1120            icont = icont + 1
1121            kproc(ii,ij) = icont
1122            kipos(icont+1) = ii
1123            kjpos(icont+1) = ij
1124            i2add = i2add - 1
1125         ENDIF
1126      END DO
1127      !
1128   END SUBROUTINE mpp_getnum
1129
1130
1131   SUBROUTINE init_excl_landpt
1132      !!----------------------------------------------------------------------
1133      !!                  ***  ROUTINE   ***
1134      !!
1135      !! ** Purpose : exclude exchanges which contain only land points
1136      !!
1137      !! ** Method  : if a send or receive buffer constains only land point we
1138      !!              flag off the corresponding communication
1139      !!              Warning: this selection depend on the halo size -> loop on halo size
1140      !!
1141      !!----------------------------------------------------------------------
1142      INTEGER ::   inumsave
1143      INTEGER ::   jh
1144      INTEGER ::   ipi, ipj
1145      INTEGER ::   iiwe, iiea, iist, iisz 
1146      INTEGER ::   ijso, ijno, ijst, ijsz 
1147      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zmsk
1148      LOGICAL , DIMENSION(Ni_0,Nj_0,1)      ::   lloce
1149      !!----------------------------------------------------------------------
1150      !
1151      ! read the land-sea mask on the inner domain
1152      CALL read_mask( nimpp, njmpp, Ni_0, Nj_0, lloce(:,:,1) )
1153      !
1154      ! Here we look only at communications excluding the NP folding.
1155      !   --> we switch off lbcnfd at this stage (init_nfdcom called after init_excl_landpt)...
1156      l_IdoNFold = .FALSE.
1157      !
1158      DO jh = 1, n_hlsmax    ! different halo size
1159         !
1160         ipi = Ni_0 + 2*jh   ! local domain size
1161         ipj = Nj_0 + 2*jh
1162         !
1163         ALLOCATE( zmsk(ipi,ipj) )
1164         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
1165         CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh)                 ! fill halos
1166         !       
1167         iiwe = jh   ;   iiea = Ni_0   ! bottom-left corfer - 1 of the sent data
1168         ijso = jh   ;   ijno = Nj_0
1169         IF( nn_comm == 1 ) THEN
1170            iist =  0   ;   iisz = ipi
1171            ijst =  0   ;   ijsz = ipj
1172         ELSE
1173            iist = jh   ;   iisz = Ni_0
1174            ijst = jh   ;   ijsz = Nj_0
1175         ENDIF
1176IF( nn_comm == 1 ) THEN       ! SM: NOT WORKING FOR NEIGHBOURHOOD COLLECTIVE COMMUNICATIONS, I DON'T KNOW WHY...
1177         ! do not send if we send only land points
1178         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiSnei(jh,jpwe) = -1
1179         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiSnei(jh,jpea) = -1
1180         IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh  ) )) == 0 )   mpiSnei(jh,jpso) = -1
1181         IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh  ) )) == 0 )   mpiSnei(jh,jpno) = -1
1182         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijso+1:ijso+jh  ) )) == 0 )   mpiSnei(jh,jpsw) = -1
1183         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijso+1:ijso+jh  ) )) == 0 )   mpiSnei(jh,jpse) = -1
1184         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijno+1:ijno+jh  ) )) == 0 )   mpiSnei(jh,jpnw) = -1
1185         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijno+1:ijno+jh  ) )) == 0 )   mpiSnei(jh,jpne) = -1
1186         !
1187         iiwe = iiwe-jh   ;   iiea = iiea+jh   ! bottom-left corfer - 1 of the received data
1188         ijso = ijso-jh   ;   ijno = ijno+jh
1189         ! do not send if we send only land points
1190         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiRnei(jh,jpwe) = -1
1191         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijst+1:ijst+ijsz) )) == 0 )   mpiRnei(jh,jpea) = -1
1192         IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh  ) )) == 0 )   mpiRnei(jh,jpso) = -1
1193         IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh  ) )) == 0 )   mpiRnei(jh,jpno) = -1
1194         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijso+1:ijso+jh  ) )) == 0 )   mpiRnei(jh,jpsw) = -1
1195         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijso+1:ijso+jh  ) )) == 0 )   mpiRnei(jh,jpse) = -1
1196         IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh  ,ijno+1:ijno+jh  ) )) == 0 )   mpiRnei(jh,jpnw) = -1
1197         IF( NINT(SUM( zmsk(iiea+1:iiea+jh  ,ijno+1:ijno+jh  ) )) == 0 )   mpiRnei(jh,jpne) = -1
1198ENDIF
1199         !
1200         ! Specific (and rare) problem in corner treatment because we do 1st West-East comm, next South-North comm
1201         IF( nn_comm == 1 ) THEN
1202            IF( mpiSnei(jh,jpwe) > -1 )   mpiSnei(jh, (/jpsw,jpnw/) ) = -1   ! SW and NW corners already sent through West nei
1203            IF( mpiSnei(jh,jpea) > -1 )   mpiSnei(jh, (/jpse,jpne/) ) = -1   ! SE and NE corners already sent through East nei
1204            IF( mpiRnei(jh,jpso) > -1 )   mpiRnei(jh, (/jpsw,jpse/) ) = -1   ! SW and SE corners will be received through South nei
1205            IF( mpiRnei(jh,jpno) > -1 )   mpiRnei(jh, (/jpnw,jpne/) ) = -1   ! NW and NE corners will be received through North nei
1206        ENDIF
1207         !
1208         DEALLOCATE( zmsk )
1209         !
1210         CALL mpp_ini_nc(jh)    ! Initialize/Update communicator for neighbourhood collective communications
1211         !
1212      END DO
1213
1214   END SUBROUTINE init_excl_landpt
1215
1216
1217   SUBROUTINE init_ioipsl
1218      !!----------------------------------------------------------------------
1219      !!                  ***  ROUTINE init_ioipsl  ***
1220      !!
1221      !! ** Purpose :
1222      !!
1223      !! ** Method  :
1224      !!
1225      !! History :
1226      !!   9.0  !  04-03  (G. Madec )  MPP-IOIPSL
1227      !!   " "  !  08-12  (A. Coward)  addition in case of jpni*jpnj < jpnij
1228      !!----------------------------------------------------------------------
1229      INTEGER, DIMENSION(2) ::   iglo, iloc, iabsf, iabsl, ihals, ihale, idid
1230      !!----------------------------------------------------------------------
1231
1232      ! The domain is split only horizontally along i- or/and j- direction
1233      ! So we need at the most only 1D arrays with 2 elements.
1234      ! Set idompar values equivalent to the jpdom_local_noextra definition
1235      ! used in IOM. This works even if jpnij .ne. jpni*jpnj.
1236      iglo( :) = (/ Ni0glo, Nj0glo /)
1237      iloc( :) = (/ Ni_0  , Nj_0   /)
1238      iabsf(:) = (/ Nis0  , Njs0   /) + (/ nimpp, njmpp /) - 1 - nn_hls   ! corresponds to mig0(Nis0) but mig0 is not yet defined!
1239      iabsl(:) = iabsf(:) + iloc(:) - 1
1240      ihals(:) = (/ 0     , 0      /)
1241      ihale(:) = (/ 0     , 0      /)
1242      idid( :) = (/ 1     , 2      /)
1243
1244      IF(lwp) THEN
1245          WRITE(numout,*)
1246          WRITE(numout,*) 'mpp init_ioipsl :   iloc  = ', iloc
1247          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf
1248          WRITE(numout,*) '                    ihals = ', ihals
1249          WRITE(numout,*) '                    ihale = ', ihale
1250      ENDIF
1251      !
1252      CALL flio_dom_set ( jpnij, narea-1, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
1253      !
1254   END SUBROUTINE init_ioipsl
1255
1256
1257   SUBROUTINE init_nfdcom( ldwrtlay, knum )
1258      !!----------------------------------------------------------------------
1259      !!                     ***  ROUTINE  init_nfdcom  ***
1260      !! ** Purpose :   Setup for north fold exchanges with explicit
1261      !!                point-to-point messaging
1262      !!
1263      !! ** Method :   Initialization of the northern neighbours lists.
1264      !!----------------------------------------------------------------------
1265      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
1266      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
1267      !!    3.0  ! 2021-09 complete rewrite using informations from gather north fold
1268      !!----------------------------------------------------------------------
1269      LOGICAL, INTENT(in   ) ::   ldwrtlay   ! true if additional prints in layout.dat
1270      INTEGER, INTENT(in   ) ::   knum       ! layout.dat unit
1271      !
1272      REAL(wp), DIMENSION(jpi,jpj,2,4) ::   zinfo
1273      INTEGER , DIMENSION(10) ::   irknei ! too many elements but safe...
1274      INTEGER                 ::   ji, jj, jg, jn   ! dummy loop indices
1275      LOGICAL                 ::   lnew
1276      !!----------------------------------------------------------------------
1277      !
1278      IF (lwp) THEN
1279         WRITE(numout,*)
1280         WRITE(numout,*) '   ==>>>   North fold boundary prepared for jpni >1'
1281      ENDIF
1282      !
1283      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
1284      !
1285      IF(ldwrtlay) THEN      ! additional prints in layout.dat
1286         WRITE(knum,*)
1287         WRITE(knum,*)
1288         WRITE(knum,*) 'Number of subdomains located along the north fold : ', ndim_rank_north
1289         WRITE(knum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
1290         DO jn = 1, ndim_rank_north, 5
1291            WRITE(knum,*) nrank_north( jn:MINVAL( (/jn+4,ndim_rank_north/) ) )
1292         END DO
1293      ENDIF
1294     
1295      nfd_nbnei = 0   ! defaul def (useless?)
1296      IF( ln_nnogather ) THEN
1297         !
1298         ! Use the "gather nfd" to know how to do the nfd: for ji point, which process send data from which of its ji-index?
1299         ! Note that nfd is perfectly symetric: I receive data from X <=> I send data to X  (-> no deadlock)
1300         !
1301         zinfo(:,:,:,:) = HUGE(0._wp)   ! default def to make sure we don't use the halos
1302         DO jg = 1, 4   ! grid type: T, U, V, F
1303            DO jj = nn_hls+1, jpj-nn_hls                ! inner domain (warning do_loop_substitute not yet defined)
1304               DO ji = nn_hls+1, jpi-nn_hls             ! inner domain (warning do_loop_substitute not yet defined)
1305                  zinfo(ji,jj,1,jg) = REAL(narea, wp)   ! mpi_rank + 1 (as default lbc_lnk fill is 0
1306                  zinfo(ji,jj,2,jg) = REAL(ji, wp)      ! ji of this proc
1307               END DO
1308            END DO
1309         END DO
1310         !
1311         ln_nnogather = .FALSE.   ! force "classical" North pole folding to fill all halos -> should be no more HUGE values...
1312         CALL lbc_lnk( 'mppini', zinfo(:,:,:,1), 'T', 1._wp )   ! Do 4 calls instead of 1 to save memory as the nogather version
1313         CALL lbc_lnk( 'mppini', zinfo(:,:,:,2), 'U', 1._wp )   ! creates buffer arrays with jpiglo as the first dimension
1314         CALL lbc_lnk( 'mppini', zinfo(:,:,:,3), 'V', 1._wp )   !
1315         CALL lbc_lnk( 'mppini', zinfo(:,:,:,4), 'F', 1._wp )   !
1316         ln_nnogather = .TRUE.
1317         
1318         IF( l_IdoNFold ) THEN   ! only the procs involed in the NFD must take care of this
1319           
1320            ALLOCATE( nfd_rksnd(jpi,4), nfd_jisnd(jpi,4) )          ! neighbour rand and remote ji-index for each grid (T, U, V, F)
1321            nfd_rksnd(:,:) = NINT( zinfo(:, jpj, 1, :) ) - 1        ! neighbour MPI rank
1322            nfd_jisnd(:,:) = NINT( zinfo(:, jpj, 2, :) ) - nn_hls   ! neighbour ji index (shifted as we don't send the halos)
1323            WHERE( nfd_rksnd == -1 )   nfd_jisnd = 1                ! use ji=1 if no neighbour, see mpp_nfd_generic.h90
1324           
1325            nfd_nbnei = 1                ! Number of neighbour sending data for the nfd. We have at least 1 neighbour!
1326            irknei(1) = nfd_rksnd(1,1)   ! which is the 1st one (I can be neighbour of myself, exclude land-proc are also ok)
1327            DO jg = 1, 4
1328               DO ji = 1, jpi     ! we must be able to fill the full line including halos
1329                  lnew = .TRUE.   ! new neighbour?
1330                  DO jn = 1, nfd_nbnei
1331                     IF( irknei(jn) == nfd_rksnd(ji,jg) )   lnew = .FALSE.   ! already found
1332                  END DO
1333                  IF( lnew ) THEN
1334                     jn = nfd_nbnei + 1
1335                     nfd_nbnei = jn
1336                     irknei(jn) = nfd_rksnd(ji,jg)
1337                  ENDIF
1338               END DO
1339            END DO
1340           
1341            ALLOCATE( nfd_rknei(nfd_nbnei) )
1342            nfd_rknei(:) = irknei(1:nfd_nbnei)
1343            ! re-number nfd_rksnd according to the indexes of nfd_rknei
1344            DO jn = 1, nfd_nbnei
1345               WHERE( nfd_rksnd == nfd_rknei(jn) )   nfd_rksnd = jn
1346            END DO
1347           
1348            IF( ldwrtlay ) THEN
1349               WRITE(knum,*)
1350               WRITE(knum,*) 'north fold exchanges with explicit point-to-point messaging :'
1351               WRITE(knum,*) '   number of neighbours for the NF: nfd_nbnei : ', nfd_nbnei
1352               IF( nfd_nbnei > 0 )   WRITE(knum,*) '   neighbours MPI ranks                       : ', nfd_rknei
1353            ENDIF
1354           
1355         ENDIF   ! l_IdoNFold
1356         !
1357      ENDIF   ! ln_nnogather
1358      !
1359   END SUBROUTINE init_nfdcom
1360
1361
1362   SUBROUTINE init_doloop
1363      !!----------------------------------------------------------------------
1364      !!                  ***  ROUTINE init_doloop  ***
1365      !!
1366      !! ** Purpose :   set the starting/ending indices of DO-loop
1367      !!              These indices are used in do_loop_substitute.h90
1368      !!----------------------------------------------------------------------
1369      !
1370      Nis0 =   1+nn_hls
1371      Njs0 =   1+nn_hls
1372      Nie0 = jpi-nn_hls
1373      Nje0 = jpj-nn_hls
1374      !
1375      Ni_0 = Nie0 - Nis0 + 1
1376      Nj_0 = Nje0 - Njs0 + 1
1377      !
1378      jpkm1 = jpk-1                             !   "           "
1379      !
1380   END SUBROUTINE init_doloop
1381
1382   
1383   SUBROUTINE init_locglo
1384      !!----------------------------------------------------------------------
1385      !!                     ***  ROUTINE init_locglo  ***
1386      !!
1387      !! ** Purpose :   initialization of global domain <--> local domain indices
1388      !!
1389      !! ** Method  :
1390      !!
1391      !! ** Action  : - mig , mjg : local  domain indices ==> global domain, including halos, indices
1392      !!              - mig0, mjg0: local  domain indices ==> global domain, excluding halos, indices
1393      !!              - mi0 , mi1 : global domain indices ==> local  domain indices
1394      !!              - mj0 , mj1   (if global point not in the local domain ==> mi0>mi1 and/or mj0>mj1)
1395      !!----------------------------------------------------------------------
1396      INTEGER ::   ji, jj   ! dummy loop argument
1397      !!----------------------------------------------------------------------
1398      !
1399      ALLOCATE( mig(jpi), mjg(jpj), mig0(jpi), mjg0(jpj) )
1400      ALLOCATE( mi0(jpiglo), mi1(jpiglo), mj0(jpjglo), mj1(jpjglo) )
1401      !
1402      DO ji = 1, jpi                 ! local domain indices ==> global domain indices, including halos
1403        mig(ji) = ji + nimpp - 1
1404      END DO
1405      DO jj = 1, jpj
1406        mjg(jj) = jj + njmpp - 1
1407      END DO
1408      !                              ! local domain indices ==> global domain indices, excluding halos
1409      !
1410      mig0(:) = mig(:) - nn_hls
1411      mjg0(:) = mjg(:) - nn_hls
1412      !                              ! global domain, including halos, indices ==> local domain indices
1413      !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the
1414      !                                   ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.
1415      DO ji = 1, jpiglo
1416        mi0(ji) = MAX( 1 , MIN( ji - nimpp + 1, jpi+1 ) )
1417        mi1(ji) = MAX( 0 , MIN( ji - nimpp + 1, jpi   ) )
1418      END DO
1419      DO jj = 1, jpjglo
1420        mj0(jj) = MAX( 1 , MIN( jj - njmpp + 1, jpj+1 ) )
1421        mj1(jj) = MAX( 0 , MIN( jj - njmpp + 1, jpj   ) )
1422      END DO
1423      !
1424   END SUBROUTINE init_locglo
1425   
1426   !!======================================================================
1427END MODULE mppini
Note: See TracBrowser for help on using the repository browser.