New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
mppini.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90 @ 9436

Last change on this file since 9436 was 9436, checked in by acc, 6 years ago

Branch 2017/dev_merge_2017. Reorganisation of nemogcm.F90 and mppini.F90 to better separate domain decomposition functions from the rest of the initialisation. Stage 1: Main reorganisation; see ticket #2070

  • Property svn:keywords set to Id
File size: 37.0 KB
Line 
1MODULE mppini
2   !!======================================================================
3   !!                       ***  MODULE mppini   ***
4   !! Ocean initialization : distributed memory computing initialization
5   !!======================================================================
6   !! History :  6.0  !  1994-11  (M. Guyon)  Original code
7   !!  OPA       7.0  !  1995-04  (J. Escobar, M. Imbard)
8   !!            8.0  !  1998-05  (M. Imbard, J. Escobar, L. Colombet )  SHMEM and MPI versions
9   !!  NEMO      1.0  !  2004-01  (G. Madec, J.M Molines)  F90 : free form , north fold jpni > 1
10   !!            3.4  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE) add mpp_init_nfdcom
11   !!            3.   ! 2013-06  (I. Epicoco, S. Mocavero, CMCC) mpp_init_nfdcom: setup avoiding MPI communication
12   !!            4.0  !  2016-06  (G. Madec)  use domain configuration file instead of bathymetry file
13   !!            4.0  !  2017-06  (J.M. Molines, T. Lovato) merge of mppini and mppini_2
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!  mpp_init          : Lay out the global domain over processors with/without land processor elimination
18   !!  mpp_init_mask     : Read global bathymetric information to facilitate land suppression
19   !!  mpp_init_ioipsl   : IOIPSL initialization in mpp
20   !!  mpp_init_partition: Calculate MPP domain decomposition
21   !!  factorise         : Calculate the factors of the no. of MPI processes
22   !!  mpp_init_nfdcom   : Setup for north fold exchanges with explicit point-to-point messaging
23   !!----------------------------------------------------------------------
24   USE dom_oce        ! ocean space and time domain
25   USE bdy_oce        ! open BounDarY 
26   !
27   USE lbcnfd  , ONLY : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges
28   USE lib_mpp        ! distribued memory computing library
29   USE iom            ! nemo I/O library
30   USE ioipsl         ! I/O IPSL library
31   USE in_out_manager ! I/O Manager
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC mpp_init       ! called by opa.F90
37
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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      nperio = jperio
81      nbondi = 2
82      nbondj = 2
83      nidom  = FLIO_DOM_NONE
84      npolj = jperio
85      !
86      IF(lwp) THEN
87         WRITE(numout,*)
88         WRITE(numout,*) 'mpp_init : NO massively parallel processing'
89         WRITE(numout,*) '~~~~~~~~ '
90         WRITE(numout,*) '   nperio = ', nperio, '   nimpp  = ', nimpp
91         WRITE(numout,*) '   npolj  = ', npolj , '   njmpp  = ', njmpp
92      ENDIF
93      !
94      IF(  jpni /= 1 .OR. jpnj /= 1 .OR. jpnij /= 1 )                                     &
95         CALL ctl_stop( 'mpp_init: equality  jpni = jpnj = jpnij = 1 is not satisfied',   &
96            &           'the domain is lay out for distributed memory computing!' )
97         !
98      IF( jperio == 7 ) CALL ctl_stop( 'mpp_init: jperio = 7 needs distributed memory computing ',       &
99         &                             'with 1 process. Add key_mpp_mpi in the list of active cpp keys ' )
100         !
101   END SUBROUTINE mpp_init
102
103#else
104   !!----------------------------------------------------------------------
105   !!   'key_mpp_mpi'                     MPI massively parallel processing
106   !!----------------------------------------------------------------------
107
108   SUBROUTINE mpp_init
109      !!----------------------------------------------------------------------
110      !!                  ***  ROUTINE mpp_init  ***
111      !!                   
112      !! ** Purpose :   Lay out the global domain over processors.
113      !!      If land processors are to be eliminated, this program requires the
114      !!      presence of the domain configuration file. Land processors elimination
115      !!      is performed if jpni x jpnj /= jpnij. In this case, using the MPP_PREP
116      !!      preprocessing tool, help for defining the best cutting out.
117      !!
118      !! ** Method  :   Global domain is distributed in smaller local domains.
119      !!      Periodic condition is a function of the local domain position
120      !!      (global boundary or neighbouring domain) and of the global
121      !!      periodic
122      !!      Type :         jperio global periodic condition
123      !!                     nperio local  periodic condition
124      !!
125      !! ** Action : - set domain parameters
126      !!                    nimpp     : longitudinal index
127      !!                    njmpp     : latitudinal  index
128      !!                    nperio    : lateral condition type
129      !!                    narea     : number for local area
130      !!                    nlci      : first dimension
131      !!                    nlcj      : second dimension
132      !!                    nbondi    : mark for "east-west local boundary"
133      !!                    nbondj    : mark for "north-south local boundary"
134      !!                    nproc     : number for local processor
135      !!                    noea      : number for local neighboring processor
136      !!                    nowe      : number for local neighboring processor
137      !!                    noso      : number for local neighboring processor
138      !!                    nono      : number for local neighboring processor
139      !!----------------------------------------------------------------------
140      INTEGER ::   ji, jj, jn, jproc, jarea   ! dummy loop indices
141      INTEGER ::   inum                       ! local logical unit
142      INTEGER ::   idir, ifreq, icont, isurf  ! local integers
143      INTEGER ::   ii, il1, ili, imil         !   -       -
144      INTEGER ::   ij, il2, ilj, ijm1         !   -       -
145      INTEGER ::   iino, ijno, iiso, ijso     !   -       -
146      INTEGER ::   iiea, ijea, iiwe, ijwe     !   -       -
147      INTEGER ::   iresti, irestj, iproc      !   -       -
148      INTEGER ::   ierr                       ! local logical unit
149      REAL(wp)::   zidom, zjdom               ! local scalars
150      INTEGER, DIMENSION(jpnij)     ::   iin, ii_nono, ii_noea          ! 1D workspace
151      INTEGER, DIMENSION(jpnij)     ::   ijn, ii_noso, ii_nowe          !  -     -
152      INTEGER, DIMENSION(jpni,jpnj) ::   iimppt, ilci, ibondi, ipproc   ! 2D workspace
153      INTEGER, DIMENSION(jpni,jpnj) ::   ijmppt, ilcj, ibondj, ipolj    !  -     -
154      INTEGER, DIMENSION(jpni,jpnj) ::   ilei, ildi, iono, ioea         !  -     -
155      INTEGER, DIMENSION(jpni,jpnj) ::   ilej, ildj, ioso, iowe         !  -     -
156      INTEGER, DIMENSION(jpiglo,jpjglo) ::   imask   ! 2D global domain workspace
157      !!----------------------------------------------------------------------
158
159      ! If dimensions of processor grid weren't specified in the namelist file
160      ! then we calculate them here now that we have our communicator size
161      IF( jpni < 1 .OR. jpnj < 1 ) THEN
162         IF( Agrif_Root() )   CALL mpp_init_partition( mppsize )
163      ENDIF
164      !
165#if defined key_agrif
166      IF( jpnij /= jpni*jpnj ) CALL ctl_stop( 'STOP', 'Cannot remove land proc with AGRIF' )
167#endif
168
169      !
170      ALLOCATE(  nfiimpp(jpni,jpnj), nfipproc(jpni,jpnj), nfilcit(jpni,jpnj) ,    &
171         &       nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) ,    &
172         &       njmppt(jpnij) , ibonjt(jpnij) , nldit(jpnij) , nldjt(jpnij) ,    &
173         &                                       nleit(jpnij) , nlejt(jpnij) , STAT=ierr )
174      CALL mpp_sum( ierr )
175      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' )
176     
177      !
178#if defined key_agrif
179      IF( .NOT. Agrif_Root() ) THEN       ! AGRIF children: specific setting (cf. agrif_user.F90)
180         jpiglo  = nbcellsx + 2 + 2*nbghostcells
181         jpjglo  = nbcellsy + 2 + 2*nbghostcells
182         jpimax  = ( jpiglo-2*nn_hls + (jpni-1+0) ) / jpni + 2*nn_hls
183         jpjmax  = ( jpjglo-2*nn_hls + (jpnj-1+0) ) / jpnj + 2*nn_hls
184         nperio  = 0
185         jperio  = 0
186         ln_use_jattr = .false.
187      ENDIF
188#endif
189
190      IF( Agrif_Root() ) THEN       ! AGRIF mother: specific setting from jpni and jpnj
191#if defined key_nemocice_decomp
192         jpimax = ( nx_global+2-2*nn_hls + (jpni-1) ) / jpni + 2*nn_hls    ! first  dim.
193         jpjmax = ( ny_global+2-2*nn_hls + (jpnj-1) ) / jpnj + 2*nn_hls    ! second dim.
194#else
195         jpimax = ( jpiglo - 2*nn_hls + (jpni-1) ) / jpni + 2*nn_hls    ! first  dim.
196         jpjmax = ( jpjglo - 2*nn_hls + (jpnj-1) ) / jpnj + 2*nn_hls    ! second dim.
197#endif
198      ENDIF
199
200      !
201      IF ( jpni * jpnj == jpnij ) THEN    ! regular domain lay out over processors
202         imask(:,:) = 1               
203      ELSEIF ( jpni*jpnj > jpnij ) THEN   ! remove land-only processor (i.e. where imask(:,:)=0)
204         CALL mpp_init_mask( imask )   
205      ELSE                                ! error
206         CALL ctl_stop( 'mpp_init: jpnij > jpni x jpnj. Check namelist setting!' )
207      ENDIF
208      !
209      !  1. Dimension arrays for subdomains
210      ! -----------------------------------
211      !  Computation of local domain sizes ilci() ilcj()
212      !  These dimensions depend on global sizes jpni,jpnj and jpiglo,jpjglo
213      !  The subdomains are squares lesser than or equal to the global
214      !  dimensions divided by the number of processors minus the overlap array.
215      !
216      nreci = 2 * nn_hls
217      nrecj = 2 * nn_hls
218      iresti = 1 + MOD( jpiglo - nreci -1 , jpni )
219      irestj = 1 + MOD( jpjglo - nrecj -1 , jpnj )
220      !
221      !  Need to use jpimax and jpjmax here since jpi and jpj not yet defined
222#if defined key_nemocice_decomp
223      ! Change padding to be consistent with CICE
224      ilci(1:jpni-1      ,:) = jpimax
225      ilci(jpni          ,:) = jpiglo - (jpni - 1) * (jpimax - nreci)
226      !
227      ilcj(:,      1:jpnj-1) = jpjmax
228      ilcj(:,          jpnj) = jpjglo - (jpnj - 1) * (jpjmax - nrecj)
229#else
230      ilci(1:iresti      ,:) = jpimax
231      ilci(iresti+1:jpni ,:) = jpimax-1
232
233      ilcj(:,      1:irestj) = jpjmax
234      ilcj(:, irestj+1:jpnj) = jpjmax-1
235#endif
236      !
237      zidom = nreci + sum(ilci(:,1) - nreci )
238      zjdom = nrecj + sum(ilcj(1,:) - nrecj )
239      !
240      IF(lwp) THEN
241         WRITE(numout,*)
242         WRITE(numout,*) 'mpp_init : MPI Message Passing MPI - domain lay out over processors'
243         WRITE(numout,*) '~~~~~~~~ '
244         WRITE(numout,*) '   defines mpp subdomains'
245         WRITE(numout,*) '      iresti = ', iresti, ' jpni = ', jpni 
246         WRITE(numout,*) '      irestj = ', irestj, ' jpnj = ', jpnj
247         WRITE(numout,*)
248         WRITE(numout,*) '      sum ilci(i,1) = ', zidom, ' jpiglo = ', jpiglo
249         WRITE(numout,*) '      sum ilcj(1,j) = ', zjdom, ' jpjglo = ', jpjglo
250      ENDIF
251
252      !  2. Index arrays for subdomains
253      ! -------------------------------
254      iimppt(:,:) =  1
255      ijmppt(:,:) =  1
256      ipproc(:,:) = -1
257      !
258      IF( jpni > 1 ) THEN
259         DO jj = 1, jpnj
260            DO ji = 2, jpni
261               iimppt(ji,jj) = iimppt(ji-1,jj) + ilci(ji-1,jj) - nreci
262            END DO
263         END DO
264      ENDIF
265      nfiimpp(:,:) = iimppt(:,:)
266      !
267      IF( jpnj > 1 )THEN
268         DO jj = 2, jpnj
269            DO ji = 1, jpni
270               ijmppt(ji,jj) = ijmppt(ji,jj-1) + ilcj(ji,jj-1) - nrecj
271            END DO
272         END DO
273      ENDIF
274
275      ! 3. Subdomain description in the Regular Case
276      ! --------------------------------------------
277      nperio = 0
278      icont = -1
279      DO jarea = 1, jpni*jpnj
280         ii = 1 + MOD(jarea-1,jpni)
281         ij = 1 +    (jarea-1)/jpni
282         ili = ilci(ii,ij)
283         ilj = ilcj(ii,ij)
284         ibondj(ii,ij) = -1
285         IF( jarea >  jpni          )   ibondj(ii,ij) = 0
286         IF( jarea >  (jpnj-1)*jpni )   ibondj(ii,ij) = 1
287         IF( jpnj  == 1             )   ibondj(ii,ij) = 2
288         ibondi(ii,ij) = 0
289         IF( MOD(jarea,jpni) ==  1  )   ibondi(ii,ij) = -1
290         IF( MOD(jarea,jpni) ==  0  )   ibondi(ii,ij) =  1
291         IF( jpni            ==  1  )   ibondi(ii,ij) =  2
292
293         ! Subdomain neighbors (get their zone number): default definition
294         iproc = jarea - 1
295         ioso(ii,ij) = iproc - jpni
296         iowe(ii,ij) = iproc - 1
297         ioea(ii,ij) = iproc + 1
298         iono(ii,ij) = iproc + jpni
299         ildi(ii,ij) =  1  + nn_hls
300         ilei(ii,ij) = ili - nn_hls
301         ildj(ii,ij) =  1  + nn_hls
302         ilej(ii,ij) = ilj - nn_hls
303
304         ! East-West periodicity: change ibondi, ioea, iowe
305         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7 ) THEN
306            IF( jpni == 1 )THEN
307               ibondi(ii,ij) = 2
308               nperio = 1
309            ELSE
310               ibondi(ii,ij) = 0
311            ENDIF
312            IF( MOD(jarea,jpni) == 0 ) THEN
313               ioea(ii,ij) = iproc - (jpni-1)
314            ENDIF
315            IF( MOD(jarea,jpni) == 1 ) THEN
316               iowe(ii,ij) = iproc + jpni - 1
317            ENDIF
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         isurf = 0
339         DO jj = 1, ilj
340            DO  ji = 1, ili
341               IF( imask(ji+iimppt(ii,ij)-1, jj+ijmppt(ii,ij)-1) == 1)   isurf = isurf+1
342            END DO
343         END DO
344         !
345         IF( isurf /= 0 ) THEN
346            icont = icont + 1
347            ipproc(ii,ij) = icont
348            iin(icont+1) = ii
349            ijn(icont+1) = ij
350         ENDIF
351      END DO
352      !
353      nfipproc(:,:) = ipproc(:,:)
354
355      ! Check potential error
356      IF( icont+1 /= jpnij ) THEN
357         WRITE(ctmp1,*) ' jpni =',jpni,' jpnj =',jpnj
358         WRITE(ctmp2,*) ' jpnij =',jpnij, '< jpni x jpnj' 
359         WRITE(ctmp3,*) ' ***********, mpp_init2 finds jpnij=',icont+1
360         CALL ctl_stop( 'STOP', 'mpp_init: Eliminate land processors algorithm', '', ctmp1, ctmp2, '', ctmp3 )
361      ENDIF
362
363      ! 4. Subdomain print
364      ! ------------------
365      IF(lwp) THEN
366         ifreq = 4
367         il1 = 1
368         DO jn = 1, (jpni-1)/ifreq+1
369            il2 = MIN(jpni,il1+ifreq-1)
370            WRITE(numout,*)
371            WRITE(numout,9400) ('***',ji=il1,il2-1)
372            DO jj = jpnj, 1, -1
373               WRITE(numout,9403) ('   ',ji=il1,il2-1)
374               WRITE(numout,9402) jj, (ilci(ji,jj),ilcj(ji,jj),ji=il1,il2)
375               WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2)
376               WRITE(numout,9403) ('   ',ji=il1,il2-1)
377               WRITE(numout,9400) ('***',ji=il1,il2-1)
378            END DO
379            WRITE(numout,9401) (ji,ji=il1,il2)
380            il1 = il1+ifreq
381         END DO
382 9400    FORMAT('           ***'   ,20('*************',a3)    )
383 9403    FORMAT('           *     ',20('         *   ',a3)    )
384 9401    FORMAT('              '   ,20('   ',i3,'          ') )
385 9402    FORMAT('       ',i3,' *  ',20(i3,'  x',i3,'   *   ') )
386 9404    FORMAT('           *  '   ,20('      ',i3,'   *   ') )
387      ENDIF
388
389      ! 5. neighbour treatment: change ibondi, ibondj if next to a land zone
390      ! ----------------------
391      DO jarea = 1, jpni*jpnj
392         iproc = jarea-1
393         ii = 1 + MOD( jarea-1  , jpni )
394         ij = 1 +     (jarea-1) / jpni
395         IF ( ipproc(ii,ij) == -1 .AND. 0 <= iono(ii,ij) .AND. iono(ii,ij) <= jpni*jpnj-1 ) THEN
396            iino = 1 + MOD( iono(ii,ij) , jpni )
397            ijno = 1 +      iono(ii,ij) / jpni
398            ! Need to reverse the logical direction of communication
399            ! for northern neighbours of northern row processors (north-fold)
400            ! i.e. need to check that the northern neighbour only communicates
401            ! to the SOUTH (or not at all) if this area is land-only (#1057)
402            idir = 1
403            IF( ij == jpnj .AND. ijno == jpnj )   idir = -1   
404            IF( ibondj(iino,ijno) == idir     )   ibondj(iino,ijno) =   2
405            IF( ibondj(iino,ijno) == 0        )   ibondj(iino,ijno) = -idir
406         ENDIF
407         IF( ipproc(ii,ij) == -1 .AND. 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= jpni*jpnj-1 ) THEN
408            iiso = 1 + MOD( ioso(ii,ij) , jpni )
409            ijso = 1 +      ioso(ii,ij) / jpni
410            IF( ibondj(iiso,ijso) == -1 )   ibondj(iiso,ijso) = 2
411            IF( ibondj(iiso,ijso) ==  0 )   ibondj(iiso,ijso) = 1
412         ENDIF
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 )
415            ijea = 1 +      ioea(ii,ij) / jpni
416            IF( ibondi(iiea,ijea) == 1 )   ibondi(iiea,ijea) =  2
417            IF( ibondi(iiea,ijea) == 0 )   ibondi(iiea,ijea) = -1
418         ENDIF
419         IF( ipproc(ii,ij) == -1 .AND. 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= jpni*jpnj-1) THEN
420            iiwe = 1 + MOD( iowe(ii,ij) , jpni )
421            ijwe = 1 +      iowe(ii,ij) / jpni
422            IF( ibondi(iiwe,ijwe) == -1 )   ibondi(iiwe,ijwe) = 2
423            IF( ibondi(iiwe,ijwe) ==  0 )   ibondi(iiwe,ijwe) = 1
424         ENDIF
425      END DO
426
427      ! Update il[de][ij] according to modified ibond[ij]
428      ! ----------------------
429      DO jarea = 1, jpni*jpnj
430         ii = iin(narea)
431         ij = ijn(narea)
432         IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) =  1
433         IF( ibondi(ii,ij) ==  1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij)
434         IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) =  1
435         IF( ibondj(ii,ij) ==  1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij)
436      END DO
437         
438      ! just to save nono etc for all proc
439      ! warning ii*ij (zone) /= nproc (processors)!
440      ! ioso = zone number, ii_noso = proc number
441      ii_noso(:) = -1
442      ii_nono(:) = -1
443      ii_noea(:) = -1
444      ii_nowe(:) = -1 
445      DO jarea = 1, jpnij
446         ii = iin(jarea)
447         ij = ijn(jarea)
448         IF( 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= (jpni*jpnj-1) ) THEN
449            iiso = 1 + MOD( ioso(ii,ij) , jpni )
450            ijso = 1 +      ioso(ii,ij) / jpni
451            ii_noso(jarea) = ipproc(iiso,ijso)
452         ENDIF
453         IF( 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= (jpni*jpnj-1) ) THEN
454          iiwe = 1 + MOD( iowe(ii,ij) , jpni )
455          ijwe = 1 +      iowe(ii,ij) / jpni
456          ii_nowe(jarea) = ipproc(iiwe,ijwe)
457         ENDIF
458         IF( 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= (jpni*jpnj-1) ) THEN
459            iiea = 1 + MOD( ioea(ii,ij) , jpni )
460            ijea = 1 +      ioea(ii,ij) / jpni
461            ii_noea(jarea)= ipproc(iiea,ijea)
462         ENDIF
463         IF( 0 <= iono(ii,ij) .AND. iono(ii,ij) <= (jpni*jpnj-1) ) THEN
464            iino = 1 + MOD( iono(ii,ij) , jpni )
465            ijno = 1 +      iono(ii,ij) / jpni
466            ii_nono(jarea)= ipproc(iino,ijno)
467         ENDIF
468      END DO
469   
470      ! 6. Change processor name
471      ! ------------------------
472      ii = iin(narea)
473      ij = ijn(narea)
474      !
475      ! set default neighbours
476      noso = ii_noso(narea)
477      nowe = ii_nowe(narea)
478      noea = ii_noea(narea)
479      nono = ii_nono(narea)
480      nlci = ilci(ii,ij) 
481      nldi = ildi(ii,ij)
482      nlei = ilei(ii,ij)
483      nlcj = ilcj(ii,ij) 
484      nldj = ildj(ii,ij)
485      nlej = ilej(ii,ij)
486      nbondi = ibondi(ii,ij)
487      nbondj = ibondj(ii,ij)
488      nimpp = iimppt(ii,ij) 
489      njmpp = ijmppt(ii,ij)
490      jpi = nlci
491      jpj = nlcj
492      jpk = jpkglo                                             ! third dim
493#if defined key_agrif
494      ! simple trick to use same vertical grid as parent but different number of levels:
495      ! Save maximum number of levels in jpkglo, then define all vertical grids with this number.
496      ! Suppress once vertical online interpolation is ok
497      IF(.NOT.Agrif_Root())   jpkglo = Agrif_Parent( jpkglo )
498#endif
499      jpim1 = jpi-1                                            ! inner domain indices
500      jpjm1 = jpj-1                                            !   "           "
501      jpkm1 = MAX( 1, jpk-1 )                                  !   "           "
502      jpij  = jpi*jpj                                          !  jpi x j
503      DO jproc = 1, jpnij
504         ii = iin(jproc)
505         ij = ijn(jproc)
506         nlcit(jproc) = ilci(ii,ij)
507         nldit(jproc) = ildi(ii,ij)
508         nleit(jproc) = ilei(ii,ij)
509         nlcjt(jproc) = ilcj(ii,ij)
510         nldjt(jproc) = ildj(ii,ij)
511         nlejt(jproc) = ilej(ii,ij)
512         ibonit(jproc) = ibondi(ii,ij)
513         ibonjt(jproc) = ibondj(ii,ij)
514         nimppt(jproc) = iimppt(ii,ij) 
515         njmppt(jproc) = ijmppt(ii,ij) 
516         nfilcit(ii,ij) = ilci(ii,ij)
517      END DO
518
519      ! Save processor layout in ascii file
520      IF (lwp) THEN
521         CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
522         WRITE(inum,'(a)') '   jpnij   jpimax  jpjmax    jpk  jpiglo  jpjglo'//&
523   &           ' ( local:    narea     jpi     jpj )'
524         WRITE(inum,'(6i8,a,3i8,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,&
525   &           ' ( local: ',narea,jpi,jpj,' )'
526         WRITE(inum,'(a)') 'nproc nlci nlcj nldi nldj nlei nlej nimp njmp nono noso nowe noea nbondi nbondj '
527
528         DO jproc = 1, jpnij
529            ii = iin(jproc)
530            ij = ijn(jproc)
531            WRITE(inum,'(13i5,2i7)')   jproc-1, nlcit  (jproc), nlcjt  (jproc),   &
532               &                                nldit  (jproc), nldjt  (jproc),   &
533               &                                nleit  (jproc), nlejt  (jproc),   &
534               &                                nimppt (jproc), njmppt (jproc),   & 
535               &                                ii_nono(jproc), ii_noso(jproc),   &
536               &                                ii_nowe(jproc), ii_noea(jproc),   &
537               &                                ibondi (ii,ij), ibondj (ii,ij) 
538         END DO
539         CLOSE(inum)   
540      END IF
541
542      !                          ! north fold parameter
543      ! Defined npolj, either 0, 3 , 4 , 5 , 6
544      ! In this case the important thing is that npolj /= 0
545      ! Because if we go through these line it is because jpni >1 and thus
546      ! we must use lbcnorthmpp, which tests only npolj =0 or npolj /= 0
547      npolj = 0
548      ij = ijn(narea)
549      IF( jperio == 3 .OR. jperio == 4 ) THEN
550         IF( ij == jpnj )   npolj = 3
551      ENDIF
552      IF( jperio == 5 .OR. jperio == 6 ) THEN
553         IF( ij == jpnj )   npolj = 5
554      ENDIF
555      !
556      nproc = narea-1
557      IF(lwp) THEN
558         WRITE(numout,*)
559         WRITE(numout,*) '   resulting internal parameters : '
560         WRITE(numout,*) '      nproc  = ', nproc
561         WRITE(numout,*) '      nowe   = ', nowe  , '   noea  =  ', noea
562         WRITE(numout,*) '      nono   = ', nono  , '   noso  =  ', noso
563         WRITE(numout,*) '      nbondi = ', nbondi
564         WRITE(numout,*) '      nbondj = ', nbondj
565         WRITE(numout,*) '      npolj  = ', npolj
566         WRITE(numout,*) '      nperio = ', nperio
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      IF( nperio == 1 .AND. jpni /= 1 )   CALL ctl_stop( 'mpp_init: error on cyclicity' )
577
578      IF( jperio == 7 .AND. ( jpni /= 1 .OR. jpnj /= 1 ) )   &
579         &                  CALL ctl_stop( ' mpp_init: error jperio = 7 works only with jpni = jpnj = 1' )
580
581      !                          ! Prepare mpp north fold
582      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
583         CALL mpp_ini_north
584         IF(lwp) WRITE(numout,*)
585         IF(lwp) WRITE(numout,*) '   ==>>>   North fold boundary prepared for jpni >1'
586      ENDIF
587      !
588      CALL mpp_init_ioipsl       ! Prepare NetCDF output file (if necessary)
589      !
590      IF( ln_nnogather )   CALL mpp_init_nfdcom     ! northfold neighbour lists
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_lim, nn_ice_lim_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.