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 @ 11398

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

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

Download in other formats: