Home

Context Navigation

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

Last change on this file since 2809 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 81.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |