Home

Context Navigation

poismg.f90 @ 141

Last change on this file since 141 was 139, checked in by raasch, 17 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:

Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:

Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: