Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |