Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |