Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |