Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |