source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/LBC/mppini.F90 @ 13130

Last change on this file since 13130 was 13130, checked in by smasson, 3 months ago

Extra_Halo: supress halos from outputs and coupling, see #2366

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