Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |