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

Last change on this file since 9796 was 9796, checked in by mathiot, 2 years ago

fix #2099

  • Property svn:keywords set to Id
File size: 39.5 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   !!----------------------------------------------------------------------
39   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (./LICENSE)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45#if ! defined key_mpp_mpi
46   !!----------------------------------------------------------------------
47   !!   Default option :                            shared memory computing
48   !!----------------------------------------------------------------------
49
50   SUBROUTINE mpp_init
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE mpp_init  ***
53      !!
54      !! ** Purpose :   Lay out the global domain over processors.
55      !!
56      !! ** Method  :   Shared memory computing, set the local processor
57      !!              variables to the value of the global domain
58      !!----------------------------------------------------------------------
59      !
60      jpimax = jpiglo
61      jpjmax = jpjglo
62      jpi    = jpiglo
63      jpj    = jpjglo
64      jpk    = jpjglo
65      jpim1  = jpi-1                                            ! inner domain indices
66      jpjm1  = jpj-1                                            !   "           "
67      jpkm1  = MAX( 1, jpk-1 )                                  !   "           "
68      jpij   = jpi*jpj
69      jpni   = 1
70      jpnj   = 1
71      jpnij  = jpni*jpnj
72      nimpp  = 1           !
73      njmpp  = 1
74      nlci   = jpi
75      nlcj   = jpj
76      nldi   = 1
77      nldj   = 1
78      nlei   = jpi
79      nlej   = jpj
80      nbondi = 2
81      nbondj = 2
82      nidom  = FLIO_DOM_NONE
83      npolj = jperio
84      l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
85      l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
86      !
87      IF(lwp) THEN
88         WRITE(numout,*)
89         WRITE(numout,*) 'mpp_init : NO massively parallel processing'
90         WRITE(numout,*) '~~~~~~~~ '
91         WRITE(numout,*) '   l_Iperio = ', l_Iperio, '    l_Jperio = ', l_Jperio 
92         WRITE(numout,*) '     npolj  = ',   npolj , '      njmpp  = ', njmpp
93      ENDIF
94      !
95      IF(  jpni /= 1 .OR. jpnj /= 1 .OR. jpnij /= 1 )                                     &
96         CALL ctl_stop( 'mpp_init: equality  jpni = jpnj = jpnij = 1 is not satisfied',   &
97            &           'the domain is lay out for distributed memory computing!' )
98         !
99   END SUBROUTINE mpp_init
100
101#else
102   !!----------------------------------------------------------------------
103   !!   'key_mpp_mpi'                     MPI massively parallel processing
104   !!----------------------------------------------------------------------
105
106
107   SUBROUTINE mpp_init
108      !!----------------------------------------------------------------------
109      !!                  ***  ROUTINE mpp_init  ***
110      !!                   
111      !! ** Purpose :   Lay out the global domain over processors.
112      !!      If land processors are to be eliminated, this program requires the
113      !!      presence of the domain configuration file. Land processors elimination
114      !!      is performed if jpni x jpnj /= jpnij. In this case, using the MPP_PREP
115      !!      preprocessing tool, help for defining the best cutting out.
116      !!
117      !! ** Method  :   Global domain is distributed in smaller local domains.
118      !!      Periodic condition is a function of the local domain position
119      !!      (global boundary or neighbouring domain) and of the global
120      !!      periodic
121      !!      Type :         jperio global periodic condition
122      !!
123      !! ** Action : - set domain parameters
124      !!                    nimpp     : longitudinal index
125      !!                    njmpp     : latitudinal  index
126      !!                    narea     : number for local area
127      !!                    nlci      : first dimension
128      !!                    nlcj      : second dimension
129      !!                    nbondi    : mark for "east-west local boundary"
130      !!                    nbondj    : mark for "north-south local boundary"
131      !!                    nproc     : number for local processor
132      !!                    noea      : number for local neighboring processor
133      !!                    nowe      : number for local neighboring processor
134      !!                    noso      : number for local neighboring processor
135      !!                    nono      : number for local neighboring processor
136      !!----------------------------------------------------------------------
137      INTEGER ::   ji, jj, jn, jproc, jarea   ! dummy loop indices
138      INTEGER ::   inum                       ! local logical unit
139      INTEGER ::   idir, ifreq, icont, isurf  ! local integers
140      INTEGER ::   ii, il1, ili, imil         !   -       -
141      INTEGER ::   ij, il2, ilj, ijm1         !   -       -
142      INTEGER ::   iino, ijno, iiso, ijso     !   -       -
143      INTEGER ::   iiea, ijea, iiwe, ijwe     !   -       -
144      INTEGER ::   iresti, irestj, iarea0     !   -       -
145      INTEGER ::   ierr                       ! local logical unit
146      REAL(wp)::   zidom, zjdom               ! local scalars
147      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   iin, ii_nono, ii_noea          ! 1D workspace
148      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ijn, ii_noso, ii_nowe          !  -     -
149      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ilci, ibondi, ipproc   ! 2D workspace
150      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ilcj, ibondj, ipolj    !  -     -
151      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilei, ildi, iono, ioea         !  -     -
152      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilej, ildj, ioso, iowe         !  -     -
153      INTEGER, DIMENSION(jpiglo,jpjglo) ::   imask   ! 2D global domain workspace
154      !!----------------------------------------------------------------------
155
156      ! If dimensions of processor grid weren't specified in the namelist file
157      ! then we calculate them here now that we have our communicator size
158      IF( jpni < 1 .OR. jpnj < 1 )   CALL mpp_init_partition( mppsize )
159      !
160#if defined key_agrif
161      IF( jpnij /= jpni*jpnj ) CALL ctl_stop( 'STOP', 'Cannot remove land proc with AGRIF' )
162#endif
163      !
164      ALLOCATE(  nfiimpp(jpni,jpnj), nfipproc(jpni,jpnj), nfilcit(jpni,jpnj) ,    &
165         &       nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) ,    &
166         &       njmppt(jpnij) , ibonjt(jpnij) , nldit(jpnij) , nldjt(jpnij) ,    &
167         &                                       nleit(jpnij) , nlejt(jpnij) ,    &
168         &       iin(jpnij), ii_nono(jpnij), ii_noea(jpnij),   &
169         &       ijn(jpnij), ii_noso(jpnij), ii_nowe(jpnij),   &
170         &       iimppt(jpni,jpnj), ilci(jpni,jpnj), ibondi(jpni,jpnj), ipproc(jpni,jpnj),   &
171         &       ijmppt(jpni,jpnj), ilcj(jpni,jpnj), ibondj(jpni,jpnj), ipolj(jpni,jpnj),   &
172         &       ilei(jpni,jpnj), ildi(jpni,jpnj), iono(jpni,jpnj), ioea(jpni,jpnj),   &
173         &       ilej(jpni,jpnj), ildj(jpni,jpnj), ioso(jpni,jpnj), iowe(jpni,jpnj),   &
174         &       STAT=ierr )
175      CALL mpp_sum( ierr )
176      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' )
177     
178      !
179#if defined key_agrif
180      IF( .NOT. Agrif_Root() ) THEN       ! AGRIF children: specific setting (cf. agrif_user.F90)
181         IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells )   &
182            CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells' )
183         IF( jpjglo /= nbcellsy + 2 + 2*nbghostcells )   &
184            CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpjglo == nbcellsy + 2 + 2*nbghostcells' )
185         IF( ln_use_jattr )   CALL ctl_stop( 'STOP', 'mpp_init:Agrif children requires ln_use_jattr = .false. ' )
186      ENDIF
187#endif
188
189#if defined key_nemocice_decomp
190      jpimax = ( nx_global+2-2*nn_hls + (jpni-1) ) / jpni + 2*nn_hls    ! first  dim.
191      jpjmax = ( ny_global+2-2*nn_hls + (jpnj-1) ) / jpnj + 2*nn_hls    ! second dim.
192#else
193      jpimax = ( jpiglo - 2*nn_hls + (jpni-1) ) / jpni + 2*nn_hls    ! first  dim.
194      jpjmax = ( jpjglo - 2*nn_hls + (jpnj-1) ) / jpnj + 2*nn_hls    ! second dim.
195#endif
196
197      !
198      IF ( jpni * jpnj == jpnij ) THEN    ! regular domain lay out over processors
199         imask(:,:) = 1               
200      ELSEIF ( jpni*jpnj > jpnij ) THEN   ! remove land-only processor (i.e. where imask(:,:)=0)
201         CALL mpp_init_mask( imask )   
202      ELSE                                ! error
203         CALL ctl_stop( 'mpp_init: jpnij > jpni x jpnj. Check namelist setting!' )
204      ENDIF
205      !
206      !  1. Dimension arrays for subdomains
207      ! -----------------------------------
208      !  Computation of local domain sizes ilci() ilcj()
209      !  These dimensions depend on global sizes jpni,jpnj and jpiglo,jpjglo
210      !  The subdomains are squares lesser than or equal to the global
211      !  dimensions divided by the number of processors minus the overlap array.
212      !
213      nreci = 2 * nn_hls
214      nrecj = 2 * nn_hls
215      iresti = 1 + MOD( jpiglo - nreci -1 , jpni )
216      irestj = 1 + MOD( jpjglo - nrecj -1 , jpnj )
217      !
218      !  Need to use jpimax and jpjmax here since jpi and jpj not yet defined
219#if defined key_nemocice_decomp
220      ! Change padding to be consistent with CICE
221      ilci(1:jpni-1      ,:) = jpimax
222      ilci(jpni          ,:) = jpiglo - (jpni - 1) * (jpimax - nreci)
223      !
224      ilcj(:,      1:jpnj-1) = jpjmax
225      ilcj(:,          jpnj) = jpjglo - (jpnj - 1) * (jpjmax - nrecj)
226#else
227      ilci(1:iresti      ,:) = jpimax
228      ilci(iresti+1:jpni ,:) = jpimax-1
229
230      ilcj(:,      1:irestj) = jpjmax
231      ilcj(:, irestj+1:jpnj) = jpjmax-1
232#endif
233      !
234      zidom = nreci + sum(ilci(:,1) - nreci )
235      zjdom = nrecj + sum(ilcj(1,:) - nrecj )
236      !
237      IF(lwp) THEN
238         WRITE(numout,*)
239         WRITE(numout,*) 'mpp_init : MPI Message Passing MPI - domain lay out over processors'
240         WRITE(numout,*) '~~~~~~~~ '
241         WRITE(numout,*) '   defines mpp subdomains'
242         WRITE(numout,*) '      iresti = ', iresti, ' jpni = ', jpni 
243         WRITE(numout,*) '      irestj = ', irestj, ' jpnj = ', jpnj
244         WRITE(numout,*)
245         WRITE(numout,*) '      sum ilci(i,1) = ', zidom, ' jpiglo = ', jpiglo
246         WRITE(numout,*) '      sum ilcj(1,j) = ', zjdom, ' jpjglo = ', jpjglo
247      ENDIF
248
249      !  2. Index arrays for subdomains
250      ! -------------------------------
251      iimppt(:,:) =  1
252      ijmppt(:,:) =  1
253      ipproc(:,:) = -1
254      !
255      IF( jpni > 1 ) THEN
256         DO jj = 1, jpnj
257            DO ji = 2, jpni
258               iimppt(ji,jj) = iimppt(ji-1,jj) + ilci(ji-1,jj) - nreci
259            END DO
260         END DO
261      ENDIF
262      nfiimpp(:,:) = iimppt(:,:)
263      !
264      IF( jpnj > 1 )THEN
265         DO jj = 2, jpnj
266            DO ji = 1, jpni
267               ijmppt(ji,jj) = ijmppt(ji,jj-1) + ilcj(ji,jj-1) - nrecj
268            END DO
269         END DO
270      ENDIF
271
272      ! 3. Subdomain description in the Regular Case
273      ! --------------------------------------------
274      ! specific cases where there is no communication -> must do the periodicity by itself
275      ! Warning: because of potential land-area suppression, do not use nbond[ij] == 2 
276      l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
277      l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
278     
279      icont = -1
280      DO jarea = 1, jpni*jpnj
281         iarea0 = jarea - 1
282         ii = 1 + MOD(iarea0,jpni)
283         ij = 1 +     iarea0/jpni
284         ili = ilci(ii,ij)
285         ilj = ilcj(ii,ij)
286         ibondi(ii,ij) = 0                         ! default: has e-w neighbours
287         IF( ii   ==    1 )   ibondi(ii,ij) = -1   ! first column, has only e neighbour
288         IF( ii   == jpni )   ibondi(ii,ij) =  1   ! last column,  has only w neighbour
289         IF( jpni ==    1 )   ibondi(ii,ij) =  2   ! has no e-w neighbour
290         ibondj(ii,ij) = 0                         ! default: has n-s neighbours
291         IF( ij   ==    1 )   ibondj(ii,ij) = -1   ! first row, has only n neighbour
292         IF( ij   == jpnj )   ibondj(ii,ij) =  1   ! last row,  has only s neighbour
293         IF( jpnj ==    1 )   ibondj(ii,ij) =  2   ! has no n-s neighbour
294
295         ! Subdomain neighbors (get their zone number): default definition
296         ioso(ii,ij) = iarea0 - jpni
297         iowe(ii,ij) = iarea0 - 1
298         ioea(ii,ij) = iarea0 + 1
299         iono(ii,ij) = iarea0 + jpni
300         ildi(ii,ij) =  1  + nn_hls
301         ilei(ii,ij) = ili - nn_hls
302         ildj(ii,ij) =  1  + nn_hls
303         ilej(ii,ij) = ilj - nn_hls
304
305         ! East-West periodicity: change ibondi, ioea, iowe
306         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7 ) THEN
307            IF( jpni  /= 1 )   ibondi(ii,ij) = 0                        ! redefine: all have e-w neighbours
308            IF( ii ==    1 )   iowe(ii,ij) = iarea0 +        (jpni-1)   ! redefine: first column, address of w neighbour
309            IF( ii == jpni )   ioea(ii,ij) = iarea0 -        (jpni-1)   ! redefine: last column,  address of e neighbour
310         ENDIF
311
312         ! Simple North-South periodicity: change ibondj, ioso, iono
313         IF( jperio == 2 .OR. jperio == 7 ) THEN
314            IF( jpnj  /= 1 )   ibondj(ii,ij) = 0                        ! redefine: all have n-s neighbours
315            IF( ij ==    1 )   ioso(ii,ij) = iarea0 + jpni * (jpnj-1)   ! redefine: first row, address of s neighbour
316            IF( ij == jpnj )   iono(ii,ij) = iarea0 - jpni * (jpnj-1)   ! redefine: last row,  address of n neighbour
317         ENDIF
318
319         ! North fold: define ipolj, change iono. Warning: we do not change ibondj...
320         ipolj(ii,ij) = 0
321         IF( jperio == 3 .OR. jperio == 4 ) THEN
322            ijm1 = jpni*(jpnj-1)
323            imil = ijm1+(jpni+1)/2
324            IF( jarea > ijm1 ) ipolj(ii,ij) = 3
325            IF( MOD(jpni,2) == 1 .AND. jarea == imil ) ipolj(ii,ij) = 4
326            IF( ipolj(ii,ij) == 3 ) iono(ii,ij) = jpni*jpnj-jarea+ijm1   ! MPI rank of northern neighbour
327         ENDIF
328         IF( jperio == 5 .OR. jperio == 6 ) THEN
329            ijm1 = jpni*(jpnj-1)
330            imil = ijm1+(jpni+1)/2
331            IF( jarea > ijm1) ipolj(ii,ij) = 5
332            IF( MOD(jpni,2) == 1 .AND. jarea == imil ) ipolj(ii,ij) = 6
333            IF( ipolj(ii,ij) == 5) iono(ii,ij) = jpni*jpnj-jarea+ijm1    ! MPI rank of northern neighbour
334         ENDIF
335         !
336         ! Check wet points over the entire domain to preserve the MPI communication stencil
337         isurf = 0
338         DO jj = 1, ilj
339            DO  ji = 1, ili
340               IF( imask(ji+iimppt(ii,ij)-1, jj+ijmppt(ii,ij)-1) == 1)   isurf = isurf+1
341            END DO
342         END DO
343         !
344         IF( isurf /= 0 ) THEN
345            icont = icont + 1
346            ipproc(ii,ij) = icont
347            iin(icont+1) = ii
348            ijn(icont+1) = ij
349         ENDIF
350      END DO
351      !
352      nfipproc(:,:) = ipproc(:,:)
353
354      ! Check potential error
355      IF( icont+1 /= jpnij ) THEN
356         WRITE(ctmp1,*) ' jpni =',jpni,' jpnj =',jpnj
357         WRITE(ctmp2,*) ' jpnij =',jpnij, '< jpni x jpnj' 
358         WRITE(ctmp3,*) ' ***********, mpp_init2 finds jpnij=',icont+1
359         CALL ctl_stop( 'STOP', 'mpp_init: Eliminate land processors algorithm', '', ctmp1, ctmp2, '', ctmp3 )
360      ENDIF
361
362      ! 4. Subdomain print
363      ! ------------------
364      IF(lwp) THEN
365         ifreq = 4
366         il1 = 1
367         DO jn = 1, (jpni-1)/ifreq+1
368            il2 = MIN(jpni,il1+ifreq-1)
369            WRITE(numout,*)
370            WRITE(numout,9400) ('***',ji=il1,il2-1)
371            DO jj = jpnj, 1, -1
372               WRITE(numout,9403) ('   ',ji=il1,il2-1)
373               WRITE(numout,9402) jj, (ilci(ji,jj),ilcj(ji,jj),ji=il1,il2)
374               WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2)
375               WRITE(numout,9403) ('   ',ji=il1,il2-1)
376               WRITE(numout,9400) ('***',ji=il1,il2-1)
377            END DO
378            WRITE(numout,9401) (ji,ji=il1,il2)
379            il1 = il1+ifreq
380         END DO
381 9400    FORMAT('           ***'   ,20('*************',a3)    )
382 9403    FORMAT('           *     ',20('         *   ',a3)    )
383 9401    FORMAT('              '   ,20('   ',i3,'          ') )
384 9402    FORMAT('       ',i3,' *  ',20(i3,'  x',i3,'   *   ') )
385 9404    FORMAT('           *  '   ,20('      ',i3,'   *   ') )
386      ENDIF
387
388      ! 5. neighbour treatment: change ibondi, ibondj if next to a land zone
389      ! ----------------------
390      DO jarea = 1, jpni*jpnj
391         ii = 1 + MOD( jarea-1  , jpni )
392         ij = 1 +     (jarea-1) / jpni
393         ! land-only area with an active n neigbour
394         IF ( ipproc(ii,ij) == -1 .AND. 0 <= iono(ii,ij) .AND. iono(ii,ij) <= jpni*jpnj-1 ) THEN
395            iino = 1 + MOD( iono(ii,ij) , jpni )                    ! ii index of this n neigbour
396            ijno = 1 +      iono(ii,ij) / jpni                      ! ij index of this n neigbour
397            ! In case of north fold exchange: I am the n neigbour of my n neigbour!! (#1057)
398            ! --> for northern neighbours of northern row processors (in case of north-fold)
399            !     need to reverse the LOGICAL direction of communication
400            idir = 1                                           ! we are indeed the s neigbour of this n neigbour
401            IF( ij == jpnj .AND. ijno == jpnj )   idir = -1    ! both are on the last row, we are in fact the n neigbour
402            IF( ibondj(iino,ijno) == idir     )   ibondj(iino,ijno) =   2     ! this n neigbour had only a s/n neigbour -> no more
403            IF( ibondj(iino,ijno) == 0        )   ibondj(iino,ijno) = -idir   ! this n neigbour had both, n-s neighbours -> keep 1
404         ENDIF
405         ! land-only area with an active s neigbour
406         IF( ipproc(ii,ij) == -1 .AND. 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= jpni*jpnj-1 ) THEN
407            iiso = 1 + MOD( ioso(ii,ij) , jpni )                    ! ii index of this s neigbour
408            ijso = 1 +      ioso(ii,ij) / jpni                      ! ij index of this s neigbour
409            IF( ibondj(iiso,ijso) == -1 )   ibondj(iiso,ijso) = 2   ! this s neigbour had only a n neigbour    -> no more neigbour
410            IF( ibondj(iiso,ijso) ==  0 )   ibondj(iiso,ijso) = 1   ! this s neigbour had both, n-s neighbours -> keep s neigbour
411         ENDIF
412         ! land-only area with an active e neigbour
413         IF( ipproc(ii,ij) == -1 .AND. 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= jpni*jpnj-1 ) THEN
414            iiea = 1 + MOD( ioea(ii,ij) , jpni )                    ! ii index of this e neigbour
415            ijea = 1 +      ioea(ii,ij) / jpni                      ! ij index of this e neigbour
416            IF( ibondi(iiea,ijea) == 1 )   ibondi(iiea,ijea) =  2   ! this e neigbour had only a w neigbour    -> no more neigbour
417            IF( ibondi(iiea,ijea) == 0 )   ibondi(iiea,ijea) = -1   ! this e neigbour had both, e-w neighbours -> keep e neigbour
418         ENDIF
419         ! land-only area with an active w neigbour
420         IF( ipproc(ii,ij) == -1 .AND. 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= jpni*jpnj-1) THEN
421            iiwe = 1 + MOD( iowe(ii,ij) , jpni )                    ! ii index of this w neigbour
422            ijwe = 1 +      iowe(ii,ij) / jpni                      ! ij index of this w neigbour
423            IF( ibondi(iiwe,ijwe) == -1 )   ibondi(iiwe,ijwe) = 2   ! this w neigbour had only a e neigbour    -> no more neigbour
424            IF( ibondi(iiwe,ijwe) ==  0 )   ibondi(iiwe,ijwe) = 1   ! this w neigbour had both, e-w neighbours -> keep w neigbour
425         ENDIF
426      END DO
427
428      ! Update il[de][ij] according to modified ibond[ij]
429      ! ----------------------
430      DO jproc = 1, jpnij
431         ii = iin(jproc)
432         ij = ijn(jproc)
433         IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) =  1
434         IF( ibondi(ii,ij) ==  1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij)
435         IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) =  1
436         IF( ibondj(ii,ij) ==  1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij)
437      END DO
438         
439      ! just to save nono etc for all proc
440      ! warning ii*ij (zone) /= nproc (processors)!
441      ! ioso = zone number, ii_noso = proc number
442      ii_noso(:) = -1
443      ii_nono(:) = -1
444      ii_noea(:) = -1
445      ii_nowe(:) = -1 
446      DO jproc = 1, jpnij
447         ii = iin(jproc)
448         ij = ijn(jproc)
449         IF( 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= (jpni*jpnj-1) ) THEN
450            iiso = 1 + MOD( ioso(ii,ij) , jpni )
451            ijso = 1 +      ioso(ii,ij) / jpni
452            ii_noso(jproc) = ipproc(iiso,ijso)
453         ENDIF
454         IF( 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= (jpni*jpnj-1) ) THEN
455          iiwe = 1 + MOD( iowe(ii,ij) , jpni )
456          ijwe = 1 +      iowe(ii,ij) / jpni
457          ii_nowe(jproc) = ipproc(iiwe,ijwe)
458         ENDIF
459         IF( 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= (jpni*jpnj-1) ) THEN
460            iiea = 1 + MOD( ioea(ii,ij) , jpni )
461            ijea = 1 +      ioea(ii,ij) / jpni
462            ii_noea(jproc)= ipproc(iiea,ijea)
463         ENDIF
464         IF( 0 <= iono(ii,ij) .AND. iono(ii,ij) <= (jpni*jpnj-1) ) THEN
465            iino = 1 + MOD( iono(ii,ij) , jpni )
466            ijno = 1 +      iono(ii,ij) / jpni
467            ii_nono(jproc)= ipproc(iino,ijno)
468         ENDIF
469      END DO
470   
471      ! 6. Change processor name
472      ! ------------------------
473      ii = iin(narea)
474      ij = ijn(narea)
475      !
476      ! set default neighbours
477      noso = ii_noso(narea)
478      nowe = ii_nowe(narea)
479      noea = ii_noea(narea)
480      nono = ii_nono(narea)
481      nlci = ilci(ii,ij) 
482      nldi = ildi(ii,ij)
483      nlei = ilei(ii,ij)
484      nlcj = ilcj(ii,ij) 
485      nldj = ildj(ii,ij)
486      nlej = ilej(ii,ij)
487      nbondi = ibondi(ii,ij)
488      nbondj = ibondj(ii,ij)
489      nimpp = iimppt(ii,ij) 
490      njmpp = ijmppt(ii,ij)
491      jpi = nlci
492      jpj = nlcj
493      jpk = jpkglo                                             ! third dim
494#if defined key_agrif
495      ! simple trick to use same vertical grid as parent but different number of levels:
496      ! Save maximum number of levels in jpkglo, then define all vertical grids with this number.
497      ! Suppress once vertical online interpolation is ok
498!!$      IF(.NOT.Agrif_Root())   jpkglo = Agrif_Parent( jpkglo )
499#endif
500      jpim1 = jpi-1                                            ! inner domain indices
501      jpjm1 = jpj-1                                            !   "           "
502      jpkm1 = MAX( 1, jpk-1 )                                  !   "           "
503      jpij  = jpi*jpj                                          !  jpi x j
504      DO jproc = 1, jpnij
505         ii = iin(jproc)
506         ij = ijn(jproc)
507         nlcit(jproc) = ilci(ii,ij)
508         nldit(jproc) = ildi(ii,ij)
509         nleit(jproc) = ilei(ii,ij)
510         nlcjt(jproc) = ilcj(ii,ij)
511         nldjt(jproc) = ildj(ii,ij)
512         nlejt(jproc) = ilej(ii,ij)
513         ibonit(jproc) = ibondi(ii,ij)
514         ibonjt(jproc) = ibondj(ii,ij)
515         nimppt(jproc) = iimppt(ii,ij) 
516         njmppt(jproc) = ijmppt(ii,ij) 
517      END DO
518      nfilcit(:,:) = ilci(:,:)
519
520      ! Save processor layout in ascii file
521      IF (lwp) THEN
522         CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
523         WRITE(inum,'(a)') '   jpnij   jpimax  jpjmax    jpk  jpiglo  jpjglo'//&
524   &           ' ( local:    narea     jpi     jpj )'
525         WRITE(inum,'(6i8,a,3i8,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,&
526   &           ' ( local: ',narea,jpi,jpj,' )'
527         WRITE(inum,'(a)') 'nproc nlci nlcj nldi nldj nlei nlej nimp njmp nono noso nowe noea nbondi nbondj '
528
529         DO jproc = 1, jpnij
530            WRITE(inum,'(13i5,2i7)')   jproc-1, nlcit  (jproc), nlcjt  (jproc),   &
531               &                                nldit  (jproc), nldjt  (jproc),   &
532               &                                nleit  (jproc), nlejt  (jproc),   &
533               &                                nimppt (jproc), njmppt (jproc),   & 
534               &                                ii_nono(jproc), ii_noso(jproc),   &
535               &                                ii_nowe(jproc), ii_noea(jproc),   &
536               &                                ibonit (jproc), ibonjt (jproc) 
537         END DO
538         CLOSE(inum)   
539      END IF
540
541      !                          ! north fold parameter
542      ! Defined npolj, either 0, 3 , 4 , 5 , 6
543      ! In this case the important thing is that npolj /= 0
544      ! Because if we go through these line it is because jpni >1 and thus
545      ! we must use lbcnorthmpp, which tests only npolj =0 or npolj /= 0
546      npolj = 0
547      ij = ijn(narea)
548      IF( jperio == 3 .OR. jperio == 4 ) THEN
549         IF( ij == jpnj )   npolj = 3
550      ENDIF
551      IF( jperio == 5 .OR. jperio == 6 ) THEN
552         IF( ij == jpnj )   npolj = 5
553      ENDIF
554      !
555      nproc = narea-1
556      IF(lwp) THEN
557         WRITE(numout,*)
558         WRITE(numout,*) '   resulting internal parameters : '
559         WRITE(numout,*) '      nproc  = ', nproc
560         WRITE(numout,*) '      nowe   = ', nowe  , '   noea  =  ', noea
561         WRITE(numout,*) '      nono   = ', nono  , '   noso  =  ', noso
562         WRITE(numout,*) '      nbondi = ', nbondi
563         WRITE(numout,*) '      nbondj = ', nbondj
564         WRITE(numout,*) '      npolj  = ', npolj
565         WRITE(numout,*) '    l_Iperio = ', l_Iperio
566         WRITE(numout,*) '    l_Jperio = ', l_Jperio
567         WRITE(numout,*) '      nlci   = ', nlci
568         WRITE(numout,*) '      nlcj   = ', nlcj
569         WRITE(numout,*) '      nimpp  = ', nimpp
570         WRITE(numout,*) '      njmpp  = ', njmpp
571         WRITE(numout,*) '      nreci  = ', nreci 
572         WRITE(numout,*) '      nrecj  = ', nrecj 
573         WRITE(numout,*) '      nn_hls = ', nn_hls 
574      ENDIF
575
576      !                          ! Prepare mpp north fold
577      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
578         CALL mpp_ini_north
579         IF(lwp) WRITE(numout,*)
580         IF(lwp) WRITE(numout,*) '   ==>>>   North fold boundary prepared for jpni >1'
581      ENDIF
582      !
583      CALL mpp_init_ioipsl       ! Prepare NetCDF output file (if necessary)
584      !
585      IF( ln_nnogather )   CALL mpp_init_nfdcom     ! northfold neighbour lists
586      !
587      DEALLOCATE(iin, ijn, ii_nono, ii_noea, ii_noso, ii_nowe,    &
588         &       iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj,   &
589         &       ilci, ilcj, ilei, ilej, ildi, ildj,              &
590         &       iono, ioea, ioso, iowe)
591      !
592    END SUBROUTINE mpp_init
593
594
595    SUBROUTINE mpp_init_mask( kmask )
596      !!----------------------------------------------------------------------
597      !!                  ***  ROUTINE mpp_init_mask  ***
598      !!
599      !! ** Purpose : Read relevant bathymetric information in a global array
600      !!              in order to provide a land/sea mask used for the elimination
601      !!              of land domains, in an mpp computation.
602      !!
603      !! ** Method  : Read the namelist ln_zco and ln_isfcav in namelist namzgr
604      !!              in order to choose the correct bathymetric information
605      !!              (file and variables) 
606      !!----------------------------------------------------------------------
607      INTEGER, DIMENSION(jpiglo,jpjglo), INTENT(out) ::   kmask   ! global domain
608 
609      INTEGER :: inum   !: logical unit for configuration file
610      INTEGER :: ios    !: iostat error flag
611      INTEGER ::  ijstartrow                   ! temporary integers
612      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zbot, zbdy          ! global workspace
613      REAL(wp) ::   zidom , zjdom          ! local scalars
614      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           &
615           &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
616           &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             & 
617           &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
618           &             cn_ice, nn_ice_dta,                                     &
619           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
620           &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
621      !!----------------------------------------------------------------------
622      ! 0. initialisation
623      ! -----------------
624      CALL iom_open( cn_domcfg, inum )
625      !
626      ! ocean bottom level
627      CALL iom_get( inum, jpdom_unknown, 'bottom_level' , zbot , lrowattr=ln_use_jattr )  ! nb of ocean T-points
628      !
629      CALL iom_close( inum )
630      !
631      ! 2D ocean mask (=1 if at least one level of the water column is ocean, =0 otherwise)
632      WHERE( zbot(:,:) > 0 )   ;   kmask(:,:) = 1
633      ELSEWHERE                ;   kmask(:,:) = 0
634      END WHERE
635 
636      ! Adjust kmask with bdy_msk if it exists
637 
638      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY
639      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
640903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp )
641      !
642      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY
643      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
644904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp )
645
646      IF( ln_bdy .AND. ln_mask_file ) THEN
647         CALL iom_open( cn_mask_file, inum )
648         CALL iom_get ( inum, jpdom_unknown, 'bdy_msk', zbdy )
649         CALL iom_close( inum )
650         WHERE ( zbdy(:,:) <= 0. ) kmask = 0
651      ENDIF
652      !
653   END SUBROUTINE mpp_init_mask
654
655
656   SUBROUTINE mpp_init_ioipsl
657      !!----------------------------------------------------------------------
658      !!                  ***  ROUTINE mpp_init_ioipsl  ***
659      !!
660      !! ** Purpose :   
661      !!
662      !! ** Method  :   
663      !!
664      !! History :
665      !!   9.0  !  04-03  (G. Madec )  MPP-IOIPSL
666      !!   " "  !  08-12  (A. Coward)  addition in case of jpni*jpnj < jpnij
667      !!----------------------------------------------------------------------
668      INTEGER, DIMENSION(2) ::   iglo, iloc, iabsf, iabsl, ihals, ihale, idid
669      !!----------------------------------------------------------------------
670
671      ! The domain is split only horizontally along i- or/and j- direction
672      ! So we need at the most only 1D arrays with 2 elements.
673      ! Set idompar values equivalent to the jpdom_local_noextra definition
674      ! used in IOM. This works even if jpnij .ne. jpni*jpnj.
675      iglo(1) = jpiglo
676      iglo(2) = jpjglo
677      iloc(1) = nlci
678      iloc(2) = nlcj
679      iabsf(1) = nimppt(narea)
680      iabsf(2) = njmppt(narea)
681      iabsl(:) = iabsf(:) + iloc(:) - 1
682      ihals(1) = nldi - 1
683      ihals(2) = nldj - 1
684      ihale(1) = nlci - nlei
685      ihale(2) = nlcj - nlej
686      idid(1) = 1
687      idid(2) = 2
688
689      IF(lwp) THEN
690          WRITE(numout,*)
691          WRITE(numout,*) 'mpp_init_ioipsl :   iloc  = ', iloc (1), iloc (2)
692          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf(1), iabsf(2)
693          WRITE(numout,*) '                    ihals = ', ihals(1), ihals(2)
694          WRITE(numout,*) '                    ihale = ', ihale(1), ihale(2)
695      ENDIF
696      !
697      CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
698      !
699   END SUBROUTINE mpp_init_ioipsl 
700
701
702   SUBROUTINE mpp_init_partition( num_pes )
703      !!----------------------------------------------------------------------
704      !!                 ***  ROUTINE mpp_init_partition  ***
705      !!
706      !! ** Purpose :
707      !!
708      !! ** Method  :
709      !!----------------------------------------------------------------------
710      INTEGER, INTENT(in) ::   num_pes   ! The number of MPI processes we have
711      !
712      INTEGER, PARAMETER :: nfactmax = 20
713      INTEGER :: nfact ! The no. of factors returned
714      INTEGER :: ierr  ! Error flag
715      INTEGER :: ji
716      INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value
717      INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors
718      !!----------------------------------------------------------------------
719      !
720      ierr = 0
721      !
722      CALL factorise( ifact, nfactmax, nfact, num_pes, ierr )
723      !
724      IF( nfact <= 1 ) THEN
725         WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
726         WRITE (numout, *) '       : using grid of ',num_pes,' x 1'
727         jpnj = 1
728         jpni = num_pes
729      ELSE
730         ! Search through factors for the pair that are closest in value
731         mindiff = 1000000
732         imin    = 1
733         DO ji = 1, nfact-1, 2
734            idiff = ABS( ifact(ji) - ifact(ji+1) )
735            IF( idiff < mindiff ) THEN
736               mindiff = idiff
737               imin = ji
738            ENDIF
739         END DO
740         jpnj = ifact(imin)
741         jpni = ifact(imin + 1)
742      ENDIF
743      !
744      jpnij = jpni*jpnj
745      !
746   END SUBROUTINE mpp_init_partition
747
748
749   SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr )
750      !!----------------------------------------------------------------------
751      !!                     ***  ROUTINE factorise  ***
752      !!
753      !! ** Purpose :   return the prime factors of n.
754      !!                knfax factors are returned in array kfax which is of
755      !!                maximum dimension kmaxfax.
756      !! ** Method  :
757      !!----------------------------------------------------------------------
758      INTEGER                    , INTENT(in   ) ::   kn, kmaxfax
759      INTEGER                    , INTENT(  out) ::   kerr, knfax
760      INTEGER, DIMENSION(kmaxfax), INTENT(  out) ::   kfax
761      !
762      INTEGER :: ifac, jl, inu
763      INTEGER, PARAMETER :: ntest = 14
764      INTEGER, DIMENSION(ntest) ::   ilfax
765      !!----------------------------------------------------------------------
766      !
767      ! lfax contains the set of allowed factors.
768      ilfax(:) = (/(2**jl,jl=ntest,1,-1)/)
769      !
770      ! Clear the error flag and initialise output vars
771      kerr  = 0
772      kfax  = 1
773      knfax = 0
774      !
775      IF( kn /= 1 ) THEN      ! Find the factors of n
776         !
777         ! nu holds the unfactorised part of the number.
778         ! knfax holds the number of factors found.
779         ! l points to the allowed factor list.
780         ! ifac holds the current factor.
781         !
782         inu   = kn
783         knfax = 0
784         !
785         DO jl = ntest, 1, -1
786            !
787            ifac = ilfax(jl)
788            IF( ifac > inu )   CYCLE
789            !
790            ! Test whether the factor will divide.
791            !
792            IF( MOD(inu,ifac) == 0 ) THEN
793               !
794               knfax = knfax + 1            ! Add the factor to the list
795               IF( knfax > kmaxfax ) THEN
796                  kerr = 6
797                  write (*,*) 'FACTOR: insufficient space in factor array ', knfax
798                  return
799               ENDIF
800               kfax(knfax) = ifac
801               ! Store the other factor that goes with this one
802               knfax = knfax + 1
803               kfax(knfax) = inu / ifac
804               !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax)
805            ENDIF
806            !
807         END DO
808         !
809      ENDIF
810      !
811   END SUBROUTINE factorise
812
813
814   SUBROUTINE mpp_init_nfdcom
815      !!----------------------------------------------------------------------
816      !!                     ***  ROUTINE  mpp_init_nfdcom  ***
817      !! ** Purpose :   Setup for north fold exchanges with explicit
818      !!                point-to-point messaging
819      !!
820      !! ** Method :   Initialization of the northern neighbours lists.
821      !!----------------------------------------------------------------------
822      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
823      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
824      !!----------------------------------------------------------------------
825      INTEGER  ::   sxM, dxM, sxT, dxT, jn
826      INTEGER  ::   njmppmax
827      !!----------------------------------------------------------------------
828      !
829      njmppmax = MAXVAL( njmppt )
830      !
831      !initializes the north-fold communication variables
832      isendto(:) = 0
833      nsndto     = 0
834      !
835      IF ( njmpp == njmppmax ) THEN      ! if I am a process in the north
836         !
837         !sxM is the first point (in the global domain) needed to compute the north-fold for the current process
838         sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
839         !dxM is the last point (in the global domain) needed to compute the north-fold for the current process
840         dxM = jpiglo - nimppt(narea) + 2
841         !
842         ! loop over the other north-fold processes to find the processes
843         ! managing the points belonging to the sxT-dxT range
844         !
845         DO jn = 1, jpni
846            !
847            sxT = nfiimpp(jn, jpnj)                            ! sxT = 1st  point (in the global domain) of the jn process
848            dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1    ! dxT = last point (in the global domain) of the jn process
849            !
850            IF    ( sxT < sxM  .AND.  sxM < dxT ) THEN
851               nsndto          = nsndto + 1
852               isendto(nsndto) = jn
853            ELSEIF( sxM <= sxT  .AND.  dxM >= dxT ) THEN
854               nsndto          = nsndto + 1
855               isendto(nsndto) = jn
856            ELSEIF( dxM <  dxT  .AND.  sxT <  dxM ) THEN
857               nsndto          = nsndto + 1
858               isendto(nsndto) = jn
859            ENDIF
860            !
861         END DO
862         nfsloop = 1
863         nfeloop = nlci
864         DO jn = 2,jpni-1
865            IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
866               IF( nfipproc(jn-1,jpnj) == -1 )   nfsloop = nldi
867               IF( nfipproc(jn+1,jpnj) == -1 )   nfeloop = nlei
868            ENDIF
869         END DO
870         !
871      ENDIF
872      l_north_nogather = .TRUE.
873      !
874   END SUBROUTINE mpp_init_nfdcom
875
876   
877#endif
878
879   !!======================================================================
880END MODULE mppini
Note: See TracBrowser for help on using the repository browser.