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/branches/2021/dev_r14312_MPI_Interface/src/OCE/LBC – NEMO

source: NEMO/branches/2021/dev_r14312_MPI_Interface/src/OCE/LBC/mppini.F90 @ 14314

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

dev_r14312_MPI_Interface: first implementation, #2598

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