source: NEMO/branches/UKMO/2018_NEMO4_beta_mirror_10037_Benchmark/src/OCE/LBC/mppini.F90 @ 10284

Last change on this file since 10284 was 10284, checked in by andmirek, 3 years ago

GMED 425 dev_r10037_mppmask

File size: 40.4 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    = jpkglo
65      jpim1  = jpi-1                                            ! inner domain indices
66      jpjm1  = jpj-1                                            !   "           "
67      jpkm1  = MAX( 1, jpk-1 )                                  !   "           "
68      jpij   = jpi*jpj
69      jpni   = 1
70      jpnj   = 1
71      jpnij  = jpni*jpnj
72      nimpp  = 1           !
73      njmpp  = 1
74      nlci   = jpi
75      nlcj   = jpj
76      nldi   = 1
77      nldj   = 1
78      nlei   = jpi
79      nlej   = jpj
80      nbondi = 2
81      nbondj = 2
82      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      END IF
539
540      !                          ! north fold parameter
541      ! Defined npolj, either 0, 3 , 4 , 5 , 6
542      ! In this case the important thing is that npolj /= 0
543      ! Because if we go through these line it is because jpni >1 and thus
544      ! we must use lbcnorthmpp, which tests only npolj =0 or npolj /= 0
545      npolj = 0
546      ij = ijn(narea)
547      IF( jperio == 3 .OR. jperio == 4 ) THEN
548         IF( ij == jpnj )   npolj = 3
549      ENDIF
550      IF( jperio == 5 .OR. jperio == 6 ) THEN
551         IF( ij == jpnj )   npolj = 5
552      ENDIF
553      !
554      nproc = narea-1
555      IF(lwp) THEN
556         WRITE(numout,*)
557         WRITE(numout,*) '   resulting internal parameters : '
558         WRITE(numout,*) '      nproc  = ', nproc
559         WRITE(numout,*) '      nowe   = ', nowe  , '   noea  =  ', noea
560         WRITE(numout,*) '      nono   = ', nono  , '   noso  =  ', noso
561         WRITE(numout,*) '      nbondi = ', nbondi
562         WRITE(numout,*) '      nbondj = ', nbondj
563         WRITE(numout,*) '      npolj  = ', npolj
564         WRITE(numout,*) '    l_Iperio = ', l_Iperio
565         WRITE(numout,*) '    l_Jperio = ', l_Jperio
566         WRITE(numout,*) '      nlci   = ', nlci
567         WRITE(numout,*) '      nlcj   = ', nlcj
568         WRITE(numout,*) '      nimpp  = ', nimpp
569         WRITE(numout,*) '      njmpp  = ', njmpp
570         WRITE(numout,*) '      nreci  = ', nreci 
571         WRITE(numout,*) '      nrecj  = ', nrecj 
572         WRITE(numout,*) '      nn_hls = ', nn_hls 
573      ENDIF
574
575      !                          ! Prepare mpp north fold
576      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
577         CALL mpp_ini_north
578         IF (lwp) THEN
579            WRITE(numout,*)
580            WRITE(numout,*) '   ==>>>   North fold boundary prepared for jpni >1'
581            ! additional prints in layout.dat
582            WRITE(inum,*)
583            WRITE(inum,*)
584            WRITE(inum,*) 'number of subdomains located along the north fold : ', ndim_rank_north
585            WRITE(inum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
586            DO jproc = 1, ndim_rank_north, 5
587               WRITE(inum,*) nrank_north( jproc:MINVAL( (/jproc+4,ndim_rank_north/) ) )
588            END DO
589         ENDIF
590      ENDIF
591      !
592      CALL mpp_init_ioipsl       ! Prepare NetCDF output file (if necessary)
593      !
594      IF( ln_nnogather ) THEN
595         CALL mpp_init_nfdcom     ! northfold neighbour lists
596         IF (lwp) THEN
597            WRITE(inum,*)
598            WRITE(inum,*)
599            WRITE(inum,*) 'north fold exchanges with explicit point-to-point messaging :'
600            WRITE(inum,*) 'nfsloop : ', nfsloop
601            WRITE(inum,*) 'nfeloop : ', nfeloop
602            WRITE(inum,*) 'nsndto : ', nsndto
603            WRITE(inum,*) 'isendto : ', isendto
604         ENDIF
605      ENDIF
606      !
607      IF (lwp) CLOSE(inum)   
608      !
609      DEALLOCATE(iin, ijn, ii_nono, ii_noea, ii_noso, ii_nowe,    &
610         &       iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj,   &
611         &       ilci, ilcj, ilei, ilej, ildi, ildj,              &
612         &       iono, ioea, ioso, iowe)
613      !
614    END SUBROUTINE mpp_init
615
616
617    SUBROUTINE mpp_init_mask( kmask )
618      !!----------------------------------------------------------------------
619      !!                  ***  ROUTINE mpp_init_mask  ***
620      !!
621      !! ** Purpose : Read relevant bathymetric information in a global array
622      !!              in order to provide a land/sea mask used for the elimination
623      !!              of land domains, in an mpp computation.
624      !!
625      !! ** Method  : Read the namelist ln_zco and ln_isfcav in namelist namzgr
626      !!              in order to choose the correct bathymetric information
627      !!              (file and variables) 
628      !!----------------------------------------------------------------------
629      INTEGER, DIMENSION(jpiglo,jpjglo), INTENT(out) ::   kmask   ! global domain
630 
631      INTEGER :: inum   !: logical unit for configuration file
632      INTEGER :: ios    !: iostat error flag
633      INTEGER ::  ijstartrow                   ! temporary integers
634      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zmsk, zbdy          ! global workspace
635      REAL(wp) ::   zidom , zjdom          ! local scalars
636      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           &
637           &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
638           &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             & 
639           &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
640           &             cn_ice, nn_ice_dta,                                     &
641           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
642           &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
643      !!----------------------------------------------------------------------
644      ! 0. initialisation
645      ! -----------------
646      CALL iom_open( cn_domcfg, inum )
647      !
648      ! mpp mask (used to define land processor)
649      CALL iom_get( inum, jpdom_unknown, 'mppmask' , zmsk , lrowattr=ln_use_jattr )  ! nb of ocean T-points
650      !
651      CALL iom_close( inum )
652      !
653      ! 2D ocean mask (=1 if at least one level of the water column is ocean, =0 otherwise)
654      WHERE( zmsk(:,:) > 0 )   ;   kmask(:,:) = 1
655      ELSEWHERE                ;   kmask(:,:) = 0
656      END WHERE
657 
658      ! Adjust kmask with bdy_msk if it exists
659 
660      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY
661      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
662903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp )
663      !
664      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY
665      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
666904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp )
667
668      IF( ln_bdy .AND. ln_mask_file ) THEN
669         CALL iom_open( cn_mask_file, inum )
670         CALL iom_get ( inum, jpdom_unknown, 'bdy_msk', zbdy )
671         CALL iom_close( inum )
672         WHERE ( zbdy(:,:) <= 0. ) kmask = 0
673      ENDIF
674      !
675   END SUBROUTINE mpp_init_mask
676
677
678   SUBROUTINE mpp_init_ioipsl
679      !!----------------------------------------------------------------------
680      !!                  ***  ROUTINE mpp_init_ioipsl  ***
681      !!
682      !! ** Purpose :   
683      !!
684      !! ** Method  :   
685      !!
686      !! History :
687      !!   9.0  !  04-03  (G. Madec )  MPP-IOIPSL
688      !!   " "  !  08-12  (A. Coward)  addition in case of jpni*jpnj < jpnij
689      !!----------------------------------------------------------------------
690      INTEGER, DIMENSION(2) ::   iglo, iloc, iabsf, iabsl, ihals, ihale, idid
691      !!----------------------------------------------------------------------
692
693      ! The domain is split only horizontally along i- or/and j- direction
694      ! So we need at the most only 1D arrays with 2 elements.
695      ! Set idompar values equivalent to the jpdom_local_noextra definition
696      ! used in IOM. This works even if jpnij .ne. jpni*jpnj.
697      iglo(1) = jpiglo
698      iglo(2) = jpjglo
699      iloc(1) = nlci
700      iloc(2) = nlcj
701      iabsf(1) = nimppt(narea)
702      iabsf(2) = njmppt(narea)
703      iabsl(:) = iabsf(:) + iloc(:) - 1
704      ihals(1) = nldi - 1
705      ihals(2) = nldj - 1
706      ihale(1) = nlci - nlei
707      ihale(2) = nlcj - nlej
708      idid(1) = 1
709      idid(2) = 2
710
711      IF(lwp) THEN
712          WRITE(numout,*)
713          WRITE(numout,*) 'mpp_init_ioipsl :   iloc  = ', iloc (1), iloc (2)
714          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf(1), iabsf(2)
715          WRITE(numout,*) '                    ihals = ', ihals(1), ihals(2)
716          WRITE(numout,*) '                    ihale = ', ihale(1), ihale(2)
717      ENDIF
718      !
719      CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
720      !
721   END SUBROUTINE mpp_init_ioipsl 
722
723
724   SUBROUTINE mpp_init_partition( num_pes )
725      !!----------------------------------------------------------------------
726      !!                 ***  ROUTINE mpp_init_partition  ***
727      !!
728      !! ** Purpose :
729      !!
730      !! ** Method  :
731      !!----------------------------------------------------------------------
732      INTEGER, INTENT(in) ::   num_pes   ! The number of MPI processes we have
733      !
734      INTEGER, PARAMETER :: nfactmax = 20
735      INTEGER :: nfact ! The no. of factors returned
736      INTEGER :: ierr  ! Error flag
737      INTEGER :: ji
738      INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value
739      INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors
740      !!----------------------------------------------------------------------
741      !
742      ierr = 0
743      !
744      CALL factorise( ifact, nfactmax, nfact, num_pes, ierr )
745      !
746      IF( nfact <= 1 ) THEN
747         WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
748         WRITE (numout, *) '       : using grid of ',num_pes,' x 1'
749         jpnj = 1
750         jpni = num_pes
751      ELSE
752         ! Search through factors for the pair that are closest in value
753         mindiff = 1000000
754         imin    = 1
755         DO ji = 1, nfact-1, 2
756            idiff = ABS( ifact(ji) - ifact(ji+1) )
757            IF( idiff < mindiff ) THEN
758               mindiff = idiff
759               imin = ji
760            ENDIF
761         END DO
762         jpnj = ifact(imin)
763         jpni = ifact(imin + 1)
764      ENDIF
765      !
766      jpnij = jpni*jpnj
767      !
768   END SUBROUTINE mpp_init_partition
769
770
771   SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr )
772      !!----------------------------------------------------------------------
773      !!                     ***  ROUTINE factorise  ***
774      !!
775      !! ** Purpose :   return the prime factors of n.
776      !!                knfax factors are returned in array kfax which is of
777      !!                maximum dimension kmaxfax.
778      !! ** Method  :
779      !!----------------------------------------------------------------------
780      INTEGER                    , INTENT(in   ) ::   kn, kmaxfax
781      INTEGER                    , INTENT(  out) ::   kerr, knfax
782      INTEGER, DIMENSION(kmaxfax), INTENT(  out) ::   kfax
783      !
784      INTEGER :: ifac, jl, inu
785      INTEGER, PARAMETER :: ntest = 14
786      INTEGER, DIMENSION(ntest) ::   ilfax
787      !!----------------------------------------------------------------------
788      !
789      ! lfax contains the set of allowed factors.
790      ilfax(:) = (/(2**jl,jl=ntest,1,-1)/)
791      !
792      ! Clear the error flag and initialise output vars
793      kerr  = 0
794      kfax  = 1
795      knfax = 0
796      !
797      IF( kn /= 1 ) THEN      ! Find the factors of n
798         !
799         ! nu holds the unfactorised part of the number.
800         ! knfax holds the number of factors found.
801         ! l points to the allowed factor list.
802         ! ifac holds the current factor.
803         !
804         inu   = kn
805         knfax = 0
806         !
807         DO jl = ntest, 1, -1
808            !
809            ifac = ilfax(jl)
810            IF( ifac > inu )   CYCLE
811            !
812            ! Test whether the factor will divide.
813            !
814            IF( MOD(inu,ifac) == 0 ) THEN
815               !
816               knfax = knfax + 1            ! Add the factor to the list
817               IF( knfax > kmaxfax ) THEN
818                  kerr = 6
819                  write (*,*) 'FACTOR: insufficient space in factor array ', knfax
820                  return
821               ENDIF
822               kfax(knfax) = ifac
823               ! Store the other factor that goes with this one
824               knfax = knfax + 1
825               kfax(knfax) = inu / ifac
826               !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax)
827            ENDIF
828            !
829         END DO
830         !
831      ENDIF
832      !
833   END SUBROUTINE factorise
834
835
836   SUBROUTINE mpp_init_nfdcom
837      !!----------------------------------------------------------------------
838      !!                     ***  ROUTINE  mpp_init_nfdcom  ***
839      !! ** Purpose :   Setup for north fold exchanges with explicit
840      !!                point-to-point messaging
841      !!
842      !! ** Method :   Initialization of the northern neighbours lists.
843      !!----------------------------------------------------------------------
844      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
845      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
846      !!----------------------------------------------------------------------
847      INTEGER  ::   sxM, dxM, sxT, dxT, jn
848      INTEGER  ::   njmppmax
849      !!----------------------------------------------------------------------
850      !
851      njmppmax = MAXVAL( njmppt )
852      !
853      !initializes the north-fold communication variables
854      isendto(:) = 0
855      nsndto     = 0
856      !
857      IF ( njmpp == njmppmax ) THEN      ! if I am a process in the north
858         !
859         !sxM is the first point (in the global domain) needed to compute the north-fold for the current process
860         sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
861         !dxM is the last point (in the global domain) needed to compute the north-fold for the current process
862         dxM = jpiglo - nimppt(narea) + 2
863         !
864         ! loop over the other north-fold processes to find the processes
865         ! managing the points belonging to the sxT-dxT range
866         !
867         DO jn = 1, jpni
868            !
869            sxT = nfiimpp(jn, jpnj)                            ! sxT = 1st  point (in the global domain) of the jn process
870            dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1    ! dxT = last point (in the global domain) of the jn process
871            !
872            IF    ( sxT < sxM  .AND.  sxM < dxT ) THEN
873               nsndto          = nsndto + 1
874               isendto(nsndto) = jn
875            ELSEIF( sxM <= sxT  .AND.  dxM >= dxT ) THEN
876               nsndto          = nsndto + 1
877               isendto(nsndto) = jn
878            ELSEIF( dxM <  dxT  .AND.  sxT <  dxM ) THEN
879               nsndto          = nsndto + 1
880               isendto(nsndto) = jn
881            ENDIF
882            !
883         END DO
884         nfsloop = 1
885         nfeloop = nlci
886         DO jn = 2,jpni-1
887            IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
888               IF( nfipproc(jn-1,jpnj) == -1 )   nfsloop = nldi
889               IF( nfipproc(jn+1,jpnj) == -1 )   nfeloop = nlei
890            ENDIF
891         END DO
892         !
893      ENDIF
894      l_north_nogather = .TRUE.
895      !
896   END SUBROUTINE mpp_init_nfdcom
897
898   
899#endif
900
901   !!======================================================================
902END MODULE mppini
Note: See TracBrowser for help on using the repository browser.