source: NEMO/trunk/src/OCE/LBC/mppini.F90 @ 9657

Last change on this file since 9657 was 9657, checked in by clem, 4 years ago

change nn_ice_lim and cn_ice_lim to nn_ice and cn_ice for BDY

  • Property svn:keywords set to Id
File size: 37.7 KB
Line 
1MODULE mppini
2   !!======================================================================
3   !!                       ***  MODULE mppini   ***
4   !! Ocean initialization : distributed memory computing initialization
5   !!======================================================================
6   !! History :  6.0  !  1994-11  (M. Guyon)  Original code
7   !!  OPA       7.0  !  1995-04  (J. Escobar, M. Imbard)
8   !!            8.0  !  1998-05  (M. Imbard, J. Escobar, L. Colombet )  SHMEM and MPI versions
9   !!  NEMO      1.0  !  2004-01  (G. Madec, J.M Molines)  F90 : free form , north fold jpni > 1
10   !!            3.4  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE) add mpp_init_nfdcom
11   !!            3.   ! 2013-06  (I. Epicoco, S. Mocavero, CMCC) mpp_init_nfdcom: setup avoiding MPI communication
12   !!            4.0  !  2016-06  (G. Madec)  use domain configuration file instead of bathymetry file
13   !!            4.0  !  2017-06  (J.M. Molines, T. Lovato) merge of mppini and mppini_2
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!  mpp_init          : Lay out the global domain over processors with/without land processor elimination
18   !!  mpp_init_mask     : Read global bathymetric information to facilitate land suppression
19   !!  mpp_init_ioipsl   : IOIPSL initialization in mpp
20   !!  mpp_init_partition: Calculate MPP domain decomposition
21   !!  factorise         : Calculate the factors of the no. of MPI processes
22   !!  mpp_init_nfdcom   : Setup for north fold exchanges with explicit point-to-point messaging
23   !!----------------------------------------------------------------------
24   USE dom_oce        ! ocean space and time domain
25   USE bdy_oce        ! open BounDarY 
26   !
27   USE lbcnfd  , ONLY : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges
28   USE lib_mpp        ! distribued memory computing library
29   USE iom            ! nemo I/O library
30   USE ioipsl         ! I/O IPSL library
31   USE in_out_manager ! I/O Manager
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC mpp_init       ! called by opa.F90
37
38   !!----------------------------------------------------------------------
39   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (./LICENSE)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45#if ! defined key_mpp_mpi
46   !!----------------------------------------------------------------------
47   !!   Default option :                            shared memory computing
48   !!----------------------------------------------------------------------
49
50   SUBROUTINE mpp_init
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE mpp_init  ***
53      !!
54      !! ** Purpose :   Lay out the global domain over processors.
55      !!
56      !! ** Method  :   Shared memory computing, set the local processor
57      !!              variables to the value of the global domain
58      !!----------------------------------------------------------------------
59      !
60      jpimax = jpiglo
61      jpjmax = jpjglo
62      jpi    = jpiglo
63      jpj    = jpjglo
64      jpk    = jpjglo
65      jpim1  = jpi-1                                            ! inner domain indices
66      jpjm1  = jpj-1                                            !   "           "
67      jpkm1  = MAX( 1, jpk-1 )                                  !   "           "
68      jpij   = jpi*jpj
69      jpni   = 1
70      jpnj   = 1
71      jpnij  = jpni*jpnj
72      nimpp  = 1           !
73      njmpp  = 1
74      nlci   = jpi
75      nlcj   = jpj
76      nldi   = 1
77      nldj   = 1
78      nlei   = jpi
79      nlej   = jpj
80      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, iarea0     !   -       -
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         iarea0 = jarea - 1
282         ii = 1 + MOD(iarea0,jpni)
283         ij = 1 +     iarea0/jpni
284         ili = ilci(ii,ij)
285         ilj = ilcj(ii,ij)
286         ibondj(ii,ij) = -1
287         IF( jarea >  jpni          )   ibondj(ii,ij) = 0
288         IF( jarea >  (jpnj-1)*jpni )   ibondj(ii,ij) = 1
289         IF( jpnj  == 1             )   ibondj(ii,ij) = 2
290         ibondi(ii,ij) = 0
291         IF( MOD(jarea,jpni) ==  1  )   ibondi(ii,ij) = -1
292         IF( MOD(jarea,jpni) ==  0  )   ibondi(ii,ij) =  1
293         IF( jpni            ==  1  )   ibondi(ii,ij) =  2
294
295         ! Subdomain neighbors (get their zone number): default definition
296         ioso(ii,ij) = iarea0 - jpni
297         iowe(ii,ij) = iarea0 - 1
298         ioea(ii,ij) = iarea0 + 1
299         iono(ii,ij) = iarea0 + jpni
300         ildi(ii,ij) =  1  + nn_hls
301         ilei(ii,ij) = ili - nn_hls
302         ildj(ii,ij) =  1  + nn_hls
303         ilej(ii,ij) = ilj - nn_hls
304
305         ! East-West periodicity: change ibondi, ioea, iowe
306         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7 ) THEN
307            IF( jpni == 1 )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) = iarea0 - (jpni-1)
315            ENDIF
316            IF( MOD(jarea,jpni) == 1 ) THEN
317               iowe(ii,ij) = iarea0 + 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         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 jproc = 1, jpnij
430         ii = iin(jproc)
431         ij = ijn(jproc)
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 jproc = 1, jpnij
446         ii = iin(jproc)
447         ij = ijn(jproc)
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(jproc) = 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(jproc) = 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(jproc)= 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(jproc)= 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            WRITE(inum,'(13i5,2i7)')   jproc-1, nlcit  (jproc), nlcjt  (jproc),   &
530               &                                nldit  (jproc), nldjt  (jproc),   &
531               &                                nleit  (jproc), nlejt  (jproc),   &
532               &                                nimppt (jproc), njmppt (jproc),   & 
533               &                                ii_nono(jproc), ii_noso(jproc),   &
534               &                                ii_nowe(jproc), ii_noea(jproc),   &
535               &                                ibonit (jproc), ibonjt (jproc) 
536         END DO
537         CLOSE(inum)   
538      END IF
539
540      !                          ! north fold parameter
541      ! Defined npolj, either 0, 3 , 4 , 5 , 6
542      ! In this case the important thing is that npolj /= 0
543      ! Because if we go through these line it is because jpni >1 and thus
544      ! we must use lbcnorthmpp, which tests only npolj =0 or npolj /= 0
545      npolj = 0
546      ij = ijn(narea)
547      IF( jperio == 3 .OR. jperio == 4 ) THEN
548         IF( ij == jpnj )   npolj = 3
549      ENDIF
550      IF( jperio == 5 .OR. jperio == 6 ) THEN
551         IF( ij == jpnj )   npolj = 5
552      ENDIF
553      !
554      nproc = narea-1
555      IF(lwp) THEN
556         WRITE(numout,*)
557         WRITE(numout,*) '   resulting internal parameters : '
558         WRITE(numout,*) '      nproc  = ', nproc
559         WRITE(numout,*) '      nowe   = ', nowe  , '   noea  =  ', noea
560         WRITE(numout,*) '      nono   = ', nono  , '   noso  =  ', noso
561         WRITE(numout,*) '      nbondi = ', nbondi
562         WRITE(numout,*) '      nbondj = ', nbondj
563         WRITE(numout,*) '      npolj  = ', npolj
564         WRITE(numout,*) '      nperio = ', nperio
565         WRITE(numout,*) '      nlci   = ', nlci
566         WRITE(numout,*) '      nlcj   = ', nlcj
567         WRITE(numout,*) '      nimpp  = ', nimpp
568         WRITE(numout,*) '      njmpp  = ', njmpp
569         WRITE(numout,*) '      nreci  = ', nreci 
570         WRITE(numout,*) '      nrecj  = ', nrecj 
571         WRITE(numout,*) '      nn_hls = ', nn_hls 
572      ENDIF
573 
574      IF( nperio == 1 .AND. jpni /= 1 )   CALL ctl_stop( 'mpp_init: error on cyclicity' )
575
576      IF( jperio == 7 .AND. ( jpni /= 1 .OR. jpnj /= 1 ) )   &
577         &                  CALL ctl_stop( ' mpp_init: error jperio = 7 works only with jpni = jpnj = 1' )
578
579      !                          ! Prepare mpp north fold
580      IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
581         CALL mpp_ini_north
582         IF(lwp) WRITE(numout,*)
583         IF(lwp) WRITE(numout,*) '   ==>>>   North fold boundary prepared for jpni >1'
584      ENDIF
585      !
586      CALL mpp_init_ioipsl       ! Prepare NetCDF output file (if necessary)
587      !
588      IF( ln_nnogather )   CALL mpp_init_nfdcom     ! northfold neighbour lists
589      !
590      DEALLOCATE(iin, ijn, ii_nono, ii_noea, ii_noso, ii_nowe,    &
591         &       iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj,   &
592         &       ilci, ilcj, ilei, ilej, ildi, ildj,              &
593         &       iono, ioea, ioso, iowe)
594      !
595    END SUBROUTINE mpp_init
596
597
598    SUBROUTINE mpp_init_mask( kmask )
599      !!----------------------------------------------------------------------
600      !!                  ***  ROUTINE mpp_init_mask  ***
601      !!
602      !! ** Purpose : Read relevant bathymetric information in a global array
603      !!              in order to provide a land/sea mask used for the elimination
604      !!              of land domains, in an mpp computation.
605      !!
606      !! ** Method  : Read the namelist ln_zco and ln_isfcav in namelist namzgr
607      !!              in order to choose the correct bathymetric information
608      !!              (file and variables) 
609      !!----------------------------------------------------------------------
610      INTEGER, DIMENSION(jpiglo,jpjglo), INTENT(out) ::   kmask   ! global domain
611 
612      INTEGER :: inum   !: logical unit for configuration file
613      INTEGER :: ios    !: iostat error flag
614      INTEGER ::  ijstartrow                   ! temporary integers
615      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zbot, zbdy          ! global workspace
616      REAL(wp) ::   zidom , zjdom          ! local scalars
617      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           &
618           &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
619           &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             & 
620           &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
621           &             cn_ice, nn_ice_dta,                                     &
622           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
623           &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
624      !!----------------------------------------------------------------------
625      ! 0. initialisation
626      ! -----------------
627      CALL iom_open( cn_domcfg, inum )
628      !
629      ! ocean bottom level
630      CALL iom_get( inum, jpdom_unknown, 'bottom_level' , zbot , lrowattr=ln_use_jattr )  ! nb of ocean T-points
631      !
632      CALL iom_close( inum )
633      !
634      ! 2D ocean mask (=1 if at least one level of the water column is ocean, =0 otherwise)
635      WHERE( zbot(:,:) > 0 )   ;   kmask(:,:) = 1
636      ELSEWHERE                ;   kmask(:,:) = 0
637      END WHERE
638 
639      ! Adjust kmask with bdy_msk if it exists
640 
641      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY
642      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
643903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp )
644      !
645      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY
646      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
647904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp )
648
649      IF( ln_bdy .AND. ln_mask_file ) THEN
650         CALL iom_open( cn_mask_file, inum )
651         CALL iom_get ( inum, jpdom_unknown, 'bdy_msk', zbdy )
652         CALL iom_close( inum )
653         WHERE ( zbdy(:,:) <= 0. ) kmask = 0
654      ENDIF
655      !
656   END SUBROUTINE mpp_init_mask
657
658
659   SUBROUTINE mpp_init_ioipsl
660      !!----------------------------------------------------------------------
661      !!                  ***  ROUTINE mpp_init_ioipsl  ***
662      !!
663      !! ** Purpose :   
664      !!
665      !! ** Method  :   
666      !!
667      !! History :
668      !!   9.0  !  04-03  (G. Madec )  MPP-IOIPSL
669      !!   " "  !  08-12  (A. Coward)  addition in case of jpni*jpnj < jpnij
670      !!----------------------------------------------------------------------
671      INTEGER, DIMENSION(2) ::   iglo, iloc, iabsf, iabsl, ihals, ihale, idid
672      !!----------------------------------------------------------------------
673
674      ! The domain is split only horizontally along i- or/and j- direction
675      ! So we need at the most only 1D arrays with 2 elements.
676      ! Set idompar values equivalent to the jpdom_local_noextra definition
677      ! used in IOM. This works even if jpnij .ne. jpni*jpnj.
678      iglo(1) = jpiglo
679      iglo(2) = jpjglo
680      iloc(1) = nlci
681      iloc(2) = nlcj
682      iabsf(1) = nimppt(narea)
683      iabsf(2) = njmppt(narea)
684      iabsl(:) = iabsf(:) + iloc(:) - 1
685      ihals(1) = nldi - 1
686      ihals(2) = nldj - 1
687      ihale(1) = nlci - nlei
688      ihale(2) = nlcj - nlej
689      idid(1) = 1
690      idid(2) = 2
691
692      IF(lwp) THEN
693          WRITE(numout,*)
694          WRITE(numout,*) 'mpp_init_ioipsl :   iloc  = ', iloc (1), iloc (2)
695          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf(1), iabsf(2)
696          WRITE(numout,*) '                    ihals = ', ihals(1), ihals(2)
697          WRITE(numout,*) '                    ihale = ', ihale(1), ihale(2)
698      ENDIF
699      !
700      CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
701      !
702   END SUBROUTINE mpp_init_ioipsl 
703
704
705   SUBROUTINE mpp_init_partition( num_pes )
706      !!----------------------------------------------------------------------
707      !!                 ***  ROUTINE mpp_init_partition  ***
708      !!
709      !! ** Purpose :
710      !!
711      !! ** Method  :
712      !!----------------------------------------------------------------------
713      INTEGER, INTENT(in) ::   num_pes   ! The number of MPI processes we have
714      !
715      INTEGER, PARAMETER :: nfactmax = 20
716      INTEGER :: nfact ! The no. of factors returned
717      INTEGER :: ierr  ! Error flag
718      INTEGER :: ji
719      INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value
720      INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors
721      !!----------------------------------------------------------------------
722      !
723      ierr = 0
724      !
725      CALL factorise( ifact, nfactmax, nfact, num_pes, ierr )
726      !
727      IF( nfact <= 1 ) THEN
728         WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
729         WRITE (numout, *) '       : using grid of ',num_pes,' x 1'
730         jpnj = 1
731         jpni = num_pes
732      ELSE
733         ! Search through factors for the pair that are closest in value
734         mindiff = 1000000
735         imin    = 1
736         DO ji = 1, nfact-1, 2
737            idiff = ABS( ifact(ji) - ifact(ji+1) )
738            IF( idiff < mindiff ) THEN
739               mindiff = idiff
740               imin = ji
741            ENDIF
742         END DO
743         jpnj = ifact(imin)
744         jpni = ifact(imin + 1)
745      ENDIF
746      !
747      jpnij = jpni*jpnj
748      !
749   END SUBROUTINE mpp_init_partition
750
751
752   SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr )
753      !!----------------------------------------------------------------------
754      !!                     ***  ROUTINE factorise  ***
755      !!
756      !! ** Purpose :   return the prime factors of n.
757      !!                knfax factors are returned in array kfax which is of
758      !!                maximum dimension kmaxfax.
759      !! ** Method  :
760      !!----------------------------------------------------------------------
761      INTEGER                    , INTENT(in   ) ::   kn, kmaxfax
762      INTEGER                    , INTENT(  out) ::   kerr, knfax
763      INTEGER, DIMENSION(kmaxfax), INTENT(  out) ::   kfax
764      !
765      INTEGER :: ifac, jl, inu
766      INTEGER, PARAMETER :: ntest = 14
767      INTEGER, DIMENSION(ntest) ::   ilfax
768      !!----------------------------------------------------------------------
769      !
770      ! lfax contains the set of allowed factors.
771      ilfax(:) = (/(2**jl,jl=ntest,1,-1)/)
772      !
773      ! Clear the error flag and initialise output vars
774      kerr  = 0
775      kfax  = 1
776      knfax = 0
777      !
778      IF( kn /= 1 ) THEN      ! Find the factors of n
779         !
780         ! nu holds the unfactorised part of the number.
781         ! knfax holds the number of factors found.
782         ! l points to the allowed factor list.
783         ! ifac holds the current factor.
784         !
785         inu   = kn
786         knfax = 0
787         !
788         DO jl = ntest, 1, -1
789            !
790            ifac = ilfax(jl)
791            IF( ifac > inu )   CYCLE
792            !
793            ! Test whether the factor will divide.
794            !
795            IF( MOD(inu,ifac) == 0 ) THEN
796               !
797               knfax = knfax + 1            ! Add the factor to the list
798               IF( knfax > kmaxfax ) THEN
799                  kerr = 6
800                  write (*,*) 'FACTOR: insufficient space in factor array ', knfax
801                  return
802               ENDIF
803               kfax(knfax) = ifac
804               ! Store the other factor that goes with this one
805               knfax = knfax + 1
806               kfax(knfax) = inu / ifac
807               !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax)
808            ENDIF
809            !
810         END DO
811         !
812      ENDIF
813      !
814   END SUBROUTINE factorise
815
816
817   SUBROUTINE mpp_init_nfdcom
818      !!----------------------------------------------------------------------
819      !!                     ***  ROUTINE  mpp_init_nfdcom  ***
820      !! ** Purpose :   Setup for north fold exchanges with explicit
821      !!                point-to-point messaging
822      !!
823      !! ** Method :   Initialization of the northern neighbours lists.
824      !!----------------------------------------------------------------------
825      !!    1.0  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)
826      !!    2.0  ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
827      !!----------------------------------------------------------------------
828      INTEGER  ::   sxM, dxM, sxT, dxT, jn
829      INTEGER  ::   njmppmax
830      !!----------------------------------------------------------------------
831      !
832      njmppmax = MAXVAL( njmppt )
833      !
834      !initializes the north-fold communication variables
835      isendto(:) = 0
836      nsndto     = 0
837      !
838      IF ( njmpp == njmppmax ) THEN      ! if I am a process in the north
839         !
840         !sxM is the first point (in the global domain) needed to compute the north-fold for the current process
841         sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
842         !dxM is the last point (in the global domain) needed to compute the north-fold for the current process
843         dxM = jpiglo - nimppt(narea) + 2
844         !
845         ! loop over the other north-fold processes to find the processes
846         ! managing the points belonging to the sxT-dxT range
847         !
848         DO jn = 1, jpni
849            !
850            sxT = nfiimpp(jn, jpnj)                            ! sxT = 1st  point (in the global domain) of the jn process
851            dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1    ! dxT = last point (in the global domain) of the jn process
852            !
853            IF    ( sxT < sxM  .AND.  sxM < dxT ) THEN
854               nsndto          = nsndto + 1
855               isendto(nsndto) = jn
856            ELSEIF( sxM <= sxT  .AND.  dxM >= dxT ) THEN
857               nsndto          = nsndto + 1
858               isendto(nsndto) = jn
859            ELSEIF( dxM <  dxT  .AND.  sxT <  dxM ) THEN
860               nsndto          = nsndto + 1
861               isendto(nsndto) = jn
862            ENDIF
863            !
864         END DO
865         nfsloop = 1
866         nfeloop = nlci
867         DO jn = 2,jpni-1
868            IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
869               IF( nfipproc(jn-1,jpnj) == -1 )   nfsloop = nldi
870               IF( nfipproc(jn+1,jpnj) == -1 )   nfeloop = nlei
871            ENDIF
872         END DO
873         !
874      ENDIF
875      l_north_nogather = .TRUE.
876      !
877   END SUBROUTINE mpp_init_nfdcom
878
879   
880#endif
881
882   !!======================================================================
883END MODULE mppini
Note: See TracBrowser for help on using the repository browser.