source: NEMO/branches/2018/dev_r9759_HPC09_ESIWACE/src/OCE/LBC/mppini.F90 @ 9772

Last change on this file since 9772 was 9772, checked in by smasson, 3 years ago

dev_r9759_HPC09_ESIWACE: add benchmark features

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