Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |