Home

Context Navigation

source: palm/trunk/SOURCE/poismg.f90 @ 1328

Last change on this file since 1328 was 1323, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 64.5 KB

Line
1	SUBROUTINE poismg( r )
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Attention: Loop unrolling and cache optimization in SOR-Red/Black method
21	! still does not give the expected speedup! Further work required.
22	!
23	! Current revisions:
24	! -----------------
25	!
26	!
27	! Former revisions:
28	! -----------------
29	! $Id: poismg.f90 1323 2014-03-20 17:09:54Z maronga $
30	!
31	! 1322 2014-03-20 16:38:49Z raasch
32	! REAL constants defined as wp-kind
33	!
34	! 1320 2014-03-20 08:40:49Z raasch
35	! ONLY-attribute added to USE-statements,
36	! kind-parameters added to all INTEGER and REAL declaration statements,
37	! kinds are defined in new module kinds,
38	! old module precision_kind is removed,
39	! revision history before 2012 removed,
40	! comment fields (!:) to be used for variable explanations added to
41	! all variable declaration statements
42	!
43	! 1318 2014-03-17 13:35:16Z raasch
44	! module interfaces removed
45	!
46	! 1159 2013-05-21 11:58:22Z fricke
47	! bc_lr/ns_dirneu/neudir removed
48	!
49	! 1092 2013-02-02 11:24:22Z raasch
50	! unused variables removed
51	!
52	! 1056 2012-11-16 15:28:04Z raasch
53	! Bugfix: all ghost points have to be used for allocating p3
54	! arrays p2, f2, and f2_l changed from allocatable to automatic
55	!
56	! 1036 2012-10-22 13:43:42Z raasch
57	! code put under GPL (PALM 3.9)
58	!
59	! 996 2012-09-07 10:41:47Z raasch
60	! little reformatting
61	!
62	! 978 2012-08-09 08:28:32Z fricke
63	! bc_lr/ns_dirneu/neudir added
64	!
65	! 880 2012-04-13 06:28:59Z raasch
66	! Bugfix: preprocessor statements for parallel execution added
67	!
68	! 778 2011-11-07 14:18:25Z fricke
69	! Allocation of p3 changes when multigrid is used and the collected field on PE0
70	! has more grid points than the subdomain of an PE.
71	!
72	! 707 2011-03-29 11:39:40Z raasch
73	! p_loc is used instead of p in the main routine (poismg).
74	! On coarse grid levels, gathered data are identically processed on all PEs
75	! (before, on PE0 only), so that the subsequent scattering of data is not
76	! neccessary any more.
77	! bc_lr/ns replaced by bc_lr/ns_cyc/dirrad/raddir
78	! Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
79	! resid and restrict. They were missed before which may have led to
80	! unpredictable results.
81	!
82	! 667 2010-12-23 12:06:00Z suehring/gryschka
83	! Calls of exchange_horiz are modified.
84	!
85	! 622 2010-12-10 08:08:13Z raasch
86	! optional barriers included in order to speed up collective operations
87	!
88	! 257 2009-03-11 15:17:42Z heinze
89	! Output of messages replaced by message handling routine.
90	!
91	! 181 2008-07-30 07:07:47Z raasch
92	! Bugfix: grid_level+1 has to be used in restrict for flags-array
93	!
94	! 114 2007-10-10 00:03:15Z raasch
95	! Boundary conditions at walls are implicitly set using flag arrays. Only
96	! Neumann BC is allowed. Upper walls are still not realized.
97	! Bottom and top BCs for array f_mg in restrict removed because boundary
98	! values are not needed (right hand side of SOR iteration).
99	!
100	! 75 2007-03-22 09:54:05Z raasch
101	! 2nd+3rd argument removed from exchange horiz
102	!
103	! RCS Log replace by Id keyword, revision history cleaned up
104	!
105	! Revision 1.6 2005/03/26 20:55:54 raasch
106	! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
107	! routine prolong simplified (one call of exchange_horiz spared)
108	!
109	! Revision 1.1 2001/07/20 13:10:51 raasch
110	! Initial revision
111	!
112	!
113	! Description:
114	! ------------
115	! Solves the Poisson equation for the perturbation pressure with a multigrid
116	! V- or W-Cycle scheme.
117	!
118	! This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
119	! September 2000 - July 2001.
120	!------------------------------------------------------------------------------!
121
122	USE arrays_3d, &
123	ONLY: d, p_loc
124
125	USE control_parameters, &
126	ONLY: gathered_size, grid_level, grid_level_count, &
127	maximum_grid_level, message_string, mgcycles, mg_cycles, &
128	mg_switch_to_pe0_level, residual_limit, subdomain_size
129
130	USE cpulog, &
131	ONLY: cpu_log, log_point_s
132
133	USE indices, &
134	ONLY: nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, nys, nysg, nys_mg, nyn, &
135	nyng, nyn_mg, nzb, nzt, nzt_mg
136
137	USE kinds
138
139	USE pegrid
140
141	IMPLICIT NONE
142
143	REAL(wp) :: maxerror !:
144	REAL(wp) :: maximum_mgcycles !:
145	REAL(wp) :: residual_norm !:
146
147	REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r !:
148
149	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p3 !:
150
151
152	CALL cpu_log( log_point_s(29), 'poismg', 'start' )
153	!
154	!-- Initialize arrays and variables used in this subroutine
155
156	!-- If the number of grid points of the gathered grid, which is collected
157	!-- on PE0, is larger than the number of grid points of an PE, than array
158	!-- p3 will be enlarged.
159	IF ( gathered_size > subdomain_size ) THEN
160	ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
161	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
162	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( &
163	mg_switch_to_pe0_level)+1) )
164	ELSE
165	ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
166	ENDIF
167
168	p3 = 0.0
169
170	!
171	!-- Ghost boundaries have to be added to divergence array.
172	!-- Exchange routine needs to know the grid level!
173	grid_level = maximum_grid_level
174	CALL exchange_horiz( d, 1)
175	d(nzb,:,:) = d(nzb+1,:,:)
176
177	!
178	!-- Initiation of the multigrid scheme. Does n cycles until the
179	!-- residual is smaller than the given limit. The accuracy of the solution
180	!-- of the poisson equation will increase with the number of cycles.
181	!-- If the number of cycles is preset by the user, this number will be
182	!-- carried out regardless of the accuracy.
183	grid_level_count = 0
184	mgcycles = 0
185	IF ( mg_cycles == -1 ) THEN
186	maximum_mgcycles = 0
187	residual_norm = 1.0
188	ELSE
189	maximum_mgcycles = mg_cycles
190	residual_norm = 0.0
191	ENDIF
192
193	DO WHILE ( residual_norm > residual_limit .OR. &
194	mgcycles < maximum_mgcycles )
195
196	CALL next_mg_level( d, p_loc, p3, r)
197
198	!
199	!-- Calculate the residual if the user has not preset the number of
200	!-- cycles to be performed
201	IF ( maximum_mgcycles == 0 ) THEN
202	CALL resid( d, p_loc, r )
203	maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
204
205	#if defined( __parallel )
206	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
207	CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, &
208	comm2d, ierr)
209	#else
210	residual_norm = maxerror
211	#endif
212	residual_norm = SQRT( residual_norm )
213	ENDIF
214
215	mgcycles = mgcycles + 1
216
217	!
218	!-- If the user has not limited the number of cycles, stop the run in case
219	!-- of insufficient convergence
220	IF ( mgcycles > 1000 .AND. mg_cycles == -1 ) THEN
221	message_string = 'no sufficient convergence within 1000 cycles'
222	CALL message( 'poismg', 'PA0283', 1, 2, 0, 6, 0 )
223	ENDIF
224
225	ENDDO
226
227	DEALLOCATE( p3 )
228
229	!
230	!-- Unset the grid level. Variable is used to determine the MPI datatypes for
231	!-- ghost point exchange
232	grid_level = 0
233
234	CALL cpu_log( log_point_s(29), 'poismg', 'stop' )
235
236	END SUBROUTINE poismg
237
238
239
240	SUBROUTINE resid( f_mg, p_mg, r )
241
242	!------------------------------------------------------------------------------!
243	! Description:
244	! ------------
245	! Computes the residual of the perturbation pressure.
246	!------------------------------------------------------------------------------!
247
248	USE arrays_3d, &
249	ONLY: f1_mg, f2_mg, f3_mg
250
251	USE control_parameters, &
252	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l, &
253	inflow_n, inflow_r, inflow_s, outflow_l, outflow_n, outflow_r, &
254	outflow_s
255
256	USE grid_variables, &
257	ONLY: ddx2_mg, ddy2_mg
258
259	USE indices, &
260	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4, &
261	wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8, &
262	wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg, &
263	nzb, nzt_mg
264
265	USE kinds
266
267	IMPLICIT NONE
268
269	INTEGER(iwp) :: i
270	INTEGER(iwp) :: j
271	INTEGER(iwp) :: k
272	INTEGER(iwp) :: l
273
274	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
275	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
276	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !:
277	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
278	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
279	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !:
280	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
281	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
282	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !:
283
284	!
285	!-- Calculate the residual
286	l = grid_level
287
288	!
289	!-- Choose flag array of this level
290	SELECT CASE ( l )
291	CASE ( 1 )
292	flags => wall_flags_1
293	CASE ( 2 )
294	flags => wall_flags_2
295	CASE ( 3 )
296	flags => wall_flags_3
297	CASE ( 4 )
298	flags => wall_flags_4
299	CASE ( 5 )
300	flags => wall_flags_5
301	CASE ( 6 )
302	flags => wall_flags_6
303	CASE ( 7 )
304	flags => wall_flags_7
305	CASE ( 8 )
306	flags => wall_flags_8
307	CASE ( 9 )
308	flags => wall_flags_9
309	CASE ( 10 )
310	flags => wall_flags_10
311	END SELECT
312
313	!$OMP PARALLEL PRIVATE (i,j,k)
314	!$OMP DO
315	DO i = nxl_mg(l), nxr_mg(l)
316	DO j = nys_mg(l), nyn_mg(l)
317	DO k = nzb+1, nzt_mg(l)
318	r(k,j,i) = f_mg(k,j,i) &
319	- ddx2_mg(l) * &
320	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
321	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
322	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
323	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
324	- ddy2_mg(l) * &
325	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
326	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
327	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
328	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
329	- f2_mg(k,l) * p_mg(k+1,j,i) &
330	- f3_mg(k,l) * &
331	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
332	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
333	+ f1_mg(k,l) * p_mg(k,j,i)
334	!
335	!-- Residual within topography should be zero
336	r(k,j,i) = r(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) )
337	ENDDO
338	ENDDO
339	ENDDO
340	!$OMP END PARALLEL
341
342	!
343	!-- Horizontal boundary conditions
344	CALL exchange_horiz( r, 1)
345
346	IF ( .NOT. bc_lr_cyc ) THEN
347	IF ( inflow_l .OR. outflow_l ) r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
348	IF ( inflow_r .OR. outflow_r ) r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
349	ENDIF
350
351	IF ( .NOT. bc_ns_cyc ) THEN
352	IF ( inflow_n .OR. outflow_n ) r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
353	IF ( inflow_s .OR. outflow_s ) r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
354	ENDIF
355
356	!
357	!-- Boundary conditions at bottom and top of the domain.
358	!-- These points are not handled by the above loop. Points may be within
359	!-- buildings, but that doesn't matter.
360	IF ( ibc_p_b == 1 ) THEN
361	r(nzb,:,: ) = r(nzb+1,:,:)
362	ELSE
363	r(nzb,:,: ) = 0.0
364	ENDIF
365
366	IF ( ibc_p_t == 1 ) THEN
367	r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
368	ELSE
369	r(nzt_mg(l)+1,:,: ) = 0.0
370	ENDIF
371
372
373	END SUBROUTINE resid
374
375
376
377	SUBROUTINE restrict( f_mg, r )
378
379	!------------------------------------------------------------------------------!
380	! Description:
381	! ------------
382	! Interpolates the residual on the next coarser grid with "full weighting"
383	! scheme
384	!------------------------------------------------------------------------------!
385
386	USE control_parameters, &
387	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l, &
388	inflow_n, inflow_r, inflow_s, outflow_l, outflow_n, outflow_r, &
389	outflow_s
390
391	USE indices, &
392	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4, &
393	wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8, &
394	wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg, &
395	nzb, nzt_mg
396
397	USE kinds
398
399	IMPLICIT NONE
400
401	INTEGER(iwp) :: i !:
402	INTEGER(iwp) :: ic !:
403	INTEGER(iwp) :: j !:
404	INTEGER(iwp) :: jc !:
405	INTEGER(iwp) :: k !:
406	INTEGER(iwp) :: kc !:
407	INTEGER(iwp) :: l !:
408
409	REAL(wp) :: rkjim !:
410	REAL(wp) :: rkjip !:
411	REAL(wp) :: rkjmi !:
412	REAL(wp) :: rkjmim !:
413	REAL(wp) :: rkjmip !:
414	REAL(wp) :: rkjpi !:
415	REAL(wp) :: rkjpim !:
416	REAL(wp) :: rkjpip !:
417	REAL(wp) :: rkmji !:
418	REAL(wp) :: rkmjim !:
419	REAL(wp) :: rkmjip !:
420	REAL(wp) :: rkmjmi !:
421	REAL(wp) :: rkmjmim !:
422	REAL(wp) :: rkmjmip !:
423	REAL(wp) :: rkmjpi !:
424	REAL(wp) :: rkmjpim !:
425	REAL(wp) :: rkmjpip !:
426
427	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
428	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
429	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !:
430
431	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1, &
432	nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, &
433	nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: r !:
434
435	!
436	!-- Interpolate the residual
437	l = grid_level
438
439	!
440	!-- Choose flag array of the upper level
441	SELECT CASE ( l+1 )
442	CASE ( 1 )
443	flags => wall_flags_1
444	CASE ( 2 )
445	flags => wall_flags_2
446	CASE ( 3 )
447	flags => wall_flags_3
448	CASE ( 4 )
449	flags => wall_flags_4
450	CASE ( 5 )
451	flags => wall_flags_5
452	CASE ( 6 )
453	flags => wall_flags_6
454	CASE ( 7 )
455	flags => wall_flags_7
456	CASE ( 8 )
457	flags => wall_flags_8
458	CASE ( 9 )
459	flags => wall_flags_9
460	CASE ( 10 )
461	flags => wall_flags_10
462	END SELECT
463
464	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc)
465	!$OMP DO
466	DO ic = nxl_mg(l), nxr_mg(l)
467	i = 2*ic
468	DO jc = nys_mg(l), nyn_mg(l)
469	j = 2*jc
470	DO kc = nzb+1, nzt_mg(l)
471	k = 2*kc-1
472	!
473	!-- Use implicit Neumann BCs if the respective gridpoint is inside
474	!-- the building
475	rkjim = r(k,j,i-1) + IBITS( flags(k,j,i-1), 6, 1 ) * &
476	( r(k,j,i) - r(k,j,i-1) )
477	rkjip = r(k,j,i+1) + IBITS( flags(k,j,i+1), 6, 1 ) * &
478	( r(k,j,i) - r(k,j,i+1) )
479	rkjpi = r(k,j+1,i) + IBITS( flags(k,j+1,i), 6, 1 ) * &
480	( r(k,j,i) - r(k,j+1,i) )
481	rkjmi = r(k,j-1,i) + IBITS( flags(k,j-1,i), 6, 1 ) * &
482	( r(k,j,i) - r(k,j-1,i) )
483	rkjmim = r(k,j-1,i-1) + IBITS( flags(k,j-1,i-1), 6, 1 ) * &
484	( r(k,j,i) - r(k,j-1,i-1) )
485	rkjpim = r(k,j+1,i-1) + IBITS( flags(k,j+1,i-1), 6, 1 ) * &
486	( r(k,j,i) - r(k,j+1,i-1) )
487	rkjmip = r(k,j-1,i+1) + IBITS( flags(k,j-1,i+1), 6, 1 ) * &
488	( r(k,j,i) - r(k,j-1,i+1) )
489	rkjpip = r(k,j+1,i+1) + IBITS( flags(k,j+1,i+1), 6, 1 ) * &
490	( r(k,j,i) - r(k,j+1,i+1) )
491	rkmji = r(k-1,j,i) + IBITS( flags(k-1,j,i), 6, 1 ) * &
492	( r(k,j,i) - r(k-1,j,i) )
493	rkmjim = r(k-1,j,i-1) + IBITS( flags(k-1,j,i-1), 6, 1 ) * &
494	( r(k,j,i) - r(k-1,j,i-1) )
495	rkmjip = r(k-1,j,i+1) + IBITS( flags(k-1,j,i+1), 6, 1 ) * &
496	( r(k,j,i) - r(k-1,j,i+1) )
497	rkmjpi = r(k-1,j+1,i) + IBITS( flags(k-1,j+1,i), 6, 1 ) * &
498	( r(k,j,i) - r(k-1,j+1,i) )
499	rkmjmi = r(k-1,j-1,i) + IBITS( flags(k-1,j-1,i), 6, 1 ) * &
500	( r(k,j,i) - r(k-1,j-1,i) )
501	rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) * &
502	( r(k,j,i) - r(k-1,j-1,i-1) )
503	rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) * &
504	( r(k,j,i) - r(k-1,j+1,i-1) )
505	rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) * &
506	( r(k,j,i) - r(k-1,j-1,i+1) )
507	rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) * &
508	( r(k,j,i) - r(k-1,j+1,i+1) )
509
510	f_mg(kc,jc,ic) = 1.0 / 64.0_wp * ( &
511	8.0 * r(k,j,i) &
512	+ 4.0 * ( rkjim + rkjip + &
513	rkjpi + rkjmi ) &
514	+ 2.0 * ( rkjmim + rkjpim + &
515	rkjmip + rkjpip ) &
516	+ 4.0 * rkmji &
517	+ 2.0 * ( rkmjim + rkmjim + &
518	rkmjpi + rkmjmi ) &
519	+ ( rkmjmim + rkmjpim + &
520	rkmjmip + rkmjpip ) &
521	+ 4.0 * r(k+1,j,i) &
522	+ 2.0 * ( r(k+1,j,i-1) + r(k+1,j,i+1) + &
523	r(k+1,j+1,i) + r(k+1,j-1,i) ) &
524	+ ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
525	r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
526	)
527
528	! f_mg(kc,jc,ic) = 1.0 / 64.0_wp * ( &
529	! 8.0 * r(k,j,i) &
530	! + 4.0 * ( r(k,j,i-1) + r(k,j,i+1) + &
531	! r(k,j+1,i) + r(k,j-1,i) ) &
532	! + 2.0 * ( r(k,j-1,i-1) + r(k,j+1,i-1) + &
533	! r(k,j-1,i+1) + r(k,j+1,i+1) ) &
534	! + 4.0 * r(k-1,j,i) &
535	! + 2.0 * ( r(k-1,j,i-1) + r(k-1,j,i+1) + &
536	! r(k-1,j+1,i) + r(k-1,j-1,i) ) &
537	! + ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
538	! r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
539	! + 4.0 * r(k+1,j,i) &
540	! + 2.0 * ( r(k+1,j,i-1) + r(k+1,j,i+1) + &
541	! r(k+1,j+1,i) + r(k+1,j-1,i) ) &
542	! + ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
543	! r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
544	! )
545	ENDDO
546	ENDDO
547	ENDDO
548	!$OMP END PARALLEL
549
550	!
551	!-- Horizontal boundary conditions
552	CALL exchange_horiz( f_mg, 1)
553
554	IF ( .NOT. bc_lr_cyc ) THEN
555	IF (inflow_l .OR. outflow_l) f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
556	IF (inflow_r .OR. outflow_r) f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
557	ENDIF
558
559	IF ( .NOT. bc_ns_cyc ) THEN
560	IF (inflow_n .OR. outflow_n) f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
561	IF (inflow_s .OR. outflow_s) f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
562	ENDIF
563
564	!
565	!-- Boundary conditions at bottom and top of the domain.
566	!-- These points are not handled by the above loop. Points may be within
567	!-- buildings, but that doesn't matter.
568	IF ( ibc_p_b == 1 ) THEN
569	f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
570	ELSE
571	f_mg(nzb,:,: ) = 0.0
572	ENDIF
573
574	IF ( ibc_p_t == 1 ) THEN
575	f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
576	ELSE
577	f_mg(nzt_mg(l)+1,:,: ) = 0.0
578	ENDIF
579
580
581	END SUBROUTINE restrict
582
583
584
585	SUBROUTINE prolong( p, temp )
586
587	!------------------------------------------------------------------------------!
588	! Description:
589	! ------------
590	! Interpolates the correction of the perturbation pressure
591	! to the next finer grid.
592	!------------------------------------------------------------------------------!
593
594	USE control_parameters, &
595	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l, &
596	inflow_n, inflow_r, inflow_s, outflow_l, outflow_n, outflow_r, &
597	outflow_s
598
599	USE indices, &
600	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
601
602	USE kinds
603
604	IMPLICIT NONE
605
606	INTEGER(iwp) :: i !:
607	INTEGER(iwp) :: j !:
608	INTEGER(iwp) :: k !:
609	INTEGER(iwp) :: l !:
610
611	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
612	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
613	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: p !:
614
615	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
616	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
617	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: temp !:
618
619
620	!
621	!-- First, store elements of the coarser grid on the next finer grid
622	l = grid_level
623
624	!$OMP PARALLEL PRIVATE (i,j,k)
625	!$OMP DO
626	DO i = nxl_mg(l-1), nxr_mg(l-1)
627	DO j = nys_mg(l-1), nyn_mg(l-1)
628	!CDIR NODEP
629	DO k = nzb+1, nzt_mg(l-1)
630	!
631	!-- Points of the coarse grid are directly stored on the next finer
632	!-- grid
633	temp(2k-1,2j,2*i) = p(k,j,i)
634	!
635	!-- Points between two coarse-grid points
636	temp(2k-1,2j,2i+1) = 0.5 ( p(k,j,i) + p(k,j,i+1) )
637	temp(2k-1,2j+1,2i) = 0.5 ( p(k,j,i) + p(k,j+1,i) )
638	temp(2k,2j,2i) = 0.5 ( p(k,j,i) + p(k+1,j,i) )
639	!
640	!-- Points in the center of the planes stretched by four points
641	!-- of the coarse grid cube
642	temp(2k-1,2j+1,2i+1) = 0.25 ( p(k,j,i) + p(k,j,i+1) + &
643	p(k,j+1,i) + p(k,j+1,i+1) )
644	temp(2k,2j,2i+1) = 0.25 ( p(k,j,i) + p(k,j,i+1) + &
645	p(k+1,j,i) + p(k+1,j,i+1) )
646	temp(2k,2j+1,2i) = 0.25 ( p(k,j,i) + p(k,j+1,i) + &
647	p(k+1,j,i) + p(k+1,j+1,i) )
648	!
649	!-- Points in the middle of coarse grid cube
650	temp(2k,2j+1,2i+1) = 0.125 ( p(k,j,i) + p(k,j,i+1) + &
651	p(k,j+1,i) + p(k,j+1,i+1) + &
652	p(k+1,j,i) + p(k+1,j,i+1) + &
653	p(k+1,j+1,i) + p(k+1,j+1,i+1) )
654	ENDDO
655	ENDDO
656	ENDDO
657	!$OMP END PARALLEL
658
659	!
660	!-- Horizontal boundary conditions
661	CALL exchange_horiz( temp, 1)
662
663	IF ( .NOT. bc_lr_cyc ) THEN
664	IF (inflow_l .OR. outflow_l) temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
665	IF (inflow_r .OR. outflow_r) temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
666	ENDIF
667
668	IF ( .NOT. bc_ns_cyc ) THEN
669	IF (inflow_n .OR. outflow_n) temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
670	IF (inflow_s .OR. outflow_s) temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
671	ENDIF
672
673	!
674	!-- Bottom and top boundary conditions
675	IF ( ibc_p_b == 1 ) THEN
676	temp(nzb,:,: ) = temp(nzb+1,:,:)
677	ELSE
678	temp(nzb,:,: ) = 0.0
679	ENDIF
680
681	IF ( ibc_p_t == 1 ) THEN
682	temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
683	ELSE
684	temp(nzt_mg(l)+1,:,: ) = 0.0
685	ENDIF
686
687
688	END SUBROUTINE prolong
689
690
691	SUBROUTINE redblack( f_mg, p_mg )
692
693	!------------------------------------------------------------------------------!
694	! Description:
695	! ------------
696	! Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
697	! 3D-Red-Black decomposition (GS-RB) is used.
698	!------------------------------------------------------------------------------!
699
700	USE arrays_3d, &
701	ONLY: f1_mg, f2_mg, f3_mg
702
703	USE control_parameters, &
704	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l, &
705	inflow_n, inflow_r, inflow_s, ngsrb, outflow_l, outflow_n, &
706	outflow_r, outflow_s
707
708	USE cpulog, &
709	ONLY: cpu_log, log_point_s
710
711	USE grid_variables, &
712	ONLY: ddx2_mg, ddy2_mg
713
714	USE indices, &
715	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4, &
716	wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8, &
717	wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg, &
718	nzb, nzt_mg
719
720	USE kinds
721
722	IMPLICIT NONE
723
724	INTEGER(iwp) :: color !:
725	INTEGER(iwp) :: i !:
726	INTEGER(iwp) :: ic !:
727	INTEGER(iwp) :: j !:
728	INTEGER(iwp) :: jc !:
729	INTEGER(iwp) :: jj !:
730	INTEGER(iwp) :: k !:
731	INTEGER(iwp) :: l !:
732	INTEGER(iwp) :: n !:
733
734	LOGICAL :: unroll !:
735
736	REAL(wp) :: wall_left !:
737	REAL(wp) :: wall_north !:
738	REAL(wp) :: wall_right !:
739	REAL(wp) :: wall_south !:
740	REAL(wp) :: wall_total !:
741	REAL(wp) :: wall_top !:
742
743	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
744	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
745	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !:
746	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
747	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
748	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !:
749
750	l = grid_level
751
752	!
753	!-- Choose flag array of this level
754	SELECT CASE ( l )
755	CASE ( 1 )
756	flags => wall_flags_1
757	CASE ( 2 )
758	flags => wall_flags_2
759	CASE ( 3 )
760	flags => wall_flags_3
761	CASE ( 4 )
762	flags => wall_flags_4
763	CASE ( 5 )
764	flags => wall_flags_5
765	CASE ( 6 )
766	flags => wall_flags_6
767	CASE ( 7 )
768	flags => wall_flags_7
769	CASE ( 8 )
770	flags => wall_flags_8
771	CASE ( 9 )
772	flags => wall_flags_9
773	CASE ( 10 )
774	flags => wall_flags_10
775	END SELECT
776
777	unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. &
778	MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
779
780	DO n = 1, ngsrb
781
782	DO color = 1, 2
783
784	IF ( .NOT. unroll ) THEN
785
786	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
787
788	!
789	!-- Without unrolling of loops, no cache optimization
790	DO i = nxl_mg(l), nxr_mg(l), 2
791	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
792	DO k = nzb+1, nzt_mg(l), 2
793	! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
794	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
795	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
796	! + f2_mg(k,l) * p_mg(k+1,j,i) &
797	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
798	! )
799
800	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
801	ddx2_mg(l) * &
802	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
803	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
804	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
805	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
806	+ ddy2_mg(l) * &
807	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
808	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
809	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
810	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
811	+ f2_mg(k,l) * p_mg(k+1,j,i) &
812	+ f3_mg(k,l) * &
813	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
814	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
815	- f_mg(k,j,i) )
816	ENDDO
817	ENDDO
818	ENDDO
819
820	DO i = nxl_mg(l)+1, nxr_mg(l), 2
821	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
822	DO k = nzb+1, nzt_mg(l), 2
823	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
824	ddx2_mg(l) * &
825	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
826	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
827	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
828	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
829	+ ddy2_mg(l) * &
830	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
831	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
832	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
833	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
834	+ f2_mg(k,l) * p_mg(k+1,j,i) &
835	+ f3_mg(k,l) * &
836	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
837	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
838	- f_mg(k,j,i) )
839	ENDDO
840	ENDDO
841	ENDDO
842
843	DO i = nxl_mg(l), nxr_mg(l), 2
844	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
845	DO k = nzb+2, nzt_mg(l), 2
846	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
847	ddx2_mg(l) * &
848	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
849	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
850	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
851	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
852	+ ddy2_mg(l) * &
853	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
854	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
855	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
856	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
857	+ f2_mg(k,l) * p_mg(k+1,j,i) &
858	+ f3_mg(k,l) * &
859	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
860	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
861	- f_mg(k,j,i) )
862	ENDDO
863	ENDDO
864	ENDDO
865
866	DO i = nxl_mg(l)+1, nxr_mg(l), 2
867	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
868	DO k = nzb+2, nzt_mg(l), 2
869	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
870	ddx2_mg(l) * &
871	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
872	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
873	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
874	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
875	+ ddy2_mg(l) * &
876	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
877	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
878	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
879	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
880	+ f2_mg(k,l) * p_mg(k+1,j,i) &
881	+ f3_mg(k,l) * &
882	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
883	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
884	- f_mg(k,j,i) )
885	ENDDO
886	ENDDO
887	ENDDO
888	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
889
890	ELSE
891
892	!
893	!-- Loop unrolling along y, only one i loop for better cache use
894	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
895	DO ic = nxl_mg(l), nxr_mg(l), 2
896	DO jc = nys_mg(l), nyn_mg(l), 4
897	i = ic
898	jj = jc+2-color
899	DO k = nzb+1, nzt_mg(l), 2
900	j = jj
901	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
902	ddx2_mg(l) * &
903	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
904	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
905	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
906	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
907	+ ddy2_mg(l) * &
908	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
909	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
910	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
911	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
912	+ f2_mg(k,l) * p_mg(k+1,j,i) &
913	+ f3_mg(k,l) * &
914	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
915	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
916	- f_mg(k,j,i) )
917	j = jj+2
918	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
919	ddx2_mg(l) * &
920	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
921	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
922	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
923	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
924	+ ddy2_mg(l) * &
925	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
926	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
927	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
928	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
929	+ f2_mg(k,l) * p_mg(k+1,j,i) &
930	+ f3_mg(k,l) * &
931	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
932	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
933	- f_mg(k,j,i) )
934	ENDDO
935
936	i = ic+1
937	jj = jc+color-1
938	DO k = nzb+1, nzt_mg(l), 2
939	j =jj
940	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
941	ddx2_mg(l) * &
942	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
943	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
944	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
945	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
946	+ ddy2_mg(l) * &
947	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
948	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
949	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
950	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
951	+ f2_mg(k,l) * p_mg(k+1,j,i) &
952	+ f3_mg(k,l) * &
953	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
954	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
955	- f_mg(k,j,i) )
956	j = jj+2
957	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
958	ddx2_mg(l) * &
959	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
960	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
961	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
962	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
963	+ ddy2_mg(l) * &
964	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
965	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
966	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
967	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
968	+ f2_mg(k,l) * p_mg(k+1,j,i) &
969	+ f3_mg(k,l) * &
970	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
971	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
972	- f_mg(k,j,i) )
973	ENDDO
974
975	i = ic
976	jj = jc+color-1
977	DO k = nzb+2, nzt_mg(l), 2
978	j =jj
979	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
980	ddx2_mg(l) * &
981	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
982	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
983	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
984	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
985	+ ddy2_mg(l) * &
986	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
987	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
988	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
989	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
990	+ f2_mg(k,l) * p_mg(k+1,j,i) &
991	+ f3_mg(k,l) * &
992	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
993	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
994	- f_mg(k,j,i) )
995	j = jj+2
996	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
997	ddx2_mg(l) * &
998	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
999	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1000	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1001	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1002	+ ddy2_mg(l) * &
1003	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1004	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1005	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1006	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1007	+ f2_mg(k,l) * p_mg(k+1,j,i) &
1008	+ f3_mg(k,l) * &
1009	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1010	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1011	- f_mg(k,j,i) )
1012	ENDDO
1013
1014	i = ic+1
1015	jj = jc+2-color
1016	DO k = nzb+2, nzt_mg(l), 2
1017	j =jj
1018	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1019	ddx2_mg(l) * &
1020	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1021	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1022	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1023	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1024	+ ddy2_mg(l) * &
1025	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1026	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1027	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1028	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1029	+ f2_mg(k,l) * p_mg(k+1,j,i) &
1030	+ f3_mg(k,l) * &
1031	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1032	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1033	- f_mg(k,j,i) )
1034	j = jj+2
1035	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1036	ddx2_mg(l) * &
1037	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1038	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1039	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1040	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1041	+ ddy2_mg(l) * &
1042	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1043	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1044	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1045	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1046	+ f2_mg(k,l) * p_mg(k+1,j,i) &
1047	+ f3_mg(k,l) * &
1048	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1049	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1050	- f_mg(k,j,i) )
1051	ENDDO
1052
1053	ENDDO
1054	ENDDO
1055	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
1056
1057	ENDIF
1058
1059	!
1060	!-- Horizontal boundary conditions
1061	CALL exchange_horiz( p_mg, 1 )
1062
1063	IF ( .NOT. bc_lr_cyc ) THEN
1064	IF ( inflow_l .OR. outflow_l ) THEN
1065	p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
1066	ENDIF
1067	IF ( inflow_r .OR. outflow_r ) THEN
1068	p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
1069	ENDIF
1070	ENDIF
1071
1072	IF ( .NOT. bc_ns_cyc ) THEN
1073	IF ( inflow_n .OR. outflow_n ) THEN
1074	p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
1075	ENDIF
1076	IF ( inflow_s .OR. outflow_s ) THEN
1077	p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
1078	ENDIF
1079	ENDIF
1080
1081	!
1082	!-- Bottom and top boundary conditions
1083	IF ( ibc_p_b == 1 ) THEN
1084	p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
1085	ELSE
1086	p_mg(nzb,:,: ) = 0.0
1087	ENDIF
1088
1089	IF ( ibc_p_t == 1 ) THEN
1090	p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
1091	ELSE
1092	p_mg(nzt_mg(l)+1,:,: ) = 0.0
1093	ENDIF
1094
1095	ENDDO
1096
1097	ENDDO
1098
1099	!
1100	!-- Set pressure within topography and at the topography surfaces
1101	!$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total)
1102	!$OMP DO
1103	DO i = nxl_mg(l), nxr_mg(l)
1104	DO j = nys_mg(l), nyn_mg(l)
1105	DO k = nzb, nzt_mg(l)
1106	!
1107	!-- First, set pressure inside topography to zero
1108	p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) )
1109	!
1110	!-- Second, determine if the gridpoint inside topography is adjacent
1111	!-- to a wall and set its value to a value given by the average of
1112	!-- those values obtained from Neumann boundary condition
1113	wall_left = IBITS( flags(k,j,i-1), 5, 1 )
1114	wall_right = IBITS( flags(k,j,i+1), 4, 1 )
1115	wall_south = IBITS( flags(k,j-1,i), 3, 1 )
1116	wall_north = IBITS( flags(k,j+1,i), 2, 1 )
1117	wall_top = IBITS( flags(k+1,j,i), 0, 1 )
1118	wall_total = wall_left + wall_right + wall_south + wall_north + &
1119	wall_top
1120
1121	IF ( wall_total > 0.0 ) THEN
1122	p_mg(k,j,i) = 1.0 / wall_total * &
1123	( wall_left * p_mg(k,j,i-1) + &
1124	wall_right * p_mg(k,j,i+1) + &
1125	wall_south * p_mg(k,j-1,i) + &
1126	wall_north * p_mg(k,j+1,i) + &
1127	wall_top * p_mg(k+1,j,i) )
1128	ENDIF
1129	ENDDO
1130	ENDDO
1131	ENDDO
1132	!$OMP END PARALLEL
1133
1134	!
1135	!-- One more time horizontal boundary conditions
1136	CALL exchange_horiz( p_mg, 1)
1137
1138
1139	END SUBROUTINE redblack
1140
1141
1142
1143	SUBROUTINE mg_gather( f2, f2_sub )
1144
1145	USE control_parameters, &
1146	ONLY: grid_level
1147
1148	USE cpulog, &
1149	ONLY: cpu_log, log_point_s
1150
1151	USE indices, &
1152	ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1153
1154	USE kinds
1155
1156	USE pegrid
1157
1158	IMPLICIT NONE
1159
1160	INTEGER(iwp) :: i !:
1161	INTEGER(iwp) :: il !:
1162	INTEGER(iwp) :: ir !:
1163	INTEGER(iwp) :: j !:
1164	INTEGER(iwp) :: jn !:
1165	INTEGER(iwp) :: js !:
1166	INTEGER(iwp) :: k !:
1167	INTEGER(iwp) :: nwords !:
1168
1169	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1170	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1171	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2 !:
1172	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1173	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1174	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2_l !:
1175
1176	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
1177	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1178	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub !:
1179
1180
1181	#if defined( __parallel )
1182	CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
1183
1184	f2_l = 0.0
1185
1186	!
1187	!-- Store the local subdomain array on the total array
1188	js = mg_loc_ind(3,myid)
1189	IF ( south_border_pe ) js = js - 1
1190	jn = mg_loc_ind(4,myid)
1191	IF ( north_border_pe ) jn = jn + 1
1192	il = mg_loc_ind(1,myid)
1193	IF ( left_border_pe ) il = il - 1
1194	ir = mg_loc_ind(2,myid)
1195	IF ( right_border_pe ) ir = ir + 1
1196	DO i = il, ir
1197	DO j = js, jn
1198	DO k = nzb, nzt_mg(grid_level)+1
1199	f2_l(k,j,i) = f2_sub(k,j,i)
1200	ENDDO
1201	ENDDO
1202	ENDDO
1203
1204	!
1205	!-- Find out the number of array elements of the total array
1206	nwords = SIZE( f2 )
1207
1208	!
1209	!-- Gather subdomain data from all PEs
1210	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1211	CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
1212	f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
1213	nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
1214
1215	CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
1216	#endif
1217
1218	END SUBROUTINE mg_gather
1219
1220
1221
1222	SUBROUTINE mg_scatter( p2, p2_sub )
1223	!
1224	!-- TODO: It may be possible to improve the speed of this routine by using
1225	!-- non-blocking communication
1226
1227	USE control_parameters, &
1228	ONLY: grid_level
1229
1230	USE cpulog, &
1231	ONLY: cpu_log, log_point_s
1232
1233	USE indices, &
1234	ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1235
1236	USE kinds
1237
1238	USE pegrid
1239
1240	IMPLICIT NONE
1241
1242	INTEGER(iwp) :: nwords !:
1243
1244	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
1245	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
1246	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !:
1247
1248	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
1249	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1250	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub !:
1251
1252	!
1253	!-- Find out the number of array elements of the subdomain array
1254	nwords = SIZE( p2_sub )
1255
1256	#if defined( __parallel )
1257	CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
1258
1259	p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1260	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
1261
1262	CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
1263	#endif
1264
1265	END SUBROUTINE mg_scatter
1266
1267
1268
1269	RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r )
1270
1271	!------------------------------------------------------------------------------!
1272	! Description:
1273	! ------------
1274	! This is where the multigrid technique takes place. V- and W- Cycle are
1275	! implemented and steered by the parameter "gamma". Parameter "nue" determines
1276	! the convergence of the multigrid iterative solution. There are nue times
1277	! RB-GS iterations. It should be set to "1" or "2", considering the time effort
1278	! one would like to invest. Last choice shows a very good converging factor,
1279	! but leads to an increase in computing time.
1280	!------------------------------------------------------------------------------!
1281
1282	USE control_parameters, &
1283	ONLY: bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir, &
1284	gamma_mg, grid_level, grid_level_count, ibc_p_b, ibc_p_t, &
1285	inflow_l, inflow_n, inflow_r, inflow_s, maximum_grid_level, &
1286	mg_switch_to_pe0_level, mg_switch_to_pe0, ngsrb, outflow_l, &
1287	outflow_n, outflow_r, outflow_s
1288
1289
1290	USE indices, &
1291	ONLY: mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn, &
1292	nyn_mg, nzb, nzt, nzt_mg
1293
1294	USE kinds
1295
1296	USE pegrid
1297
1298	IMPLICIT NONE
1299
1300	INTEGER(iwp) :: i !:
1301	INTEGER(iwp) :: j !:
1302	INTEGER(iwp) :: k !:
1303	INTEGER(iwp) :: nxl_mg_save !:
1304	INTEGER(iwp) :: nxr_mg_save !:
1305	INTEGER(iwp) :: nyn_mg_save !:
1306	INTEGER(iwp) :: nys_mg_save !:
1307	INTEGER(iwp) :: nzt_mg_save !:
1308
1309	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1310	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1311	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !:
1312	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1313	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1314	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !:
1315	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1316	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1317	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3 !:
1318	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1319	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1320	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !:
1321
1322	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
1323	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
1324	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: f2 !:
1325	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
1326	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
1327	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !:
1328
1329	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f2_sub !:
1330	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p2_sub !:
1331
1332	!
1333	!-- Restriction to the coarsest grid
1334	10 IF ( grid_level == 1 ) THEN
1335
1336	!
1337	!-- Solution on the coarsest grid. Double the number of Gauss-Seidel
1338	!-- iterations in order to get a more accurate solution.
1339	ngsrb = 2 * ngsrb
1340
1341	CALL redblack( f_mg, p_mg )
1342
1343	ngsrb = ngsrb / 2
1344
1345
1346	ELSEIF ( grid_level /= 1 ) THEN
1347
1348	grid_level_count(grid_level) = grid_level_count(grid_level) + 1
1349
1350	!
1351	!-- Solution on the actual grid level
1352	CALL redblack( f_mg, p_mg )
1353
1354	!
1355	!-- Determination of the actual residual
1356	CALL resid( f_mg, p_mg, r )
1357
1358	!
1359	!-- Restriction of the residual (finer grid values!) to the next coarser
1360	!-- grid. Therefore, the grid level has to be decremented now. nxl..nzt have
1361	!-- to be set to the coarse grid values, because these variables are needed
1362	!-- for the exchange of ghost points in routine exchange_horiz
1363	grid_level = grid_level - 1
1364	nxl = nxl_mg(grid_level)
1365	nys = nys_mg(grid_level)
1366	nxr = nxr_mg(grid_level)
1367	nyn = nyn_mg(grid_level)
1368	nzt = nzt_mg(grid_level)
1369
1370	IF ( grid_level == mg_switch_to_pe0_level ) THEN
1371
1372	!
1373	!-- From this level on, calculations are done on PE0 only.
1374	!-- First, carry out restriction on the subdomain.
1375	!-- Therefore, indices of the level have to be changed to subdomain values
1376	!-- in between (otherwise, the restrict routine would expect
1377	!-- the gathered array)
1378
1379	nxl_mg_save = nxl_mg(grid_level)
1380	nxr_mg_save = nxr_mg(grid_level)
1381	nys_mg_save = nys_mg(grid_level)
1382	nyn_mg_save = nyn_mg(grid_level)
1383	nzt_mg_save = nzt_mg(grid_level)
1384	nxl_mg(grid_level) = mg_loc_ind(1,myid)
1385	nxr_mg(grid_level) = mg_loc_ind(2,myid)
1386	nys_mg(grid_level) = mg_loc_ind(3,myid)
1387	nyn_mg(grid_level) = mg_loc_ind(4,myid)
1388	nzt_mg(grid_level) = mg_loc_ind(5,myid)
1389	nxl = mg_loc_ind(1,myid)
1390	nxr = mg_loc_ind(2,myid)
1391	nys = mg_loc_ind(3,myid)
1392	nyn = mg_loc_ind(4,myid)
1393	nzt = mg_loc_ind(5,myid)
1394
1395	ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1, &
1396	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1397	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
1398
1399	CALL restrict( f2_sub, r )
1400
1401	!
1402	!-- Restore the correct indices of this level
1403	nxl_mg(grid_level) = nxl_mg_save
1404	nxr_mg(grid_level) = nxr_mg_save
1405	nys_mg(grid_level) = nys_mg_save
1406	nyn_mg(grid_level) = nyn_mg_save
1407	nzt_mg(grid_level) = nzt_mg_save
1408	nxl = nxl_mg(grid_level)
1409	nxr = nxr_mg(grid_level)
1410	nys = nys_mg(grid_level)
1411	nyn = nyn_mg(grid_level)
1412	nzt = nzt_mg(grid_level)
1413	!
1414	!-- Gather all arrays from the subdomains on PE0
1415	CALL mg_gather( f2, f2_sub )
1416
1417	!
1418	!-- Set switch for routine exchange_horiz, that no ghostpoint exchange
1419	!-- has to be carried out from now on
1420	mg_switch_to_pe0 = .TRUE.
1421
1422	!
1423	!-- In case of non-cyclic lateral boundary conditions, both in- and
1424	!-- outflow conditions have to be used on all PEs after the switch,
1425	!-- because then they have the total domain.
1426	IF ( bc_lr_dirrad ) THEN
1427	inflow_l = .TRUE.
1428	inflow_r = .FALSE.
1429	outflow_l = .FALSE.
1430	outflow_r = .TRUE.
1431	ELSEIF ( bc_lr_raddir ) THEN
1432	inflow_l = .FALSE.
1433	inflow_r = .TRUE.
1434	outflow_l = .TRUE.
1435	outflow_r = .FALSE.
1436	ENDIF
1437
1438	IF ( bc_ns_dirrad ) THEN
1439	inflow_n = .TRUE.
1440	inflow_s = .FALSE.
1441	outflow_n = .FALSE.
1442	outflow_s = .TRUE.
1443	ELSEIF ( bc_ns_raddir ) THEN
1444	inflow_n = .FALSE.
1445	inflow_s = .TRUE.
1446	outflow_n = .TRUE.
1447	outflow_s = .FALSE.
1448	ENDIF
1449
1450	DEALLOCATE( f2_sub )
1451
1452	ELSE
1453
1454	CALL restrict( f2, r )
1455
1456	ENDIF
1457
1458	p2 = 0.0
1459
1460	!
1461	!-- Repeat the same procedure till the coarsest grid is reached
1462	CALL next_mg_level( f2, p2, p3, r )
1463
1464	ENDIF
1465
1466	!
1467	!-- Now follows the prolongation
1468	IF ( grid_level >= 2 ) THEN
1469
1470	!
1471	!-- Prolongation of the new residual. The values are transferred
1472	!-- from the coarse to the next finer grid.
1473	IF ( grid_level == mg_switch_to_pe0_level+1 ) THEN
1474
1475	#if defined( __parallel )
1476	!
1477	!-- At this level, the new residual first has to be scattered from
1478	!-- PE0 to the other PEs
1479	ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1, &
1480	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1481	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
1482
1483	CALL mg_scatter( p2, p2_sub )
1484
1485	!
1486	!-- Therefore, indices of the previous level have to be changed to
1487	!-- subdomain values in between (otherwise, the prolong routine would
1488	!-- expect the gathered array)
1489	nxl_mg_save = nxl_mg(grid_level-1)
1490	nxr_mg_save = nxr_mg(grid_level-1)
1491	nys_mg_save = nys_mg(grid_level-1)
1492	nyn_mg_save = nyn_mg(grid_level-1)
1493	nzt_mg_save = nzt_mg(grid_level-1)
1494	nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
1495	nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
1496	nys_mg(grid_level-1) = mg_loc_ind(3,myid)
1497	nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
1498	nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
1499
1500	!
1501	!-- Set switch for routine exchange_horiz, that ghostpoint exchange
1502	!-- has to be carried again out from now on
1503	mg_switch_to_pe0 = .FALSE.
1504
1505	!
1506	!-- For non-cyclic lateral boundary conditions, restore the
1507	!-- in-/outflow conditions
1508	inflow_l = .FALSE.; inflow_r = .FALSE.
1509	inflow_n = .FALSE.; inflow_s = .FALSE.
1510	outflow_l = .FALSE.; outflow_r = .FALSE.
1511	outflow_n = .FALSE.; outflow_s = .FALSE.
1512
1513	IF ( pleft == MPI_PROC_NULL ) THEN
1514	IF ( bc_lr_dirrad ) THEN
1515	inflow_l = .TRUE.
1516	ELSEIF ( bc_lr_raddir ) THEN
1517	outflow_l = .TRUE.
1518	ENDIF
1519	ENDIF
1520
1521	IF ( pright == MPI_PROC_NULL ) THEN
1522	IF ( bc_lr_dirrad ) THEN
1523	outflow_r = .TRUE.
1524	ELSEIF ( bc_lr_raddir ) THEN
1525	inflow_r = .TRUE.
1526	ENDIF
1527	ENDIF
1528
1529	IF ( psouth == MPI_PROC_NULL ) THEN
1530	IF ( bc_ns_dirrad ) THEN
1531	outflow_s = .TRUE.
1532	ELSEIF ( bc_ns_raddir ) THEN
1533	inflow_s = .TRUE.
1534	ENDIF
1535	ENDIF
1536
1537	IF ( pnorth == MPI_PROC_NULL ) THEN
1538	IF ( bc_ns_dirrad ) THEN
1539	inflow_n = .TRUE.
1540	ELSEIF ( bc_ns_raddir ) THEN
1541	outflow_n = .TRUE.
1542	ENDIF
1543	ENDIF
1544
1545	CALL prolong( p2_sub, p3 )
1546
1547	!
1548	!-- Restore the correct indices of the previous level
1549	nxl_mg(grid_level-1) = nxl_mg_save
1550	nxr_mg(grid_level-1) = nxr_mg_save
1551	nys_mg(grid_level-1) = nys_mg_save
1552	nyn_mg(grid_level-1) = nyn_mg_save
1553	nzt_mg(grid_level-1) = nzt_mg_save
1554
1555	DEALLOCATE( p2_sub )
1556	#endif
1557
1558	ELSE
1559
1560	CALL prolong( p2, p3 )
1561
1562	ENDIF
1563
1564	!
1565	!-- Computation of the new pressure correction. Therefore,
1566	!-- values from prior grids are added up automatically stage by stage.
1567	DO i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
1568	DO j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
1569	DO k = nzb, nzt_mg(grid_level)+1
1570	p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
1571	ENDDO
1572	ENDDO
1573	ENDDO
1574
1575	!
1576	!-- Relaxation of the new solution
1577	CALL redblack( f_mg, p_mg )
1578
1579	ENDIF
1580
1581
1582	!
1583	!-- The following few lines serve the steering of the multigrid scheme
1584	IF ( grid_level == maximum_grid_level ) THEN
1585
1586	GOTO 20
1587
1588	ELSEIF ( grid_level /= maximum_grid_level .AND. grid_level /= 1 .AND. &
1589	grid_level_count(grid_level) /= gamma_mg ) THEN
1590
1591	GOTO 10
1592
1593	ENDIF
1594
1595	!
1596	!-- Reset counter for the next call of poismg
1597	grid_level_count(grid_level) = 0
1598
1599	!
1600	!-- Continue with the next finer level. nxl..nzt have to be
1601	!-- set to the finer grid values, because these variables are needed for the
1602	!-- exchange of ghost points in routine exchange_horiz
1603	grid_level = grid_level + 1
1604	nxl = nxl_mg(grid_level)
1605	nxr = nxr_mg(grid_level)
1606	nys = nys_mg(grid_level)
1607	nyn = nyn_mg(grid_level)
1608	nzt = nzt_mg(grid_level)
1609
1610	20 CONTINUE
1611
1612	END SUBROUTINE next_mg_level

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |