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 utils/tools/DOMAINcfg/src – NEMO

Context Navigation

source: utils/tools/DOMAINcfg/src/mppini.F90 @ 13204

Last change on this file since 13204 was 13204, checked in by smasson, 4 years ago
tools: update with tools_dev_r12970_AGRIF_CMEMS
File size: 55.2 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_partition: Calculate MPP domain decomposition
20	!! factorise : Calculate the factors of the no. of MPI processes
21	!! mpp_init_nfdcom : Setup for north fold exchanges with explicit point-to-point messaging
22	!!----------------------------------------------------------------------
23	USE dom_oce ! ocean space and time domain
24	!
25	USE lbcnfd , ONLY : isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges
26	USE lib_mpp ! distribued memory computing library
27	USE iom ! nemo I/O library
28	USE ioipsl ! I/O IPSL library
29	USE in_out_manager ! I/O Manager
30
31	IMPLICIT NONE
32	PRIVATE
33
34	PUBLIC mpp_init ! called by opa.F90
35
36	INTEGER :: numbot = -1 ! 'bottom_level' local logical unit
37
38	!!----------------------------------------------------------------------
39	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
40	!! $Id: mppini.F90 10570 2019-01-24 15:14:49Z acc $
41	!! Software governed by the CeCILL license (see ./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 = jpkglo
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	npolj = jperio
83	l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
84	l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
85	!
86	IF(lwp) THEN
87	WRITE(numout,*)
88	WRITE(numout,*) 'mpp_init : NO massively parallel processing'
89	WRITE(numout,*) '~~~~~~~~ '
90	WRITE(numout,*) ' l_Iperio = ', l_Iperio, ' l_Jperio = ', l_Jperio
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 defined key_agrif
99	IF (.not.agrif_root()) THEN
100	CALL agrif_nemo_init
101	ENDIF
102
103	IF( .NOT. Agrif_Root() ) THEN ! AGRIF children: specific setting (cf. agrif_user.F90)
104	print *,'nbcellsx = ',nbcellsx,nbghostcells_x
105	print *,'nbcellsy = ',nbcellsy,nbghostcells_y_s,nbghostcells_y_n
106	IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells_x ) THEN
107	IF(lwp) THEN
108	WRITE(numout,*)
109	WRITE(numout,) 'jpiglo should be: ', nbcellsx + 2 + 2nbghostcells_x
110	ENDIF
111	CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells_x' )
112	ENDIF
113	IF( jpjglo /= nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n ) THEN
114	IF(lwp) THEN
115	WRITE(numout,*)
116	WRITE(numout,*) 'jpjglo shoud be: ', nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n
117	ENDIF
118	CALL ctl_stop( 'STOP', &
119	'mpp_init: Agrif children requires jpjglo == nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n' )
120	ENDIF
121	IF( ln_use_jattr ) CALL ctl_stop( 'STOP', 'mpp_init:Agrif children requires ln_use_jattr = .false. ' )
122	ENDIF
123	#endif
124	!
125	END SUBROUTINE mpp_init
126
127	#else
128	!!----------------------------------------------------------------------
129	!! 'key_mpp_mpi' MPI massively parallel processing
130	!!----------------------------------------------------------------------
131
132
133	SUBROUTINE mpp_init
134	!!----------------------------------------------------------------------
135	!! * ROUTINE mpp_init *
136	!!
137	!! ** Purpose : Lay out the global domain over processors.
138	!! If land processors are to be eliminated, this program requires the
139	!! presence of the domain configuration file. Land processors elimination
140	!! is performed if jpni x jpnj /= jpnij. In this case, using the MPP_PREP
141	!! preprocessing tool, help for defining the best cutting out.
142	!!
143	!! ** Method : Global domain is distributed in smaller local domains.
144	!! Periodic condition is a function of the local domain position
145	!! (global boundary or neighbouring domain) and of the global
146	!! periodic
147	!! Type : jperio global periodic condition
148	!!
149	!! ** Action : - set domain parameters
150	!! nimpp : longitudinal index
151	!! njmpp : latitudinal index
152	!! narea : number for local area
153	!! nlci : first dimension
154	!! nlcj : second dimension
155	!! nbondi : mark for "east-west local boundary"
156	!! nbondj : mark for "north-south local boundary"
157	!! nproc : number for local processor
158	!! noea : number for local neighboring processor
159	!! nowe : number for local neighboring processor
160	!! noso : number for local neighboring processor
161	!! nono : number for local neighboring processor
162	!!----------------------------------------------------------------------
163	INTEGER :: ji, jj, jn, jproc, jarea ! dummy loop indices
164	INTEGER :: inijmin
165	INTEGER :: i2add
166	INTEGER :: inum ! local logical unit
167	INTEGER :: idir, ifreq, icont ! local integers
168	INTEGER :: ii, il1, ili, imil ! - -
169	INTEGER :: ij, il2, ilj, ijm1 ! - -
170	INTEGER :: iino, ijno, iiso, ijso ! - -
171	INTEGER :: iiea, ijea, iiwe, ijwe ! - -
172	INTEGER :: iarea0 ! - -
173	INTEGER :: ierr, ios !
174	INTEGER :: inbi, inbj, iimax, ijmax, icnt1, icnt2
175	LOGICAL :: llbest
176	LOGICAL :: llwrtlay
177	INTEGER, ALLOCATABLE, DIMENSION(:) :: iin, ii_nono, ii_noea ! 1D workspace
178	INTEGER, ALLOCATABLE, DIMENSION(:) :: ijn, ii_noso, ii_nowe ! - -
179	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ilci, ibondi, ipproc ! 2D workspace
180	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ilcj, ibondj, ipolj ! - -
181	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilei, ildi, iono, ioea ! - -
182	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilej, ildj, ioso, iowe ! - -
183	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: llisoce ! - -
184	!!----------------------------------------------------------------------
185
186	llwrtlay = lwp
187	!
188	IF( ln_read_cfg ) CALL iom_open( cn_domcfg, numbot )
189	!
190	! 1. Dimension arrays for subdomains
191	! -----------------------------------
192	!
193	! If dimensions of processor grid weren't specified in the namelist file
194	! then we calculate them here now that we have our communicator size
195	IF( jpni < 1 .OR. jpnj < 1 ) THEN
196	CALL mpp_init_bestpartition( mppsize, jpni, jpnj )
197	llbest = .TRUE.
198	ELSE
199	CALL mpp_init_bestpartition( mppsize, inbi, inbj, icnt2 )
200	CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax )
201	CALL mpp_basic_decomposition( inbi, inbj, iimax, ijmax )
202	IF( iimaxijmax < jpimaxjpjmax ) THEN
203	llbest = .FALSE.
204	icnt1 = jpni*jpnj - mppsize
205	WRITE(ctmp1,9000) ' The chosen domain decomposition ', jpni, ' x ', jpnj, ' with ', icnt1, ' land sub-domains'
206	WRITE(ctmp2,9000) ' has larger MPI subdomains (jpi = ', jpimax, ', jpj = ', jpjmax, ', jpijpj = ', jpimaxjpjmax, ')'
207	WRITE(ctmp3,9000) ' than the following domain decompostion ', inbi, ' x ', inbj, ' with ', icnt2, ' land sub-domains'
208	WRITE(ctmp4,9000) ' which MPI subdomains size is jpi = ', iimax, ', jpj = ', ijmax, ', jpijpj = ', iimaxijmax, ' '
209	CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4, ' ', ' --- YOU ARE WASTING CPU... ---', ' ' )
210	ELSE
211	llbest = .TRUE.
212	ENDIF
213	ENDIF
214
215	! look for land mpi subdomains...
216	ALLOCATE( llisoce(jpni,jpnj) )
217	CALL mpp_init_isoce( jpni, jpnj, llisoce )
218	inijmin = COUNT( llisoce ) ! number of oce subdomains
219
220	IF( mppsize < inijmin ) THEN
221	WRITE(ctmp1,9001) ' With this specified domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj
222	WRITE(ctmp2,9002) ' we can eliminate only ', jpni*jpnj - inijmin, ' land mpi subdomains therefore '
223	WRITE(ctmp3,9001) ' the number of ocean mpi subdomains (', inijmin,') exceed the number of MPI processes:', mppsize
224	WRITE(ctmp4,*) ' ==>>> There is the list of best domain decompositions you should use: '
225	CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 )
226	CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. ) ! must be done by all core
227	CALL ctl_stop( 'STOP' )
228	ENDIF
229
230	IF( mppsize > jpni*jpnj ) THEN
231	WRITE(ctmp1,9003) ' The number of mpi processes: ', mppsize
232	WRITE(ctmp2,9003) ' exceeds the maximum number of subdomains (ocean+land) = ', jpni*jpnj
233	WRITE(ctmp3,9001) ' defined by the following domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj
234	WRITE(ctmp4,*) ' ==>>> There is the list of best domain decompositions you should use: '
235	CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 )
236	CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. ) ! must be done by all core
237	CALL ctl_stop( 'STOP' )
238	ENDIF
239
240	jpnij = mppsize ! force jpnij definition <-- remove as much land subdomains as needed to reach this condition
241	IF( mppsize > inijmin ) THEN
242	WRITE(ctmp1,9003) ' The number of mpi processes: ', mppsize
243	WRITE(ctmp2,9003) ' exceeds the maximum number of ocean subdomains = ', inijmin
244	WRITE(ctmp3,9002) ' we suppressed ', jpni*jpnj - mppsize, ' land subdomains '
245	WRITE(ctmp4,9002) ' BUT we had to keep ', mppsize - inijmin, ' land subdomains that are useless...'
246	CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4, ' ', ' --- YOU ARE WASTING CPU... ---', ' ' )
247	ELSE ! mppsize = inijmin
248	IF(lwp) THEN
249	IF(llbest) WRITE(numout,*) 'mpp_init: You use an optimal domain decomposition'
250	WRITE(numout,*) '~~~~~~~~ '
251	WRITE(numout,9003) ' Number of mpi processes: ', mppsize
252	WRITE(numout,9003) ' Number of ocean subdomains = ', inijmin
253	WRITE(numout,9003) ' Number of suppressed land subdomains = ', jpni*jpnj - inijmin
254	WRITE(numout,*)
255	ENDIF
256	ENDIF
257	9000 FORMAT (a, i4, a, i4, a, i7, a)
258	9001 FORMAT (a, i4, a, i4)
259	9002 FORMAT (a, i4, a)
260	9003 FORMAT (a, i5)
261
262	IF( numbot /= -1 ) CALL iom_close( numbot )
263
264	ALLOCATE( nfiimpp(jpni,jpnj), nfipproc(jpni,jpnj), nfilcit(jpni,jpnj) , &
265	& nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) , &
266	& njmppt(jpnij) , ibonjt(jpnij) , nldit(jpnij) , nldjt(jpnij) , &
267	& nleit(jpnij) , nlejt(jpnij) , &
268	& iin(jpnij), ii_nono(jpnij), ii_noea(jpnij), &
269	& ijn(jpnij), ii_noso(jpnij), ii_nowe(jpnij), &
270	& iimppt(jpni,jpnj), ilci(jpni,jpnj), ibondi(jpni,jpnj), ipproc(jpni,jpnj), &
271	& ijmppt(jpni,jpnj), ilcj(jpni,jpnj), ibondj(jpni,jpnj), ipolj(jpni,jpnj), &
272	& ilei(jpni,jpnj), ildi(jpni,jpnj), iono(jpni,jpnj), ioea(jpni,jpnj), &
273	& ilej(jpni,jpnj), ildj(jpni,jpnj), ioso(jpni,jpnj), iowe(jpni,jpnj), &
274	& STAT=ierr )
275	CALL mpp_sum( 'mppini', ierr )
276	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' )
277
278	#if defined key_agrif
279	IF( .NOT. Agrif_Root() ) THEN ! AGRIF children: specific setting (cf. agrif_user.F90)
280	CALL agrif_nemo_init
281	IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells_x ) THEN
282	IF(lwp) THEN
283	WRITE(numout,*)
284	WRITE(numout,) 'jpiglo should be: ', nbcellsx + 2 + 2nbghostcells_x
285	ENDIF
286	CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells' )
287	ENDIF
288	IF( jpjglo /= nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n ) THEN
289	IF(lwp) THEN
290	WRITE(numout,*)
291	WRITE(numout,*) 'jpjglo shoud be: ', nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n
292	ENDIF
293	CALL ctl_stop( 'STOP', &
294	'mpp_init: Agrif children requires jpjglo == nbcellsy + 2 + nbghostcells_y_s + nbghostcells_y_n' )
295	ENDIF
296	IF( ln_use_jattr ) CALL ctl_stop( 'STOP', 'mpp_init:Agrif children requires ln_use_jattr = .false. ' )
297	ENDIF
298	#endif
299	!
300	! 2. Index arrays for subdomains
301	! -----------------------------------
302	!
303	nreci = 2 * nn_hls
304	nrecj = 2 * nn_hls
305	CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ilci, ilcj )
306	nfiimpp(:,:) = iimppt(:,:)
307	nfilcit(:,:) = ilci(:,:)
308	!
309	IF(lwp) THEN
310	WRITE(numout,*)
311	WRITE(numout,*) 'MPI Message Passing MPI - domain lay out over processors'
312	WRITE(numout,*)
313	WRITE(numout,*) ' defines mpp subdomains'
314	WRITE(numout,*) ' jpni = ', jpni
315	WRITE(numout,*) ' jpnj = ', jpnj
316	WRITE(numout,*)
317	WRITE(numout,*) ' sum ilci(i,1) = ', sum(ilci(:,1)), ' jpiglo = ', jpiglo
318	WRITE(numout,*) ' sum ilcj(1,j) = ', sum(ilcj(1,:)), ' jpjglo = ', jpjglo
319	ENDIF
320
321	! 3. Subdomain description in the Regular Case
322	! --------------------------------------------
323	! specific cases where there is no communication -> must do the periodicity by itself
324	! Warning: because of potential land-area suppression, do not use nbond[ij] == 2
325	l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
326	l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
327
328	DO jarea = 1, jpni*jpnj
329	!
330	iarea0 = jarea - 1
331	ii = 1 + MOD(iarea0,jpni)
332	ij = 1 + iarea0/jpni
333	ili = ilci(ii,ij)
334	ilj = ilcj(ii,ij)
335	ibondi(ii,ij) = 0 ! default: has e-w neighbours
336	IF( ii == 1 ) ibondi(ii,ij) = -1 ! first column, has only e neighbour
337	IF( ii == jpni ) ibondi(ii,ij) = 1 ! last column, has only w neighbour
338	IF( jpni == 1 ) ibondi(ii,ij) = 2 ! has no e-w neighbour
339	ibondj(ii,ij) = 0 ! default: has n-s neighbours
340	IF( ij == 1 ) ibondj(ii,ij) = -1 ! first row, has only n neighbour
341	IF( ij == jpnj ) ibondj(ii,ij) = 1 ! last row, has only s neighbour
342	IF( jpnj == 1 ) ibondj(ii,ij) = 2 ! has no n-s neighbour
343
344	! Subdomain neighbors (get their zone number): default definition
345	ioso(ii,ij) = iarea0 - jpni
346	iowe(ii,ij) = iarea0 - 1
347	ioea(ii,ij) = iarea0 + 1
348	iono(ii,ij) = iarea0 + jpni
349	ildi(ii,ij) = 1 + nn_hls
350	ilei(ii,ij) = ili - nn_hls
351	ildj(ii,ij) = 1 + nn_hls
352	ilej(ii,ij) = ilj - nn_hls
353
354	! East-West periodicity: change ibondi, ioea, iowe
355	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7 ) THEN
356	IF( jpni /= 1 ) ibondi(ii,ij) = 0 ! redefine: all have e-w neighbours
357	IF( ii == 1 ) iowe(ii,ij) = iarea0 + (jpni-1) ! redefine: first column, address of w neighbour
358	IF( ii == jpni ) ioea(ii,ij) = iarea0 - (jpni-1) ! redefine: last column, address of e neighbour
359	ENDIF
360
361	! Simple North-South periodicity: change ibondj, ioso, iono
362	IF( jperio == 2 .OR. jperio == 7 ) THEN
363	IF( jpnj /= 1 ) ibondj(ii,ij) = 0 ! redefine: all have n-s neighbours
364	IF( ij == 1 ) ioso(ii,ij) = iarea0 + jpni * (jpnj-1) ! redefine: first row, address of s neighbour
365	IF( ij == jpnj ) iono(ii,ij) = iarea0 - jpni * (jpnj-1) ! redefine: last row, address of n neighbour
366	ENDIF
367
368	! North fold: define ipolj, change iono. Warning: we do not change ibondj...
369	ipolj(ii,ij) = 0
370	IF( jperio == 3 .OR. jperio == 4 ) THEN
371	ijm1 = jpni*(jpnj-1)
372	imil = ijm1+(jpni+1)/2
373	IF( jarea > ijm1 ) ipolj(ii,ij) = 3
374	IF( MOD(jpni,2) == 1 .AND. jarea == imil ) ipolj(ii,ij) = 4
375	IF( ipolj(ii,ij) == 3 ) iono(ii,ij) = jpni*jpnj-jarea+ijm1 ! MPI rank of northern neighbour
376	ENDIF
377	IF( jperio == 5 .OR. jperio == 6 ) THEN
378	ijm1 = jpni*(jpnj-1)
379	imil = ijm1+(jpni+1)/2
380	IF( jarea > ijm1) ipolj(ii,ij) = 5
381	IF( MOD(jpni,2) == 1 .AND. jarea == imil ) ipolj(ii,ij) = 6
382	IF( ipolj(ii,ij) == 5) iono(ii,ij) = jpni*jpnj-jarea+ijm1 ! MPI rank of northern neighbour
383	ENDIF
384	!
385	END DO
386
387	! 4. deal with land subdomains
388	! ----------------------------
389	!
390	! specify which subdomains are oce subdomains; other are land subdomains
391	ipproc(:,:) = -1
392	icont = -1
393	DO jarea = 1, jpni*jpnj
394	iarea0 = jarea - 1
395	ii = 1 + MOD(iarea0,jpni)
396	ij = 1 + iarea0/jpni
397	IF( llisoce(ii,ij) ) THEN
398	icont = icont + 1
399	ipproc(ii,ij) = icont
400	iin(icont+1) = ii
401	ijn(icont+1) = ij
402	ENDIF
403	END DO
404	! if needed add some land subdomains to reach jpnij active subdomains
405	i2add = jpnij - inijmin
406	DO jarea = 1, jpni*jpnj
407	iarea0 = jarea - 1
408	ii = 1 + MOD(iarea0,jpni)
409	ij = 1 + iarea0/jpni
410	IF( .NOT. llisoce(ii,ij) .AND. i2add > 0 ) THEN
411	icont = icont + 1
412	ipproc(ii,ij) = icont
413	iin(icont+1) = ii
414	ijn(icont+1) = ij
415	i2add = i2add - 1
416	ENDIF
417	END DO
418	nfipproc(:,:) = ipproc(:,:)
419
420	! neighbour treatment: change ibondi, ibondj if next to a land zone
421	DO jarea = 1, jpni*jpnj
422	ii = 1 + MOD( jarea-1 , jpni )
423	ij = 1 + (jarea-1) / jpni
424	! land-only area with an active n neigbour
425	IF ( ipproc(ii,ij) == -1 .AND. 0 <= iono(ii,ij) .AND. iono(ii,ij) <= jpni*jpnj-1 ) THEN
426	iino = 1 + MOD( iono(ii,ij) , jpni ) ! ii index of this n neigbour
427	ijno = 1 + iono(ii,ij) / jpni ! ij index of this n neigbour
428	! In case of north fold exchange: I am the n neigbour of my n neigbour!! (#1057)
429	! --> for northern neighbours of northern row processors (in case of north-fold)
430	! need to reverse the LOGICAL direction of communication
431	idir = 1 ! we are indeed the s neigbour of this n neigbour
432	IF( ij == jpnj .AND. ijno == jpnj ) idir = -1 ! both are on the last row, we are in fact the n neigbour
433	IF( ibondj(iino,ijno) == idir ) ibondj(iino,ijno) = 2 ! this n neigbour had only a s/n neigbour -> no more
434	IF( ibondj(iino,ijno) == 0 ) ibondj(iino,ijno) = -idir ! this n neigbour had both, n-s neighbours -> keep 1
435	ENDIF
436	! land-only area with an active s neigbour
437	IF( ipproc(ii,ij) == -1 .AND. 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= jpni*jpnj-1 ) THEN
438	iiso = 1 + MOD( ioso(ii,ij) , jpni ) ! ii index of this s neigbour
439	ijso = 1 + ioso(ii,ij) / jpni ! ij index of this s neigbour
440	IF( ibondj(iiso,ijso) == -1 ) ibondj(iiso,ijso) = 2 ! this s neigbour had only a n neigbour -> no more neigbour
441	IF( ibondj(iiso,ijso) == 0 ) ibondj(iiso,ijso) = 1 ! this s neigbour had both, n-s neighbours -> keep s neigbour
442	ENDIF
443	! land-only area with an active e neigbour
444	IF( ipproc(ii,ij) == -1 .AND. 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= jpni*jpnj-1 ) THEN
445	iiea = 1 + MOD( ioea(ii,ij) , jpni ) ! ii index of this e neigbour
446	ijea = 1 + ioea(ii,ij) / jpni ! ij index of this e neigbour
447	IF( ibondi(iiea,ijea) == 1 ) ibondi(iiea,ijea) = 2 ! this e neigbour had only a w neigbour -> no more neigbour
448	IF( ibondi(iiea,ijea) == 0 ) ibondi(iiea,ijea) = -1 ! this e neigbour had both, e-w neighbours -> keep e neigbour
449	ENDIF
450	! land-only area with an active w neigbour
451	IF( ipproc(ii,ij) == -1 .AND. 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= jpni*jpnj-1) THEN
452	iiwe = 1 + MOD( iowe(ii,ij) , jpni ) ! ii index of this w neigbour
453	ijwe = 1 + iowe(ii,ij) / jpni ! ij index of this w neigbour
454	IF( ibondi(iiwe,ijwe) == -1 ) ibondi(iiwe,ijwe) = 2 ! this w neigbour had only a e neigbour -> no more neigbour
455	IF( ibondi(iiwe,ijwe) == 0 ) ibondi(iiwe,ijwe) = 1 ! this w neigbour had both, e-w neighbours -> keep w neigbour
456	ENDIF
457	END DO
458
459	! Update il[de][ij] according to modified ibond[ij]
460	! ----------------------
461	DO jproc = 1, jpnij
462	ii = iin(jproc)
463	ij = ijn(jproc)
464	IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) = 1
465	IF( ibondi(ii,ij) == 1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij)
466	IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) = 1
467	IF( ibondj(ii,ij) == 1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij)
468	END DO
469
470	! 5. Subdomain print
471	! ------------------
472	IF(lwp) THEN
473	ifreq = 4
474	il1 = 1
475	DO jn = 1, (jpni-1)/ifreq+1
476	il2 = MIN(jpni,il1+ifreq-1)
477	WRITE(numout,*)
478	WRITE(numout,9400) ('***',ji=il1,il2-1)
479	DO jj = jpnj, 1, -1
480	WRITE(numout,9403) (' ',ji=il1,il2-1)
481	WRITE(numout,9402) jj, (ilci(ji,jj),ilcj(ji,jj),ji=il1,il2)
482	WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2)
483	WRITE(numout,9403) (' ',ji=il1,il2-1)
484	WRITE(numout,9400) ('***',ji=il1,il2-1)
485	END DO
486	WRITE(numout,9401) (ji,ji=il1,il2)
487	il1 = il1+ifreq
488	END DO
489	9400 FORMAT(' *' ,20('***********',a3) )
490	9403 FORMAT(' * ',20(' * ',a3) )
491	9401 FORMAT(' ' ,20(' ',i3,' ') )
492	9402 FORMAT(' ',i3,' * ',20(i3,' x',i3,' * ') )
493	9404 FORMAT(' * ' ,20(' ',i3,' * ') )
494	ENDIF
495
496	! just to save nono etc for all proc
497	! warning ii*ij (zone) /= nproc (processors)!
498	! ioso = zone number, ii_noso = proc number
499	ii_noso(:) = -1
500	ii_nono(:) = -1
501	ii_noea(:) = -1
502	ii_nowe(:) = -1
503	DO jproc = 1, jpnij
504	ii = iin(jproc)
505	ij = ijn(jproc)
506	IF( 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= (jpni*jpnj-1) ) THEN
507	iiso = 1 + MOD( ioso(ii,ij) , jpni )
508	ijso = 1 + ioso(ii,ij) / jpni
509	ii_noso(jproc) = ipproc(iiso,ijso)
510	ENDIF
511	IF( 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= (jpni*jpnj-1) ) THEN
512	iiwe = 1 + MOD( iowe(ii,ij) , jpni )
513	ijwe = 1 + iowe(ii,ij) / jpni
514	ii_nowe(jproc) = ipproc(iiwe,ijwe)
515	ENDIF
516	IF( 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= (jpni*jpnj-1) ) THEN
517	iiea = 1 + MOD( ioea(ii,ij) , jpni )
518	ijea = 1 + ioea(ii,ij) / jpni
519	ii_noea(jproc)= ipproc(iiea,ijea)
520	ENDIF
521	IF( 0 <= iono(ii,ij) .AND. iono(ii,ij) <= (jpni*jpnj-1) ) THEN
522	iino = 1 + MOD( iono(ii,ij) , jpni )
523	ijno = 1 + iono(ii,ij) / jpni
524	ii_nono(jproc)= ipproc(iino,ijno)
525	ENDIF
526	END DO
527
528	! 6. Change processor name
529	! ------------------------
530	ii = iin(narea)
531	ij = ijn(narea)
532	!
533	! set default neighbours
534	noso = ii_noso(narea)
535	nowe = ii_nowe(narea)
536	noea = ii_noea(narea)
537	nono = ii_nono(narea)
538	nlci = ilci(ii,ij)
539	nldi = ildi(ii,ij)
540	nlei = ilei(ii,ij)
541	nlcj = ilcj(ii,ij)
542	nldj = ildj(ii,ij)
543	nlej = ilej(ii,ij)
544	nbondi = ibondi(ii,ij)
545	nbondj = ibondj(ii,ij)
546	nimpp = iimppt(ii,ij)
547	njmpp = ijmppt(ii,ij)
548	jpi = nlci
549	jpj = nlcj
550	jpk = jpkglo ! third dim
551	#if defined key_agrif
552	! simple trick to use same vertical grid as parent but different number of levels:
553	! Save maximum number of levels in jpkglo, then define all vertical grids with this number.
554	! Suppress once vertical online interpolation is ok
555	!!$ IF(.NOT.Agrif_Root()) jpkglo = Agrif_Parent( jpkglo )
556	#endif
557	jpim1 = jpi-1 ! inner domain indices
558	jpjm1 = jpj-1 ! " "
559	jpkm1 = MAX( 1, jpk-1 ) ! " "
560	jpij = jpi*jpj ! jpi x j
561	DO jproc = 1, jpnij
562	ii = iin(jproc)
563	ij = ijn(jproc)
564	nlcit(jproc) = ilci(ii,ij)
565	nldit(jproc) = ildi(ii,ij)
566	nleit(jproc) = ilei(ii,ij)
567	nlcjt(jproc) = ilcj(ii,ij)
568	nldjt(jproc) = ildj(ii,ij)
569	nlejt(jproc) = ilej(ii,ij)
570	ibonit(jproc) = ibondi(ii,ij)
571	ibonjt(jproc) = ibondj(ii,ij)
572	nimppt(jproc) = iimppt(ii,ij)
573	njmppt(jproc) = ijmppt(ii,ij)
574	END DO
575
576	! Save processor layout in ascii file
577	IF (llwrtlay) THEN
578	CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
579	WRITE(inum,'(a)') ' jpnij jpimax jpjmax jpk jpiglo jpjglo'//&
580	& ' ( local: narea jpi jpj )'
581	WRITE(inum,'(6i8,a,3i8,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,&
582	& ' ( local: ',narea,jpi,jpj,' )'
583	WRITE(inum,'(a)') 'nproc nlci nlcj nldi nldj nlei nlej nimp njmp nono noso nowe noea nbondi nbondj '
584
585	DO jproc = 1, jpnij
586	WRITE(inum,'(13i5,2i7)') jproc-1, nlcit (jproc), nlcjt (jproc), &
587	& nldit (jproc), nldjt (jproc), &
588	& nleit (jproc), nlejt (jproc), &
589	& nimppt (jproc), njmppt (jproc), &
590	& ii_nono(jproc), ii_noso(jproc), &
591	& ii_nowe(jproc), ii_noea(jproc), &
592	& ibonit (jproc), ibonjt (jproc)
593	END DO
594	END IF
595
596	! ! north fold parameter
597	! Defined npolj, either 0, 3 , 4 , 5 , 6
598	! In this case the important thing is that npolj /= 0
599	! Because if we go through these line it is because jpni >1 and thus
600	! we must use lbcnorthmpp, which tests only npolj =0 or npolj /= 0
601	npolj = 0
602	ij = ijn(narea)
603	IF( jperio == 3 .OR. jperio == 4 ) THEN
604	IF( ij == jpnj ) npolj = 3
605	ENDIF
606	IF( jperio == 5 .OR. jperio == 6 ) THEN
607	IF( ij == jpnj ) npolj = 5
608	ENDIF
609	!
610	nproc = narea-1
611	IF(lwp) THEN
612	WRITE(numout,*)
613	WRITE(numout,*) ' resulting internal parameters : '
614	WRITE(numout,*) ' nproc = ', nproc
615	WRITE(numout,*) ' nowe = ', nowe , ' noea = ', noea
616	WRITE(numout,*) ' nono = ', nono , ' noso = ', noso
617	WRITE(numout,*) ' nbondi = ', nbondi
618	WRITE(numout,*) ' nbondj = ', nbondj
619	WRITE(numout,*) ' npolj = ', npolj
620	WRITE(numout,*) ' l_Iperio = ', l_Iperio
621	WRITE(numout,*) ' l_Jperio = ', l_Jperio
622	WRITE(numout,*) ' nlci = ', nlci
623	WRITE(numout,*) ' nlcj = ', nlcj
624	WRITE(numout,*) ' nimpp = ', nimpp
625	WRITE(numout,*) ' njmpp = ', njmpp
626	WRITE(numout,*) ' nreci = ', nreci
627	WRITE(numout,*) ' nrecj = ', nrecj
628	WRITE(numout,*) ' nn_hls = ', nn_hls
629	ENDIF
630
631	! ! Prepare mpp north fold
632	IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
633	CALL mpp_ini_north
634	IF (lwp) THEN
635	WRITE(numout,*)
636	WRITE(numout,*) ' ==>>> North fold boundary prepared for jpni >1'
637	! additional prints in layout.dat
638	ENDIF
639	IF (llwrtlay) THEN
640	WRITE(inum,*)
641	WRITE(inum,*)
642	WRITE(inum,*) 'number of subdomains located along the north fold : ', ndim_rank_north
643	WRITE(inum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
644	DO jproc = 1, ndim_rank_north, 5
645	WRITE(inum,*) nrank_north( jproc:MINVAL( (/jproc+4,ndim_rank_north/) ) )
646	END DO
647	ENDIF
648	ENDIF
649	!
650	IF( ln_nnogather ) THEN
651	CALL mpp_init_nfdcom ! northfold neighbour lists
652	IF (llwrtlay) THEN
653	WRITE(inum,*)
654	WRITE(inum,*)
655	WRITE(inum,*) 'north fold exchanges with explicit point-to-point messaging :'
656	WRITE(inum,*) 'nfsloop : ', nfsloop
657	WRITE(inum,*) 'nfeloop : ', nfeloop
658	WRITE(inum,*) 'nsndto : ', nsndto
659	WRITE(inum,*) 'isendto : ', isendto
660	ENDIF
661	ENDIF
662	!
663	IF (llwrtlay) CLOSE(inum)
664	!
665	DEALLOCATE(iin, ijn, ii_nono, ii_noea, ii_noso, ii_nowe, &
666	& iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj, &
667	& ilci, ilcj, ilei, ilej, ildi, ildj, &
668	& iono, ioea, ioso, iowe, llisoce)
669	!
670	END SUBROUTINE mpp_init
671
672
673	SUBROUTINE mpp_basic_decomposition( knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj)
674	!!----------------------------------------------------------------------
675	!! * ROUTINE mpp_basic_decomposition *
676	!!
677	!! ** Purpose : Lay out the global domain over processors.
678	!!
679	!! ** Method : Global domain is distributed in smaller local domains.
680	!!
681	!! ** Action : - set for all knbi*knbj domains:
682	!! kimppt : longitudinal index
683	!! kjmppt : latitudinal index
684	!! klci : first dimension
685	!! klcj : second dimension
686	!!----------------------------------------------------------------------
687	INTEGER, INTENT(in ) :: knbi, knbj
688	INTEGER, INTENT( out) :: kimax, kjmax
689	INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: kimppt, kjmppt
690	INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: klci, klcj
691	!
692	INTEGER :: ji, jj
693	INTEGER :: iresti, irestj, irm, ijpjmin
694	INTEGER :: ireci, irecj
695	!!----------------------------------------------------------------------
696	!
697	#if defined key_nemocice_decomp
698	kimax = ( nx_global+2-2nn_hls + (knbi-1) ) / knbi + 2nn_hls ! first dim.
699	kjmax = ( ny_global+2-2nn_hls + (knbj-1) ) / knbj + 2nn_hls ! second dim.
700	#else
701	kimax = ( jpiglo - 2nn_hls + (knbi-1) ) / knbi + 2nn_hls ! first dim.
702	kjmax = ( jpjglo - 2nn_hls + (knbj-1) ) / knbj + 2nn_hls ! second dim.
703	#endif
704	IF( .NOT. PRESENT(kimppt) ) RETURN
705	!
706	! 1. Dimension arrays for subdomains
707	! -----------------------------------
708	! Computation of local domain sizes klci() klcj()
709	! These dimensions depend on global sizes knbi,knbj and jpiglo,jpjglo
710	! The subdomains are squares lesser than or equal to the global
711	! dimensions divided by the number of processors minus the overlap array.
712	!
713	ireci = 2 * nn_hls
714	irecj = 2 * nn_hls
715	iresti = 1 + MOD( jpiglo - ireci -1 , knbi )
716	irestj = 1 + MOD( jpjglo - irecj -1 , knbj )
717	!
718	! Need to use kimax and kjmax here since jpi and jpj not yet defined
719	#if defined key_nemocice_decomp
720	! Change padding to be consistent with CICE
721	klci(1:knbi-1 ,:) = kimax
722	klci(knbi ,:) = jpiglo - (knbi - 1) * (kimax - nreci)
723	klcj(:, 1:knbj-1) = kjmax
724	klcj(:, knbj) = jpjglo - (knbj - 1) * (kjmax - nrecj)
725	#else
726	klci(1:iresti ,:) = kimax
727	klci(iresti+1:knbi ,:) = kimax-1
728	IF( MINVAL(klci) < 3 ) THEN
729	WRITE(ctmp1,*) ' mpp_basic_decomposition: minimum value of jpi must be >= 3'
730	WRITE(ctmp2,*) ' We have ', MINVAL(klci)
731	CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
732	ENDIF
733	IF( jperio == 3 .OR. jperio == 4 .OR. jperio == 5 .OR. jperio == 6 ) THEN
734	! minimize the size of the last row to compensate for the north pole folding coast
735	IF( jperio == 3 .OR. jperio == 4 ) ijpjmin = 5 ! V and F folding involves line jpj-3 that must not be south boundary
736	IF( jperio == 5 .OR. jperio == 6 ) ijpjmin = 4 ! V and F folding involves line jpj-2 that must not be south boundary
737	irm = knbj - irestj ! total number of lines to be removed
738	klcj(:, knbj) = MAX( ijpjmin, kjmax-irm ) ! we must have jpj >= ijpjmin in the last row
739	irm = irm - ( kjmax - klcj(1,knbj) ) ! remaining number of lines to remove
740	irestj = knbj - 1 - irm
741	klcj(:, 1:irestj) = kjmax
742	klcj(:, irestj+1:knbj-1) = kjmax-1
743	ELSE
744	ijpjmin = 3
745	klcj(:, 1:irestj) = kjmax
746	klcj(:, irestj+1:knbj) = kjmax-1
747	ENDIF
748	IF( MINVAL(klcj) < ijpjmin ) THEN
749	WRITE(ctmp1,*) ' mpp_basic_decomposition: minimum value of jpj must be >= ', ijpjmin
750	WRITE(ctmp2,*) ' We have ', MINVAL(klcj)
751	CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
752	ENDIF
753	#endif
754
755	! 2. Index arrays for subdomains
756	! -------------------------------
757	kimppt(:,:) = 1
758	kjmppt(:,:) = 1
759	!
760	IF( knbi > 1 ) THEN
761	DO jj = 1, knbj
762	DO ji = 2, knbi
763	kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - ireci
764	END DO
765	END DO
766	ENDIF
767	!
768	IF( knbj > 1 )THEN
769	DO jj = 2, knbj
770	DO ji = 1, knbi
771	kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - irecj
772	END DO
773	END DO
774	ENDIF
775
776	END SUBROUTINE mpp_basic_decomposition
777
778
779	SUBROUTINE mpp_init_bestpartition( knbij, knbi, knbj, knbcnt, ldlist )
780	!!----------------------------------------------------------------------
781	!! * ROUTINE mpp_init_bestpartition *
782	!!
783	!! ** Purpose :
784	!!
785	!! ** Method :
786	!!----------------------------------------------------------------------
787	INTEGER, INTENT(in ) :: knbij ! total number if subdomains (knbi*knbj)
788	INTEGER, OPTIONAL, INTENT( out) :: knbi, knbj ! number if subdomains along i and j (knbi and knbj)
789	INTEGER, OPTIONAL, INTENT( out) :: knbcnt ! number of land subdomains
790	LOGICAL, OPTIONAL, INTENT(in ) :: ldlist ! .true.: print the list the best domain decompositions (with land)
791	!
792	INTEGER :: ji, jj, ii, iitarget
793	INTEGER :: iszitst, iszjtst
794	INTEGER :: isziref, iszjref
795	INTEGER :: inbij, iszij
796	INTEGER :: inbimax, inbjmax, inbijmax
797	INTEGER :: isz0, isz1
798	INTEGER, DIMENSION( :), ALLOCATABLE :: indexok
799	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi0, inbj0, inbij0 ! number of subdomains along i,j
800	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi0, iszj0, iszij0 ! max size of the subdomains along i,j
801	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi1, inbj1, inbij1 ! number of subdomains along i,j
802	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi1, iszj1, iszij1 ! max size of the subdomains along i,j
803	LOGICAL :: llist
804	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d ! max size of the subdomains along i,j
805	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisoce ! - -
806	REAL(wp):: zpropland
807	!!----------------------------------------------------------------------
808	!
809	llist = .FALSE.
810	IF( PRESENT(ldlist) ) llist = ldlist
811
812	CALL mpp_init_landprop( zpropland ) ! get the proportion of land point over the gloal domain
813	inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) ) ! define the largest possible value for jpni*jpnj
814	!
815	IF( llist ) THEN ; inbijmax = inbij*2
816	ELSE ; inbijmax = inbij
817	ENDIF
818	!
819	ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax))
820	!
821	inbimax = 0
822	inbjmax = 0
823	isziref = jpiglo*jpjglo+1
824	iszjref = jpiglo*jpjglo+1
825	!
826	! get the list of knbi that gives a smaller jpimax than knbi-1
827	! get the list of knbj that gives a smaller jpjmax than knbj-1
828	DO ji = 1, inbijmax
829	#if defined key_nemocice_decomp
830	iszitst = ( nx_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
831	#else
832	iszitst = ( jpiglo - 2nn_hls + (ji-1) ) / ji + 2nn_hls
833	#endif
834	IF( iszitst < isziref ) THEN
835	isziref = iszitst
836	inbimax = inbimax + 1
837	inbi0(inbimax) = ji
838	iszi0(inbimax) = isziref
839	ENDIF
840	#if defined key_nemocice_decomp
841	iszjtst = ( ny_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
842	#else
843	iszjtst = ( jpjglo - 2nn_hls + (ji-1) ) / ji + 2nn_hls
844	#endif
845	IF( iszjtst < iszjref ) THEN
846	iszjref = iszjtst
847	inbjmax = inbjmax + 1
848	inbj0(inbjmax) = ji
849	iszj0(inbjmax) = iszjref
850	ENDIF
851	END DO
852
853	! combine these 2 lists to get all possible knbi*knbj < inbijmax
854	ALLOCATE( llmsk2d(inbimax,inbjmax) )
855	DO jj = 1, inbjmax
856	DO ji = 1, inbimax
857	IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN ; llmsk2d(ji,jj) = .TRUE.
858	ELSE ; llmsk2d(ji,jj) = .FALSE.
859	ENDIF
860	END DO
861	END DO
862	isz1 = COUNT(llmsk2d)
863	ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) )
864	ii = 0
865	DO jj = 1, inbjmax
866	DO ji = 1, inbimax
867	IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN
868	ii = ii + 1
869	inbi1(ii) = inbi0(ji)
870	inbj1(ii) = inbj0(jj)
871	iszi1(ii) = iszi0(ji)
872	iszj1(ii) = iszj0(jj)
873	END IF
874	END DO
875	END DO
876	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
877	DEALLOCATE( llmsk2d )
878
879	ALLOCATE( inbij1(isz1), iszij1(isz1) )
880	inbij1(:) = inbi1(:) * inbj1(:)
881	iszij1(:) = iszi1(:) * iszj1(:)
882
883	! if therr is no land and no print
884	IF( .NOT. llist .AND. numbot == -1 ) THEN
885	! get the smaller partition which gives the smallest subdomain size
886	ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1)
887	knbi = inbi1(ii)
888	knbj = inbj1(ii)
889	IF(PRESENT(knbcnt)) knbcnt = 0
890	DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 )
891	RETURN
892	ENDIF
893
894	! extract only the partitions which reduce the subdomain size in comparison with smaller partitions
895	ALLOCATE( indexok(isz1) ) ! to store indices of the best partitions
896	isz0 = 0 ! number of best partitions
897	inbij = 1 ! start with the min value of inbij1 => 1
898	iszij = jpiglo*jpjglo+1 ! default: larger than global domain
899	DO WHILE( inbij <= inbijmax ) ! if we did not reach the max of inbij1
900	ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1) ! warning: send back the first occurence if multiple results
901	IF ( iszij1(ii) < iszij ) THEN
902	isz0 = isz0 + 1
903	indexok(isz0) = ii
904	iszij = iszij1(ii)
905	ENDIF
906	inbij = MINVAL(inbij1, mask = inbij1 > inbij) ! warning: return largest integer value if mask = .false. everywhere
907	END DO
908	DEALLOCATE( inbij1, iszij1 )
909
910	! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size)
911	ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) )
912	DO ji = 1, isz0
913	ii = indexok(ji)
914	inbi0(ji) = inbi1(ii)
915	inbj0(ji) = inbj1(ii)
916	iszi0(ji) = iszi1(ii)
917	iszj0(ji) = iszj1(ii)
918	END DO
919	DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 )
920
921	IF( llist ) THEN ! we print about 21 best partitions
922	IF(lwp) THEN
923	WRITE(numout,*)
924	WRITE(numout, *) ' For your information:'
925	WRITE(numout,'(a,i5,a)') ' list of the best partitions around ', knbij, ' mpi processes'
926	WRITE(numout, *) ' --------------------------------------', '-----', '--------------'
927	WRITE(numout,*)
928	END IF
929	iitarget = MINLOC( inbi0(:)inbj0(:), mask = inbi0(:)inbj0(:) >= knbij, dim = 1 )
930	DO ji = MAX(1,iitarget-10), MIN(isz0,iitarget+10)
931	ALLOCATE( llisoce(inbi0(ji), inbj0(ji)) )
932	CALL mpp_init_isoce( inbi0(ji), inbj0(ji), llisoce ) ! Warning: must be call by all cores (call mpp_sum)
933	inbij = COUNT(llisoce)
934	DEALLOCATE( llisoce )
935	IF(lwp) WRITE(numout,'(a, i5, a, i5, a, i4, a, i4, a, i9, a, i5, a, i5, a)') &
936	& 'nb_cores ' , inbij,' oce + ', inbi0(ji)*inbj0(ji) - inbij &
937	& , ' land ( ', inbi0(ji),' x ', inbj0(ji), &
938	& ' ), nb_points ', iszi0(ji)*iszj0(ji),' ( ', iszi0(ji),' x ', iszj0(ji),' )'
939	END DO
940	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
941	RETURN
942	ENDIF
943
944	DEALLOCATE( iszi0, iszj0 )
945	inbij = inbijmax + 1 ! default: larger than possible
946	ii = isz0+1 ! start from the end of the list (smaller subdomains)
947	DO WHILE( inbij > knbij ) ! while the number of ocean subdomains exceed the number of procs
948	ii = ii -1
949	ALLOCATE( llisoce(inbi0(ii), inbj0(ii)) )
950	CALL mpp_init_isoce( inbi0(ii), inbj0(ii), llisoce ) ! must be done by all core
951	inbij = COUNT(llisoce)
952	DEALLOCATE( llisoce )
953	END DO
954	knbi = inbi0(ii)
955	knbj = inbj0(ii)
956	IF(PRESENT(knbcnt)) knbcnt = knbi * knbj - inbij
957	DEALLOCATE( inbi0, inbj0 )
958	!
959	END SUBROUTINE mpp_init_bestpartition
960
961
962	SUBROUTINE mpp_init_landprop( propland )
963	!!----------------------------------------------------------------------
964	!! * ROUTINE mpp_init_landprop *
965	!!
966	!! ** Purpose : the the proportion of land points in the surface land-sea mask
967	!!
968	!! ** Method : read iproc strips (of length jpiglo) of the land-sea mask
969	!!----------------------------------------------------------------------
970	REAL(wp), INTENT( out) :: propland ! proportion of land points in the global domain (between 0 and 1)
971	!
972	INTEGER, DIMENSION(jpni*jpnj) :: kusedom_1d
973	INTEGER :: inboce, iarea
974	INTEGER :: iproc, idiv, ijsz
975	INTEGER :: ijstr
976	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce
977	!!----------------------------------------------------------------------
978	! do nothing if there is no land-sea mask
979	IF( numbot == -1 ) THEN
980	propland = 0.
981	RETURN
982	ENDIF
983
984	! number of processes reading the bathymetry file
985	iproc = MINVAL( (/mppsize, jpjglo/2, 100/) ) ! read a least 2 lines, no more that 100 processes reading at the same time
986
987	! we want to read iproc strips of the land-sea mask. -> pick up iproc processes every idiv processes starting at 1
988	IF( iproc == 1 ) THEN ; idiv = mppsize
989	ELSE ; idiv = ( mppsize - 1 ) / ( iproc - 1 )
990	ENDIF
991
992	iarea = (narea-1)/idiv ! involed process number (starting counting at 0)
993	IF( MOD( narea-1, idiv ) == 0 .AND. iarea < iproc ) THEN ! beware idiv can be = to 1
994	!
995	ijsz = jpjglo / iproc ! width of the stripe to read
996	IF( iarea < MOD(jpjglo,iproc) ) ijsz = ijsz + 1
997	ijstr = iarea*(jpjglo/iproc) + MIN(iarea, MOD(jpjglo,iproc)) + 1 ! starting j position of the reading
998	!
999	ALLOCATE( lloce(jpiglo, ijsz) ) ! allocate the strip
1000	CALL mpp_init_readbot_strip( ijstr, ijsz, lloce )
1001	inboce = COUNT(lloce) ! number of ocean point in the stripe
1002	DEALLOCATE(lloce)
1003	!
1004	ELSE
1005	inboce = 0
1006	ENDIF
1007	CALL mpp_sum( 'mppini', inboce ) ! total number of ocean points over the global domain
1008	!
1009	propland = REAL( jpiglojpjglo - inboce, wp ) / REAL( jpiglojpjglo, wp )
1010	!
1011	END SUBROUTINE mpp_init_landprop
1012
1013
1014	SUBROUTINE mpp_init_isoce( knbi, knbj, ldisoce )
1015	!!----------------------------------------------------------------------
1016	!! * ROUTINE mpp_init_nboce *
1017	!!
1018	!! ** Purpose : check for a mpi domain decomposition knbi x knbj which
1019	!! subdomains contain at least 1 ocean point
1020	!!
1021	!! ** Method : read knbj strips (of length jpiglo) of the land-sea mask
1022	!!----------------------------------------------------------------------
1023	INTEGER, INTENT(in ) :: knbi, knbj ! domain decomposition
1024	LOGICAL, DIMENSION(knbi,knbj), INTENT( out) :: ldisoce ! .true. if a sub domain constains 1 ocean point
1025	!
1026	INTEGER, DIMENSION(knbi,knbj) :: inboce ! number oce oce pint in each mpi subdomain
1027	INTEGER, DIMENSION(knbi*knbj) :: inboce_1d
1028	INTEGER :: idiv, iimax, ijmax, iarea
1029	INTEGER :: ji, jn
1030	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce ! lloce(i,j) = .true. if the point (i,j) is ocean
1031	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ilci
1032	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ilcj
1033	!!----------------------------------------------------------------------
1034	! do nothing if there is no land-sea mask
1035	IF( numbot == -1 ) THEN
1036	ldisoce(:,:) = .TRUE.
1037	RETURN
1038	ENDIF
1039
1040	! we want to read knbj strips of the land-sea mask. -> pick up knbj processes every idiv processes starting at 1
1041	IF ( knbj == 1 ) THEN ; idiv = mppsize
1042	ELSE IF ( mppsize < knbj ) THEN ; idiv = 1
1043	ELSE ; idiv = ( mppsize - 1 ) / ( knbj - 1 )
1044	ENDIF
1045	inboce(:,:) = 0 ! default no ocean point found
1046
1047	DO jn = 0, (knbj-1)/mppsize ! if mppsize < knbj : more strips than mpi processes (because of potential land domains)
1048	!
1049	iarea = (narea-1)/idiv + jn * mppsize ! involed process number (starting counting at 0)
1050	IF( MOD( narea-1, idiv ) == 0 .AND. iarea < knbj ) THEN ! beware idiv can be = to 1
1051	!
1052	ALLOCATE( iimppt(knbi,knbj), ijmppt(knbi,knbj), ilci(knbi,knbj), ilcj(knbi,knbj) )
1053	CALL mpp_basic_decomposition( knbi, knbj, iimax, ijmax, iimppt, ijmppt, ilci, ilcj )
1054	!
1055	ALLOCATE( lloce(jpiglo, ilcj(1,iarea+1)) ) ! allocate the strip
1056	CALL mpp_init_readbot_strip( ijmppt(1,iarea+1), ilcj(1,iarea+1), lloce ) ! read the strip
1057	DO ji = 1, knbi
1058	inboce(ji,iarea+1) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ilci(ji,1)-1,:) ) ! number of ocean point in subdomain
1059	END DO
1060	!
1061	DEALLOCATE(lloce)
1062	DEALLOCATE(iimppt, ijmppt, ilci, ilcj)
1063	!
1064	ENDIF
1065	END DO
1066
1067	inboce_1d = RESHAPE(inboce, (/ knbi*knbj /))
1068	CALL mpp_sum( 'mppini', inboce_1d )
1069	inboce = RESHAPE(inboce_1d, (/knbi, knbj/))
1070	ldisoce(:,:) = inboce(:,:) /= 0
1071	!
1072	END SUBROUTINE mpp_init_isoce
1073
1074
1075	SUBROUTINE mpp_init_readbot_strip( kjstr, kjcnt, ldoce )
1076	!!----------------------------------------------------------------------
1077	!! * ROUTINE mpp_init_readbot_strip *
1078	!!
1079	!! ** Purpose : Read relevant bathymetric information in order to
1080	!! provide a land/sea mask used for the elimination
1081	!! of land domains, in an mpp computation.
1082	!!
1083	!! ** Method : read stipe of size (jpiglo,...)
1084	!!----------------------------------------------------------------------
1085	INTEGER , INTENT(in ) :: kjstr ! starting j position of the reading
1086	INTEGER , INTENT(in ) :: kjcnt ! number of lines to read
1087	LOGICAL, DIMENSION(jpiglo,kjcnt), INTENT( out) :: ldoce ! ldoce(i,j) = .true. if the point (i,j) is ocean
1088	!
1089	INTEGER :: inumsave ! local logical unit
1090	REAL(wp), DIMENSION(jpiglo,kjcnt) :: zbot
1091	!!----------------------------------------------------------------------
1092	!
1093	inumsave = numout ; numout = numnul ! redirect all print to /dev/null
1094	!
1095	IF( numbot /= -1 ) THEN
1096	CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
1097	ELSE
1098	zbot(:,:) = 1. ! put a non-null value
1099	ENDIF
1100
1101	!
1102	ldoce(:,:) = zbot(:,:) > 0.
1103	numout = inumsave
1104	!
1105	END SUBROUTINE mpp_init_readbot_strip
1106
1107	SUBROUTINE mpp_init_nfdcom
1108	!!----------------------------------------------------------------------
1109	!! * ROUTINE mpp_init_nfdcom *
1110	!! ** Purpose : Setup for north fold exchanges with explicit
1111	!! point-to-point messaging
1112	!!
1113	!! ** Method : Initialization of the northern neighbours lists.
1114	!!----------------------------------------------------------------------
1115	!! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE)
1116	!! 2.0 ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
1117	!!----------------------------------------------------------------------
1118	INTEGER :: sxM, dxM, sxT, dxT, jn
1119	INTEGER :: njmppmax
1120	!!----------------------------------------------------------------------
1121	!
1122	njmppmax = MAXVAL( njmppt )
1123	!
1124	!initializes the north-fold communication variables
1125	isendto(:) = 0
1126	nsndto = 0
1127	!
1128	IF ( njmpp == njmppmax ) THEN ! if I am a process in the north
1129	!
1130	!sxM is the first point (in the global domain) needed to compute the north-fold for the current process
1131	sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
1132	!dxM is the last point (in the global domain) needed to compute the north-fold for the current process
1133	dxM = jpiglo - nimppt(narea) + 2
1134	!
1135	! loop over the other north-fold processes to find the processes
1136	! managing the points belonging to the sxT-dxT range
1137	!
1138	DO jn = 1, jpni
1139	!
1140	sxT = nfiimpp(jn, jpnj) ! sxT = 1st point (in the global domain) of the jn process
1141	dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1 ! dxT = last point (in the global domain) of the jn process
1142	!
1143	IF ( sxT < sxM .AND. sxM < dxT ) THEN
1144	nsndto = nsndto + 1
1145	isendto(nsndto) = jn
1146	ELSEIF( sxM <= sxT .AND. dxM >= dxT ) THEN
1147	nsndto = nsndto + 1
1148	isendto(nsndto) = jn
1149	ELSEIF( dxM < dxT .AND. sxT < dxM ) THEN
1150	nsndto = nsndto + 1
1151	isendto(nsndto) = jn
1152	ENDIF
1153	!
1154	END DO
1155	nfsloop = 1
1156	nfeloop = nlci
1157	DO jn = 2,jpni-1
1158	IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
1159	IF( nfipproc(jn-1,jpnj) == -1 ) nfsloop = nldi
1160	IF( nfipproc(jn+1,jpnj) == -1 ) nfeloop = nlei
1161	ENDIF
1162	END DO
1163	!
1164	ENDIF
1165	l_north_nogather = .TRUE.
1166	!
1167	END SUBROUTINE mpp_init_nfdcom
1168
1169
1170	#endif
1171
1172	!!======================================================================
1173	END MODULE mppini

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

Download in other formats: