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

Last change on this file since 10570 was 10570, checked in by acc, 19 months ago

Trunk update to implement finer control over the choice of text report files generated. See ticket: #2167

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