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/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/LBC – NEMO

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90 @ 9498

Last change on this file since 9498 was 9498, checked in by davestorkey, 6 years ago

branches/UKMO/dev_merge_2017_CICE_interface : revert changes so that I can clear SVN keywords first.

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