Home

Context Navigation

init_grid.f90 @ 449

Last change on this file since 449 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: 39.6 KB

Line
1	SUBROUTINE init_grid
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: init_grid.f90 392 2009-09-24 10:39:14Z raasch $
11	!
12	! 274 2009-03-26 15:11:21Z heinze
13	! Output of messages replaced by message handling routine.
14	! new topography case 'single_street_canyon'
15	!
16	! 217 2008-12-09 18:00:48Z letzel
17	! +topography_grid_convention
18	!
19	! 134 2007-11-21 07:28:38Z letzel
20	! Redefine initial nzb_local as the actual total size of topography (later the
21	! extent of topography in nzb_local is reduced by 1dx at the E topography walls
22	! and by 1dy at the N topography walls to form the basis for nzb_s_inner);
23	! for consistency redefine 'single_building' case.
24	! Calculation of wall flag arrays
25	!
26	! 94 2007-06-01 15:25:22Z raasch
27	! Grid definition for ocean version
28	!
29	! 75 2007-03-22 09:54:05Z raasch
30	! storage of topography height arrays zu_s_inner and zw_s_inner,
31	! 2nd+3rd argument removed from exchange horiz
32	!
33	! 19 2007-02-23 04:53:48Z raasch
34	! Setting of nzt_diff
35	!
36	! RCS Log replace by Id keyword, revision history cleaned up
37	!
38	! Revision 1.17 2006/08/22 14:00:05 raasch
39	! +dz_max to limit vertical stretching,
40	! bugfix in index array initialization for line- or point-like topography
41	! structures
42	!
43	! Revision 1.1 1997/08/11 06:17:45 raasch
44	! Initial revision (Testversion)
45	!
46	!
47	! Description:
48	! ------------
49	! Creating grid depending constants
50	!------------------------------------------------------------------------------!
51
52	USE arrays_3d
53	USE control_parameters
54	USE grid_variables
55	USE indices
56	USE pegrid
57
58	IMPLICIT NONE
59
60	INTEGER :: bh, blx, bly, bxl, bxr, byn, bys, ch, cwx, cwy, cxl, cxr, cyn, &
61	cys, gls, i, inc, i_center, j, j_center, k, l, nxl_l, nxr_l, &
62	nyn_l, nys_l, nzb_si, nzt_l, vi
63
64	INTEGER, DIMENSION(:), ALLOCATABLE :: vertical_influence
65
66	INTEGER, DIMENSION(:,:), ALLOCATABLE :: corner_nl, corner_nr, corner_sl, &
67	corner_sr, wall_l, wall_n, wall_r,&
68	wall_s, nzb_local, nzb_tmp
69
70	REAL :: dx_l, dy_l, dz_stretched
71
72	REAL, DIMENSION(0:ny,0:nx) :: topo_height
73
74	REAL, DIMENSION(:,:,:), ALLOCATABLE :: distance
75
76	!
77	!-- Allocate grid arrays
78	ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &
79	dzw(1:nzt+1), l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1) )
80
81	!
82	!-- Compute height of u-levels from constant grid length and dz stretch factors
83	IF ( dz == -1.0 ) THEN
84	message_string = 'missing dz'
85	CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
86	ELSEIF ( dz <= 0.0 ) THEN
87	WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
88	CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
89	ENDIF
90
91	!
92	!-- Define the vertical grid levels
93	IF ( .NOT. ocean ) THEN
94	!
95	!-- Grid for atmosphere with surface at z=0 (k=0, w-grid).
96	!-- Since the w-level lies on the surface, the first u-level (staggered!)
97	!-- lies below the surface (used for "mirror" boundary condition).
98	!-- The first u-level above the surface corresponds to the top of the
99	!-- Prandtl-layer.
100	zu(0) = - dz * 0.5
101	zu(1) = dz * 0.5
102
103	dz_stretch_level_index = nzt+1
104	dz_stretched = dz
105	DO k = 2, nzt+1
106	IF ( dz_stretch_level <= zu(k-1) .AND. dz_stretched < dz_max ) THEN
107	dz_stretched = dz_stretched * dz_stretch_factor
108	dz_stretched = MIN( dz_stretched, dz_max )
109	IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
110	ENDIF
111	zu(k) = zu(k-1) + dz_stretched
112	ENDDO
113
114	!
115	!-- Compute the w-levels. They are always staggered half-way between the
116	!-- corresponding u-levels. The top w-level is extrapolated linearly.
117	zw(0) = 0.0
118	DO k = 1, nzt
119	zw(k) = ( zu(k) + zu(k+1) ) * 0.5
120	ENDDO
121	zw(nzt+1) = zw(nzt) + 2.0 * ( zu(nzt+1) - zw(nzt) )
122
123	ELSE
124	!
125	!-- Grid for ocean with solid surface at z=0 (k=0, w-grid). The free water
126	!-- surface is at k=nzt (w-grid).
127	!-- Since the w-level lies always on the surface, the first/last u-level
128	!-- (staggered!) lies below the bottom surface / above the free surface.
129	!-- It is used for "mirror" boundary condition.
130	!-- The first u-level above the bottom surface corresponds to the top of the
131	!-- Prandtl-layer.
132	zu(nzt+1) = dz * 0.5
133	zu(nzt) = - dz * 0.5
134
135	dz_stretch_level_index = 0
136	dz_stretched = dz
137	DO k = nzt-1, 0, -1
138	IF ( dz_stretch_level <= ABS( zu(k+1) ) .AND. &
139	dz_stretched < dz_max ) THEN
140	dz_stretched = dz_stretched * dz_stretch_factor
141	dz_stretched = MIN( dz_stretched, dz_max )
142	IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
143	ENDIF
144	zu(k) = zu(k+1) - dz_stretched
145	ENDDO
146
147	!
148	!-- Compute the w-levels. They are always staggered half-way between the
149	!-- corresponding u-levels.
150	!-- The top w-level (nzt+1) is not used but set for consistency, since
151	!-- w and all scalar variables are defined up tp nzt+1.
152	zw(nzt+1) = dz
153	zw(nzt) = 0.0
154	DO k = 0, nzt
155	zw(k) = ( zu(k) + zu(k+1) ) * 0.5
156	ENDDO
157
158	ENDIF
159
160	!
161	!-- Compute grid lengths.
162	DO k = 1, nzt+1
163	dzu(k) = zu(k) - zu(k-1)
164	ddzu(k) = 1.0 / dzu(k)
165	dzw(k) = zw(k) - zw(k-1)
166	ddzw(k) = 1.0 / dzw(k)
167	ENDDO
168
169	DO k = 1, nzt
170	dd2zu(k) = 1.0 / ( dzu(k) + dzu(k+1) )
171	ENDDO
172
173	!
174	!-- In case of multigrid method, compute grid lengths and grid factors for the
175	!-- grid levels
176	IF ( psolver == 'multigrid' ) THEN
177
178	ALLOCATE( ddx2_mg(maximum_grid_level), ddy2_mg(maximum_grid_level), &
179	dzu_mg(nzb+1:nzt+1,maximum_grid_level), &
180	dzw_mg(nzb+1:nzt+1,maximum_grid_level), &
181	f1_mg(nzb+1:nzt,maximum_grid_level), &
182	f2_mg(nzb+1:nzt,maximum_grid_level), &
183	f3_mg(nzb+1:nzt,maximum_grid_level) )
184
185	dzu_mg(:,maximum_grid_level) = dzu
186	dzw_mg(:,maximum_grid_level) = dzw
187	nzt_l = nzt
188	DO l = maximum_grid_level-1, 1, -1
189	dzu_mg(nzb+1,l) = 2.0 * dzu_mg(nzb+1,l+1)
190	dzw_mg(nzb+1,l) = 2.0 * dzw_mg(nzb+1,l+1)
191	nzt_l = nzt_l / 2
192	DO k = 2, nzt_l+1
193	dzu_mg(k,l) = dzu_mg(2k-2,l+1) + dzu_mg(2k-1,l+1)
194	dzw_mg(k,l) = dzw_mg(2k-2,l+1) + dzw_mg(2k-1,l+1)
195	ENDDO
196	ENDDO
197
198	nzt_l = nzt
199	dx_l = dx
200	dy_l = dy
201	DO l = maximum_grid_level, 1, -1
202	ddx2_mg(l) = 1.0 / dx_l**2
203	ddy2_mg(l) = 1.0 / dy_l**2
204	DO k = nzb+1, nzt_l
205	f2_mg(k,l) = 1.0 / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
206	f3_mg(k,l) = 1.0 / ( dzu_mg(k,l) * dzw_mg(k,l) )
207	f1_mg(k,l) = 2.0 * ( ddx2_mg(l) + ddy2_mg(l) ) + &
208	f2_mg(k,l) + f3_mg(k,l)
209	ENDDO
210	nzt_l = nzt_l / 2
211	dx_l = dx_l * 2.0
212	dy_l = dy_l * 2.0
213	ENDDO
214
215	ENDIF
216
217	!
218	!-- Compute the reciprocal values of the horizontal grid lengths.
219	ddx = 1.0 / dx
220	ddy = 1.0 / dy
221	dx2 = dx * dx
222	dy2 = dy * dy
223	ddx2 = 1.0 / dx2
224	ddy2 = 1.0 / dy2
225
226	!
227	!-- Compute the grid-dependent mixing length.
228	DO k = 1, nzt
229	l_grid(k) = ( dx * dy * dzw(k) )**0.33333333333333
230	ENDDO
231
232	!
233	!-- Allocate outer and inner index arrays for topography and set
234	!-- defaults.
235	!-- nzb_local has to contain additional layers of ghost points for calculating
236	!-- the flag arrays needed for the multigrid method
237	gls = 2**( maximum_grid_level )
238	ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), &
239	corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), &
240	nzb_local(-gls:ny+gls,-gls:nx+gls), nzb_tmp(-1:ny+1,-1:nx+1), &
241	wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr), &
242	wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
243	ALLOCATE( fwxm(nys-1:nyn+1,nxl-1:nxr+1), fwxp(nys-1:nyn+1,nxl-1:nxr+1), &
244	fwym(nys-1:nyn+1,nxl-1:nxr+1), fwyp(nys-1:nyn+1,nxl-1:nxr+1), &
245	fxm(nys-1:nyn+1,nxl-1:nxr+1), fxp(nys-1:nyn+1,nxl-1:nxr+1), &
246	fym(nys-1:nyn+1,nxl-1:nxr+1), fyp(nys-1:nyn+1,nxl-1:nxr+1), &
247	nzb_s_inner(nys-1:nyn+1,nxl-1:nxr+1), &
248	nzb_s_outer(nys-1:nyn+1,nxl-1:nxr+1), &
249	nzb_u_inner(nys-1:nyn+1,nxl-1:nxr+1), &
250	nzb_u_outer(nys-1:nyn+1,nxl-1:nxr+1), &
251	nzb_v_inner(nys-1:nyn+1,nxl-1:nxr+1), &
252	nzb_v_outer(nys-1:nyn+1,nxl-1:nxr+1), &
253	nzb_w_inner(nys-1:nyn+1,nxl-1:nxr+1), &
254	nzb_w_outer(nys-1:nyn+1,nxl-1:nxr+1), &
255	nzb_diff_s_inner(nys-1:nyn+1,nxl-1:nxr+1), &
256	nzb_diff_s_outer(nys-1:nyn+1,nxl-1:nxr+1), &
257	nzb_diff_u(nys-1:nyn+1,nxl-1:nxr+1), &
258	nzb_diff_v(nys-1:nyn+1,nxl-1:nxr+1), &
259	nzb_2d(nys-1:nyn+1,nxl-1:nxr+1), &
260	wall_e_x(nys-1:nyn+1,nxl-1:nxr+1), &
261	wall_e_y(nys-1:nyn+1,nxl-1:nxr+1), &
262	wall_u(nys-1:nyn+1,nxl-1:nxr+1), &
263	wall_v(nys-1:nyn+1,nxl-1:nxr+1), &
264	wall_w_x(nys-1:nyn+1,nxl-1:nxr+1), &
265	wall_w_y(nys-1:nyn+1,nxl-1:nxr+1) )
266
267	ALLOCATE( l_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
268
269	nzb_s_inner = nzb; nzb_s_outer = nzb
270	nzb_u_inner = nzb; nzb_u_outer = nzb
271	nzb_v_inner = nzb; nzb_v_outer = nzb
272	nzb_w_inner = nzb; nzb_w_outer = nzb
273
274	!
275	!-- Define vertical gridpoint from (or to) which on the usual finite difference
276	!-- form (which does not use surface fluxes) is applied
277	IF ( prandtl_layer .OR. use_surface_fluxes ) THEN
278	nzb_diff = nzb + 2
279	ELSE
280	nzb_diff = nzb + 1
281	ENDIF
282	IF ( use_top_fluxes ) THEN
283	nzt_diff = nzt - 1
284	ELSE
285	nzt_diff = nzt
286	ENDIF
287
288	nzb_diff_s_inner = nzb_diff; nzb_diff_s_outer = nzb_diff
289	nzb_diff_u = nzb_diff; nzb_diff_v = nzb_diff
290
291	wall_e_x = 0.0; wall_e_y = 0.0; wall_u = 0.0; wall_v = 0.0
292	wall_w_x = 0.0; wall_w_y = 0.0
293	fwxp = 1.0; fwxm = 1.0; fwyp = 1.0; fwym = 1.0
294	fxp = 1.0; fxm = 1.0; fyp = 1.0; fym = 1.0
295
296	!
297	!-- Initialize near-wall mixing length l_wall only in the vertical direction
298	!-- for the moment,
299	!-- multiplication with wall_adjustment_factor near the end of this routine
300	l_wall(nzb,:,:) = l_grid(1)
301	DO k = nzb+1, nzt
302	l_wall(k,:,:) = l_grid(k)
303	ENDDO
304	l_wall(nzt+1,:,:) = l_grid(nzt)
305
306	ALLOCATE ( vertical_influence(nzb:nzt) )
307	DO k = 1, nzt
308	vertical_influence(k) = MIN ( INT( l_grid(k) / &
309	( wall_adjustment_factor * dzw(k) ) + 0.5 ), nzt - k )
310	ENDDO
311
312	DO k = 1, MAXVAL( nzb_s_inner )
313	IF ( l_grid(k) > 1.5 * dx * wall_adjustment_factor .OR. &
314	l_grid(k) > 1.5 * dy * wall_adjustment_factor ) THEN
315	WRITE( message_string, * ) 'grid anisotropy exceeds ', &
316	'threshold given by only local', &
317	' &horizontal reduction of near_wall ', &
318	'mixing length l_wall', &
319	' &starting from height level k = ', k, '.'
320	CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
321	EXIT
322	ENDIF
323	ENDDO
324	vertical_influence(0) = vertical_influence(1)
325
326	DO i = nxl-1, nxr+1
327	DO j = nys-1, nyn+1
328	DO k = nzb_s_inner(j,i) + 1, &
329	nzb_s_inner(j,i) + vertical_influence(nzb_s_inner(j,i))
330	l_wall(k,j,i) = zu(k) - zw(nzb_s_inner(j,i))
331	ENDDO
332	ENDDO
333	ENDDO
334
335	!
336	!-- Set outer and inner index arrays for non-flat topography.
337	!-- Here consistency checks concerning domain size and periodicity are
338	!-- necessary.
339	!-- Within this SELECT CASE structure only nzb_local is initialized
340	!-- individually depending on the chosen topography type, all other index
341	!-- arrays are initialized further below.
342	SELECT CASE ( TRIM( topography ) )
343
344	CASE ( 'flat' )
345	!
346	!-- No actions necessary
347
348	CASE ( 'single_building' )
349	!
350	!-- Single rectangular building, by default centered in the middle of the
351	!-- total domain
352	blx = NINT( building_length_x / dx )
353	bly = NINT( building_length_y / dy )
354	bh = NINT( building_height / dz )
355
356	IF ( building_wall_left == 9999999.9 ) THEN
357	building_wall_left = ( nx + 1 - blx ) / 2 * dx
358	ENDIF
359	bxl = NINT( building_wall_left / dx )
360	bxr = bxl + blx
361
362	IF ( building_wall_south == 9999999.9 ) THEN
363	building_wall_south = ( ny + 1 - bly ) / 2 * dy
364	ENDIF
365	bys = NINT( building_wall_south / dy )
366	byn = bys + bly
367
368	!
369	!-- Building size has to meet some requirements
370	IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR. &
371	( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) ) THEN
372	WRITE( message_string, * ) 'inconsistent building parameters:', &
373	'& bxl=', bxl, 'bxr=', bxr, 'bys=', bys, &
374	'byn=', byn, 'nx=', nx, 'ny=', ny
375	CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
376	ENDIF
377
378	!
379	!-- Define the building.
380	nzb_local = 0
381	nzb_local(bys:byn,bxl:bxr) = bh
382
383	CASE ( 'single_street_canyon' )
384	!
385	!-- Single quasi-2D street canyon of infinite length in x or y direction.
386	!-- The canyon is centered in the other direction by default.
387	IF ( canyon_width_x /= 9999999.9 ) THEN
388	!
389	!-- Street canyon in y direction
390	cwx = NINT( canyon_width_x / dx )
391	IF ( canyon_wall_left == 9999999.9 ) THEN
392	canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
393	ENDIF
394	cxl = NINT( canyon_wall_left / dx )
395	cxr = cxl + cwx
396
397	ELSEIF ( canyon_width_y /= 9999999.9 ) THEN
398	!
399	!-- Street canyon in x direction
400	cwy = NINT( canyon_width_y / dy )
401	IF ( canyon_wall_south == 9999999.9 ) THEN
402	canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
403	ENDIF
404	cys = NINT( canyon_wall_south / dy )
405	cyn = cys + cwy
406
407	ELSE
408
409	message_string = 'no street canyon width given'
410	CALL message( 'init_grid', 'PA0204', 1, 2, 0, 6, 0 )
411
412	ENDIF
413
414	ch = NINT( canyon_height / dz )
415	dp_level_ind_b = ch
416	!
417	!-- Street canyon size has to meet some requirements
418	IF ( canyon_width_x /= 9999999.9 ) THEN
419	IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR. &
420	( ch < 3 ) ) THEN
421	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
422	'&cxl=', cxl, 'cxr=', cxr, &
423	'cwx=', cwx, &
424	'ch=', ch, 'nx=', nx, 'ny=', ny
425	CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 )
426	ENDIF
427	ELSEIF ( canyon_width_y /= 9999999.9 ) THEN
428	IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR. &
429	( ch < 3 ) ) THEN
430	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
431	'&cys=', cys, 'cyn=', cyn, &
432	'cwy=', cwy, &
433	'ch=', ch, 'nx=', nx, 'ny=', ny
434	CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 )
435	ENDIF
436	ENDIF
437	IF ( canyon_width_x /= 9999999.9 .AND. canyon_width_y /= 9999999.9 ) &
438	THEN
439	message_string = 'inconsistent canyon parameters:' // &
440	'&street canyon can only be oriented' // &
441	'&either in x- or in y-direction'
442	CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 )
443	ENDIF
444
445	nzb_local = ch
446	IF ( canyon_width_x /= 9999999.9 ) THEN
447	nzb_local(:,cxl+1:cxr-1) = 0
448	ELSEIF ( canyon_width_y /= 9999999.9 ) THEN
449	nzb_local(cys+1:cyn-1,:) = 0
450	ENDIF
451
452	CASE ( 'read_from_file' )
453	!
454	!-- Arbitrary irregular topography data in PALM format (exactly matching
455	!-- the grid size and total domain size)
456	OPEN( 90, FILE='TOPOGRAPHY_DATA', STATUS='OLD', FORM='FORMATTED', &
457	ERR=10 )
458	DO j = ny, 0, -1
459	READ( 90, *, ERR=11, END=11 ) ( topo_height(j,i), i = 0, nx )
460	ENDDO
461	!
462	!-- Calculate the index height of the topography
463	DO i = 0, nx
464	DO j = 0, ny
465	nzb_local(j,i) = NINT( topo_height(j,i) / dz )
466	ENDDO
467	ENDDO
468	!
469	!-- Add cyclic boundaries (additional layers are for calculating flag
470	!-- arrays needed for the multigrid sover)
471	nzb_local(-gls:-1,0:nx) = nzb_local(ny-gls+1:ny,0:nx)
472	nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,0:nx)
473	nzb_local(:,-gls:-1) = nzb_local(:,nx-gls+1:nx)
474	nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
475
476	GOTO 12
477
478	10 message_string = 'file TOPOGRAPHY_DATA does not exist'
479	CALL message( 'init_grid', 'PA0208', 1, 2, 0, 6, 0 )
480
481	11 message_string = 'errors in file TOPOGRAPHY_DATA'
482	CALL message( 'init_grid', 'PA0209', 1, 2, 0, 6, 0 )
483
484	12 CLOSE( 90 )
485
486	CASE DEFAULT
487	!
488	!-- The DEFAULT case is reached either if the parameter topography
489	!-- contains a wrong character string or if the user has defined a special
490	!-- case in the user interface. There, the subroutine user_init_grid
491	!-- checks which of these two conditions applies.
492	CALL user_init_grid( gls, nzb_local )
493
494	END SELECT
495
496	!
497	!-- Test output of nzb_local -1:ny+1,-1:nx+1
498	! WRITE (9,) ' nzb_local *'
499	! DO j = ny+1, -1, -1
500	! WRITE (9,'(194(1X,I2))') ( nzb_local(j,i), i = -1, nx+1 )
501	! ENDDO
502
503	!
504	!-- Consistency checks and index array initialization are only required for
505	!-- non-flat topography, also the initialization of topography height arrays
506	!-- zu_s_inner and zw_w_inner
507	IF ( TRIM( topography ) /= 'flat' ) THEN
508
509	!
510	!-- Consistency checks
511	IF ( MINVAL( nzb_local ) < 0 .OR. MAXVAL( nzb_local ) > nz + 1 ) THEN
512	WRITE( message_string, * ) 'nzb_local values are outside the', &
513	'model domain', &
514	'&MINVAL( nzb_local ) = ', MINVAL(nzb_local), &
515	'&MAXVAL( nzb_local ) = ', MAXVAL(nzb_local)
516	CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
517	ENDIF
518
519	IF ( bc_lr == 'cyclic' ) THEN
520	IF ( ANY( nzb_local(:,-1) /= nzb_local(:,nx) ) .OR. &
521	ANY( nzb_local(:,0) /= nzb_local(:,nx+1) ) ) THEN
522	message_string = 'nzb_local does not fulfill cyclic' // &
523	' boundary condition in x-direction'
524	CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
525	ENDIF
526	ENDIF
527	IF ( bc_ns == 'cyclic' ) THEN
528	IF ( ANY( nzb_local(-1,:) /= nzb_local(ny,:) ) .OR. &
529	ANY( nzb_local(0,:) /= nzb_local(ny+1,:) ) ) THEN
530	message_string = 'nzb_local does not fulfill cyclic' // &
531	' boundary condition in y-direction'
532	CALL message( 'init_grid', 'PA0212', 1, 2, 0, 6, 0 )
533	ENDIF
534	ENDIF
535
536	IF ( topography_grid_convention == 'cell_edge' ) THEN
537	!
538	!-- The array nzb_local as defined using the 'cell_edge' convention
539	!-- describes the actual total size of topography which is defined at the
540	!-- cell edges where u=0 on the topography walls in x-direction and v=0
541	!-- on the topography walls in y-direction. However, PALM uses individual
542	!-- arrays nzb_u\|v\|w\|s_inner\|outer that are based on nzb_s_inner.
543	!-- Therefore, the extent of topography in nzb_local is now reduced by
544	!-- 1dx at the E topography walls and by 1dy at the N topography walls
545	!-- to form the basis for nzb_s_inner.
546	DO j = -gls, ny + gls
547	DO i = -gls, nx
548	nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j,i+1) )
549	ENDDO
550	ENDDO
551	!-- apply cyclic boundary conditions in x-direction
552	!(ist das erforderlich? Ursache von Seung Bus Fehler?)
553	nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
554	DO i = -gls, nx + gls
555	DO j = -gls, ny
556	nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j+1,i) )
557	ENDDO
558	ENDDO
559	!-- apply cyclic boundary conditions in y-direction
560	!(ist das erforderlich? Ursache von Seung Bus Fehler?)
561	nzb_local(ny+1:ny+gls,:) = nzb_local(0:gls-1,:)
562	ENDIF
563
564	!
565	!-- Initialize index arrays nzb_s_inner and nzb_w_inner
566	nzb_s_inner = nzb_local(nys-1:nyn+1,nxl-1:nxr+1)
567	nzb_w_inner = nzb_local(nys-1:nyn+1,nxl-1:nxr+1)
568
569	!
570	!-- Initialize remaining index arrays:
571	!-- first pre-initialize them with nzb_s_inner...
572	nzb_u_inner = nzb_s_inner
573	nzb_u_outer = nzb_s_inner
574	nzb_v_inner = nzb_s_inner
575	nzb_v_outer = nzb_s_inner
576	nzb_w_outer = nzb_s_inner
577	nzb_s_outer = nzb_s_inner
578
579	!
580	!-- ...then extend pre-initialized arrays in their according directions
581	!-- based on nzb_local using nzb_tmp as a temporary global index array
582
583	!
584	!-- nzb_s_outer:
585	!-- extend nzb_local east-/westwards first, then north-/southwards
586	nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
587	DO j = -1, ny + 1
588	DO i = 0, nx
589	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i), &
590	nzb_local(j,i+1) )
591	ENDDO
592	ENDDO
593	DO i = nxl, nxr
594	DO j = nys, nyn
595	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
596	nzb_tmp(j+1,i) )
597	ENDDO
598	!
599	!-- non-cyclic boundary conditions (overwritten by call of
600	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
601	IF ( nys == 0 ) THEN
602	j = -1
603	nzb_s_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
604	ENDIF
605	IF ( nys == ny ) THEN
606	j = ny + 1
607	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
608	ENDIF
609	ENDDO
610	!
611	!-- nzb_w_outer:
612	!-- identical to nzb_s_outer
613	nzb_w_outer = nzb_s_outer
614
615	!
616	!-- nzb_u_inner:
617	!-- extend nzb_local rightwards only
618	nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
619	DO j = -1, ny + 1
620	DO i = 0, nx + 1
621	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )
622	ENDDO
623	ENDDO
624	nzb_u_inner = nzb_tmp(nys-1:nyn+1,nxl-1:nxr+1)
625
626	!
627	!-- nzb_u_outer:
628	!-- extend current nzb_tmp (nzb_u_inner) north-/southwards
629	DO i = nxl, nxr
630	DO j = nys, nyn
631	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
632	nzb_tmp(j+1,i) )
633	ENDDO
634	!
635	!-- non-cyclic boundary conditions (overwritten by call of
636	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
637	IF ( nys == 0 ) THEN
638	j = -1
639	nzb_u_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
640	ENDIF
641	IF ( nys == ny ) THEN
642	j = ny + 1
643	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
644	ENDIF
645	ENDDO
646
647	!
648	!-- nzb_v_inner:
649	!-- extend nzb_local northwards only
650	nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
651	DO i = -1, nx + 1
652	DO j = 0, ny + 1
653	nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )
654	ENDDO
655	ENDDO
656	nzb_v_inner = nzb_tmp(nys-1:nyn+1,nxl-1:nxr+1)
657
658	!
659	!-- nzb_v_outer:
660	!-- extend current nzb_tmp (nzb_v_inner) right-/leftwards
661	DO j = nys, nyn
662	DO i = nxl, nxr
663	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i), &
664	nzb_tmp(j,i+1) )
665	ENDDO
666	!
667	!-- non-cyclic boundary conditions (overwritten by call of
668	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
669	IF ( nxl == 0 ) THEN
670	i = -1
671	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i+1), nzb_tmp(j,i) )
672	ENDIF
673	IF ( nxr == nx ) THEN
674	i = nx + 1
675	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i) )
676	ENDIF
677	ENDDO
678
679	!
680	!-- Exchange of lateral boundary values (parallel computers) and cyclic
681	!-- boundary conditions, if applicable.
682	!-- Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local
683	!-- they do not require exchange and are not included here.
684	CALL exchange_horiz_2d_int( nzb_u_inner )
685	CALL exchange_horiz_2d_int( nzb_u_outer )
686	CALL exchange_horiz_2d_int( nzb_v_inner )
687	CALL exchange_horiz_2d_int( nzb_v_outer )
688	CALL exchange_horiz_2d_int( nzb_w_outer )
689	CALL exchange_horiz_2d_int( nzb_s_outer )
690
691	!
692	!-- Allocate and set the arrays containing the topography height
693	IF ( myid == 0 ) THEN
694
695	ALLOCATE( zu_s_inner(0:nx+1,0:ny+1), zw_w_inner(0:nx+1,0:ny+1) )
696
697	DO i = 0, nx + 1
698	DO j = 0, ny + 1
699	zu_s_inner(i,j) = zu(nzb_local(j,i))
700	zw_w_inner(i,j) = zw(nzb_local(j,i))
701	ENDDO
702	ENDDO
703
704	ENDIF
705
706	ENDIF
707
708	!
709	!-- Preliminary: to be removed after completion of the topography code!
710	!-- Set the former default k index arrays nzb_2d
711	nzb_2d = nzb
712
713	!
714	!-- Set the individual index arrays which define the k index from which on
715	!-- the usual finite difference form (which does not use surface fluxes) is
716	!-- applied
717	IF ( prandtl_layer .OR. use_surface_fluxes ) THEN
718	nzb_diff_u = nzb_u_inner + 2
719	nzb_diff_v = nzb_v_inner + 2
720	nzb_diff_s_inner = nzb_s_inner + 2
721	nzb_diff_s_outer = nzb_s_outer + 2
722	ELSE
723	nzb_diff_u = nzb_u_inner + 1
724	nzb_diff_v = nzb_v_inner + 1
725	nzb_diff_s_inner = nzb_s_inner + 1
726	nzb_diff_s_outer = nzb_s_outer + 1
727	ENDIF
728
729	!
730	!-- Calculation of wall switches and factors required by diffusion_u/v.f90 and
731	!-- for limitation of near-wall mixing length l_wall further below
732	corner_nl = 0
733	corner_nr = 0
734	corner_sl = 0
735	corner_sr = 0
736	wall_l = 0
737	wall_n = 0
738	wall_r = 0
739	wall_s = 0
740
741	DO i = nxl, nxr
742	DO j = nys, nyn
743	!
744	!-- u-component
745	IF ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) THEN
746	wall_u(j,i) = 1.0 ! north wall (location of adjacent fluid)
747	fym(j,i) = 0.0
748	fyp(j,i) = 1.0
749	ELSEIF ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) THEN
750	wall_u(j,i) = 1.0 ! south wall (location of adjacent fluid)
751	fym(j,i) = 1.0
752	fyp(j,i) = 0.0
753	ENDIF
754	!
755	!-- v-component
756	IF ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) THEN
757	wall_v(j,i) = 1.0 ! rigth wall (location of adjacent fluid)
758	fxm(j,i) = 0.0
759	fxp(j,i) = 1.0
760	ELSEIF ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) ) THEN
761	wall_v(j,i) = 1.0 ! left wall (location of adjacent fluid)
762	fxm(j,i) = 1.0
763	fxp(j,i) = 0.0
764	ENDIF
765	!
766	!-- w-component, also used for scalars, separate arrays for shear
767	!-- production of tke
768	IF ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) THEN
769	wall_e_y(j,i) = 1.0 ! north wall (location of adjacent fluid)
770	wall_w_y(j,i) = 1.0
771	fwym(j,i) = 0.0
772	fwyp(j,i) = 1.0
773	ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) THEN
774	wall_e_y(j,i) = -1.0 ! south wall (location of adjacent fluid)
775	wall_w_y(j,i) = 1.0
776	fwym(j,i) = 1.0
777	fwyp(j,i) = 0.0
778	ENDIF
779	IF ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) THEN
780	wall_e_x(j,i) = 1.0 ! right wall (location of adjacent fluid)
781	wall_w_x(j,i) = 1.0
782	fwxm(j,i) = 0.0
783	fwxp(j,i) = 1.0
784	ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) THEN
785	wall_e_x(j,i) = -1.0 ! left wall (location of adjacent fluid)
786	wall_w_x(j,i) = 1.0
787	fwxm(j,i) = 1.0
788	fwxp(j,i) = 0.0
789	ENDIF
790	!
791	!-- Wall and corner locations inside buildings for limitation of
792	!-- near-wall mixing length l_wall
793	IF ( nzb_s_inner(j,i) > nzb_s_inner(j+1,i) ) THEN
794
795	wall_n(j,i) = nzb_s_inner(j+1,i) + 1 ! North wall
796
797	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
798	corner_nl(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northleft corner
799	nzb_s_inner(j,i-1) ) + 1
800	ENDIF
801
802	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
803	corner_nr(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northright corner
804	nzb_s_inner(j,i+1) ) + 1
805	ENDIF
806
807	ENDIF
808
809	IF ( nzb_s_inner(j,i) > nzb_s_inner(j-1,i) ) THEN
810
811	wall_s(j,i) = nzb_s_inner(j-1,i) + 1 ! South wall
812	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
813	corner_sl(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southleft corner
814	nzb_s_inner(j,i-1) ) + 1
815	ENDIF
816
817	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
818	corner_sr(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southright corner
819	nzb_s_inner(j,i+1) ) + 1
820	ENDIF
821
822	ENDIF
823
824	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
825	wall_l(j,i) = nzb_s_inner(j,i-1) + 1 ! Left wall
826	ENDIF
827
828	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
829	wall_r(j,i) = nzb_s_inner(j,i+1) + 1 ! Right wall
830	ENDIF
831
832	ENDDO
833	ENDDO
834
835	!
836	!-- Calculate wall flag arrays for the multigrid method
837	IF ( psolver == 'multigrid' ) THEN
838	!
839	!-- Gridpoint increment of the current level
840	inc = 1
841
842	DO l = maximum_grid_level, 1 , -1
843
844	nxl_l = nxl_mg(l)
845	nxr_l = nxr_mg(l)
846	nys_l = nys_mg(l)
847	nyn_l = nyn_mg(l)
848	nzt_l = nzt_mg(l)
849
850	!
851	!-- Assign the flag level to be calculated
852	SELECT CASE ( l )
853	CASE ( 1 )
854	flags => wall_flags_1
855	CASE ( 2 )
856	flags => wall_flags_2
857	CASE ( 3 )
858	flags => wall_flags_3
859	CASE ( 4 )
860	flags => wall_flags_4
861	CASE ( 5 )
862	flags => wall_flags_5
863	CASE ( 6 )
864	flags => wall_flags_6
865	CASE ( 7 )
866	flags => wall_flags_7
867	CASE ( 8 )
868	flags => wall_flags_8
869	CASE ( 9 )
870	flags => wall_flags_9
871	CASE ( 10 )
872	flags => wall_flags_10
873	END SELECT
874
875	!
876	!-- Depending on the grid level, set the respective bits in case of
877	!-- neighbouring walls
878	!-- Bit 0: wall to the bottom
879	!-- Bit 1: wall to the top (not realized in remaining PALM code so far)
880	!-- Bit 2: wall to the south
881	!-- Bit 3: wall to the north
882	!-- Bit 4: wall to the left
883	!-- Bit 5: wall to the right
884	!-- Bit 6: inside building
885
886	flags = 0
887
888	DO i = nxl_l-1, nxr_l+1
889	DO j = nys_l-1, nyn_l+1
890	DO k = nzb, nzt_l+1
891
892	!
893	!-- Inside/outside building (inside building does not need
894	!-- further tests for walls)
895	IF ( kinc <= nzb_local(jinc,i*inc) ) THEN
896
897	flags(k,j,i) = IBSET( flags(k,j,i), 6 )
898
899	ELSE
900	!
901	!-- Bottom wall
902	IF ( (k-1)inc <= nzb_local(jinc,i*inc) ) THEN
903	flags(k,j,i) = IBSET( flags(k,j,i), 0 )
904	ENDIF
905	!
906	!-- South wall
907	IF ( kinc <= nzb_local((j-1)inc,i*inc) ) THEN
908	flags(k,j,i) = IBSET( flags(k,j,i), 2 )
909	ENDIF
910	!
911	!-- North wall
912	IF ( kinc <= nzb_local((j+1)inc,i*inc) ) THEN
913	flags(k,j,i) = IBSET( flags(k,j,i), 3 )
914	ENDIF
915	!
916	!-- Left wall
917	IF ( kinc <= nzb_local(jinc,(i-1)*inc) ) THEN
918	flags(k,j,i) = IBSET( flags(k,j,i), 4 )
919	ENDIF
920	!
921	!-- Right wall
922	IF ( kinc <= nzb_local(jinc,(i+1)*inc) ) THEN
923	flags(k,j,i) = IBSET( flags(k,j,i), 5 )
924	ENDIF
925
926	ENDIF
927
928	ENDDO
929	ENDDO
930	ENDDO
931
932	!
933	!-- Test output of flag arrays
934	! i = nxl_l
935	! WRITE (9,*) ' '
936	! WRITE (9,) ' mg level ', l, ' *', mg_switch_to_pe0_level
937	! WRITE (9,*) ' inc=', inc, ' i =', nxl_l
938	! WRITE (9,*) ' nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
939	! DO k = nzt_l+1, nzb, -1
940	! WRITE (9,'(194(1X,I2))') ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
941	! ENDDO
942
943	inc = inc * 2
944
945	ENDDO
946
947	ENDIF
948
949	!
950	!-- In case of topography: limit near-wall mixing length l_wall further:
951	!-- Go through all points of the subdomain one by one and look for the closest
952	!-- surface
953	IF ( TRIM(topography) /= 'flat' ) THEN
954	DO i = nxl, nxr
955	DO j = nys, nyn
956
957	nzb_si = nzb_s_inner(j,i)
958	vi = vertical_influence(nzb_si)
959
960	IF ( wall_n(j,i) > 0 ) THEN
961	!
962	!-- North wall (y distance)
963	DO k = wall_n(j,i), nzb_si
964	l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5 * dy )
965	ENDDO
966	!
967	!-- Above North wall (yz distance)
968	DO k = nzb_si + 1, nzb_si + vi
969	l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), &
970	SQRT( 0.25 * dy**2 + &
971	( zu(k) - zw(nzb_si) )**2 ) )
972	ENDDO
973	!
974	!-- Northleft corner (xy distance)
975	IF ( corner_nl(j,i) > 0 ) THEN
976	DO k = corner_nl(j,i), nzb_si
977	l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
978	0.5 * SQRT( dx2 + dy2 ) )
979	ENDDO
980	!
981	!-- Above Northleft corner (xyz distance)
982	DO k = nzb_si + 1, nzb_si + vi
983	l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
984	SQRT( 0.25 * (dx2 + dy2) + &
985	( zu(k) - zw(nzb_si) )**2 ) )
986	ENDDO
987	ENDIF
988	!
989	!-- Northright corner (xy distance)
990	IF ( corner_nr(j,i) > 0 ) THEN
991	DO k = corner_nr(j,i), nzb_si
992	l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
993	0.5 * SQRT( dx2 + dy2 ) )
994	ENDDO
995	!
996	!-- Above northright corner (xyz distance)
997	DO k = nzb_si + 1, nzb_si + vi
998	l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
999	SQRT( 0.25 * (dx2 + dy2) + &
1000	( zu(k) - zw(nzb_si) )**2 ) )
1001	ENDDO
1002	ENDIF
1003	ENDIF
1004
1005	IF ( wall_s(j,i) > 0 ) THEN
1006	!
1007	!-- South wall (y distance)
1008	DO k = wall_s(j,i), nzb_si
1009	l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), 0.5 * dy )
1010	ENDDO
1011	!
1012	!-- Above south wall (yz distance)
1013	DO k = nzb_si + 1, &
1014	nzb_si + vi
1015	l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), &
1016	SQRT( 0.25 * dy**2 + &
1017	( zu(k) - zw(nzb_si) )**2 ) )
1018	ENDDO
1019	!
1020	!-- Southleft corner (xy distance)
1021	IF ( corner_sl(j,i) > 0 ) THEN
1022	DO k = corner_sl(j,i), nzb_si
1023	l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
1024	0.5 * SQRT( dx2 + dy2 ) )
1025	ENDDO
1026	!
1027	!-- Above southleft corner (xyz distance)
1028	DO k = nzb_si + 1, nzb_si + vi
1029	l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
1030	SQRT( 0.25 * (dx2 + dy2) + &
1031	( zu(k) - zw(nzb_si) )**2 ) )
1032	ENDDO
1033	ENDIF
1034	!
1035	!-- Southright corner (xy distance)
1036	IF ( corner_sr(j,i) > 0 ) THEN
1037	DO k = corner_sr(j,i), nzb_si
1038	l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
1039	0.5 * SQRT( dx2 + dy2 ) )
1040	ENDDO
1041	!
1042	!-- Above southright corner (xyz distance)
1043	DO k = nzb_si + 1, nzb_si + vi
1044	l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
1045	SQRT( 0.25 * (dx2 + dy2) + &
1046	( zu(k) - zw(nzb_si) )**2 ) )
1047	ENDDO
1048	ENDIF
1049
1050	ENDIF
1051
1052	IF ( wall_l(j,i) > 0 ) THEN
1053	!
1054	!-- Left wall (x distance)
1055	DO k = wall_l(j,i), nzb_si
1056	l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), 0.5 * dx )
1057	ENDDO
1058	!
1059	!-- Above left wall (xz distance)
1060	DO k = nzb_si + 1, nzb_si + vi
1061	l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), &
1062	SQRT( 0.25 * dx**2 + &
1063	( zu(k) - zw(nzb_si) )**2 ) )
1064	ENDDO
1065	ENDIF
1066
1067	IF ( wall_r(j,i) > 0 ) THEN
1068	!
1069	!-- Right wall (x distance)
1070	DO k = wall_r(j,i), nzb_si
1071	l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5 * dx )
1072	ENDDO
1073	!
1074	!-- Above right wall (xz distance)
1075	DO k = nzb_si + 1, nzb_si + vi
1076	l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), &
1077	SQRT( 0.25 * dx**2 + &
1078	( zu(k) - zw(nzb_si) )**2 ) )
1079	ENDDO
1080
1081	ENDIF
1082
1083	ENDDO
1084	ENDDO
1085
1086	ENDIF
1087
1088	!
1089	!-- Multiplication with wall_adjustment_factor
1090	l_wall = wall_adjustment_factor * l_wall
1091
1092	!
1093	!-- Need to set lateral boundary conditions for l_wall
1094	CALL exchange_horiz( l_wall )
1095
1096	DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, &
1097	nzb_tmp, vertical_influence, wall_l, wall_n, wall_r, wall_s )
1098
1099
1100	END SUBROUTINE init_grid

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

Context Navigation

source: palm/trunk/SOURCE/init_grid.f90 @ 449

Download in other formats: