Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |