Home

Context Navigation

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

Last change on this file since 986 was 979, checked in by fricke, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 55.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |