Home

Context Navigation

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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago
Code annotations made doxygen readable
Property svn:keywords set to `Id`
File size: 42.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |