Home

Context Navigation

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

Last change on this file since 687 was 669, checked in by gryschka, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 40.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |