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

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

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

Context Navigation

source: utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/mppini.F90 @ 13024

Last change on this file since 13024 was 13024, checked in by rblod, 5 years ago
First version of new nesting tools merged with domaincfg, see ticket #2129
File size: 58.6 KB

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

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

Download in other formats: