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/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/LBC – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/LBC/mppini.F90 @ 12745

Last change on this file since 12745 was 12586, checked in by francesca, 4 years ago
Add extra-halo support (jperio 3,4) - ticket #2366
Property svn:keywords set to `Id`
File size: 60.8 KB

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

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

Download in other formats: