Home

Context Navigation

poismg.f90 @ 600

Last change on this file since 600 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: