Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |