Home

Context Navigation

poismg_noopt_mod.f90 @ 4875

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

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

Context Navigation

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

Download in other formats: