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_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/LBC – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/LBC/mppini.F90 @ 11192

Last change on this file since 11192 was 11192, checked in by smasson, 5 years ago
dev_r10984_HPC-13 : reorganization of lbclnk, part 1: simpler mpp_lnk_generic.h90 supress lbc_lnk_generic.h90, see #2285
Property svn:keywords set to `Id`
File size: 59.8 KB

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

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

Download in other formats: