Home

Context Navigation

source: palm/trunk/SOURCE/poismg_noopt.f90 @ 2605

Last change on this file since 2605 was 2233, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 68.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |