Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |