Home

Context Navigation

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

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