Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |