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

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

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

Last change on this file since 13383 was 13305, checked in by acc, 4 years ago

Trunk changes required to avoid issues with the outer halo in ORCA2_ICE_PISCES,
REPRO_8_4 tests with nn_hls=2. These changes ensure that tmask and output
from sbc_blk are set correctly in the outer halo. Failure to set valid
values in the outer halo can generate NaNs? which lead to OOB errors in the
XRGB lookup table used for the TRC optics.See #2366 for details. With these
changes all variants of the ORCA2_ICE_PISCES SETTE test will complete. There
are still differences between the 1 and 2 halo width runs but running with:
no land suppression; partial land suppression or full land suppression does
not alter either set of results. Likewise setting ln_nnogather either true
or false does not alter results. Differences in run.stat start after 140
timesteps and differences in tracer.stat start after 60 timesteps between
the different halo width sets. Equivalent tests with ln_icebergs = F show no
differences in run.stat but halo-width dependent differences in tracer.stat
persist (now after 64 timesteps).

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