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 @ 9446

Last change on this file since 9446 was 9446, checked in by smasson, 6 years ago

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