New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

mppini.F90 in utils/tools_ticket2457/DOMAINcfg/src – NEMO

Context Navigation

source: utils/tools_ticket2457/DOMAINcfg/src/mppini.F90 @ 12871

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

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

Download in other formats: