Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |