Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |