Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |