Home

Context Navigation

source: palm/trunk/SOURCE/poismg_noopt.f90 @ 2016

Last change on this file since 2016 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 67.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |