Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |