Home

Context Navigation

init_grid.f90 @ 546

Last change on this file since 546 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

Rev	Line
[1]	1	SUBROUTINE init_grid
	2
	3	!------------------------------------------------------------------------------!
[254]	4	! Current revisions:
[1]	5	! -----------------
[392]	6	!
[139]	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Id: init_grid.f90 392 2009-09-24 10:39:14Z maronga $
	11	!
[392]	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	!
[226]	16	! 217 2008-12-09 18:00:48Z letzel
	17	! +topography_grid_convention
	18	!
[139]	19	! 134 2007-11-21 07:28:38Z letzel
[134]	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.
[114]	24	! Calculation of wall flag arrays
[1]	25	!
[98]	26	! 94 2007-06-01 15:25:22Z raasch
	27	! Grid definition for ocean version
	28	!
[77]	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	!
[39]	33	! 19 2007-02-23 04:53:48Z raasch
	34	! Setting of nzt_diff
	35	!
[3]	36	! RCS Log replace by Id keyword, revision history cleaned up
	37	!
[1]	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
[240]	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
[1]	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
[254]	84	message_string = 'missing dz'
	85	CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
[1]	86	ELSEIF ( dz <= 0.0 ) THEN
[254]	87	WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
	88	CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
[1]	89	ENDIF
[94]	90
[1]	91	!
[94]	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
[1]	102
[94]	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
[1]	113
	114	!
[94]	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) )
[1]	122
[94]	123	ELSE
[1]	124	!
[94]	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	!
[1]	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
[114]	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), &
[1]	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	!
[19]	275	!-- Define vertical gridpoint from (or to) which on the usual finite difference
[1]	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
[19]	282	IF ( use_top_fluxes ) THEN
	283	nzt_diff = nzt - 1
	284	ELSE
	285	nzt_diff = nzt
	286	ENDIF
[1]	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
[254]	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 )
[1]	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
[274]	372	WRITE( message_string, * ) 'inconsistent building parameters:', &
	373	'& bxl=', bxl, 'bxr=', bxr, 'bys=', bys, &
	374	'byn=', byn, 'nx=', nx, 'ny=', ny
[254]	375	CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
[1]	376	ENDIF
	377
	378	!
[217]	379	!-- Define the building.
[1]	380	nzb_local = 0
[134]	381	nzb_local(bys:byn,bxl:bxr) = bh
[1]	382
[240]	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
[254]	408
	409	message_string = 'no street canyon width given'
	410	CALL message( 'init_grid', 'PA0204', 1, 2, 0, 6, 0 )
	411
[240]	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
[254]	421	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
[274]	422	'&cxl=', cxl, 'cxr=', cxr, &
	423	'cwx=', cwx, &
[254]	424	'ch=', ch, 'nx=', nx, 'ny=', ny
	425	CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 )
[240]	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
[254]	430	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
[274]	431	'&cys=', cys, 'cyn=', cyn, &
	432	'cwy=', cwy, &
[254]	433	'ch=', ch, 'nx=', nx, 'ny=', ny
	434	CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 )
[240]	435	ENDIF
	436	ENDIF
[274]	437	IF ( canyon_width_x /= 9999999.9 .AND. canyon_width_y /= 9999999.9 ) &
[240]	438	THEN
[274]	439	message_string = 'inconsistent canyon parameters:' // &
[254]	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 )
[240]	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
[1]	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
[114]	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)
[1]	475
	476	GOTO 12
[254]	477
	478	10 message_string = 'file TOPOGRAPHY_DATA does not exist'
	479	CALL message( 'init_grid', 'PA0208', 1, 2, 0, 6, 0 )
[1]	480
[254]	481	11 message_string = 'errors in file TOPOGRAPHY_DATA'
	482	CALL message( 'init_grid', 'PA0209', 1, 2, 0, 6, 0 )
[1]	483
	484	12 CLOSE( 90 )
	485
	486	CASE DEFAULT
	487	!
	488	!-- The DEFAULT case is reached either if the parameter topography
[217]	489	!-- contains a wrong character string or if the user has defined a special
[1]	490	!-- case in the user interface. There, the subroutine user_init_grid
	491	!-- checks which of these two conditions applies.
[114]	492	CALL user_init_grid( gls, nzb_local )
[1]	493
	494	END SELECT
	495
	496	!
[114]	497	!-- Test output of nzb_local -1:ny+1,-1:nx+1
[145]	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
[114]	502
	503	!
[1]	504	!-- Consistency checks and index array initialization are only required for
[217]	505	!-- non-flat topography, also the initialization of topography height arrays
[49]	506	!-- zu_s_inner and zw_w_inner
[1]	507	IF ( TRIM( topography ) /= 'flat' ) THEN
	508
	509	!
	510	!-- Consistency checks
	511	IF ( MINVAL( nzb_local ) < 0 .OR. MAXVAL( nzb_local ) > nz + 1 ) THEN
[274]	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)
[254]	516	CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
[1]	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
[254]	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 )
[1]	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
[254]	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 )
[1]	533	ENDIF
	534	ENDIF
	535
[217]	536	IF ( topography_grid_convention == 'cell_edge' ) THEN
[134]	537	!
[217]	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
[134]	550	ENDDO
[217]	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
[134]	558	ENDDO
[217]	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
[134]	563
[1]	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
[114]	586	nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
[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
[114]	618	nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
[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
[114]	650	nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
[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
[49]	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
[1]	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	!
[114]	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
[116]	884	!-- Bit 6: inside building
[114]	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
[145]	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
[114]	942
	943	inc = inc * 2
	944
	945	ENDDO
	946
	947	ENDIF
	948
	949	!
[1]	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
[75]	1094	CALL exchange_horiz( l_wall )
[1]	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 @ 546

Download in other formats: