Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |