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 NEMO/trunk/src/OCE/LBC – NEMO

Context Navigation

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

Last change on this file since 9667 was 9667, checked in by smasson, 7 years ago
trunk: cyclic north-south periodicity and nperio cleaning, see #2093
Property svn:keywords set to `Id`
File size: 39.5 KB

Line
1	MODULE 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	!!----------------------------------------------------------------------
43	CONTAINS
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	nbondi = 2
81	nbondj = 2
82	nidom = FLIO_DOM_NONE
83	npolj = jperio
84	l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
85	l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
86	!
87	IF(lwp) THEN
88	WRITE(numout,*)
89	WRITE(numout,*) 'mpp_init : NO massively parallel processing'
90	WRITE(numout,*) '~~~~~~~~ '
91	WRITE(numout,*) ' l_Iperio = ', l_Iperio, ' l_Jperio = ', l_Jperio
92	WRITE(numout,*) ' npolj = ', npolj , ' njmpp = ', njmpp
93	ENDIF
94	!
95	IF( jpni /= 1 .OR. jpnj /= 1 .OR. jpnij /= 1 ) &
96	CALL ctl_stop( 'mpp_init: equality jpni = jpnj = jpnij = 1 is not satisfied', &
97	& 'the domain is lay out for distributed memory computing!' )
98	!
99	END SUBROUTINE mpp_init
100
101	#else
102	!!----------------------------------------------------------------------
103	!! 'key_mpp_mpi' MPI massively parallel processing
104	!!----------------------------------------------------------------------
105
106
107	SUBROUTINE mpp_init
108	!!----------------------------------------------------------------------
109	!! * ROUTINE mpp_init *
110	!!
111	!! ** Purpose : Lay out the global domain over processors.
112	!! If land processors are to be eliminated, this program requires the
113	!! presence of the domain configuration file. Land processors elimination
114	!! is performed if jpni x jpnj /= jpnij. In this case, using the MPP_PREP
115	!! preprocessing tool, help for defining the best cutting out.
116	!!
117	!! ** Method : Global domain is distributed in smaller local domains.
118	!! Periodic condition is a function of the local domain position
119	!! (global boundary or neighbouring domain) and of the global
120	!! periodic
121	!! Type : jperio global periodic condition
122	!!
123	!! ** Action : - set domain parameters
124	!! nimpp : longitudinal index
125	!! njmpp : latitudinal index
126	!! narea : number for local area
127	!! nlci : first dimension
128	!! nlcj : second dimension
129	!! nbondi : mark for "east-west local boundary"
130	!! nbondj : mark for "north-south local boundary"
131	!! nproc : number for local processor
132	!! noea : number for local neighboring processor
133	!! nowe : number for local neighboring processor
134	!! noso : number for local neighboring processor
135	!! nono : number for local neighboring processor
136	!!----------------------------------------------------------------------
137	INTEGER :: ji, jj, jn, jproc, jarea ! dummy loop indices
138	INTEGER :: inum ! local logical unit
139	INTEGER :: idir, ifreq, icont, isurf ! local integers
140	INTEGER :: ii, il1, ili, imil ! - -
141	INTEGER :: ij, il2, ilj, ijm1 ! - -
142	INTEGER :: iino, ijno, iiso, ijso ! - -
143	INTEGER :: iiea, ijea, iiwe, ijwe ! - -
144	INTEGER :: iresti, irestj, iarea0 ! - -
145	INTEGER :: ierr ! local logical unit
146	REAL(wp):: zidom, zjdom ! local scalars
147	INTEGER, ALLOCATABLE, DIMENSION(:) :: iin, ii_nono, ii_noea ! 1D workspace
148	INTEGER, ALLOCATABLE, DIMENSION(:) :: ijn, ii_noso, ii_nowe ! - -
149	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ilci, ibondi, ipproc ! 2D workspace
150	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ilcj, ibondj, ipolj ! - -
151	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilei, ildi, iono, ioea ! - -
152	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilej, ildj, ioso, iowe ! - -
153	INTEGER, DIMENSION(jpiglo,jpjglo) :: imask ! 2D global domain workspace
154	!!----------------------------------------------------------------------
155
156	! If dimensions of processor grid weren't specified in the namelist file
157	! then we calculate them here now that we have our communicator size
158	IF( jpni < 1 .OR. jpnj < 1 ) CALL mpp_init_partition( mppsize )
159	!
160	#if defined key_agrif
161	IF( jpnij /= jpni*jpnj ) CALL ctl_stop( 'STOP', 'Cannot remove land proc with AGRIF' )
162	#endif
163	!
164	ALLOCATE( nfiimpp(jpni,jpnj), nfipproc(jpni,jpnj), nfilcit(jpni,jpnj) , &
165	& nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) , &
166	& njmppt(jpnij) , ibonjt(jpnij) , nldit(jpnij) , nldjt(jpnij) , &
167	& nleit(jpnij) , nlejt(jpnij) , &
168	& iin(jpnij), ii_nono(jpnij), ii_noea(jpnij), &
169	& ijn(jpnij), ii_noso(jpnij), ii_nowe(jpnij), &
170	& iimppt(jpni,jpnj), ilci(jpni,jpnj), ibondi(jpni,jpnj), ipproc(jpni,jpnj), &
171	& ijmppt(jpni,jpnj), ilcj(jpni,jpnj), ibondj(jpni,jpnj), ipolj(jpni,jpnj), &
172	& ilei(jpni,jpnj), ildi(jpni,jpnj), iono(jpni,jpnj), ioea(jpni,jpnj), &
173	& ilej(jpni,jpnj), ildj(jpni,jpnj), ioso(jpni,jpnj), iowe(jpni,jpnj), &
174	& STAT=ierr )
175	CALL mpp_sum( ierr )
176	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' )
177
178	!
179	#if defined key_agrif
180	IF( .NOT. Agrif_Root() ) THEN ! AGRIF children: specific setting (cf. agrif_user.F90)
181	IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells ) &
182	CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells' )
183	IF( jpjglo /= nbcellsy + 2 + 2*nbghostcells ) &
184	CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpjglo == nbcellsy + 2 + 2*nbghostcells' )
185	IF( ln_use_jattr ) CALL ctl_stop( 'STOP', 'mpp_init:Agrif children requires ln_use_jattr = .false. ' )
186	ENDIF
187	#endif
188
189	#if defined key_nemocice_decomp
190	jpimax = ( nx_global+2-2nn_hls + (jpni-1) ) / jpni + 2nn_hls ! first dim.
191	jpjmax = ( ny_global+2-2nn_hls + (jpnj-1) ) / jpnj + 2nn_hls ! second dim.
192	#else
193	jpimax = ( jpiglo - 2nn_hls + (jpni-1) ) / jpni + 2nn_hls ! first dim.
194	jpjmax = ( jpjglo - 2nn_hls + (jpnj-1) ) / jpnj + 2nn_hls ! second dim.
195	#endif
196
197	!
198	IF ( jpni * jpnj == jpnij ) THEN ! regular domain lay out over processors
199	imask(:,:) = 1
200	ELSEIF ( jpni*jpnj > jpnij ) THEN ! remove land-only processor (i.e. where imask(:,:)=0)
201	CALL mpp_init_mask( imask )
202	ELSE ! error
203	CALL ctl_stop( 'mpp_init: jpnij > jpni x jpnj. Check namelist setting!' )
204	ENDIF
205	!
206	! 1. Dimension arrays for subdomains
207	! -----------------------------------
208	! Computation of local domain sizes ilci() ilcj()
209	! These dimensions depend on global sizes jpni,jpnj and jpiglo,jpjglo
210	! The subdomains are squares lesser than or equal to the global
211	! dimensions divided by the number of processors minus the overlap array.
212	!
213	nreci = 2 * nn_hls
214	nrecj = 2 * nn_hls
215	iresti = 1 + MOD( jpiglo - nreci -1 , jpni )
216	irestj = 1 + MOD( jpjglo - nrecj -1 , jpnj )
217	!
218	! Need to use jpimax and jpjmax here since jpi and jpj not yet defined
219	#if defined key_nemocice_decomp
220	! Change padding to be consistent with CICE
221	ilci(1:jpni-1 ,:) = jpimax
222	ilci(jpni ,:) = jpiglo - (jpni - 1) * (jpimax - nreci)
223	!
224	ilcj(:, 1:jpnj-1) = jpjmax
225	ilcj(:, jpnj) = jpjglo - (jpnj - 1) * (jpjmax - nrecj)
226	#else
227	ilci(1:iresti ,:) = jpimax
228	ilci(iresti+1:jpni ,:) = jpimax-1
229
230	ilcj(:, 1:irestj) = jpjmax
231	ilcj(:, irestj+1:jpnj) = jpjmax-1
232	#endif
233	!
234	zidom = nreci + sum(ilci(:,1) - nreci )
235	zjdom = nrecj + sum(ilcj(1,:) - nrecj )
236	!
237	IF(lwp) THEN
238	WRITE(numout,*)
239	WRITE(numout,*) 'mpp_init : MPI Message Passing MPI - domain lay out over processors'
240	WRITE(numout,*) '~~~~~~~~ '
241	WRITE(numout,*) ' defines mpp subdomains'
242	WRITE(numout,*) ' iresti = ', iresti, ' jpni = ', jpni
243	WRITE(numout,*) ' irestj = ', irestj, ' jpnj = ', jpnj
244	WRITE(numout,*)
245	WRITE(numout,*) ' sum ilci(i,1) = ', zidom, ' jpiglo = ', jpiglo
246	WRITE(numout,*) ' sum ilcj(1,j) = ', zjdom, ' jpjglo = ', jpjglo
247	ENDIF
248
249	! 2. Index arrays for subdomains
250	! -------------------------------
251	iimppt(:,:) = 1
252	ijmppt(:,:) = 1
253	ipproc(:,:) = -1
254	!
255	IF( jpni > 1 ) THEN
256	DO jj = 1, jpnj
257	DO ji = 2, jpni
258	iimppt(ji,jj) = iimppt(ji-1,jj) + ilci(ji-1,jj) - nreci
259	END DO
260	END DO
261	ENDIF
262	nfiimpp(:,:) = iimppt(:,:)
263	!
264	IF( jpnj > 1 )THEN
265	DO jj = 2, jpnj
266	DO ji = 1, jpni
267	ijmppt(ji,jj) = ijmppt(ji,jj-1) + ilcj(ji,jj-1) - nrecj
268	END DO
269	END DO
270	ENDIF
271
272	! 3. Subdomain description in the Regular Case
273	! --------------------------------------------
274	! specific cases where there is no communication -> must do the periodicity by itself
275	! Warning: because of potential land-area suppression, do not use nbond[ij] == 2
276	l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
277	l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
278
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	ibondi(ii,ij) = 0 ! default: has e-w neighbours
287	IF( ii == 1 ) ibondi(ii,ij) = -1 ! first column, has only e neighbour
288	IF( ii == jpni ) ibondi(ii,ij) = 1 ! last column, has only w neighbour
289	IF( jpni == 1 ) ibondi(ii,ij) = 2 ! has no e-w neighbour
290	ibondj(ii,ij) = 0 ! default: has n-s neighbours
291	IF( ij == 1 ) ibondj(ii,ij) = -1 ! first row, has only n neighbour
292	IF( ij == jpnj ) ibondj(ii,ij) = 1 ! last row, has only s neighbour
293	IF( jpnj == 1 ) ibondj(ii,ij) = 2 ! has no n-s neighbour
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 ) ibondi(ii,ij) = 0 ! redefine: all have e-w neighbours
308	IF( ii == 1 ) iowe(ii,ij) = iarea0 + (jpni-1) ! redefine: first column, address of w neighbour
309	IF( ii == jpni ) ioea(ii,ij) = iarea0 - (jpni-1) ! redefine: last column, address of e neighbour
310	ENDIF
311
312	! Simple North-South periodicity: change ibondj, ioso, iono
313	IF( jperio == 2 .OR. jperio == 7 ) THEN
314	IF( jpnj /= 1 ) ibondj(ii,ij) = 0 ! redefine: all have n-s neighbours
315	IF( ij == 1 ) ioso(ii,ij) = iarea0 + jpni * (jpnj-1) ! redefine: first row, address of s neighbour
316	IF( ij == jpnj ) iono(ii,ij) = iarea0 - jpni * (jpnj-1) ! redefine: last row, address of n neighbour
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	ii = 1 + MOD( jarea-1 , jpni )
392	ij = 1 + (jarea-1) / jpni
393	! land-only area with an active n neigbour
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 ) ! ii index of this n neigbour
396	ijno = 1 + iono(ii,ij) / jpni ! ij index of this n neigbour
397	! In case of north fold exchange: I am the n neigbour of my n neigbour!! (#1057)
398	! --> for northern neighbours of northern row processors (in case of north-fold)
399	! need to reverse the LOGICAL direction of communication
400	idir = 1 ! we are indeed the s neigbour of this n neigbour
401	IF( ij == jpnj .AND. ijno == jpnj ) idir = -1 ! both are on the last row, we are in fact the n neigbour
402	IF( ibondj(iino,ijno) == idir ) ibondj(iino,ijno) = 2 ! this n neigbour had only a s/n neigbour -> no more
403	IF( ibondj(iino,ijno) == 0 ) ibondj(iino,ijno) = -idir ! this n neigbour had both, n-s neighbours -> keep 1
404	ENDIF
405	! land-only area with an active s neigbour
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 ) ! ii index of this s neigbour
408	ijso = 1 + ioso(ii,ij) / jpni ! ij index of this s neigbour
409	IF( ibondj(iiso,ijso) == -1 ) ibondj(iiso,ijso) = 2 ! this s neigbour had only a n neigbour -> no more neigbour
410	IF( ibondj(iiso,ijso) == 0 ) ibondj(iiso,ijso) = 1 ! this s neigbour had both, n-s neighbours -> keep s neigbour
411	ENDIF
412	! land-only area with an active e neigbour
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 ) ! ii index of this e neigbour
415	ijea = 1 + ioea(ii,ij) / jpni ! ij index of this e neigbour
416	IF( ibondi(iiea,ijea) == 1 ) ibondi(iiea,ijea) = 2 ! this e neigbour had only a w neigbour -> no more neigbour
417	IF( ibondi(iiea,ijea) == 0 ) ibondi(iiea,ijea) = -1 ! this e neigbour had both, e-w neighbours -> keep e neigbour
418	ENDIF
419	! land-only area with an active w neigbour
420	IF( ipproc(ii,ij) == -1 .AND. 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= jpni*jpnj-1) THEN
421	iiwe = 1 + MOD( iowe(ii,ij) , jpni ) ! ii index of this w neigbour
422	ijwe = 1 + iowe(ii,ij) / jpni ! ij index of this w neigbour
423	IF( ibondi(iiwe,ijwe) == -1 ) ibondi(iiwe,ijwe) = 2 ! this w neigbour had only a e neigbour -> no more neigbour
424	IF( ibondi(iiwe,ijwe) == 0 ) ibondi(iiwe,ijwe) = 1 ! this w neigbour had both, e-w neighbours -> keep w neigbour
425	ENDIF
426	END DO
427
428	! Update il[de][ij] according to modified ibond[ij]
429	! ----------------------
430	DO jproc = 1, jpnij
431	ii = iin(jproc)
432	ij = ijn(jproc)
433	IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) = 1
434	IF( ibondi(ii,ij) == 1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij)
435	IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) = 1
436	IF( ibondj(ii,ij) == 1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij)
437	END DO
438
439	! just to save nono etc for all proc
440	! warning ii*ij (zone) /= nproc (processors)!
441	! ioso = zone number, ii_noso = proc number
442	ii_noso(:) = -1
443	ii_nono(:) = -1
444	ii_noea(:) = -1
445	ii_nowe(:) = -1
446	DO jproc = 1, jpnij
447	ii = iin(jproc)
448	ij = ijn(jproc)
449	IF( 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= (jpni*jpnj-1) ) THEN
450	iiso = 1 + MOD( ioso(ii,ij) , jpni )
451	ijso = 1 + ioso(ii,ij) / jpni
452	ii_noso(jproc) = ipproc(iiso,ijso)
453	ENDIF
454	IF( 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= (jpni*jpnj-1) ) THEN
455	iiwe = 1 + MOD( iowe(ii,ij) , jpni )
456	ijwe = 1 + iowe(ii,ij) / jpni
457	ii_nowe(jproc) = ipproc(iiwe,ijwe)
458	ENDIF
459	IF( 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= (jpni*jpnj-1) ) THEN
460	iiea = 1 + MOD( ioea(ii,ij) , jpni )
461	ijea = 1 + ioea(ii,ij) / jpni
462	ii_noea(jproc)= ipproc(iiea,ijea)
463	ENDIF
464	IF( 0 <= iono(ii,ij) .AND. iono(ii,ij) <= (jpni*jpnj-1) ) THEN
465	iino = 1 + MOD( iono(ii,ij) , jpni )
466	ijno = 1 + iono(ii,ij) / jpni
467	ii_nono(jproc)= ipproc(iino,ijno)
468	ENDIF
469	END DO
470
471	! 6. Change processor name
472	! ------------------------
473	ii = iin(narea)
474	ij = ijn(narea)
475	!
476	! set default neighbours
477	noso = ii_noso(narea)
478	nowe = ii_nowe(narea)
479	noea = ii_noea(narea)
480	nono = ii_nono(narea)
481	nlci = ilci(ii,ij)
482	nldi = ildi(ii,ij)
483	nlei = ilei(ii,ij)
484	nlcj = ilcj(ii,ij)
485	nldj = ildj(ii,ij)
486	nlej = ilej(ii,ij)
487	nbondi = ibondi(ii,ij)
488	nbondj = ibondj(ii,ij)
489	nimpp = iimppt(ii,ij)
490	njmpp = ijmppt(ii,ij)
491	jpi = nlci
492	jpj = nlcj
493	jpk = jpkglo ! third dim
494	#if defined key_agrif
495	! simple trick to use same vertical grid as parent but different number of levels:
496	! Save maximum number of levels in jpkglo, then define all vertical grids with this number.
497	! Suppress once vertical online interpolation is ok
498	!!$ IF(.NOT.Agrif_Root()) jpkglo = Agrif_Parent( jpkglo )
499	#endif
500	jpim1 = jpi-1 ! inner domain indices
501	jpjm1 = jpj-1 ! " "
502	jpkm1 = MAX( 1, jpk-1 ) ! " "
503	jpij = jpi*jpj ! jpi x j
504	DO jproc = 1, jpnij
505	ii = iin(jproc)
506	ij = ijn(jproc)
507	nlcit(jproc) = ilci(ii,ij)
508	nldit(jproc) = ildi(ii,ij)
509	nleit(jproc) = ilei(ii,ij)
510	nlcjt(jproc) = ilcj(ii,ij)
511	nldjt(jproc) = ildj(ii,ij)
512	nlejt(jproc) = ilej(ii,ij)
513	ibonit(jproc) = ibondi(ii,ij)
514	ibonjt(jproc) = ibondj(ii,ij)
515	nimppt(jproc) = iimppt(ii,ij)
516	njmppt(jproc) = ijmppt(ii,ij)
517	nfilcit(ii,ij) = ilci(ii,ij)
518	END DO
519
520	! Save processor layout in ascii file
521	IF (lwp) THEN
522	CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
523	WRITE(inum,'(a)') ' jpnij jpimax jpjmax jpk jpiglo jpjglo'//&
524	& ' ( local: narea jpi jpj )'
525	WRITE(inum,'(6i8,a,3i8,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,&
526	& ' ( local: ',narea,jpi,jpj,' )'
527	WRITE(inum,'(a)') 'nproc nlci nlcj nldi nldj nlei nlej nimp njmp nono noso nowe noea nbondi nbondj '
528
529	DO jproc = 1, jpnij
530	WRITE(inum,'(13i5,2i7)') jproc-1, nlcit (jproc), nlcjt (jproc), &
531	& nldit (jproc), nldjt (jproc), &
532	& nleit (jproc), nlejt (jproc), &
533	& nimppt (jproc), njmppt (jproc), &
534	& ii_nono(jproc), ii_noso(jproc), &
535	& ii_nowe(jproc), ii_noea(jproc), &
536	& ibonit (jproc), ibonjt (jproc)
537	END DO
538	CLOSE(inum)
539	END IF
540
541	! ! north fold parameter
542	! Defined npolj, either 0, 3 , 4 , 5 , 6
543	! In this case the important thing is that npolj /= 0
544	! Because if we go through these line it is because jpni >1 and thus
545	! we must use lbcnorthmpp, which tests only npolj =0 or npolj /= 0
546	npolj = 0
547	ij = ijn(narea)
548	IF( jperio == 3 .OR. jperio == 4 ) THEN
549	IF( ij == jpnj ) npolj = 3
550	ENDIF
551	IF( jperio == 5 .OR. jperio == 6 ) THEN
552	IF( ij == jpnj ) npolj = 5
553	ENDIF
554	!
555	nproc = narea-1
556	IF(lwp) THEN
557	WRITE(numout,*)
558	WRITE(numout,*) ' resulting internal parameters : '
559	WRITE(numout,*) ' nproc = ', nproc
560	WRITE(numout,*) ' nowe = ', nowe , ' noea = ', noea
561	WRITE(numout,*) ' nono = ', nono , ' noso = ', noso
562	WRITE(numout,*) ' nbondi = ', nbondi
563	WRITE(numout,*) ' nbondj = ', nbondj
564	WRITE(numout,*) ' npolj = ', npolj
565	WRITE(numout,*) ' l_Iperio = ', l_Iperio
566	WRITE(numout,*) ' l_Jperio = ', l_Jperio
567	WRITE(numout,*) ' nlci = ', nlci
568	WRITE(numout,*) ' nlcj = ', nlcj
569	WRITE(numout,*) ' nimpp = ', nimpp
570	WRITE(numout,*) ' njmpp = ', njmpp
571	WRITE(numout,*) ' nreci = ', nreci
572	WRITE(numout,*) ' nrecj = ', nrecj
573	WRITE(numout,*) ' nn_hls = ', nn_hls
574	ENDIF
575
576	! ! Prepare mpp north fold
577	IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
578	CALL mpp_ini_north
579	IF(lwp) WRITE(numout,*)
580	IF(lwp) WRITE(numout,*) ' ==>>> North fold boundary prepared for jpni >1'
581	ENDIF
582	!
583	CALL mpp_init_ioipsl ! Prepare NetCDF output file (if necessary)
584	!
585	IF( ln_nnogather ) CALL mpp_init_nfdcom ! northfold neighbour lists
586	!
587	DEALLOCATE(iin, ijn, ii_nono, ii_noea, ii_noso, ii_nowe, &
588	& iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj, &
589	& ilci, ilcj, ilei, ilej, ildi, ildj, &
590	& iono, ioea, ioso, iowe)
591	!
592	END SUBROUTINE mpp_init
593
594
595	SUBROUTINE mpp_init_mask( kmask )
596	!!----------------------------------------------------------------------
597	!! * ROUTINE mpp_init_mask *
598	!!
599	!! ** Purpose : Read relevant bathymetric information in a global array
600	!! in order to provide a land/sea mask used for the elimination
601	!! of land domains, in an mpp computation.
602	!!
603	!! ** Method : Read the namelist ln_zco and ln_isfcav in namelist namzgr
604	!! in order to choose the correct bathymetric information
605	!! (file and variables)
606	!!----------------------------------------------------------------------
607	INTEGER, DIMENSION(jpiglo,jpjglo), INTENT(out) :: kmask ! global domain
608
609	INTEGER :: inum !: logical unit for configuration file
610	INTEGER :: ios !: iostat error flag
611	INTEGER :: ijstartrow ! temporary integers
612	REAL(wp), DIMENSION(jpiglo,jpjglo) :: zbot, zbdy ! global workspace
613	REAL(wp) :: zidom , zjdom ! local scalars
614	NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file, &
615	& ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, &
616	& cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, &
617	& ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
618	& cn_ice, nn_ice_dta, &
619	& rn_ice_tem, rn_ice_sal, rn_ice_age, &
620	& ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
621	!!----------------------------------------------------------------------
622	! 0. initialisation
623	! -----------------
624	CALL iom_open( cn_domcfg, inum )
625	!
626	! ocean bottom level
627	CALL iom_get( inum, jpdom_unknown, 'bottom_level' , zbot , lrowattr=ln_use_jattr ) ! nb of ocean T-points
628	!
629	CALL iom_close( inum )
630	!
631	! 2D ocean mask (=1 if at least one level of the water column is ocean, =0 otherwise)
632	WHERE( zbot(:,:) > 0 ) ; kmask(:,:) = 1
633	ELSEWHERE ; kmask(:,:) = 0
634	END WHERE
635
636	! Adjust kmask with bdy_msk if it exists
637
638	REWIND( numnam_ref ) ! Namelist nambdy in reference namelist : BDY
639	READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
640	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp )
641	!
642	REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist : BDY
643	READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
644	904 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp )
645
646	IF( ln_bdy .AND. ln_mask_file ) THEN
647	CALL iom_open( cn_mask_file, inum )
648	CALL iom_get ( inum, jpdom_unknown, 'bdy_msk', zbdy )
649	CALL iom_close( inum )
650	WHERE ( zbdy(:,:) <= 0. ) kmask = 0
651	ENDIF
652	!
653	END SUBROUTINE mpp_init_mask
654
655
656	SUBROUTINE mpp_init_ioipsl
657	!!----------------------------------------------------------------------
658	!! * ROUTINE mpp_init_ioipsl *
659	!!
660	!! ** Purpose :
661	!!
662	!! ** Method :
663	!!
664	!! History :
665	!! 9.0 ! 04-03 (G. Madec ) MPP-IOIPSL
666	!! " " ! 08-12 (A. Coward) addition in case of jpni*jpnj < jpnij
667	!!----------------------------------------------------------------------
668	INTEGER, DIMENSION(2) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
669	!!----------------------------------------------------------------------
670
671	! The domain is split only horizontally along i- or/and j- direction
672	! So we need at the most only 1D arrays with 2 elements.
673	! Set idompar values equivalent to the jpdom_local_noextra definition
674	! used in IOM. This works even if jpnij .ne. jpni*jpnj.
675	iglo(1) = jpiglo
676	iglo(2) = jpjglo
677	iloc(1) = nlci
678	iloc(2) = nlcj
679	iabsf(1) = nimppt(narea)
680	iabsf(2) = njmppt(narea)
681	iabsl(:) = iabsf(:) + iloc(:) - 1
682	ihals(1) = nldi - 1
683	ihals(2) = nldj - 1
684	ihale(1) = nlci - nlei
685	ihale(2) = nlcj - nlej
686	idid(1) = 1
687	idid(2) = 2
688
689	IF(lwp) THEN
690	WRITE(numout,*)
691	WRITE(numout,*) 'mpp_init_ioipsl : iloc = ', iloc (1), iloc (2)
692	WRITE(numout,*) '~~~~~~~~~~~~~~~ iabsf = ', iabsf(1), iabsf(2)
693	WRITE(numout,*) ' ihals = ', ihals(1), ihals(2)
694	WRITE(numout,*) ' ihale = ', ihale(1), ihale(2)
695	ENDIF
696	!
697	CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
698	!
699	END SUBROUTINE mpp_init_ioipsl
700
701
702	SUBROUTINE mpp_init_partition( num_pes )
703	!!----------------------------------------------------------------------
704	!! * ROUTINE mpp_init_partition *
705	!!
706	!! ** Purpose :
707	!!
708	!! ** Method :
709	!!----------------------------------------------------------------------
710	INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have
711	!
712	INTEGER, PARAMETER :: nfactmax = 20
713	INTEGER :: nfact ! The no. of factors returned
714	INTEGER :: ierr ! Error flag
715	INTEGER :: ji
716	INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value
717	INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors
718	!!----------------------------------------------------------------------
719	!
720	ierr = 0
721	!
722	CALL factorise( ifact, nfactmax, nfact, num_pes, ierr )
723	!
724	IF( nfact <= 1 ) THEN
725	WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
726	WRITE (numout, *) ' : using grid of ',num_pes,' x 1'
727	jpnj = 1
728	jpni = num_pes
729	ELSE
730	! Search through factors for the pair that are closest in value
731	mindiff = 1000000
732	imin = 1
733	DO ji = 1, nfact-1, 2
734	idiff = ABS( ifact(ji) - ifact(ji+1) )
735	IF( idiff < mindiff ) THEN
736	mindiff = idiff
737	imin = ji
738	ENDIF
739	END DO
740	jpnj = ifact(imin)
741	jpni = ifact(imin + 1)
742	ENDIF
743	!
744	jpnij = jpni*jpnj
745	!
746	END SUBROUTINE mpp_init_partition
747
748
749	SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr )
750	!!----------------------------------------------------------------------
751	!! * ROUTINE factorise *
752	!!
753	!! ** Purpose : return the prime factors of n.
754	!! knfax factors are returned in array kfax which is of
755	!! maximum dimension kmaxfax.
756	!! ** Method :
757	!!----------------------------------------------------------------------
758	INTEGER , INTENT(in ) :: kn, kmaxfax
759	INTEGER , INTENT( out) :: kerr, knfax
760	INTEGER, DIMENSION(kmaxfax), INTENT( out) :: kfax
761	!
762	INTEGER :: ifac, jl, inu
763	INTEGER, PARAMETER :: ntest = 14
764	INTEGER, DIMENSION(ntest) :: ilfax
765	!!----------------------------------------------------------------------
766	!
767	! lfax contains the set of allowed factors.
768	ilfax(:) = (/(2**jl,jl=ntest,1,-1)/)
769	!
770	! Clear the error flag and initialise output vars
771	kerr = 0
772	kfax = 1
773	knfax = 0
774	!
775	IF( kn /= 1 ) THEN ! Find the factors of n
776	!
777	! nu holds the unfactorised part of the number.
778	! knfax holds the number of factors found.
779	! l points to the allowed factor list.
780	! ifac holds the current factor.
781	!
782	inu = kn
783	knfax = 0
784	!
785	DO jl = ntest, 1, -1
786	!
787	ifac = ilfax(jl)
788	IF( ifac > inu ) CYCLE
789	!
790	! Test whether the factor will divide.
791	!
792	IF( MOD(inu,ifac) == 0 ) THEN
793	!
794	knfax = knfax + 1 ! Add the factor to the list
795	IF( knfax > kmaxfax ) THEN
796	kerr = 6
797	write (,) 'FACTOR: insufficient space in factor array ', knfax
798	return
799	ENDIF
800	kfax(knfax) = ifac
801	! Store the other factor that goes with this one
802	knfax = knfax + 1
803	kfax(knfax) = inu / ifac
804	!WRITE (,) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax)
805	ENDIF
806	!
807	END DO
808	!
809	ENDIF
810	!
811	END SUBROUTINE factorise
812
813
814	SUBROUTINE mpp_init_nfdcom
815	!!----------------------------------------------------------------------
816	!! * ROUTINE mpp_init_nfdcom *
817	!! ** Purpose : Setup for north fold exchanges with explicit
818	!! point-to-point messaging
819	!!
820	!! ** Method : Initialization of the northern neighbours lists.
821	!!----------------------------------------------------------------------
822	!! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE)
823	!! 2.0 ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
824	!!----------------------------------------------------------------------
825	INTEGER :: sxM, dxM, sxT, dxT, jn
826	INTEGER :: njmppmax
827	!!----------------------------------------------------------------------
828	!
829	njmppmax = MAXVAL( njmppt )
830	!
831	!initializes the north-fold communication variables
832	isendto(:) = 0
833	nsndto = 0
834	!
835	IF ( njmpp == njmppmax ) THEN ! if I am a process in the north
836	!
837	!sxM is the first point (in the global domain) needed to compute the north-fold for the current process
838	sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
839	!dxM is the last point (in the global domain) needed to compute the north-fold for the current process
840	dxM = jpiglo - nimppt(narea) + 2
841	!
842	! loop over the other north-fold processes to find the processes
843	! managing the points belonging to the sxT-dxT range
844	!
845	DO jn = 1, jpni
846	!
847	sxT = nfiimpp(jn, jpnj) ! sxT = 1st point (in the global domain) of the jn process
848	dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1 ! dxT = last point (in the global domain) of the jn process
849	!
850	IF ( sxT < sxM .AND. sxM < dxT ) THEN
851	nsndto = nsndto + 1
852	isendto(nsndto) = jn
853	ELSEIF( sxM <= sxT .AND. dxM >= dxT ) THEN
854	nsndto = nsndto + 1
855	isendto(nsndto) = jn
856	ELSEIF( dxM < dxT .AND. sxT < dxM ) THEN
857	nsndto = nsndto + 1
858	isendto(nsndto) = jn
859	ENDIF
860	!
861	END DO
862	nfsloop = 1
863	nfeloop = nlci
864	DO jn = 2,jpni-1
865	IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
866	IF( nfipproc(jn-1,jpnj) == -1 ) nfsloop = nldi
867	IF( nfipproc(jn+1,jpnj) == -1 ) nfeloop = nlei
868	ENDIF
869	END DO
870	!
871	ENDIF
872	l_north_nogather = .TRUE.
873	!
874	END SUBROUTINE mpp_init_nfdcom
875
876
877	#endif
878
879	!!======================================================================
880	END MODULE mppini

Note: See TracBrowser for help on using the repository browser.

Download in other formats: