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

Last change on this file since 9508 was 9508, checked in by clem, 2 years ago

add forgotten key_agrif

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