Home

Context Navigation

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

Last change on this file since 4558 was 4457, checked in by raasch, 5 years ago
ghost point exchange modularized, bugfix for wrong 2d-exchange
Property svn:keywords set to `Id`
File size: 77.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |