Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |