Context Navigation

mppini.F90 @ 10986

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

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/LBC/mppini.F90 @ 10986

Download in other formats: