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

Last change on this file since 392 was 392, checked in by raasch, 12 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(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
194              dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-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 ( k*inc <= nzb_local(j*inc,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(j*inc,i*inc) )  THEN
903                         flags(k,j,i) = IBSET( flags(k,j,i), 0 )
904                      ENDIF
905!
906!--                   South wall
907                      IF ( k*inc <= 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 ( k*inc <= 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 ( k*inc <= nzb_local(j*inc,(i-1)*inc) )  THEN
918                         flags(k,j,i) = IBSET( flags(k,j,i), 4 )
919                      ENDIF
920!
921!--                   Right wall
922                      IF ( k*inc <= nzb_local(j*inc,(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( dx**2 + dy**2 ) )
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 * (dx**2 + dy**2) + &
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( dx**2 + dy**2 ) )
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 * (dx**2 + dy**2) + &
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( dx**2 + dy**2 ) )
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 * (dx**2 + dy**2) + &
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( dx**2 + dy**2 ) )
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 * (dx**2 + dy**2) + &
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.