Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |