Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |