Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |