Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |