Home

Context Navigation

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

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