Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |