Home

Context Navigation

source: palm/trunk/SOURCE/poismg_noopt_mod.f90 @ 3016

Last change on this file since 3016 was 2939, checked in by suehring, 7 years ago
Set lateral boundary conditions for divergence in multigrid solver
Property svn:keywords set to `Id`
File size: 82.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |