Home

Context Navigation

source: palm/trunk/SOURCE/poismg.f90 @ 1762

Last change on this file since 1762 was 1762, checked in by hellstea, 8 years ago
Introduction of nested domain system
Property svn:keywords set to `Id`
File size: 66.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |