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

Last change on this file since 10542 was 10542, checked in by smasson, 19 months ago

trunk: improve NP folding balance, see #2214

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