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

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

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

Download in other formats: