Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |