Home

Context Navigation

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

Last change on this file since 624 was 623, checked in by raasch, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 54.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |