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

Last change on this file since 10560 was 10560, checked in by smasson, 5 years ago
trunk: bugfix in mppini for automatic domain decomposition with low number of cores
Property svn:keywords set to `Id`
File size: 56.8 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, irm, ijpjmin
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	IF( MINVAL(klci) < 3 ) THEN
708	WRITE(ctmp1,*) ' mpp_basic_decomposition: minimum value of jpi must be >= 3'
709	WRITE(ctmp2,*) ' We have ', MINVAL(klci)
710	CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
711	ENDIF
712	IF( jperio == 3 .OR. jperio == 4 .OR. jperio == 5 .OR. jperio == 6 ) THEN
713	! minimize the size of the last row to compensate for the north pole folding coast
714	IF( jperio == 3 .OR. jperio == 4 ) ijpjmin = 5 ! V and F folding involves line jpj-3 that must not be south boundary
715	IF( jperio == 5 .OR. jperio == 6 ) ijpjmin = 4 ! V and F folding involves line jpj-2 that must not be south boundary
716	irm = knbj - irestj ! total number of lines to be removed
717	klcj(:, knbj) = MAX( ijpjmin, kjmax-irm ) ! we must have jpj >= ijpjmin in the last row
718	irm = irm - ( kjmax - klcj(1,knbj) ) ! remaining number of lines to remove
719	irestj = knbj - 1 - irm
720	klcj(:, 1:irestj) = kjmax
721	klcj(:, irestj+1:knbj-1) = kjmax-1
722	ELSE
723	ijpjmin = 3
724	klcj(:, 1:irestj) = kjmax
725	klcj(:, irestj+1:knbj) = kjmax-1
726	ENDIF
727	IF( MINVAL(klcj) < ijpjmin ) THEN
728	WRITE(ctmp1,*) ' mpp_basic_decomposition: minimum value of jpj must be >= ', ijpjmin
729	WRITE(ctmp2,*) ' We have ', MINVAL(klcj)
730	CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
731	ENDIF
732	#endif
733
734	! 2. Index arrays for subdomains
735	! -------------------------------
736	kimppt(:,:) = 1
737	kjmppt(:,:) = 1
738	!
739	IF( knbi > 1 ) THEN
740	DO jj = 1, knbj
741	DO ji = 2, knbi
742	kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - ireci
743	END DO
744	END DO
745	ENDIF
746	!
747	IF( knbj > 1 )THEN
748	DO jj = 2, knbj
749	DO ji = 1, knbi
750	kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - irecj
751	END DO
752	END DO
753	ENDIF
754
755	END SUBROUTINE mpp_basic_decomposition
756
757
758	SUBROUTINE mpp_init_bestpartition( knbij, knbi, knbj, knbcnt, ldlist )
759	!!----------------------------------------------------------------------
760	!! * ROUTINE mpp_init_bestpartition *
761	!!
762	!! ** Purpose :
763	!!
764	!! ** Method :
765	!!----------------------------------------------------------------------
766	INTEGER, INTENT(in ) :: knbij ! total number if subdomains (knbi*knbj)
767	INTEGER, OPTIONAL, INTENT( out) :: knbi, knbj ! number if subdomains along i and j (knbi and knbj)
768	INTEGER, OPTIONAL, INTENT( out) :: knbcnt ! number of land subdomains
769	LOGICAL, OPTIONAL, INTENT(in ) :: ldlist ! .true.: print the list the best domain decompositions (with land)
770	!
771	INTEGER :: ji, jj, ii, iitarget
772	INTEGER :: iszitst, iszjtst
773	INTEGER :: isziref, iszjref
774	INTEGER :: inbij, iszij
775	INTEGER :: inbimax, inbjmax, inbijmax
776	INTEGER :: isz0, isz1
777	INTEGER, DIMENSION( :), ALLOCATABLE :: indexok
778	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi0, inbj0, inbij0 ! number of subdomains along i,j
779	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi0, iszj0, iszij0 ! max size of the subdomains along i,j
780	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi1, inbj1, inbij1 ! number of subdomains along i,j
781	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi1, iszj1, iszij1 ! max size of the subdomains along i,j
782	LOGICAL :: llist
783	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d ! max size of the subdomains along i,j
784	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisoce ! - -
785	REAL(wp):: zpropland
786	!!----------------------------------------------------------------------
787	!
788	llist = .FALSE.
789	IF( PRESENT(ldlist) ) llist = ldlist
790
791	CALL mpp_init_landprop( zpropland ) ! get the proportion of land point over the gloal domain
792	inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) ) ! define the largest possible value for jpni*jpnj
793	!
794	IF( llist ) THEN ; inbijmax = inbij*2
795	ELSE ; inbijmax = inbij
796	ENDIF
797	!
798	ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax))
799	!
800	inbimax = 0
801	inbjmax = 0
802	isziref = jpiglo*jpjglo+1
803	iszjref = jpiglo*jpjglo+1
804	!
805	! get the list of knbi that gives a smaller jpimax than knbi-1
806	! get the list of knbj that gives a smaller jpjmax than knbj-1
807	DO ji = 1, inbijmax
808	#if defined key_nemocice_decomp
809	iszitst = ( nx_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
810	#else
811	iszitst = ( jpiglo - 2nn_hls + (ji-1) ) / ji + 2nn_hls
812	#endif
813	IF( iszitst < isziref ) THEN
814	isziref = iszitst
815	inbimax = inbimax + 1
816	inbi0(inbimax) = ji
817	iszi0(inbimax) = isziref
818	ENDIF
819	#if defined key_nemocice_decomp
820	iszjtst = ( ny_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
821	#else
822	iszjtst = ( jpjglo - 2nn_hls + (ji-1) ) / ji + 2nn_hls
823	#endif
824	IF( iszjtst < iszjref ) THEN
825	iszjref = iszjtst
826	inbjmax = inbjmax + 1
827	inbj0(inbjmax) = ji
828	iszj0(inbjmax) = iszjref
829	ENDIF
830	END DO
831
832	! combine these 2 lists to get all possible knbi*knbj < inbijmax
833	ALLOCATE( llmsk2d(inbimax,inbjmax) )
834	DO jj = 1, inbjmax
835	DO ji = 1, inbimax
836	IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN ; llmsk2d(ji,jj) = .TRUE.
837	ELSE ; llmsk2d(ji,jj) = .FALSE.
838	ENDIF
839	END DO
840	END DO
841	isz1 = COUNT(llmsk2d)
842	ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) )
843	ii = 0
844	DO jj = 1, inbjmax
845	DO ji = 1, inbimax
846	IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN
847	ii = ii + 1
848	inbi1(ii) = inbi0(ji)
849	inbj1(ii) = inbj0(jj)
850	iszi1(ii) = iszi0(ji)
851	iszj1(ii) = iszj0(jj)
852	END IF
853	END DO
854	END DO
855	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
856	DEALLOCATE( llmsk2d )
857
858	ALLOCATE( inbij1(isz1), iszij1(isz1) )
859	inbij1(:) = inbi1(:) * inbj1(:)
860	iszij1(:) = iszi1(:) * iszj1(:)
861
862	! if therr is no land and no print
863	IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN
864	! get the smaller partition which gives the smallest subdomain size
865	ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1)
866	knbi = inbi1(ii)
867	knbj = inbj1(ii)
868	IF(PRESENT(knbcnt)) knbcnt = 0
869	DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 )
870	RETURN
871	ENDIF
872
873	! extract only the partitions which reduce the subdomain size in comparison with smaller partitions
874	ALLOCATE( indexok(isz1) ) ! to store indices of the best partitions
875	isz0 = 0 ! number of best partitions
876	inbij = 1 ! start with the min value of inbij1 => 1
877	iszij = jpiglo*jpjglo+1 ! default: larger than global domain
878	DO WHILE( inbij <= inbijmax ) ! if we did not reach the max of inbij1
879	ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1) ! warning: send back the first occurence if multiple results
880	IF ( iszij1(ii) < iszij ) THEN
881	isz0 = isz0 + 1
882	indexok(isz0) = ii
883	iszij = iszij1(ii)
884	ENDIF
885	inbij = MINVAL(inbij1, mask = inbij1 > inbij) ! warning: return largest integer value if mask = .false. everywhere
886	END DO
887	DEALLOCATE( inbij1, iszij1 )
888
889	! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size)
890	ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) )
891	DO ji = 1, isz0
892	ii = indexok(ji)
893	inbi0(ji) = inbi1(ii)
894	inbj0(ji) = inbj1(ii)
895	iszi0(ji) = iszi1(ii)
896	iszj0(ji) = iszj1(ii)
897	END DO
898	DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 )
899
900	IF( llist ) THEN ! we print about 21 best partitions
901	IF(lwp) THEN
902	WRITE(numout,*)
903	WRITE(numout, *) ' For your information:'
904	WRITE(numout,'(a,i5,a)') ' list of the best partitions around ', knbij, ' mpi processes'
905	WRITE(numout, *) ' --------------------------------------', '-----', '--------------'
906	WRITE(numout,*)
907	END IF
908	iitarget = MINLOC( inbi0(:)inbj0(:), mask = inbi0(:)inbj0(:) >= knbij, dim = 1 )
909	DO ji = MAX(1,iitarget-10), MIN(isz0,iitarget+10)
910	ALLOCATE( llisoce(inbi0(ji), inbj0(ji)) )
911	CALL mpp_init_isoce( inbi0(ji), inbj0(ji), llisoce ) ! Warning: must be call by all cores (call mpp_sum)
912	inbij = COUNT(llisoce)
913	DEALLOCATE( llisoce )
914	IF(lwp) WRITE(numout,'(a, i5, a, i5, a, i4, a, i4, a, i9, a, i5, a, i5, a)') &
915	& 'nb_cores ' , inbij,' oce + ', inbi0(ji)*inbj0(ji) - inbij &
916	& , ' land ( ', inbi0(ji),' x ', inbj0(ji), &
917	& ' ), nb_points ', iszi0(ji)*iszj0(ji),' ( ', iszi0(ji),' x ', iszj0(ji),' )'
918	END DO
919	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
920	RETURN
921	ENDIF
922
923	DEALLOCATE( iszi0, iszj0 )
924	inbij = inbijmax + 1 ! default: larger than possible
925	ii = isz0+1 ! start from the end of the list (smaller subdomains)
926	DO WHILE( inbij > knbij ) ! while the number of ocean subdomains exceed the number of procs
927	ii = ii -1
928	ALLOCATE( llisoce(inbi0(ii), inbj0(ii)) )
929	CALL mpp_init_isoce( inbi0(ii), inbj0(ii), llisoce ) ! must be done by all core
930	inbij = COUNT(llisoce)
931	DEALLOCATE( llisoce )
932	END DO
933	knbi = inbi0(ii)
934	knbj = inbj0(ii)
935	IF(PRESENT(knbcnt)) knbcnt = knbi * knbj - inbij
936	DEALLOCATE( inbi0, inbj0 )
937	!
938	END SUBROUTINE mpp_init_bestpartition
939
940
941	SUBROUTINE mpp_init_landprop( propland )
942	!!----------------------------------------------------------------------
943	!! * ROUTINE mpp_init_landprop *
944	!!
945	!! ** Purpose : the the proportion of land points in the surface land-sea mask
946	!!
947	!! ** Method : read iproc strips (of length jpiglo) of the land-sea mask
948	!!----------------------------------------------------------------------
949	REAL(wp), INTENT( out) :: propland ! proportion of land points in the global domain (between 0 and 1)
950	!
951	INTEGER, DIMENSION(jpni*jpnj) :: kusedom_1d
952	INTEGER :: inboce, iarea
953	INTEGER :: iproc, idiv, ijsz
954	INTEGER :: ijstr
955	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce
956	!!----------------------------------------------------------------------
957	! do nothing if there is no land-sea mask
958	IF( numbot == -1 .and. numbdy == -1 ) THEN
959	propland = 0.
960	RETURN
961	ENDIF
962
963	! number of processes reading the bathymetry file
964	iproc = MINVAL( (/mppsize, jpjglo/2, 100/) ) ! read a least 2 lines, no more that 100 processes reading at the same time
965
966	! we want to read iproc strips of the land-sea mask. -> pick up iproc processes every idiv processes starting at 1
967	IF( iproc == 1 ) THEN ; idiv = mppsize
968	ELSE ; idiv = ( mppsize - 1 ) / ( iproc - 1 )
969	ENDIF
970
971	iarea = (narea-1)/idiv ! involed process number (starting counting at 0)
972	IF( MOD( narea-1, idiv ) == 0 .AND. iarea < iproc ) THEN ! beware idiv can be = to 1
973	!
974	ijsz = jpjglo / iproc ! width of the stripe to read
975	IF( iarea < MOD(jpjglo,iproc) ) ijsz = ijsz + 1
976	ijstr = iarea*(jpjglo/iproc) + MIN(iarea, MOD(jpjglo,iproc)) + 1 ! starting j position of the reading
977	!
978	ALLOCATE( lloce(jpiglo, ijsz) ) ! allocate the strip
979	CALL mpp_init_readbot_strip( ijstr, ijsz, lloce )
980	inboce = COUNT(lloce) ! number of ocean point in the stripe
981	DEALLOCATE(lloce)
982	!
983	ELSE
984	inboce = 0
985	ENDIF
986	CALL mpp_sum( 'mppini', inboce ) ! total number of ocean points over the global domain
987	!
988	propland = REAL( jpiglojpjglo - inboce, wp ) / REAL( jpiglojpjglo, wp )
989	!
990	END SUBROUTINE mpp_init_landprop
991
992
993	SUBROUTINE mpp_init_isoce( knbi, knbj, ldisoce )
994	!!----------------------------------------------------------------------
995	!! * ROUTINE mpp_init_nboce *
996	!!
997	!! ** Purpose : check for a mpi domain decomposition knbi x knbj which
998	!! subdomains contain at least 1 ocean point
999	!!
1000	!! ** Method : read knbj strips (of length jpiglo) of the land-sea mask
1001	!!----------------------------------------------------------------------
1002	INTEGER, INTENT(in ) :: knbi, knbj ! domain decomposition
1003	LOGICAL, DIMENSION(knbi,knbj), INTENT( out) :: ldisoce ! .true. if a sub domain constains 1 ocean point
1004	!
1005	INTEGER, DIMENSION(knbi,knbj) :: inboce ! number oce oce pint in each mpi subdomain
1006	INTEGER, DIMENSION(knbi*knbj) :: inboce_1d
1007	INTEGER :: idiv, iimax, ijmax, iarea
1008	INTEGER :: ji, jn
1009	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce ! lloce(i,j) = .true. if the point (i,j) is ocean
1010	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ilci
1011	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ilcj
1012	!!----------------------------------------------------------------------
1013	! do nothing if there is no land-sea mask
1014	IF( numbot == -1 .AND. numbdy == -1 ) THEN
1015	ldisoce(:,:) = .TRUE.
1016	RETURN
1017	ENDIF
1018
1019	! we want to read knbj strips of the land-sea mask. -> pick up knbj processes every idiv processes starting at 1
1020	IF ( knbj == 1 ) THEN ; idiv = mppsize
1021	ELSE IF ( mppsize < knbj ) THEN ; idiv = 1
1022	ELSE ; idiv = ( mppsize - 1 ) / ( knbj - 1 )
1023	ENDIF
1024	inboce(:,:) = 0 ! default no ocean point found
1025
1026	DO jn = 0, (knbj-1)/mppsize ! if mppsize < knbj : more strips than mpi processes (because of potential land domains)
1027	!
1028	iarea = (narea-1)/idiv + jn * mppsize ! involed process number (starting counting at 0)
1029	IF( MOD( narea-1, idiv ) == 0 .AND. iarea < knbj ) THEN ! beware idiv can be = to 1
1030	!
1031	ALLOCATE( iimppt(knbi,knbj), ijmppt(knbi,knbj), ilci(knbi,knbj), ilcj(knbi,knbj) )
1032	CALL mpp_basic_decomposition( knbi, knbj, iimax, ijmax, iimppt, ijmppt, ilci, ilcj )
1033	!
1034	ALLOCATE( lloce(jpiglo, ilcj(1,iarea+1)) ) ! allocate the strip
1035	CALL mpp_init_readbot_strip( ijmppt(1,iarea+1), ilcj(1,iarea+1), lloce ) ! read the strip
1036	DO ji = 1, knbi
1037	inboce(ji,iarea+1) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ilci(ji,1)-1,:) ) ! number of ocean point in subdomain
1038	END DO
1039	!
1040	DEALLOCATE(lloce)
1041	DEALLOCATE(iimppt, ijmppt, ilci, ilcj)
1042	!
1043	ENDIF
1044	END DO
1045
1046	inboce_1d = RESHAPE(inboce, (/ knbi*knbj /))
1047	CALL mpp_sum( 'mppini', inboce_1d )
1048	inboce = RESHAPE(inboce_1d, (/knbi, knbj/))
1049	ldisoce(:,:) = inboce(:,:) /= 0
1050	!
1051	END SUBROUTINE mpp_init_isoce
1052
1053
1054	SUBROUTINE mpp_init_readbot_strip( kjstr, kjcnt, ldoce )
1055	!!----------------------------------------------------------------------
1056	!! * ROUTINE mpp_init_readbot_strip *
1057	!!
1058	!! ** Purpose : Read relevant bathymetric information in order to
1059	!! provide a land/sea mask used for the elimination
1060	!! of land domains, in an mpp computation.
1061	!!
1062	!! ** Method : read stipe of size (jpiglo,...)
1063	!!----------------------------------------------------------------------
1064	INTEGER , INTENT(in ) :: kjstr ! starting j position of the reading
1065	INTEGER , INTENT(in ) :: kjcnt ! number of lines to read
1066	LOGICAL, DIMENSION(jpiglo,kjcnt), INTENT( out) :: ldoce ! ldoce(i,j) = .true. if the point (i,j) is ocean
1067	!
1068	INTEGER :: inumsave ! local logical unit
1069	REAL(wp), DIMENSION(jpiglo,kjcnt) :: zbot, zbdy
1070	!!----------------------------------------------------------------------
1071	!
1072	inumsave = numout ; numout = numnul ! redirect all print to /dev/null
1073	!
1074	IF( numbot /= -1 ) THEN
1075	CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
1076	ELSE
1077	zbot(:,:) = 1. ! put a non-null value
1078	ENDIF
1079
1080	IF( numbdy /= -1 ) THEN ! Adjust with bdy_msk if it exists
1081	CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
1082	zbot(:,:) = zbot(:,:) * zbdy(:,:)
1083	ENDIF
1084	!
1085	ldoce(:,:) = zbot(:,:) > 0.
1086	numout = inumsave
1087	!
1088	END SUBROUTINE mpp_init_readbot_strip
1089
1090
1091	SUBROUTINE mpp_init_ioipsl
1092	!!----------------------------------------------------------------------
1093	!! * ROUTINE mpp_init_ioipsl *
1094	!!
1095	!! ** Purpose :
1096	!!
1097	!! ** Method :
1098	!!
1099	!! History :
1100	!! 9.0 ! 04-03 (G. Madec ) MPP-IOIPSL
1101	!! " " ! 08-12 (A. Coward) addition in case of jpni*jpnj < jpnij
1102	!!----------------------------------------------------------------------
1103	INTEGER, DIMENSION(2) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
1104	!!----------------------------------------------------------------------
1105
1106	! The domain is split only horizontally along i- or/and j- direction
1107	! So we need at the most only 1D arrays with 2 elements.
1108	! Set idompar values equivalent to the jpdom_local_noextra definition
1109	! used in IOM. This works even if jpnij .ne. jpni*jpnj.
1110	iglo(1) = jpiglo
1111	iglo(2) = jpjglo
1112	iloc(1) = nlci
1113	iloc(2) = nlcj
1114	iabsf(1) = nimppt(narea)
1115	iabsf(2) = njmppt(narea)
1116	iabsl(:) = iabsf(:) + iloc(:) - 1
1117	ihals(1) = nldi - 1
1118	ihals(2) = nldj - 1
1119	ihale(1) = nlci - nlei
1120	ihale(2) = nlcj - nlej
1121	idid(1) = 1
1122	idid(2) = 2
1123
1124	IF(lwp) THEN
1125	WRITE(numout,*)
1126	WRITE(numout,*) 'mpp_init_ioipsl : iloc = ', iloc (1), iloc (2)
1127	WRITE(numout,*) '~~~~~~~~~~~~~~~ iabsf = ', iabsf(1), iabsf(2)
1128	WRITE(numout,*) ' ihals = ', ihals(1), ihals(2)
1129	WRITE(numout,*) ' ihale = ', ihale(1), ihale(2)
1130	ENDIF
1131	!
1132	CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
1133	!
1134	END SUBROUTINE mpp_init_ioipsl
1135
1136
1137	SUBROUTINE mpp_init_nfdcom
1138	!!----------------------------------------------------------------------
1139	!! * ROUTINE mpp_init_nfdcom *
1140	!! ** Purpose : Setup for north fold exchanges with explicit
1141	!! point-to-point messaging
1142	!!
1143	!! ** Method : Initialization of the northern neighbours lists.
1144	!!----------------------------------------------------------------------
1145	!! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE)
1146	!! 2.0 ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
1147	!!----------------------------------------------------------------------
1148	INTEGER :: sxM, dxM, sxT, dxT, jn
1149	INTEGER :: njmppmax
1150	!!----------------------------------------------------------------------
1151	!
1152	njmppmax = MAXVAL( njmppt )
1153	!
1154	!initializes the north-fold communication variables
1155	isendto(:) = 0
1156	nsndto = 0
1157	!
1158	IF ( njmpp == njmppmax ) THEN ! if I am a process in the north
1159	!
1160	!sxM is the first point (in the global domain) needed to compute the north-fold for the current process
1161	sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
1162	!dxM is the last point (in the global domain) needed to compute the north-fold for the current process
1163	dxM = jpiglo - nimppt(narea) + 2
1164	!
1165	! loop over the other north-fold processes to find the processes
1166	! managing the points belonging to the sxT-dxT range
1167	!
1168	DO jn = 1, jpni
1169	!
1170	sxT = nfiimpp(jn, jpnj) ! sxT = 1st point (in the global domain) of the jn process
1171	dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1 ! dxT = last point (in the global domain) of the jn process
1172	!
1173	IF ( sxT < sxM .AND. sxM < dxT ) THEN
1174	nsndto = nsndto + 1
1175	isendto(nsndto) = jn
1176	ELSEIF( sxM <= sxT .AND. dxM >= dxT ) THEN
1177	nsndto = nsndto + 1
1178	isendto(nsndto) = jn
1179	ELSEIF( dxM < dxT .AND. sxT < dxM ) THEN
1180	nsndto = nsndto + 1
1181	isendto(nsndto) = jn
1182	ENDIF
1183	!
1184	END DO
1185	nfsloop = 1
1186	nfeloop = nlci
1187	DO jn = 2,jpni-1
1188	IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
1189	IF( nfipproc(jn-1,jpnj) == -1 ) nfsloop = nldi
1190	IF( nfipproc(jn+1,jpnj) == -1 ) nfeloop = nlei
1191	ENDIF
1192	END DO
1193	!
1194	ENDIF
1195	l_north_nogather = .TRUE.
1196	!
1197	END SUBROUTINE mpp_init_nfdcom
1198
1199
1200	#endif
1201
1202	!!======================================================================
1203	END MODULE mppini

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

Download in other formats: