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

Last change on this file was 13490, checked in by smasson, 2 months ago

trunk: fix bestpartition in mppini, see #2526

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