Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |