Home

Context Navigation

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

Last change on this file since 1238 was 1160, checked in by fricke, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 56.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |