Home

Context Navigation

poismg.f90 @ 622

Last change on this file since 622 was 622, checked in by raasch, 13 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:

Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: