Home

Context Navigation

source: palm/trunk/SOURCE/init_pegrid.f90 @ 709

Last change on this file since 709 was 709, checked in by raasch, 14 years ago
formatting adjustments
Property svn:keywords set to `Id`
File size: 40.0 KB

Line
1	SUBROUTINE init_pegrid
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	! formatting adjustments
7	!
8	! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
9	! TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values
10	!
11	! Former revisions:
12	! -----------------
13	! $Id: init_pegrid.f90 709 2011-03-30 09:31:40Z raasch $
14	!
15	! 707 2011-03-29 11:39:40Z raasch
16	! bc_lr/ns replaced by bc_lr/ns_cyc/dirrad/raddir
17	!
18	! 667 2010-12-23 12:06:00Z suehring/gryschka
19	! Moved determination of target_id's from init_coupling
20	! Determination of parameters needed for coupling (coupling_topology, ngp_a,
21	! ngp_o) with different grid/processor-topology in ocean and atmosphere
22	! Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
23	! The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
24	! maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
25	! This distinction is due to reasons of data exchange and performance for the
26	! normal grid and grids in poismg.
27	! The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
28	! New MPI-Vectors for data exchange between left and right boundaries added.
29	! This is due to reasons of performance (10% faster).
30	!
31	! 646 2010-12-15 13:03:52Z raasch
32	! lctit is now using a 2d decomposition by default
33	!
34	! 622 2010-12-10 08:08:13Z raasch
35	! optional barriers included in order to speed up collective operations
36	!
37	! 438 2010-02-01 04:32:43Z raasch
38	! 2d-decomposition is default for Cray-XT machines
39	!
40	! 274 2009-03-26 15:11:21Z heinze
41	! Output of messages replaced by message handling routine.
42	!
43	! 206 2008-10-13 14:59:11Z raasch
44	! Implementation of a MPI-1 coupling: added __parallel within the __mpi2 part
45	! 2d-decomposition is default on SGI-ICE systems
46	!
47	! 197 2008-09-16 15:29:03Z raasch
48	! multigrid levels are limited by subdomains if mg_switch_to_pe0_level = -1,
49	! nz is used instead nnz for calculating mg-levels
50	! Collect on PE0 horizontal index bounds from all other PEs,
51	! broadcast the id of the inflow PE (using the respective communicator)
52	!
53	! 114 2007-10-10 00:03:15Z raasch
54	! Allocation of wall flag arrays for multigrid solver
55	!
56	! 108 2007-08-24 15:10:38Z letzel
57	! Intercommunicator (comm_inter) and derived data type (type_xy) for
58	! coupled model runs created, assign coupling_mode_remote,
59	! indices nxlu and nysv are calculated (needed for non-cyclic boundary
60	! conditions)
61	!
62	! 82 2007-04-16 15:40:52Z raasch
63	! Cpp-directive lcmuk changed to intel_openmp_bug, setting of host on lcmuk by
64	! cpp-directive removed
65	!
66	! 75 2007-03-22 09:54:05Z raasch
67	! uxrp, vynp eliminated,
68	! dirichlet/neumann changed to dirichlet/radiation, etc.,
69	! poisfft_init is only called if fft-solver is switched on
70	!
71	! RCS Log replace by Id keyword, revision history cleaned up
72	!
73	! Revision 1.28 2006/04/26 13:23:32 raasch
74	! lcmuk does not understand the !$ comment so a cpp-directive is required
75	!
76	! Revision 1.1 1997/07/24 11:15:09 raasch
77	! Initial revision
78	!
79	!
80	! Description:
81	! ------------
82	! Determination of the virtual processor topology (if not prescribed by the
83	! user)and computation of the grid point number and array bounds of the local
84	! domains.
85	!------------------------------------------------------------------------------!
86
87	USE control_parameters
88	USE fft_xy
89	USE grid_variables
90	USE indices
91	USE pegrid
92	USE poisfft_mod
93	USE poisfft_hybrid_mod
94	USE statistics
95	USE transpose_indices
96
97
98
99	IMPLICIT NONE
100
101	INTEGER :: gathered_size, i, id_inflow_l, id_recycling_l, ind(5), j, k, &
102	maximum_grid_level_l, mg_switch_to_pe0_level_l, mg_levels_x, &
103	mg_levels_y, mg_levels_z, nnx_y, nnx_z, nny_x, nny_z, nnz_x, &
104	nnz_y, numproc_sqr, nx_total, nxl_l, nxr_l, nyn_l, nys_l, &
105	nzb_l, nzt_l, omp_get_num_threads, subdomain_size
106
107	INTEGER, DIMENSION(:), ALLOCATABLE :: ind_all, nxlf, nxrf, nynf, nysf
108
109	INTEGER, DIMENSION(2) :: pdims_remote
110
111	LOGICAL :: found
112
113	!
114	!-- Get the number of OpenMP threads
115	!$OMP PARALLEL
116	#if defined( __intel_openmp_bug )
117	threads_per_task = omp_get_num_threads()
118	#else
119	!$ threads_per_task = omp_get_num_threads()
120	#endif
121	!$OMP END PARALLEL
122
123
124	#if defined( __parallel )
125
126	!
127	!-- Determine the processor topology or check it, if prescribed by the user
128	IF ( npex == -1 .AND. npey == -1 ) THEN
129
130	!
131	!-- Automatic determination of the topology
132	!-- The default on SMP- and cluster-hosts is a 1d-decomposition along x
133	IF ( host(1:3) == 'ibm' .OR. host(1:3) == 'nec' .OR. &
134	( host(1:2) == 'lc' .AND. host(3:5) /= 'sgi' .AND. &
135	host(3:4) /= 'xt' .AND. host(3:5) /= 'tit' ) .OR. &
136	host(1:3) == 'dec' ) THEN
137
138	pdims(1) = numprocs
139	pdims(2) = 1
140
141	ELSE
142
143	numproc_sqr = SQRT( REAL( numprocs ) )
144	pdims(1) = MAX( numproc_sqr , 1 )
145	DO WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
146	pdims(1) = pdims(1) - 1
147	ENDDO
148	pdims(2) = numprocs / pdims(1)
149
150	ENDIF
151
152	ELSEIF ( npex /= -1 .AND. npey /= -1 ) THEN
153
154	!
155	!-- Prescribed by user. Number of processors on the prescribed topology
156	!-- must be equal to the number of PEs available to the job
157	IF ( ( npex * npey ) /= numprocs ) THEN
158	WRITE( message_string, * ) 'number of PEs of the prescribed ', &
159	'topology (', npex*npey,') does not match & the number of ', &
160	'PEs available to the job (', numprocs, ')'
161	CALL message( 'init_pegrid', 'PA0221', 1, 2, 0, 6, 0 )
162	ENDIF
163	pdims(1) = npex
164	pdims(2) = npey
165
166	ELSE
167	!
168	!-- If the processor topology is prescribed by the user, the number of
169	!-- PEs must be given in both directions
170	message_string = 'if the processor topology is prescribed by the, ' // &
171	' user& both values of "npex" and "npey" must be given ' // &
172	'in the &NAMELIST-parameter file'
173	CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
174
175	ENDIF
176
177	!
178	!-- The hybrid solver can only be used in case of a 1d-decomposition along x
179	IF ( pdims(2) /= 1 .AND. psolver == 'poisfft_hybrid' ) THEN
180	message_string = 'psolver = "poisfft_hybrid" can only be' // &
181	'& used in case of a 1d-decomposition along x'
182	CALL message( 'init_pegrid', 'PA0223', 1, 2, 0, 6, 0 )
183	ENDIF
184
185	!
186	!-- For communication speedup, set barriers in front of collective
187	!-- communications by default on SGI-type systems
188	IF ( host(3:5) == 'sgi' ) collective_wait = .TRUE.
189
190	!
191	!-- If necessary, set horizontal boundary conditions to non-cyclic
192	IF ( .NOT. bc_lr_cyc ) cyclic(1) = .FALSE.
193	IF ( .NOT. bc_ns_cyc ) cyclic(2) = .FALSE.
194
195	!
196	!-- Create the virtual processor grid
197	CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
198	comm2d, ierr )
199	CALL MPI_COMM_RANK( comm2d, myid, ierr )
200	WRITE (myid_char,'(''_'',I4.4)') myid
201
202	CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
203	CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
204	CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
205
206	!
207	!-- Determine sub-topologies for transpositions
208	!-- Transposition from z to x:
209	remain_dims(1) = .TRUE.
210	remain_dims(2) = .FALSE.
211	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
212	CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
213	!
214	!-- Transposition from x to y
215	remain_dims(1) = .FALSE.
216	remain_dims(2) = .TRUE.
217	CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
218	CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
219
220
221	!
222	!-- Find a grid (used for array d) which will match the transposition demands
223	IF ( grid_matching == 'strict' ) THEN
224
225	nxa = nx; nya = ny; nza = nz
226
227	ELSE
228
229	found = .FALSE.
230	xn: DO nxa = nx, 2*nx
231	!
232	!-- Meet conditions for nx
233	IF ( MOD( nxa+1, pdims(1) ) /= 0 .OR. &
234	MOD( nxa+1, pdims(2) ) /= 0 ) CYCLE xn
235
236	yn: DO nya = ny, 2*ny
237	!
238	!-- Meet conditions for ny
239	IF ( MOD( nya+1, pdims(2) ) /= 0 .OR. &
240	MOD( nya+1, pdims(1) ) /= 0 ) CYCLE yn
241
242
243	zn: DO nza = nz, 2*nz
244	!
245	!-- Meet conditions for nz
246	IF ( ( MOD( nza, pdims(1) ) /= 0 .AND. pdims(1) /= 1 .AND. &
247	pdims(2) /= 1 ) .OR. &
248	( MOD( nza, pdims(2) ) /= 0 .AND. dt_dosp /= 9999999.9 &
249	) ) THEN
250	CYCLE zn
251	ELSE
252	found = .TRUE.
253	EXIT xn
254	ENDIF
255
256	ENDDO zn
257
258	ENDDO yn
259
260	ENDDO xn
261
262	IF ( .NOT. found ) THEN
263	message_string = 'no matching grid for transpositions found'
264	CALL message( 'init_pegrid', 'PA0224', 1, 2, 0, 6, 0 )
265	ENDIF
266
267	ENDIF
268
269	!
270	!-- Calculate array bounds in x-direction for every PE.
271	!-- The last PE along x may get less grid points than the others
272	ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
273	nysf(0:pdims(2)-1), nnx_pe(0:pdims(1)-1), nny_pe(0:pdims(2)-1) )
274
275	IF ( MOD( nxa+1 , pdims(1) ) /= 0 ) THEN
276	WRITE( message_string, * ) 'x-direction: gridpoint number (',nx+1,') ',&
277	'is not an& integral divisor of the number ', &
278	'processors (', pdims(1),')'
279	CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
280	ELSE
281	nnx = ( nxa + 1 ) / pdims(1)
282	IF ( nnx*pdims(1) - ( nx + 1) > nnx ) THEN
283	WRITE( message_string, * ) 'x-direction: nx does not match the', &
284	'requirements given by the number of PEs &used', &
285	'& please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
286	- ( nx + 1 ) ) ), ' instead of nx =', nx
287	CALL message( 'init_pegrid', 'PA0226', 1, 2, 0, 6, 0 )
288	ENDIF
289	ENDIF
290
291	!
292	!-- Left and right array bounds, number of gridpoints
293	DO i = 0, pdims(1)-1
294	nxlf(i) = i * nnx
295	nxrf(i) = ( i + 1 ) * nnx - 1
296	nnx_pe(i) = MIN( nx, nxrf(i) ) - nxlf(i) + 1
297	ENDDO
298
299	!
300	!-- Calculate array bounds in y-direction for every PE.
301	IF ( MOD( nya+1 , pdims(2) ) /= 0 ) THEN
302	WRITE( message_string, * ) 'y-direction: gridpoint number (',ny+1,') ', &
303	'is not an& integral divisor of the number of', &
304	'processors (', pdims(2),')'
305	CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
306	ELSE
307	nny = ( nya + 1 ) / pdims(2)
308	IF ( nny*pdims(2) - ( ny + 1) > nny ) THEN
309	WRITE( message_string, * ) 'y-direction: ny does not match the', &
310	'requirements given by the number of PEs &used ', &
311	'& please use ny = ', ny - ( pdims(2) - ( nnx*pdims(2) &
312	- ( ny + 1 ) ) ), ' instead of ny =', ny
313	CALL message( 'init_pegrid', 'PA0228', 1, 2, 0, 6, 0 )
314	ENDIF
315	ENDIF
316
317	!
318	!-- South and north array bounds
319	DO j = 0, pdims(2)-1
320	nysf(j) = j * nny
321	nynf(j) = ( j + 1 ) * nny - 1
322	nny_pe(j) = MIN( ny, nynf(j) ) - nysf(j) + 1
323	ENDDO
324
325	!
326	!-- Local array bounds of the respective PEs
327	nxl = nxlf(pcoord(1))
328	nxra = nxrf(pcoord(1))
329	nxr = MIN( nx, nxra )
330	nys = nysf(pcoord(2))
331	nyna = nynf(pcoord(2))
332	nyn = MIN( ny, nyna )
333	nzb = 0
334	nzta = nza
335	nzt = MIN( nz, nzta )
336	nnz = nza
337
338	!
339	!-- Set switches to define if the PE is situated at the border of the virtual
340	!-- processor grid
341	IF ( nxl == 0 ) left_border_pe = .TRUE.
342	IF ( nxr == nx ) right_border_pe = .TRUE.
343	IF ( nys == 0 ) south_border_pe = .TRUE.
344	IF ( nyn == ny ) north_border_pe = .TRUE.
345
346	!
347	!-- Calculate array bounds and gridpoint numbers for the transposed arrays
348	!-- (needed in the pressure solver)
349	!-- For the transposed arrays, cyclic boundaries as well as top and bottom
350	!-- boundaries are omitted, because they are obstructive to the transposition
351
352	!
353	!-- 1. transposition z --> x
354	!-- This transposition is not neccessary in case of a 1d-decomposition along x,
355	!-- except that the uptream-spline method is switched on
356	IF ( pdims(2) /= 1 .OR. momentum_advec == 'ups-scheme' .OR. &
357	scalar_advec == 'ups-scheme' ) THEN
358
359	IF ( pdims(2) == 1 .AND. ( momentum_advec == 'ups-scheme' .OR. &
360	scalar_advec == 'ups-scheme' ) ) THEN
361	message_string = '1d-decomposition along x ' // &
362	'chosen but nz restrictions may occur' // &
363	'& since ups-scheme is activated'
364	CALL message( 'init_pegrid', 'PA0229', 0, 1, 0, 6, 0 )
365	ENDIF
366	nys_x = nys
367	nyn_xa = nyna
368	nyn_x = nyn
369	nny_x = nny
370	IF ( MOD( nza , pdims(1) ) /= 0 ) THEN
371	WRITE( message_string, * ) 'transposition z --> x:', &
372	'&nz=',nz,' is not an integral divisior of pdims(1)=', &
373	pdims(1)
374	CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
375	ENDIF
376	nnz_x = nza / pdims(1)
377	nzb_x = 1 + myidx * nnz_x
378	nzt_xa = ( myidx + 1 ) * nnz_x
379	nzt_x = MIN( nzt, nzt_xa )
380
381	sendrecvcount_zx = nnx * nny * nnz_x
382
383	ELSE
384	!
385	!--- Setting of dummy values because otherwise variables are undefined in
386	!--- the next step x --> y
387	!--- WARNING: This case has still to be clarified!!!!!!!!!!!!
388	nnz_x = 1
389	nzb_x = 1
390	nzt_xa = 1
391	nzt_x = 1
392	nny_x = nny
393
394	ENDIF
395
396	!
397	!-- 2. transposition x --> y
398	nnz_y = nnz_x
399	nzb_y = nzb_x
400	nzt_ya = nzt_xa
401	nzt_y = nzt_x
402	IF ( MOD( nxa+1 , pdims(2) ) /= 0 ) THEN
403	WRITE( message_string, * ) 'transposition x --> y:', &
404	'&nx+1=',nx+1,' is not an integral divisor of ',&
405	'pdims(2)=',pdims(2)
406	CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
407	ENDIF
408	nnx_y = (nxa+1) / pdims(2)
409	nxl_y = myidy * nnx_y
410	nxr_ya = ( myidy + 1 ) * nnx_y - 1
411	nxr_y = MIN( nx, nxr_ya )
412
413	sendrecvcount_xy = nnx_y * nny_x * nnz_y
414
415	!
416	!-- 3. transposition y --> z (ELSE: x --> y in case of 1D-decomposition
417	!-- along x)
418	IF ( pdims(2) /= 1 .OR. momentum_advec == 'ups-scheme' .OR. &
419	scalar_advec == 'ups-scheme' ) THEN
420	!
421	!-- y --> z
422	!-- This transposition is not neccessary in case of a 1d-decomposition
423	!-- along x, except that the uptream-spline method is switched on
424	nnx_z = nnx_y
425	nxl_z = nxl_y
426	nxr_za = nxr_ya
427	nxr_z = nxr_y
428	IF ( MOD( nya+1 , pdims(1) ) /= 0 ) THEN
429	WRITE( message_string, * ) 'transposition y --> z:', &
430	'& ny+1=',ny+1,' is not an integral divisor of ',&
431	'pdims(1)=',pdims(1)
432	CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
433	ENDIF
434	nny_z = (nya+1) / pdims(1)
435	nys_z = myidx * nny_z
436	nyn_za = ( myidx + 1 ) * nny_z - 1
437	nyn_z = MIN( ny, nyn_za )
438
439	sendrecvcount_yz = nnx_y * nny_z * nnz_y
440
441	ELSE
442	!
443	!-- x --> y. This condition must be fulfilled for a 1D-decomposition along x
444	IF ( MOD( nya+1 , pdims(1) ) /= 0 ) THEN
445	WRITE( message_string, * ) 'transposition x --> y:', &
446	'& ny+1=',ny+1,' is not an integral divisor of ',&
447	'pdims(1)=',pdims(1)
448	CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
449	ENDIF
450
451	ENDIF
452
453	!
454	!-- Indices for direct transpositions z --> y (used for calculating spectra)
455	IF ( dt_dosp /= 9999999.9 ) THEN
456	IF ( MOD( nza, pdims(2) ) /= 0 ) THEN
457	WRITE( message_string, * ) 'direct transposition z --> y (needed ', &
458	'for spectra):& nz=',nz,' is not an integral divisor of ',&
459	'pdims(2)=',pdims(2)
460	CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
461	ELSE
462	nxl_yd = nxl
463	nxr_yda = nxra
464	nxr_yd = nxr
465	nzb_yd = 1 + myidy * ( nza / pdims(2) )
466	nzt_yda = ( myidy + 1 ) * ( nza / pdims(2) )
467	nzt_yd = MIN( nzt, nzt_yda )
468
469	sendrecvcount_zyd = nnx * nny * ( nza / pdims(2) )
470	ENDIF
471	ENDIF
472
473	!
474	!-- Indices for direct transpositions y --> x (they are only possible in case
475	!-- of a 1d-decomposition along x)
476	IF ( pdims(2) == 1 ) THEN
477	nny_x = nny / pdims(1)
478	nys_x = myid * nny_x
479	nyn_xa = ( myid + 1 ) * nny_x - 1
480	nyn_x = MIN( ny, nyn_xa )
481	nzb_x = 1
482	nzt_xa = nza
483	nzt_x = nz
484	sendrecvcount_xy = nnx * nny_x * nza
485	ENDIF
486
487	!
488	!-- Indices for direct transpositions x --> y (they are only possible in case
489	!-- of a 1d-decomposition along y)
490	IF ( pdims(1) == 1 ) THEN
491	nnx_y = nnx / pdims(2)
492	nxl_y = myid * nnx_y
493	nxr_ya = ( myid + 1 ) * nnx_y - 1
494	nxr_y = MIN( nx, nxr_ya )
495	nzb_y = 1
496	nzt_ya = nza
497	nzt_y = nz
498	sendrecvcount_xy = nnx_y * nny * nza
499	ENDIF
500
501	!
502	!-- Arrays for storing the array bounds are needed any more
503	DEALLOCATE( nxlf , nxrf , nynf , nysf )
504
505	!
506	!-- Collect index bounds from other PEs (to be written to restart file later)
507	ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
508
509	IF ( myid == 0 ) THEN
510
511	hor_index_bounds(1,0) = nxl
512	hor_index_bounds(2,0) = nxr
513	hor_index_bounds(3,0) = nys
514	hor_index_bounds(4,0) = nyn
515
516	!
517	!-- Receive data from all other PEs
518	DO i = 1, numprocs-1
519	CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
520	ierr )
521	hor_index_bounds(:,i) = ibuf(1:4)
522	ENDDO
523
524	ELSE
525	!
526	!-- Send index bounds to PE0
527	ibuf(1) = nxl
528	ibuf(2) = nxr
529	ibuf(3) = nys
530	ibuf(4) = nyn
531	CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
532
533	ENDIF
534
535	#if defined( __print )
536	!
537	!-- Control output
538	IF ( myid == 0 ) THEN
539	PRINT, ' processor topology *'
540	PRINT*, ' '
541	PRINT*, 'myid pcoord left right south north idx idy nxl: nxr',&
542	&' nys: nyn'
543	PRINT*, '------------------------------------------------------------',&
544	&'-----------'
545	WRITE (*,1000) 0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
546	myidx, myidy, nxl, nxr, nys, nyn
547	1000 FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
548	2(2X,I4,':',I4))
549
550	!
551	!-- Receive data from the other PEs
552	DO i = 1,numprocs-1
553	CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
554	ierr )
555	WRITE (*,1000) i, ( ibuf(j) , j = 1,12 )
556	ENDDO
557	ELSE
558
559	!
560	!-- Send data to PE0
561	ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
562	ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
563	ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
564	ibuf(12) = nyn
565	CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )
566	ENDIF
567	#endif
568
569	#if defined( __parallel )
570	#if defined( __mpi2 )
571	!
572	!-- In case of coupled runs, get the port name on PE0 of the atmosphere model
573	!-- and pass it to PE0 of the ocean model
574	IF ( myid == 0 ) THEN
575
576	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
577
578	CALL MPI_OPEN_PORT( MPI_INFO_NULL, port_name, ierr )
579
580	CALL MPI_PUBLISH_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, &
581	ierr )
582
583	!
584	!-- Write a flag file for the ocean model and the other atmosphere
585	!-- processes.
586	!-- There seems to be a bug in MPICH2 which causes hanging processes
587	!-- in case that execution of LOOKUP_NAME is continued too early
588	!-- (i.e. before the port has been created)
589	OPEN( 90, FILE='COUPLING_PORT_OPENED', FORM='FORMATTED' )
590	WRITE ( 90, '(''TRUE'')' )
591	CLOSE ( 90 )
592
593	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
594
595	!
596	!-- Continue only if the atmosphere model has created the port.
597	!-- There seems to be a bug in MPICH2 which causes hanging processes
598	!-- in case that execution of LOOKUP_NAME is continued too early
599	!-- (i.e. before the port has been created)
600	INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
601	DO WHILE ( .NOT. found )
602	INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
603	ENDDO
604
605	CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
606
607	ENDIF
608
609	ENDIF
610
611	!
612	!-- In case of coupled runs, establish the connection between the atmosphere
613	!-- and the ocean model and define the intercommunicator (comm_inter)
614	CALL MPI_BARRIER( comm2d, ierr )
615	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
616
617	CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
618	comm_inter, ierr )
619	coupling_mode_remote = 'ocean_to_atmosphere'
620
621	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
622
623	CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
624	comm_inter, ierr )
625	coupling_mode_remote = 'atmosphere_to_ocean'
626
627	ENDIF
628	#endif
629
630	!
631	!-- Determine the number of ghost point layers
632	IF ( scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' ) THEN
633	nbgp = 3
634	ELSE
635	nbgp = 1
636	ENDIF
637
638	!
639	!-- Create a new MPI derived datatype for the exchange of surface (xy) data,
640	!-- which is needed for coupled atmosphere-ocean runs.
641	!-- First, calculate number of grid points of an xy-plane.
642	ngp_xy = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
643	CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
644	CALL MPI_TYPE_COMMIT( type_xy, ierr )
645
646	IF ( TRIM( coupling_mode ) /= 'uncoupled' ) THEN
647
648	!
649	!-- Pass the number of grid points of the atmosphere model to
650	!-- the ocean model and vice versa
651	IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN
652
653	nx_a = nx
654	ny_a = ny
655
656	IF ( myid == 0 ) THEN
657
658	CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter, &
659	ierr )
660	CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter, &
661	ierr )
662	CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, &
663	ierr )
664	CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter, &
665	status, ierr )
666	CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter, &
667	status, ierr )
668	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &
669	comm_inter, status, ierr )
670	ENDIF
671
672	CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
673	CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr )
674	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
675
676	ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN
677
678	nx_o = nx
679	ny_o = ny
680
681	IF ( myid == 0 ) THEN
682
683	CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
684	ierr )
685	CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
686	ierr )
687	CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, &
688	status, ierr )
689	CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
690	CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
691	CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
692	ENDIF
693
694	CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
695	CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr)
696	CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
697
698	ENDIF
699
700	ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
701	ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
702
703	!
704	!-- Determine if the horizontal grid and the number of PEs in ocean and
705	!-- atmosphere is same or not
706	IF ( nx_o == nx_a .AND. ny_o == ny_a .AND. &
707	pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
708	THEN
709	coupling_topology = 0
710	ELSE
711	coupling_topology = 1
712	ENDIF
713
714	!
715	!-- Determine the target PEs for the exchange between ocean and
716	!-- atmosphere (comm2d)
717	IF ( coupling_topology == 0 ) THEN
718	!
719	!-- In case of identical topologies, every atmosphere PE has exactly one
720	!-- ocean PE counterpart and vice versa
721	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
722	target_id = myid + numprocs
723	ELSE
724	target_id = myid
725	ENDIF
726
727	ELSE
728	!
729	!-- In case of nonequivalent topology in ocean and atmosphere only for
730	!-- PE0 in ocean and PE0 in atmosphere a target_id is needed, since
731	!-- data echxchange between ocean and atmosphere will be done only
732	!-- between these PEs.
733	IF ( myid == 0 ) THEN
734
735	IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
736	target_id = numprocs
737	ELSE
738	target_id = 0
739	ENDIF
740
741	ENDIF
742
743	ENDIF
744
745	ENDIF
746
747
748	#endif
749
750	#else
751
752	!
753	!-- Array bounds when running on a single PE (respectively a non-parallel
754	!-- machine)
755	nxl = 0
756	nxr = nx
757	nxra = nx
758	nnx = nxr - nxl + 1
759	nys = 0
760	nyn = ny
761	nyna = ny
762	nny = nyn - nys + 1
763	nzb = 0
764	nzt = nz
765	nzta = nz
766	nnz = nz
767
768	ALLOCATE( hor_index_bounds(4,0:0) )
769	hor_index_bounds(1,0) = nxl
770	hor_index_bounds(2,0) = nxr
771	hor_index_bounds(3,0) = nys
772	hor_index_bounds(4,0) = nyn
773
774	!
775	!-- Array bounds for the pressure solver (in the parallel code, these bounds
776	!-- are the ones for the transposed arrays)
777	nys_x = nys
778	nyn_x = nyn
779	nyn_xa = nyn
780	nzb_x = nzb + 1
781	nzt_x = nzt
782	nzt_xa = nzt
783
784	nxl_y = nxl
785	nxr_y = nxr
786	nxr_ya = nxr
787	nzb_y = nzb + 1
788	nzt_y = nzt
789	nzt_ya = nzt
790
791	nxl_z = nxl
792	nxr_z = nxr
793	nxr_za = nxr
794	nys_z = nys
795	nyn_z = nyn
796	nyn_za = nyn
797
798	#endif
799
800	!
801	!-- Calculate number of grid levels necessary for the multigrid poisson solver
802	!-- as well as the gridpoint indices on each level
803	IF ( psolver == 'multigrid' ) THEN
804
805	!
806	!-- First calculate number of possible grid levels for the subdomains
807	mg_levels_x = 1
808	mg_levels_y = 1
809	mg_levels_z = 1
810
811	i = nnx
812	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
813	i = i / 2
814	mg_levels_x = mg_levels_x + 1
815	ENDDO
816
817	j = nny
818	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
819	j = j / 2
820	mg_levels_y = mg_levels_y + 1
821	ENDDO
822
823	k = nz ! do not use nnz because it might be > nz due to transposition
824	! requirements
825	DO WHILE ( MOD( k, 2 ) == 0 .AND. k /= 2 )
826	k = k / 2
827	mg_levels_z = mg_levels_z + 1
828	ENDDO
829
830	maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
831
832	!
833	!-- Find out, if the total domain allows more levels. These additional
834	!-- levels are identically processed on all PEs.
835	IF ( numprocs > 1 .AND. mg_switch_to_pe0_level /= -1 ) THEN
836
837	IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) ) THEN
838
839	mg_switch_to_pe0_level_l = maximum_grid_level
840
841	mg_levels_x = 1
842	mg_levels_y = 1
843
844	i = nx+1
845	DO WHILE ( MOD( i, 2 ) == 0 .AND. i /= 2 )
846	i = i / 2
847	mg_levels_x = mg_levels_x + 1
848	ENDDO
849
850	j = ny+1
851	DO WHILE ( MOD( j, 2 ) == 0 .AND. j /= 2 )
852	j = j / 2
853	mg_levels_y = mg_levels_y + 1
854	ENDDO
855
856	maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
857
858	IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l ) THEN
859	mg_switch_to_pe0_level_l = maximum_grid_level_l - &
860	mg_switch_to_pe0_level_l + 1
861	ELSE
862	mg_switch_to_pe0_level_l = 0
863	ENDIF
864
865	ELSE
866
867	mg_switch_to_pe0_level_l = 0
868	maximum_grid_level_l = maximum_grid_level
869
870	ENDIF
871
872	!
873	!-- Use switch level calculated above only if it is not pre-defined
874	!-- by user
875	IF ( mg_switch_to_pe0_level == 0 ) THEN
876
877	IF ( mg_switch_to_pe0_level_l /= 0 ) THEN
878	mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
879	maximum_grid_level = maximum_grid_level_l
880	ENDIF
881
882	ELSE
883	!
884	!-- Check pre-defined value and reset to default, if neccessary
885	IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l .OR. &
886	mg_switch_to_pe0_level >= maximum_grid_level_l ) THEN
887	message_string = 'mg_switch_to_pe0_level ' // &
888	'out of range and reset to default (=0)'
889	CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
890	mg_switch_to_pe0_level = 0
891	ELSE
892	!
893	!-- Use the largest number of possible levels anyway and recalculate
894	!-- the switch level to this largest number of possible values
895	maximum_grid_level = maximum_grid_level_l
896
897	ENDIF
898
899	ENDIF
900
901	ENDIF
902
903	ALLOCATE( grid_level_count(maximum_grid_level), &
904	nxl_mg(maximum_grid_level), nxr_mg(maximum_grid_level), &
905	nyn_mg(maximum_grid_level), nys_mg(maximum_grid_level), &
906	nzt_mg(maximum_grid_level) )
907
908	grid_level_count = 0
909	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
910
911	DO i = maximum_grid_level, 1 , -1
912
913	IF ( i == mg_switch_to_pe0_level ) THEN
914	#if defined( __parallel )
915	!
916	!-- Save the grid size of the subdomain at the switch level, because
917	!-- it is needed in poismg.
918	ind(1) = nxl_l; ind(2) = nxr_l
919	ind(3) = nys_l; ind(4) = nyn_l
920	ind(5) = nzt_l
921	ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
922	CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
923	MPI_INTEGER, comm2d, ierr )
924	DO j = 0, numprocs-1
925	DO k = 1, 5
926	mg_loc_ind(k,j) = ind_all(k+j*5)
927	ENDDO
928	ENDDO
929	DEALLOCATE( ind_all )
930	!
931	!-- Calculate the grid size of the total domain
932	nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
933	nxl_l = 0
934	nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
935	nys_l = 0
936	!
937	!-- The size of this gathered array must not be larger than the
938	!-- array tend, which is used in the multigrid scheme as a temporary
939	!-- array
940	subdomain_size = ( nxr - nxl + 3 ) * ( nyn - nys + 3 ) * &
941	( nzt - nzb + 2 )
942	gathered_size = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
943	( nzt_l - nzb + 2 )
944
945	IF ( gathered_size > subdomain_size ) THEN
946	message_string = 'not enough memory for storing ' // &
947	'gathered multigrid data on PE0'
948	CALL message( 'init_pegrid', 'PA0236', 1, 2, 0, 6, 0 )
949	ENDIF
950	#else
951	message_string = 'multigrid gather/scatter impossible ' // &
952	'in non parallel mode'
953	CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
954	#endif
955	ENDIF
956
957	nxl_mg(i) = nxl_l
958	nxr_mg(i) = nxr_l
959	nys_mg(i) = nys_l
960	nyn_mg(i) = nyn_l
961	nzt_mg(i) = nzt_l
962
963	nxl_l = nxl_l / 2
964	nxr_l = nxr_l / 2
965	nys_l = nys_l / 2
966	nyn_l = nyn_l / 2
967	nzt_l = nzt_l / 2
968	ENDDO
969
970	ELSE
971
972	maximum_grid_level = 0
973
974	ENDIF
975
976	grid_level = maximum_grid_level
977
978	#if defined( __parallel )
979	!
980	!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
981	ngp_y = nyn - nys + 1 + 2 * nbgp
982
983	!
984	!-- Define new MPI derived datatypes for the exchange of ghost points in
985	!-- x- and y-direction for 2D-arrays (line)
986	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, &
987	ierr )
988	CALL MPI_TYPE_COMMIT( type_x, ierr )
989	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, &
990	type_x_int, ierr )
991	CALL MPI_TYPE_COMMIT( type_x_int, ierr )
992
993	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
994	CALL MPI_TYPE_COMMIT( type_y, ierr )
995	CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int, ierr )
996	CALL MPI_TYPE_COMMIT( type_y_int, ierr )
997
998
999	!
1000	!-- Calculate gridpoint numbers for the exchange of ghost points along x
1001	!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
1002	!-- exchange of ghost points in y-direction (xz-plane).
1003	!-- Do these calculations for the model grid and (if necessary) also
1004	!-- for the coarser grid levels used in the multigrid method
1005	ALLOCATE ( ngp_yz(0:maximum_grid_level), type_xz(0:maximum_grid_level),&
1006	type_yz(0:maximum_grid_level) )
1007
1008	nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
1009
1010	!
1011	!-- Discern between the model grid, which needs nbgp ghost points and
1012	!-- grid levels for the multigrid scheme. In the latter case only one
1013	!-- ghost point is necessary.
1014	!-- First definition of MPI-datatypes for exchange of ghost layers on normal
1015	!-- grid. The following loop is needed for data exchange in poismg.f90.
1016	!
1017	!-- Determine number of grid points of yz-layer for exchange
1018	ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
1019
1020	!
1021	!-- Define an MPI-datatype for the exchange of left/right boundaries.
1022	!-- Although data are contiguous in physical memory (which does not
1023	!-- necessarily require an MPI-derived datatype), the data exchange between
1024	!-- left and right PE's using the MPI-derived type is 10% faster than without.
1025	CALL MPI_TYPE_VECTOR( nxr-nxl+1+2nbgp, nbgp(nzt-nzb+2), ngp_yz(0), &
1026	MPI_REAL, type_xz(0), ierr )
1027	CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
1028
1029	CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), &
1030	ierr )
1031	CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
1032
1033	!
1034	!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
1035	IF ( psolver == 'multigrid' ) THEN
1036	!
1037	!-- Definition of MPI-datatyoe as above, but only 1 ghost level is used
1038	DO i = maximum_grid_level, 1 , -1
1039
1040	ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
1041
1042	CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
1043	MPI_REAL, type_xz(i), ierr )
1044	CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
1045
1046	CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), &
1047	ierr )
1048	CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
1049
1050	nxl_l = nxl_l / 2
1051	nxr_l = nxr_l / 2
1052	nys_l = nys_l / 2
1053	nyn_l = nyn_l / 2
1054	nzt_l = nzt_l / 2
1055
1056	ENDDO
1057
1058	ENDIF
1059	#endif
1060
1061	#if defined( __parallel )
1062	!
1063	!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
1064	!-- horizontal boundary conditions.
1065	IF ( pleft == MPI_PROC_NULL ) THEN
1066	IF ( bc_lr_dirrad ) THEN
1067	inflow_l = .TRUE.
1068	ELSEIF ( bc_lr_raddir ) THEN
1069	outflow_l = .TRUE.
1070	ENDIF
1071	ENDIF
1072
1073	IF ( pright == MPI_PROC_NULL ) THEN
1074	IF ( bc_lr_dirrad ) THEN
1075	outflow_r = .TRUE.
1076	ELSEIF ( bc_lr_raddir ) THEN
1077	inflow_r = .TRUE.
1078	ENDIF
1079	ENDIF
1080
1081	IF ( psouth == MPI_PROC_NULL ) THEN
1082	IF ( bc_ns_dirrad ) THEN
1083	outflow_s = .TRUE.
1084	ELSEIF ( bc_ns_raddir ) THEN
1085	inflow_s = .TRUE.
1086	ENDIF
1087	ENDIF
1088
1089	IF ( pnorth == MPI_PROC_NULL ) THEN
1090	IF ( bc_ns_dirrad ) THEN
1091	inflow_n = .TRUE.
1092	ELSEIF ( bc_ns_raddir ) THEN
1093	outflow_n = .TRUE.
1094	ENDIF
1095	ENDIF
1096
1097	!
1098	!-- Broadcast the id of the inflow PE
1099	IF ( inflow_l ) THEN
1100	id_inflow_l = myidx
1101	ELSE
1102	id_inflow_l = 0
1103	ENDIF
1104	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1105	CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, &
1106	comm1dx, ierr )
1107
1108	!
1109	!-- Broadcast the id of the recycling plane
1110	!-- WARNING: needs to be adjusted in case of inflows other than from left side!
1111	IF ( ( recycling_width / dx ) >= nxl .AND. &
1112	( recycling_width / dx ) <= nxr ) THEN
1113	id_recycling_l = myidx
1114	ELSE
1115	id_recycling_l = 0
1116	ENDIF
1117	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1118	CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, &
1119	comm1dx, ierr )
1120
1121	#else
1122	IF ( bc_lr_dirrad ) THEN
1123	inflow_l = .TRUE.
1124	outflow_r = .TRUE.
1125	ELSEIF ( bc_lr_raddir ) THEN
1126	outflow_l = .TRUE.
1127	inflow_r = .TRUE.
1128	ENDIF
1129
1130	IF ( bc_ns_dirrad ) THEN
1131	inflow_n = .TRUE.
1132	outflow_s = .TRUE.
1133	ELSEIF ( bc_ns_raddir ) THEN
1134	outflow_n = .TRUE.
1135	inflow_s = .TRUE.
1136	ENDIF
1137	#endif
1138	!
1139	!-- At the outflow, u or v, respectively, have to be calculated for one more
1140	!-- grid point.
1141	IF ( outflow_l ) THEN
1142	nxlu = nxl + 1
1143	ELSE
1144	nxlu = nxl
1145	ENDIF
1146	IF ( outflow_s ) THEN
1147	nysv = nys + 1
1148	ELSE
1149	nysv = nys
1150	ENDIF
1151
1152	IF ( psolver == 'poisfft_hybrid' ) THEN
1153	CALL poisfft_hybrid_ini
1154	ELSEIF ( psolver == 'poisfft' ) THEN
1155	CALL poisfft_init
1156	ENDIF
1157
1158	!
1159	!-- Allocate wall flag arrays used in the multigrid solver
1160	IF ( psolver == 'multigrid' ) THEN
1161
1162	DO i = maximum_grid_level, 1, -1
1163
1164	SELECT CASE ( i )
1165
1166	CASE ( 1 )
1167	ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1, &
1168	nys_mg(i)-1:nyn_mg(i)+1, &
1169	nxl_mg(i)-1:nxr_mg(i)+1) )
1170
1171	CASE ( 2 )
1172	ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1, &
1173	nys_mg(i)-1:nyn_mg(i)+1, &
1174	nxl_mg(i)-1:nxr_mg(i)+1) )
1175
1176	CASE ( 3 )
1177	ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1, &
1178	nys_mg(i)-1:nyn_mg(i)+1, &
1179	nxl_mg(i)-1:nxr_mg(i)+1) )
1180
1181	CASE ( 4 )
1182	ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1, &
1183	nys_mg(i)-1:nyn_mg(i)+1, &
1184	nxl_mg(i)-1:nxr_mg(i)+1) )
1185
1186	CASE ( 5 )
1187	ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1, &
1188	nys_mg(i)-1:nyn_mg(i)+1, &
1189	nxl_mg(i)-1:nxr_mg(i)+1) )
1190
1191	CASE ( 6 )
1192	ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1, &
1193	nys_mg(i)-1:nyn_mg(i)+1, &
1194	nxl_mg(i)-1:nxr_mg(i)+1) )
1195
1196	CASE ( 7 )
1197	ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1, &
1198	nys_mg(i)-1:nyn_mg(i)+1, &
1199	nxl_mg(i)-1:nxr_mg(i)+1) )
1200
1201	CASE ( 8 )
1202	ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1, &
1203	nys_mg(i)-1:nyn_mg(i)+1, &
1204	nxl_mg(i)-1:nxr_mg(i)+1) )
1205
1206	CASE ( 9 )
1207	ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1, &
1208	nys_mg(i)-1:nyn_mg(i)+1, &
1209	nxl_mg(i)-1:nxr_mg(i)+1) )
1210
1211	CASE ( 10 )
1212	ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1, &
1213	nys_mg(i)-1:nyn_mg(i)+1, &
1214	nxl_mg(i)-1:nxr_mg(i)+1) )
1215
1216	CASE DEFAULT
1217	message_string = 'more than 10 multigrid levels'
1218	CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
1219
1220	END SELECT
1221
1222	ENDDO
1223
1224	ENDIF
1225
1226	END SUBROUTINE init_pegrid

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |