source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/LBC/mppini.F90 @ 10345

Last change on this file since 10345 was 10345, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10344, see #2133

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