Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |