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/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/OCE/LBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/OCE/LBC/mppini.F90 @ 11719

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

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

Download in other formats: