source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LBC/mppini.F90 @ 13219

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

better e3: update with trunk@13218 see #2385

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