Home

Context Navigation

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

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