Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |