Home

Context Navigation

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

Last change on this file since 1745 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 65.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |