Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |