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

Last change on this file since 10425 was 10425, checked in by smasson, 5 years ago
trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk
Property svn:keywords set to `Id`
File size: 54.7 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 (between 0 and 1)
927	!
928	INTEGER, DIMENSION(jpni*jpnj) :: kusedom_1d
929	INTEGER :: inboce
930	INTEGER :: iproc, idiv, ijsz
931	INTEGER :: ijstr, ijend, ijcnt
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 among mppsize processes
944	IF( iproc == 1 ) THEN ; idiv = mppsize
945	ELSE ; idiv = ( mppsize - 1 ) / ( iproc - 1 )
946	ENDIF
947	ijsz = jpjglo / iproc
948	IF( narea <= MOD(jpjglo,iproc) ) ijsz = ijsz + 1
949
950	IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= iproc ) ) THEN
951	!
952	ijstr = (narea-1)*(jpjglo/iproc) + MIN(narea-1, MOD(jpjglo,iproc)) + 1
953	ijend = ijstr + ijsz - 1
954	ijcnt = ijend - ijstr + 1
955	!
956	ALLOCATE( lloce(jpiglo, ijcnt) ) ! allocate the strip
957	CALL mpp_init_readbot_strip( ijstr, ijcnt, lloce )
958	inboce = COUNT(lloce)
959	DEALLOCATE(lloce)
960	!
961	ELSE
962	inboce = 0
963	ENDIF
964	CALL mpp_sum( 'mppini', inboce )
965	!
966	propland = REAL( jpiglojpjglo - inboce, wp ) / REAL( jpiglojpjglo, wp )
967	!
968	END SUBROUTINE mpp_init_landprop
969
970
971	SUBROUTINE mpp_init_isoce( knbi, knbj, ldisoce )
972	!!----------------------------------------------------------------------
973	!! * ROUTINE mpp_init_nboce *
974	!!
975	!! ** Purpose : check for a mpi domain decomposition knbi x knbj which
976	!! subdomains contain at least 1 ocean point
977	!!
978	!! ** Method : read knbj strips (of length jpiglo) of the land-sea mask
979	!!----------------------------------------------------------------------
980	INTEGER, INTENT(in ) :: knbi, knbj
981	LOGICAL, DIMENSION(knbi,knbj), INTENT( out) :: ldisoce !
982	!
983	INTEGER, DIMENSION(knbi,knbj) :: inboce
984	INTEGER, DIMENSION(knbi*knbj) :: inboce_1d
985	INTEGER :: idiv, i2read, inj
986	INTEGER :: iimax, ijmax
987	INTEGER :: ji,jj
988	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce
989	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ilci
990	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ilcj
991	!!----------------------------------------------------------------------
992	! do nothing if there is no land-sea mask
993	IF( numbot == -1 .AND. numbdy == -1 ) THEN
994	ldisoce(:,:) = .TRUE.
995	RETURN
996	ENDIF
997
998	! we want to read knbj strips of the land-sea mask. -> pick up knbj processes among mppsize processes
999	IF( knbj == 1 ) THEN ; idiv = mppsize
1000	ELSE ; idiv = ( mppsize - 1 ) / ( knbj - 1 )
1001	ENDIF
1002	inboce(:,:) = 0
1003	IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= knbj ) ) THEN
1004	!
1005	ALLOCATE( iimppt(knbi,knbj), ijmppt(knbi,knbj), ilci(knbi,knbj), ilcj(knbi,knbj) )
1006	CALL mpp_basic_decomposition( knbi, knbj, iimax, ijmax, iimppt, ijmppt, ilci, ilcj )
1007	!
1008	i2read = knbj / mppsize ! strip number to be read by this process
1009	IF( ( narea - 1 ) / idiv < MOD(knbj,mppsize) ) i2read = i2read + 1
1010	DO jj = 1, i2read
1011	! strip number to be read (from 1 to knbj)
1012	inj = ( narea - 1 ) * ( knbj / mppsize ) + MIN( MOD(knbj,mppsize), ( narea - 1 ) / idiv ) + jj
1013	ALLOCATE( lloce(jpiglo, ilcj(1,inj)) ) ! allocate the strip
1014	CALL mpp_init_readbot_strip( ijmppt(1,inj), ilcj(1,inj), lloce ) ! read the strip
1015	DO ji = 1, knbi
1016	inboce(ji,inj) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ilci(ji,1)-1,:) )
1017	END DO
1018	DEALLOCATE(lloce)
1019	END DO
1020	!
1021	DEALLOCATE(iimppt, ijmppt, ilci, ilcj)
1022	!
1023	ENDIF
1024
1025	inboce_1d = RESHAPE(inboce, (/ knbi*knbj /))
1026	CALL mpp_sum( 'mppini', inboce_1d )
1027	inboce = RESHAPE(inboce_1d, (/knbi, knbj/))
1028	ldisoce = inboce /= 0
1029	!
1030	END SUBROUTINE mpp_init_isoce
1031
1032
1033	SUBROUTINE mpp_init_readbot_strip( kjstr, kjcnt, ldoce )
1034	!!----------------------------------------------------------------------
1035	!! * ROUTINE mpp_init_readbot_strip *
1036	!!
1037	!! ** Purpose : Read relevant bathymetric information in order to
1038	!! provide a land/sea mask used for the elimination
1039	!! of land domains, in an mpp computation.
1040	!!
1041	!! ** Method : read stipe of size (jpiglo,...)
1042	!!----------------------------------------------------------------------
1043	INTEGER , INTENT(in ) :: kjstr !
1044	INTEGER , INTENT(in ) :: kjcnt !
1045	LOGICAL, DIMENSION(jpiglo,kjcnt), INTENT( out) :: ldoce !
1046	!
1047	INTEGER :: inumsave ! local logical unit
1048	REAL(wp), DIMENSION(jpiglo,kjcnt) :: zbot, zbdy
1049	!!----------------------------------------------------------------------
1050	!
1051	inumsave = numout ; numout = numnul ! redirect all print to /dev/null
1052	!
1053	IF( numbot /= -1 ) THEN
1054	CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
1055	ELSE
1056	zbot(:,:) = 1. ! put a non-null value
1057	ENDIF
1058
1059	IF( numbdy /= -1 ) THEN ! Adjust with bdy_msk if it exists
1060	CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
1061	zbot(:,:) = zbot(:,:) * zbdy(:,:)
1062	ENDIF
1063	!
1064	ldoce = zbot > 0.
1065	numout = inumsave
1066	!
1067	END SUBROUTINE mpp_init_readbot_strip
1068
1069
1070	SUBROUTINE mpp_init_ioipsl
1071	!!----------------------------------------------------------------------
1072	!! * ROUTINE mpp_init_ioipsl *
1073	!!
1074	!! ** Purpose :
1075	!!
1076	!! ** Method :
1077	!!
1078	!! History :
1079	!! 9.0 ! 04-03 (G. Madec ) MPP-IOIPSL
1080	!! " " ! 08-12 (A. Coward) addition in case of jpni*jpnj < jpnij
1081	!!----------------------------------------------------------------------
1082	INTEGER, DIMENSION(2) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
1083	!!----------------------------------------------------------------------
1084
1085	! The domain is split only horizontally along i- or/and j- direction
1086	! So we need at the most only 1D arrays with 2 elements.
1087	! Set idompar values equivalent to the jpdom_local_noextra definition
1088	! used in IOM. This works even if jpnij .ne. jpni*jpnj.
1089	iglo(1) = jpiglo
1090	iglo(2) = jpjglo
1091	iloc(1) = nlci
1092	iloc(2) = nlcj
1093	iabsf(1) = nimppt(narea)
1094	iabsf(2) = njmppt(narea)
1095	iabsl(:) = iabsf(:) + iloc(:) - 1
1096	ihals(1) = nldi - 1
1097	ihals(2) = nldj - 1
1098	ihale(1) = nlci - nlei
1099	ihale(2) = nlcj - nlej
1100	idid(1) = 1
1101	idid(2) = 2
1102
1103	IF(lwp) THEN
1104	WRITE(numout,*)
1105	WRITE(numout,*) 'mpp_init_ioipsl : iloc = ', iloc (1), iloc (2)
1106	WRITE(numout,*) '~~~~~~~~~~~~~~~ iabsf = ', iabsf(1), iabsf(2)
1107	WRITE(numout,*) ' ihals = ', ihals(1), ihals(2)
1108	WRITE(numout,*) ' ihale = ', ihale(1), ihale(2)
1109	ENDIF
1110	!
1111	CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
1112	!
1113	END SUBROUTINE mpp_init_ioipsl
1114
1115
1116	SUBROUTINE mpp_init_nfdcom
1117	!!----------------------------------------------------------------------
1118	!! * ROUTINE mpp_init_nfdcom *
1119	!! ** Purpose : Setup for north fold exchanges with explicit
1120	!! point-to-point messaging
1121	!!
1122	!! ** Method : Initialization of the northern neighbours lists.
1123	!!----------------------------------------------------------------------
1124	!! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE)
1125	!! 2.0 ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
1126	!!----------------------------------------------------------------------
1127	INTEGER :: sxM, dxM, sxT, dxT, jn
1128	INTEGER :: njmppmax
1129	!!----------------------------------------------------------------------
1130	!
1131	njmppmax = MAXVAL( njmppt )
1132	!
1133	!initializes the north-fold communication variables
1134	isendto(:) = 0
1135	nsndto = 0
1136	!
1137	IF ( njmpp == njmppmax ) THEN ! if I am a process in the north
1138	!
1139	!sxM is the first point (in the global domain) needed to compute the north-fold for the current process
1140	sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
1141	!dxM is the last point (in the global domain) needed to compute the north-fold for the current process
1142	dxM = jpiglo - nimppt(narea) + 2
1143	!
1144	! loop over the other north-fold processes to find the processes
1145	! managing the points belonging to the sxT-dxT range
1146	!
1147	DO jn = 1, jpni
1148	!
1149	sxT = nfiimpp(jn, jpnj) ! sxT = 1st point (in the global domain) of the jn process
1150	dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1 ! dxT = last point (in the global domain) of the jn process
1151	!
1152	IF ( sxT < sxM .AND. sxM < dxT ) THEN
1153	nsndto = nsndto + 1
1154	isendto(nsndto) = jn
1155	ELSEIF( sxM <= sxT .AND. dxM >= dxT ) THEN
1156	nsndto = nsndto + 1
1157	isendto(nsndto) = jn
1158	ELSEIF( dxM < dxT .AND. sxT < dxM ) THEN
1159	nsndto = nsndto + 1
1160	isendto(nsndto) = jn
1161	ENDIF
1162	!
1163	END DO
1164	nfsloop = 1
1165	nfeloop = nlci
1166	DO jn = 2,jpni-1
1167	IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
1168	IF( nfipproc(jn-1,jpnj) == -1 ) nfsloop = nldi
1169	IF( nfipproc(jn+1,jpnj) == -1 ) nfeloop = nlei
1170	ENDIF
1171	END DO
1172	!
1173	ENDIF
1174	l_north_nogather = .TRUE.
1175	!
1176	END SUBROUTINE mpp_init_nfdcom
1177
1178
1179	#endif
1180
1181	!!======================================================================
1182	END MODULE mppini

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

Download in other formats: