source: utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/mppini.F90 @ 13109

Last change on this file since 13109 was 13109, checked in by rblod, 3 months ago

ticket #2129 : major corrections in domcfg

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