Home

Context Navigation

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

Last change on this file since 1905 was 1903, checked in by suehring, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 75.4 KB

Rev	Line
[1682]	1	!> @file init_grid.f90
[1743]	2	!-------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the terms
	6	! of the GNU General Public License as published by the Free Software Foundation,
	7	! either version 3 of the License, or (at your option) any later version.
	8	!
	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	12	!
	13	! You should have received a copy of the GNU General Public License along with
	14	! PALM. If not, see <http://www.gnu.org/licenses/>.
	15	!
[1818]	16	! Copyright 1997-2016 Leibniz Universitaet Hannover
[1743]	17	!-------------------------------------------------------------------------------!
[1036]	18	!
[254]	19	! Current revisions:
[1]	20	! -----------------
[1805]	21	!
[1903]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: init_grid.f90 1903 2016-05-09 11:20:41Z suehring $
	26	!
[1903]	27	! 1902 2016-05-09 11:18:56Z suehring
	28	! Set topography flags only for multigrid solver, not for multigrid_fast
	29	!
[1887]	30	! 1886 2016-04-21 11:20:47Z suehring
	31	! Bugfix: setting advection flags near walls
	32	! reformulated index values for nzb_v_inner
	33	! variable discriptions added in declaration block
	34	!
[1846]	35	! 1845 2016-04-08 08:29:13Z raasch
	36	! nzb_2d removed
	37	!
[1805]	38	! 1804 2016-04-05 16:30:18Z maronga
	39	! Removed code for parameter file check (__check)
	40	!
[1780]	41	! 1779 2016-03-03 08:01:28Z raasch
	42	! coupling_char is trimmed at every place it occurs, because it can have
	43	! different length now
	44	!
[1763]	45	! 1762 2016-02-25 12:31:13Z hellstea
	46	! Introduction of nested domain feature
	47	!
[1744]	48	! 1743 2016-01-13 10:23:51Z raasch
	49	! Bugfix for calculation of nzb_s_outer and nzb_u_outer at north boundary of
	50	! total domain
	51	!
[1692]	52	! 1691 2015-10-26 16:17:44Z maronga
	53	! Renamed prandtl_layer to constant_flux_layer.
	54	!
[1683]	55	! 1682 2015-10-07 23:56:08Z knoop
	56	! Code annotations made doxygen readable
	57	!
[1678]	58	! 1677 2015-10-02 13:25:23Z boeske
	59	! Bugfix: Ghost points are included in wall_flags_0 and wall_flags_00
	60	!
[1676]	61	! 1675 2015-10-02 08:28:59Z gronemeier
	62	! Bugfix: Definition of topography grid levels
	63	!
[1662]	64	! 1660 2015-09-21 08:15:16Z gronemeier
	65	! Bugfix: Definition of topography grid levels if vertical grid stretching
	66	! starts below the maximum topography height.
	67	!
[1581]	68	! 1580 2015-04-10 13:43:49Z suehring
	69	! Bugfix: setting flags for 5th order scheme near buildings
	70	!
[1576]	71	! 1575 2015-03-27 09:56:27Z raasch
	72	! adjustments for psolver-queries
	73	!
[1558]	74	! 1557 2015-03-05 16:43:04Z suehring
	75	! Adjustment for monotoinic limiter
	76	!
[1419]	77	! 1418 2014-06-06 13:05:08Z fricke
	78	! Bugfix: Change if-condition for stretched grid in the ocean, with the old
	79	! condition and a negative value for dz_stretch_level the condition
	80	! was always true for the whole model domain
	81	!
[1410]	82	! 1409 2014-05-23 12:11:32Z suehring
	83	! Bugfix: set wall_flags_0 at inflow and outflow boundary also for i <= nxlu
	84	! j <= nysv
	85	!
[1354]	86	! 1353 2014-04-08 15:21:23Z heinze
	87	! REAL constants provided with KIND-attribute
	88	!
[1323]	89	! 1322 2014-03-20 16:38:49Z raasch
	90	! REAL constants defined as wp-kind
	91	!
[1321]	92	! 1320 2014-03-20 08:40:49Z raasch
[1320]	93	! ONLY-attribute added to USE-statements,
	94	! kind-parameters added to all INTEGER and REAL declaration statements,
	95	! kinds are defined in new module kinds,
	96	! revision history before 2012 removed,
	97	! comment fields (!:) to be used for variable explanations added to
	98	! all variable declaration statements
[1321]	99	!
[1222]	100	! 1221 2013-09-10 08:59:13Z raasch
	101	! wall_flags_00 introduced to hold bits 32-63,
	102	! additional 3D-flag arrays for replacing the 2D-index array nzb_s_inner in
	103	! loops optimized for openACC (pres + flow_statistics)
	104	!
[1093]	105	! 1092 2013-02-02 11:24:22Z raasch
	106	! unused variables removed
	107	!
[1070]	108	! 1069 2012-11-28 16:18:43Z maronga
[1779]	109	! bugfix: added coupling_char to TOPOGRAPHY_DATA to allow topography in the
	110	! ocean model in case of coupled runs
[1070]	111	!
[1037]	112	! 1036 2012-10-22 13:43:42Z raasch
	113	! code put under GPL (PALM 3.9)
	114	!
[1017]	115	! 1015 2012-09-27 09:23:24Z raasch
	116	! lower index for calculating wall_flags_0 set to nzb_w_inner instead of
	117	! nzb_w_inner+1
	118	!
[997]	119	! 996 2012-09-07 10:41:47Z raasch
	120	! little reformatting
	121	!
[979]	122	! 978 2012-08-09 08:28:32Z fricke
	123	! Bugfix: nzb_max is set to nzt at non-cyclic lateral boundaries
	124	! Bugfix: Set wall_flags_0 for inflow boundary
	125	!
[928]	126	! 927 2012-06-06 19:15:04Z raasch
	127	! Wall flags are not set for multigrid method in case of masking method
	128	!
[865]	129	! 864 2012-03-27 15:10:33Z gryschka
[927]	130	! In case of ocean and Dirichlet bottom bc for u and v dzu_mg and ddzu_pres
	131	! were not correctly defined for k=1.
[865]	132	!
[863]	133	! 861 2012-03-26 14:18:34Z suehring
[861]	134	! Set wall_flags_0. The array is needed for degradation in ws-scheme near walls,
	135	! inflow and outflow boundaries as well as near the bottom and the top of the
[863]	136	! model domain.!
[861]	137	! Initialization of nzb_s_inner and nzb_w_inner.
	138	! gls has to be at least nbgp to do not exceed the array bounds of nzb_local
	139	! while setting wall_flags_0
	140	!
[844]	141	! 843 2012-02-29 15:16:21Z gryschka
	142	! In case of ocean and dirichlet bc for u and v at the bottom
	143	! the first u-level ist defined at same height as the first w-level
	144	!
[819]	145	! 818 2012-02-08 16:11:23Z maronga
	146	! Bugfix: topo_height is only required if topography is used. It is thus now
	147	! allocated in the topography branch
	148	!
[810]	149	! 809 2012-01-30 13:32:58Z maronga
	150	! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
	151	!
[808]	152	! 807 2012-01-25 11:53:51Z maronga
	153	! New cpp directive "__check" implemented which is used by check_namelist_files
	154	!
[1]	155	! Revision 1.1 1997/08/11 06:17:45 raasch
	156	! Initial revision (Testversion)
	157	!
	158	!
	159	! Description:
	160	! ------------
[1682]	161	!> Creating grid depending constants
[1]	162	!------------------------------------------------------------------------------!
[1682]	163	SUBROUTINE init_grid
	164
[1]	165
[1320]	166	USE arrays_3d, &
	167	ONLY: dd2zu, ddzu, ddzu_pres, ddzw, dzu, dzu_mg, dzw, dzw_mg, f1_mg, &
	168	f2_mg, f3_mg, l_grid, l_wall, zu, zw
	169
[1353]	170	USE control_parameters, &
[1320]	171	ONLY: bc_lr, bc_ns, building_height, building_length_x, &
	172	building_length_y, building_wall_left, building_wall_south, &
	173	canyon_height, canyon_wall_left, canyon_wall_south, &
[1691]	174	canyon_width_x, canyon_width_y, constant_flux_layer, &
	175	coupling_char, dp_level_ind_b, dz, dz_max, dz_stretch_factor, &
	176	dz_stretch_level, dz_stretch_level_index, ibc_uv_b, io_blocks, &
	177	io_group, inflow_l, inflow_n, inflow_r, inflow_s, &
	178	masking_method, maximum_grid_level, message_string, &
[1762]	179	momentum_advec, nest_domain, nest_bound_l, nest_bound_n, &
	180	nest_bound_r, nest_bound_s, ocean, outflow_l, outflow_n, &
	181	outflow_r, outflow_s, psolver, scalar_advec, topography, &
[1691]	182	topography_grid_convention, use_surface_fluxes, use_top_fluxes, &
	183	wall_adjustment_factor
[1320]	184
	185	USE grid_variables, &
	186	ONLY: ddx, ddx2, ddx2_mg, ddy, ddy2, ddy2_mg, dx, dx2, dy, dy2, fwxm, &
	187	fwxp, fwym, fwyp, fxm, fxp, fym, fyp, wall_e_x, wall_e_y, &
	188	wall_u, wall_v, wall_w_x, wall_w_y, zu_s_inner, zw_w_inner
	189
	190	USE indices, &
	191	ONLY: flags, nbgp, nx, nxl, nxlg, nxlu, nxl_mg, nxr, nxrg, nxr_mg, &
	192	ny, nyn, nyng, nyn_mg, nys, nysv, nys_mg, nysg, nz, nzb, &
[1845]	193	nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u, &
	194	nzb_diff_v, nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner, &
	195	nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner, &
	196	nzb_w_outer, nzt, nzt_diff, nzt_mg, rflags_invers, &
[1320]	197	rflags_s_inner, wall_flags_0, wall_flags_00, wall_flags_1, &
	198	wall_flags_10, wall_flags_2, wall_flags_3, wall_flags_4, &
	199	wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8, &
	200	wall_flags_9
	201
	202	USE kinds
	203
[1]	204	USE pegrid
	205
	206	IMPLICIT NONE
	207
[1886]	208	INTEGER(iwp) :: bh !< temporary vertical index of building height
	209	INTEGER(iwp) :: blx !< grid point number of building size along x
	210	INTEGER(iwp) :: bly !< grid point number of building size along y
	211	INTEGER(iwp) :: bxl !< index for left building wall
	212	INTEGER(iwp) :: bxr !< index for right building wall
	213	INTEGER(iwp) :: byn !< index for north building wall
	214	INTEGER(iwp) :: bys !< index for south building wall
	215	INTEGER(iwp) :: ch !< temporary vertical index for canyon height
	216	INTEGER(iwp) :: cwx !< grid point number of canyon size along x
	217	INTEGER(iwp) :: cwy !< grid point number of canyon size along y
	218	INTEGER(iwp) :: cxl !< index for left canyon wall
	219	INTEGER(iwp) :: cxr !< index for right canyon wall
	220	INTEGER(iwp) :: cyn !< index for north canyon wall
	221	INTEGER(iwp) :: cys !< index for south canyon wall
	222	INTEGER(iwp) :: gls !< number of lateral ghost points at total model
	223	!< domain boundaries required for multigrid solver
	224	INTEGER(iwp) :: i !< index variable along x
	225	INTEGER(iwp) :: ii !< loop variable for reading topography file
	226	INTEGER(iwp) :: inc !< incremental parameter for coarsening grid level
	227	INTEGER(iwp) :: j !< index variable along y
	228	INTEGER(iwp) :: k !< index variable along z
	229	INTEGER(iwp) :: l !< loop variable
	230	INTEGER(iwp) :: nxl_l !< index of left PE boundary for multigrid level
	231	INTEGER(iwp) :: nxr_l !< index of right PE boundary for multigrid level
	232	INTEGER(iwp) :: nyn_l !< index of north PE boundary for multigrid level
	233	INTEGER(iwp) :: nys_l !< index of south PE boundary for multigrid level
	234	INTEGER(iwp) :: nzb_si !< dummy index for local nzb_s_inner
	235	INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level
	236	INTEGER(iwp) :: vi !< dummy for vertical influence
[1]	237
[1886]	238	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: &
	239	vertical_influence !< number of vertical grid points above
	240	!< obstacle where adjustment of near-
	241	!< wall mixing length is required
	242
	243	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_nl !< index of
	244	!< north-left corner location to limit
	245	!< near-wall mixing length
	246	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_nr !< north-right
	247	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_sl !< south-left
	248	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_sr !< south-right
	249	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_l !< distance to
	250	!< adjacent left-facing
	251	!< wall
	252	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_n !< north-facing
	253	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_r !< right-facing
	254	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_s !< right-facing
	255	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography
	256	!< top at cell-center
	257	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp !< dummy to calculate
	258	!< topography indices
	259	!< on u- and v-grid
[1]	260
[1886]	261	LOGICAL :: flag_set = .FALSE. !< steering variable for advection flags
[1]	262
[1886]	263	REAL(wp) :: dx_l !< grid spacing along x on different multigrid level
	264	REAL(wp) :: dy_l !< grid spacing along y on different multigrid level
	265	REAL(wp) :: dz_stretched !< stretched vertical grid spacing
[861]	266
[1886]	267	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: topo_height !< input variable for
	268	!< topography height
[1]	269
[667]	270
[1]	271	!
[709]	272	!-- Calculation of horizontal array bounds including ghost layers
[667]	273	nxlg = nxl - nbgp
	274	nxrg = nxr + nbgp
	275	nysg = nys - nbgp
	276	nyng = nyn + nbgp
[709]	277
[667]	278	!
[1]	279	!-- Allocate grid arrays
[1353]	280	ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &
[667]	281	dzw(1:nzt+1), l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) )
[1]	282
	283	!
	284	!-- Compute height of u-levels from constant grid length and dz stretch factors
[1353]	285	IF ( dz == -1.0_wp ) THEN
[254]	286	message_string = 'missing dz'
	287	CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
[1353]	288	ELSEIF ( dz <= 0.0_wp ) THEN
[254]	289	WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
	290	CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
[1]	291	ENDIF
[94]	292
[1]	293	!
[94]	294	!-- Define the vertical grid levels
	295	IF ( .NOT. ocean ) THEN
	296	!
	297	!-- Grid for atmosphere with surface at z=0 (k=0, w-grid).
[843]	298	!-- The second u-level (k=1) corresponds to the top of the
[94]	299	!-- Prandtl-layer.
[667]	300
	301	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
[1353]	302	zu(0) = 0.0_wp
	303	! zu(0) = - dz * 0.5_wp
[667]	304	ELSE
[1353]	305	zu(0) = - dz * 0.5_wp
[667]	306	ENDIF
[1353]	307	zu(1) = dz * 0.5_wp
[1]	308
[94]	309	dz_stretch_level_index = nzt+1
	310	dz_stretched = dz
	311	DO k = 2, nzt+1
	312	IF ( dz_stretch_level <= zu(k-1) .AND. dz_stretched < dz_max ) THEN
	313	dz_stretched = dz_stretched * dz_stretch_factor
	314	dz_stretched = MIN( dz_stretched, dz_max )
	315	IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
	316	ENDIF
	317	zu(k) = zu(k-1) + dz_stretched
	318	ENDDO
[1]	319
	320	!
[94]	321	!-- Compute the w-levels. They are always staggered half-way between the
[843]	322	!-- corresponding u-levels. In case of dirichlet bc for u and v at the
	323	!-- ground the first u- and w-level (k=0) are defined at same height (z=0).
	324	!-- The top w-level is extrapolated linearly.
[1353]	325	zw(0) = 0.0_wp
[94]	326	DO k = 1, nzt
[1353]	327	zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
[94]	328	ENDDO
[1353]	329	zw(nzt+1) = zw(nzt) + 2.0_wp * ( zu(nzt+1) - zw(nzt) )
[1]	330
[94]	331	ELSE
[1]	332	!
[843]	333	!-- Grid for ocean with free water surface is at k=nzt (w-grid).
	334	!-- In case of neumann bc at the ground the first first u-level (k=0) lies
	335	!-- below the first w-level (k=0). In case of dirichlet bc the first u- and
	336	!-- w-level are defined at same height, but staggered from the second level.
	337	!-- The second u-level (k=1) corresponds to the top of the Prandtl-layer.
[1353]	338	zu(nzt+1) = dz * 0.5_wp
	339	zu(nzt) = - dz * 0.5_wp
[94]	340
	341	dz_stretch_level_index = 0
	342	dz_stretched = dz
	343	DO k = nzt-1, 0, -1
[1418]	344	!
	345	!-- The default value of dz_stretch_level is positive, thus the first
	346	!-- condition is always true. Hence, the second condition is necessary.
	347	IF ( dz_stretch_level >= zu(k+1) .AND. dz_stretch_level <= 0.0 &
	348	.AND. dz_stretched < dz_max ) THEN
[94]	349	dz_stretched = dz_stretched * dz_stretch_factor
	350	dz_stretched = MIN( dz_stretched, dz_max )
	351	IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
	352	ENDIF
	353	zu(k) = zu(k+1) - dz_stretched
	354	ENDDO
	355
	356	!
	357	!-- Compute the w-levels. They are always staggered half-way between the
[843]	358	!-- corresponding u-levels, except in case of dirichlet bc for u and v
	359	!-- at the ground. In this case the first u- and w-level are defined at
	360	!-- same height. The top w-level (nzt+1) is not used but set for
	361	!-- consistency, since w and all scalar variables are defined up tp nzt+1.
[94]	362	zw(nzt+1) = dz
[1353]	363	zw(nzt) = 0.0_wp
[94]	364	DO k = 0, nzt
[1353]	365	zw(k) = ( zu(k) + zu(k+1) ) * 0.5_wp
[94]	366	ENDDO
	367
[843]	368	!
	369	!-- In case of dirichlet bc for u and v the first u- and w-level are defined
	370	!-- at same height.
	371	IF ( ibc_uv_b == 0 ) THEN
	372	zu(0) = zw(0)
	373	ENDIF
	374
[94]	375	ENDIF
	376
	377	!
[1]	378	!-- Compute grid lengths.
	379	DO k = 1, nzt+1
	380	dzu(k) = zu(k) - zu(k-1)
[1353]	381	ddzu(k) = 1.0_wp / dzu(k)
[1]	382	dzw(k) = zw(k) - zw(k-1)
[1353]	383	ddzw(k) = 1.0_wp / dzw(k)
[1]	384	ENDDO
	385
	386	DO k = 1, nzt
[1353]	387	dd2zu(k) = 1.0_wp / ( dzu(k) + dzu(k+1) )
[1]	388	ENDDO
[667]	389
	390	!
[709]	391	!-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
	392	!-- everywhere. For the actual grid, the grid spacing at the lowest level
	393	!-- is only dz/2, but should be dz. Therefore, an additional array
	394	!-- containing with appropriate grid information is created for these
	395	!-- solvers.
[1575]	396	IF ( psolver(1:9) /= 'multigrid' ) THEN
[667]	397	ALLOCATE( ddzu_pres(1:nzt+1) )
	398	ddzu_pres = ddzu
[864]	399	ddzu_pres(1) = ddzu_pres(2) ! change for lowest level
[667]	400	ENDIF
[1]	401
	402	!
	403	!-- In case of multigrid method, compute grid lengths and grid factors for the
	404	!-- grid levels
[1575]	405	IF ( psolver(1:9) == 'multigrid' ) THEN
[1]	406
	407	ALLOCATE( ddx2_mg(maximum_grid_level), ddy2_mg(maximum_grid_level), &
	408	dzu_mg(nzb+1:nzt+1,maximum_grid_level), &
	409	dzw_mg(nzb+1:nzt+1,maximum_grid_level), &
	410	f1_mg(nzb+1:nzt,maximum_grid_level), &
	411	f2_mg(nzb+1:nzt,maximum_grid_level), &
	412	f3_mg(nzb+1:nzt,maximum_grid_level) )
	413
	414	dzu_mg(:,maximum_grid_level) = dzu
[667]	415	!
[864]	416	!-- Next line to ensure an equally spaced grid.
	417	dzu_mg(1,maximum_grid_level) = dzu(2)
[709]	418
[1]	419	dzw_mg(:,maximum_grid_level) = dzw
	420	nzt_l = nzt
	421	DO l = maximum_grid_level-1, 1, -1
[1353]	422	dzu_mg(nzb+1,l) = 2.0_wp * dzu_mg(nzb+1,l+1)
	423	dzw_mg(nzb+1,l) = 2.0_wp * dzw_mg(nzb+1,l+1)
[1]	424	nzt_l = nzt_l / 2
	425	DO k = 2, nzt_l+1
	426	dzu_mg(k,l) = dzu_mg(2k-2,l+1) + dzu_mg(2k-1,l+1)
	427	dzw_mg(k,l) = dzw_mg(2k-2,l+1) + dzw_mg(2k-1,l+1)
	428	ENDDO
	429	ENDDO
	430
	431	nzt_l = nzt
	432	dx_l = dx
	433	dy_l = dy
	434	DO l = maximum_grid_level, 1, -1
[1353]	435	ddx2_mg(l) = 1.0_wp / dx_l**2
	436	ddy2_mg(l) = 1.0_wp / dy_l**2
[1]	437	DO k = nzb+1, nzt_l
[1353]	438	f2_mg(k,l) = 1.0_wp / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
	439	f3_mg(k,l) = 1.0_wp / ( dzu_mg(k,l) * dzw_mg(k,l) )
	440	f1_mg(k,l) = 2.0_wp * ( ddx2_mg(l) + ddy2_mg(l) ) + &
[1]	441	f2_mg(k,l) + f3_mg(k,l)
	442	ENDDO
	443	nzt_l = nzt_l / 2
[1353]	444	dx_l = dx_l * 2.0_wp
	445	dy_l = dy_l * 2.0_wp
[1]	446	ENDDO
	447
	448	ENDIF
	449
	450	!
	451	!-- Compute the reciprocal values of the horizontal grid lengths.
[1353]	452	ddx = 1.0_wp / dx
	453	ddy = 1.0_wp / dy
[1]	454	dx2 = dx * dx
	455	dy2 = dy * dy
[1353]	456	ddx2 = 1.0_wp / dx2
	457	ddy2 = 1.0_wp / dy2
[1]	458
	459	!
	460	!-- Compute the grid-dependent mixing length.
	461	DO k = 1, nzt
[1322]	462	l_grid(k) = ( dx * dy * dzw(k) )**0.33333333333333_wp
[1]	463	ENDDO
	464
	465	!
	466	!-- Allocate outer and inner index arrays for topography and set
[114]	467	!-- defaults.
	468	!-- nzb_local has to contain additional layers of ghost points for calculating
	469	!-- the flag arrays needed for the multigrid method
	470	gls = 2**( maximum_grid_level )
[861]	471	IF ( gls < nbgp ) gls = nbgp
[667]	472
[114]	473	ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), &
	474	corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), &
[667]	475	nzb_local(-gls:ny+gls,-gls:nx+gls), &
	476	nzb_tmp(-nbgp:ny+nbgp,-nbgp:nx+nbgp), &
[114]	477	wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr), &
[1]	478	wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
[667]	479	ALLOCATE( fwxm(nysg:nyng,nxlg:nxrg), fwxp(nysg:nyng,nxlg:nxrg), &
	480	fwym(nysg:nyng,nxlg:nxrg), fwyp(nysg:nyng,nxlg:nxrg), &
	481	fxm(nysg:nyng,nxlg:nxrg), fxp(nysg:nyng,nxlg:nxrg), &
	482	fym(nysg:nyng,nxlg:nxrg), fyp(nysg:nyng,nxlg:nxrg), &
	483	nzb_s_inner(nysg:nyng,nxlg:nxrg), &
	484	nzb_s_outer(nysg:nyng,nxlg:nxrg), &
	485	nzb_u_inner(nysg:nyng,nxlg:nxrg), &
	486	nzb_u_outer(nysg:nyng,nxlg:nxrg), &
	487	nzb_v_inner(nysg:nyng,nxlg:nxrg), &
	488	nzb_v_outer(nysg:nyng,nxlg:nxrg), &
	489	nzb_w_inner(nysg:nyng,nxlg:nxrg), &
	490	nzb_w_outer(nysg:nyng,nxlg:nxrg), &
	491	nzb_diff_s_inner(nysg:nyng,nxlg:nxrg), &
	492	nzb_diff_s_outer(nysg:nyng,nxlg:nxrg), &
	493	nzb_diff_u(nysg:nyng,nxlg:nxrg), &
	494	nzb_diff_v(nysg:nyng,nxlg:nxrg), &
[1221]	495	rflags_s_inner(nzb:nzt+2,nysg:nyng,nxlg:nxrg), &
	496	rflags_invers(nysg:nyng,nxlg:nxrg,nzb:nzt+2), &
[667]	497	wall_e_x(nysg:nyng,nxlg:nxrg), &
	498	wall_e_y(nysg:nyng,nxlg:nxrg), &
	499	wall_u(nysg:nyng,nxlg:nxrg), &
	500	wall_v(nysg:nyng,nxlg:nxrg), &
	501	wall_w_x(nysg:nyng,nxlg:nxrg), &
	502	wall_w_y(nysg:nyng,nxlg:nxrg) )
[1]	503
	504
[667]	505
	506	ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	507
[818]	508
[1]	509	nzb_s_inner = nzb; nzb_s_outer = nzb
	510	nzb_u_inner = nzb; nzb_u_outer = nzb
	511	nzb_v_inner = nzb; nzb_v_outer = nzb
	512	nzb_w_inner = nzb; nzb_w_outer = nzb
	513
[1353]	514	rflags_s_inner = 1.0_wp
	515	rflags_invers = 1.0_wp
[1221]	516
[1]	517	!
[19]	518	!-- Define vertical gridpoint from (or to) which on the usual finite difference
[1]	519	!-- form (which does not use surface fluxes) is applied
[1691]	520	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
[1]	521	nzb_diff = nzb + 2
	522	ELSE
	523	nzb_diff = nzb + 1
	524	ENDIF
[19]	525	IF ( use_top_fluxes ) THEN
	526	nzt_diff = nzt - 1
	527	ELSE
	528	nzt_diff = nzt
	529	ENDIF
[1]	530
	531	nzb_diff_s_inner = nzb_diff; nzb_diff_s_outer = nzb_diff
	532	nzb_diff_u = nzb_diff; nzb_diff_v = nzb_diff
	533
[1353]	534	wall_e_x = 0.0_wp; wall_e_y = 0.0_wp; wall_u = 0.0_wp; wall_v = 0.0_wp
	535	wall_w_x = 0.0_wp; wall_w_y = 0.0_wp
	536	fwxp = 1.0_wp; fwxm = 1.0_wp; fwyp = 1.0_wp; fwym = 1.0_wp
	537	fxp = 1.0_wp; fxm = 1.0_wp; fyp = 1.0_wp; fym = 1.0_wp
[1]	538
	539	!
	540	!-- Initialize near-wall mixing length l_wall only in the vertical direction
	541	!-- for the moment,
	542	!-- multiplication with wall_adjustment_factor near the end of this routine
	543	l_wall(nzb,:,:) = l_grid(1)
	544	DO k = nzb+1, nzt
	545	l_wall(k,:,:) = l_grid(k)
	546	ENDDO
	547	l_wall(nzt+1,:,:) = l_grid(nzt)
	548
	549	ALLOCATE ( vertical_influence(nzb:nzt) )
	550	DO k = 1, nzt
	551	vertical_influence(k) = MIN ( INT( l_grid(k) / &
[1353]	552	( wall_adjustment_factor * dzw(k) ) + 0.5_wp ), nzt - k )
[1]	553	ENDDO
	554
	555	DO k = 1, MAXVAL( nzb_s_inner )
[1353]	556	IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR. &
	557	l_grid(k) > 1.5_wp * dy * wall_adjustment_factor ) THEN
[254]	558	WRITE( message_string, * ) 'grid anisotropy exceeds ', &
	559	'threshold given by only local', &
	560	' &horizontal reduction of near_wall ', &
	561	'mixing length l_wall', &
	562	' &starting from height level k = ', k, '.'
	563	CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
[1]	564	EXIT
	565	ENDIF
	566	ENDDO
	567	vertical_influence(0) = vertical_influence(1)
	568
[667]	569	DO i = nxlg, nxrg
	570	DO j = nysg, nyng
[1]	571	DO k = nzb_s_inner(j,i) + 1, &
	572	nzb_s_inner(j,i) + vertical_influence(nzb_s_inner(j,i))
	573	l_wall(k,j,i) = zu(k) - zw(nzb_s_inner(j,i))
	574	ENDDO
	575	ENDDO
	576	ENDDO
	577
	578	!
	579	!-- Set outer and inner index arrays for non-flat topography.
	580	!-- Here consistency checks concerning domain size and periodicity are
	581	!-- necessary.
	582	!-- Within this SELECT CASE structure only nzb_local is initialized
	583	!-- individually depending on the chosen topography type, all other index
	584	!-- arrays are initialized further below.
	585	SELECT CASE ( TRIM( topography ) )
	586
	587	CASE ( 'flat' )
	588	!
[555]	589	!-- nzb_local is required for the multigrid solver
	590	nzb_local = 0
[1]	591
	592	CASE ( 'single_building' )
	593	!
	594	!-- Single rectangular building, by default centered in the middle of the
	595	!-- total domain
	596	blx = NINT( building_length_x / dx )
	597	bly = NINT( building_length_y / dy )
[1675]	598	bh = MINLOC( ABS( zw - building_height ), 1 ) - 1
	599	IF ( ABS( zw(bh ) - building_height ) == &
	600	ABS( zw(bh+1) - building_height ) ) bh = bh + 1
[1]	601
[1322]	602	IF ( building_wall_left == 9999999.9_wp ) THEN
[1]	603	building_wall_left = ( nx + 1 - blx ) / 2 * dx
	604	ENDIF
	605	bxl = NINT( building_wall_left / dx )
	606	bxr = bxl + blx
	607
[1322]	608	IF ( building_wall_south == 9999999.9_wp ) THEN
[1]	609	building_wall_south = ( ny + 1 - bly ) / 2 * dy
	610	ENDIF
	611	bys = NINT( building_wall_south / dy )
	612	byn = bys + bly
	613
	614	!
	615	!-- Building size has to meet some requirements
	616	IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR. &
	617	( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) ) THEN
[274]	618	WRITE( message_string, * ) 'inconsistent building parameters:', &
	619	'& bxl=', bxl, 'bxr=', bxr, 'bys=', bys, &
	620	'byn=', byn, 'nx=', nx, 'ny=', ny
[254]	621	CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
[1]	622	ENDIF
	623
	624	!
[217]	625	!-- Define the building.
[1]	626	nzb_local = 0
[134]	627	nzb_local(bys:byn,bxl:bxr) = bh
[1]	628
[240]	629	CASE ( 'single_street_canyon' )
	630	!
	631	!-- Single quasi-2D street canyon of infinite length in x or y direction.
	632	!-- The canyon is centered in the other direction by default.
[1322]	633	IF ( canyon_width_x /= 9999999.9_wp ) THEN
[240]	634	!
	635	!-- Street canyon in y direction
	636	cwx = NINT( canyon_width_x / dx )
[1322]	637	IF ( canyon_wall_left == 9999999.9_wp ) THEN
[240]	638	canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
	639	ENDIF
	640	cxl = NINT( canyon_wall_left / dx )
	641	cxr = cxl + cwx
	642
[1322]	643	ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN
[240]	644	!
	645	!-- Street canyon in x direction
	646	cwy = NINT( canyon_width_y / dy )
[1322]	647	IF ( canyon_wall_south == 9999999.9_wp ) THEN
[240]	648	canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
	649	ENDIF
	650	cys = NINT( canyon_wall_south / dy )
	651	cyn = cys + cwy
	652
	653	ELSE
[254]	654
	655	message_string = 'no street canyon width given'
	656	CALL message( 'init_grid', 'PA0204', 1, 2, 0, 6, 0 )
	657
[240]	658	ENDIF
	659
[1675]	660	ch = MINLOC( ABS( zw - canyon_height ), 1 ) - 1
	661	IF ( ABS( zw(ch ) - canyon_height ) == &
	662	ABS( zw(ch+1) - canyon_height ) ) ch = ch + 1
	663
[240]	664	dp_level_ind_b = ch
	665	!
	666	!-- Street canyon size has to meet some requirements
[1322]	667	IF ( canyon_width_x /= 9999999.9_wp ) THEN
[1353]	668	IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR. &
[240]	669	( ch < 3 ) ) THEN
[1353]	670	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
	671	'&cxl=', cxl, 'cxr=', cxr, &
	672	'cwx=', cwx, &
[254]	673	'ch=', ch, 'nx=', nx, 'ny=', ny
	674	CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 )
[240]	675	ENDIF
[1322]	676	ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN
[1353]	677	IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR. &
[240]	678	( ch < 3 ) ) THEN
[1353]	679	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
	680	'&cys=', cys, 'cyn=', cyn, &
	681	'cwy=', cwy, &
[254]	682	'ch=', ch, 'nx=', nx, 'ny=', ny
	683	CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 )
[240]	684	ENDIF
	685	ENDIF
[1353]	686	IF ( canyon_width_x /= 9999999.9_wp .AND. &
	687	canyon_width_y /= 9999999.9_wp ) THEN
	688	message_string = 'inconsistent canyon parameters:' // &
	689	'&street canyon can only be oriented' // &
[254]	690	'&either in x- or in y-direction'
	691	CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 )
[240]	692	ENDIF
	693
	694	nzb_local = ch
[1322]	695	IF ( canyon_width_x /= 9999999.9_wp ) THEN
[240]	696	nzb_local(:,cxl+1:cxr-1) = 0
[1322]	697	ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN
[240]	698	nzb_local(cys+1:cyn-1,:) = 0
	699	ENDIF
	700
[1]	701	CASE ( 'read_from_file' )
[759]	702
[818]	703	ALLOCATE ( topo_height(0:ny,0:nx) )
	704
[759]	705	DO ii = 0, io_blocks-1
	706	IF ( ii == io_group ) THEN
	707
[1]	708	!
[759]	709	!-- Arbitrary irregular topography data in PALM format (exactly
	710	!-- matching the grid size and total domain size)
[1779]	711	OPEN( 90, FILE='TOPOGRAPHY_DATA'//TRIM( coupling_char ), &
	712	STATUS='OLD', FORM='FORMATTED', ERR=10 )
[759]	713	DO j = ny, 0, -1
	714	READ( 90, *, ERR=11, END=11 ) ( topo_height(j,i), i = 0,nx )
	715	ENDDO
	716
	717	GOTO 12
	718
[1779]	719	10 message_string = 'file TOPOGRAPHY'//TRIM( coupling_char )// &
	720	' does not exist'
[759]	721	CALL message( 'init_grid', 'PA0208', 1, 2, 0, 6, 0 )
	722
[1779]	723	11 message_string = 'errors in file TOPOGRAPHY_DATA'// &
	724	TRIM( coupling_char )
[759]	725	CALL message( 'init_grid', 'PA0209', 1, 2, 0, 6, 0 )
	726
	727	12 CLOSE( 90 )
	728
	729	ENDIF
[1804]	730	#if defined( __parallel )
[759]	731	CALL MPI_BARRIER( comm2d, ierr )
	732	#endif
[559]	733	ENDDO
[759]	734
[1]	735	!
	736	!-- Calculate the index height of the topography
	737	DO i = 0, nx
	738	DO j = 0, ny
[1675]	739	nzb_local(j,i) = MINLOC( ABS( zw - topo_height(j,i) ), 1 ) - 1
	740	IF ( ABS( zw(nzb_local(j,i) ) - topo_height(j,i) ) == &
	741	ABS( zw(nzb_local(j,i)+1) - topo_height(j,i) ) ) &
	742	nzb_local(j,i) = nzb_local(j,i) + 1
[1]	743	ENDDO
	744	ENDDO
[818]	745
	746	DEALLOCATE ( topo_height )
[114]	747	!
[759]	748	!-- Add cyclic boundaries (additional layers are for calculating
	749	!-- flag arrays needed for the multigrid sover)
[114]	750	nzb_local(-gls:-1,0:nx) = nzb_local(ny-gls+1:ny,0:nx)
	751	nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,0:nx)
	752	nzb_local(:,-gls:-1) = nzb_local(:,nx-gls+1:nx)
	753	nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
[667]	754
[1]	755	CASE DEFAULT
	756	!
	757	!-- The DEFAULT case is reached either if the parameter topography
[217]	758	!-- contains a wrong character string or if the user has defined a special
[1]	759	!-- case in the user interface. There, the subroutine user_init_grid
	760	!-- checks which of these two conditions applies.
[114]	761	CALL user_init_grid( gls, nzb_local )
[1]	762
	763	END SELECT
	764	!
[861]	765	!-- Determine the maximum level of topography. Furthermore it is used for
	766	!-- steering the degradation of order of the applied advection scheme.
[978]	767	!-- In case of non-cyclic lateral boundaries, the order of the advection
[996]	768	!-- scheme have to be reduced up to nzt (required at the lateral boundaries).
[1677]	769	nzb_max = MAXVAL( nzb_local ) + 1
[1353]	770	IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR. &
[1762]	771	inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s .OR. &
	772	nest_domain ) &
	773	THEN
	774	nzb_max = nzt
[978]	775	ENDIF
	776
[861]	777	!
[1]	778	!-- Consistency checks and index array initialization are only required for
[217]	779	!-- non-flat topography, also the initialization of topography height arrays
[49]	780	!-- zu_s_inner and zw_w_inner
[1]	781	IF ( TRIM( topography ) /= 'flat' ) THEN
	782
	783	!
	784	!-- Consistency checks
	785	IF ( MINVAL( nzb_local ) < 0 .OR. MAXVAL( nzb_local ) > nz + 1 ) THEN
[1353]	786	WRITE( message_string, * ) 'nzb_local values are outside the', &
	787	'model domain', &
	788	'&MINVAL( nzb_local ) = ', MINVAL(nzb_local), &
[274]	789	'&MAXVAL( nzb_local ) = ', MAXVAL(nzb_local)
[254]	790	CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
[1]	791	ENDIF
	792
[722]	793	IF ( bc_lr == 'cyclic' ) THEN
[1353]	794	IF ( ANY( nzb_local(:,-1) /= nzb_local(:,nx) ) .OR. &
[1]	795	ANY( nzb_local(:,0) /= nzb_local(:,nx+1) ) ) THEN
[1353]	796	message_string = 'nzb_local does not fulfill cyclic' // &
[254]	797	' boundary condition in x-direction'
	798	CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
[1]	799	ENDIF
	800	ENDIF
[722]	801	IF ( bc_ns == 'cyclic' ) THEN
[1353]	802	IF ( ANY( nzb_local(-1,:) /= nzb_local(ny,:) ) .OR. &
[1]	803	ANY( nzb_local(0,:) /= nzb_local(ny+1,:) ) ) THEN
[1353]	804	message_string = 'nzb_local does not fulfill cyclic' // &
[254]	805	' boundary condition in y-direction'
	806	CALL message( 'init_grid', 'PA0212', 1, 2, 0, 6, 0 )
[1]	807	ENDIF
	808	ENDIF
	809
[217]	810	IF ( topography_grid_convention == 'cell_edge' ) THEN
[134]	811	!
[217]	812	!-- The array nzb_local as defined using the 'cell_edge' convention
	813	!-- describes the actual total size of topography which is defined at the
	814	!-- cell edges where u=0 on the topography walls in x-direction and v=0
	815	!-- on the topography walls in y-direction. However, PALM uses individual
	816	!-- arrays nzb_u\|v\|w\|s_inner\|outer that are based on nzb_s_inner.
	817	!-- Therefore, the extent of topography in nzb_local is now reduced by
	818	!-- 1dx at the E topography walls and by 1dy at the N topography walls
	819	!-- to form the basis for nzb_s_inner.
	820	DO j = -gls, ny + gls
	821	DO i = -gls, nx
	822	nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j,i+1) )
	823	ENDDO
[134]	824	ENDDO
[217]	825	!-- apply cyclic boundary conditions in x-direction
	826	!(ist das erforderlich? Ursache von Seung Bus Fehler?)
	827	nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
	828	DO i = -gls, nx + gls
	829	DO j = -gls, ny
	830	nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j+1,i) )
	831	ENDDO
[134]	832	ENDDO
[217]	833	!-- apply cyclic boundary conditions in y-direction
	834	!(ist das erforderlich? Ursache von Seung Bus Fehler?)
	835	nzb_local(ny+1:ny+gls,:) = nzb_local(0:gls-1,:)
	836	ENDIF
[134]	837
[1]	838	!
	839	!-- Initialize index arrays nzb_s_inner and nzb_w_inner
[861]	840	nzb_s_inner = nzb_local(nysg:nyng,nxlg:nxrg)
	841	nzb_w_inner = nzb_local(nysg:nyng,nxlg:nxrg)
[1]	842
	843	!
	844	!-- Initialize remaining index arrays:
	845	!-- first pre-initialize them with nzb_s_inner...
	846	nzb_u_inner = nzb_s_inner
	847	nzb_u_outer = nzb_s_inner
	848	nzb_v_inner = nzb_s_inner
	849	nzb_v_outer = nzb_s_inner
	850	nzb_w_outer = nzb_s_inner
	851	nzb_s_outer = nzb_s_inner
	852
	853	!
	854	!-- ...then extend pre-initialized arrays in their according directions
	855	!-- based on nzb_local using nzb_tmp as a temporary global index array
	856
	857	!
	858	!-- nzb_s_outer:
	859	!-- extend nzb_local east-/westwards first, then north-/southwards
[667]	860	nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
[1]	861	DO j = -1, ny + 1
	862	DO i = 0, nx
[1353]	863	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i), &
[1]	864	nzb_local(j,i+1) )
	865	ENDDO
	866	ENDDO
	867	DO i = nxl, nxr
	868	DO j = nys, nyn
[1353]	869	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
[1]	870	nzb_tmp(j+1,i) )
	871	ENDDO
	872	!
	873	!-- non-cyclic boundary conditions (overwritten by call of
	874	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
	875	IF ( nys == 0 ) THEN
	876	j = -1
	877	nzb_s_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
	878	ENDIF
[1743]	879	IF ( nyn == ny ) THEN
[1]	880	j = ny + 1
	881	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
	882	ENDIF
	883	ENDDO
	884	!
	885	!-- nzb_w_outer:
	886	!-- identical to nzb_s_outer
	887	nzb_w_outer = nzb_s_outer
	888
	889	!
	890	!-- nzb_u_inner:
	891	!-- extend nzb_local rightwards only
[667]	892	nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
[1]	893	DO j = -1, ny + 1
	894	DO i = 0, nx + 1
	895	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )
	896	ENDDO
	897	ENDDO
[667]	898	nzb_u_inner = nzb_tmp(nysg:nyng,nxlg:nxrg)
[1]	899
	900	!
	901	!-- nzb_u_outer:
	902	!-- extend current nzb_tmp (nzb_u_inner) north-/southwards
	903	DO i = nxl, nxr
	904	DO j = nys, nyn
[1353]	905	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
[1]	906	nzb_tmp(j+1,i) )
	907	ENDDO
	908	!
	909	!-- non-cyclic boundary conditions (overwritten by call of
	910	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
	911	IF ( nys == 0 ) THEN
	912	j = -1
	913	nzb_u_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
	914	ENDIF
[1743]	915	IF ( nyn == ny ) THEN
[1]	916	j = ny + 1
	917	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
	918	ENDIF
	919	ENDDO
	920
	921	!
	922	!-- nzb_v_inner:
	923	!-- extend nzb_local northwards only
[667]	924	nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
[1]	925	DO i = -1, nx + 1
	926	DO j = 0, ny + 1
	927	nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )
	928	ENDDO
	929	ENDDO
[1886]	930	nzb_v_inner = nzb_tmp(nysg:nyng,nxlg:nxrg)
[1]	931
	932	!
	933	!-- nzb_v_outer:
	934	!-- extend current nzb_tmp (nzb_v_inner) right-/leftwards
	935	DO j = nys, nyn
	936	DO i = nxl, nxr
[1353]	937	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i), &
[1]	938	nzb_tmp(j,i+1) )
	939	ENDDO
	940	!
	941	!-- non-cyclic boundary conditions (overwritten by call of
	942	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
	943	IF ( nxl == 0 ) THEN
	944	i = -1
	945	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i+1), nzb_tmp(j,i) )
	946	ENDIF
	947	IF ( nxr == nx ) THEN
	948	i = nx + 1
	949	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i) )
	950	ENDIF
	951	ENDDO
[1804]	952
[1]	953	!
	954	!-- Exchange of lateral boundary values (parallel computers) and cyclic
	955	!-- boundary conditions, if applicable.
	956	!-- Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local
	957	!-- they do not require exchange and are not included here.
	958	CALL exchange_horiz_2d_int( nzb_u_inner )
	959	CALL exchange_horiz_2d_int( nzb_u_outer )
	960	CALL exchange_horiz_2d_int( nzb_v_inner )
	961	CALL exchange_horiz_2d_int( nzb_v_outer )
	962	CALL exchange_horiz_2d_int( nzb_w_outer )
	963	CALL exchange_horiz_2d_int( nzb_s_outer )
	964
[49]	965	!
	966	!-- Allocate and set the arrays containing the topography height
	967	IF ( myid == 0 ) THEN
	968
	969	ALLOCATE( zu_s_inner(0:nx+1,0:ny+1), zw_w_inner(0:nx+1,0:ny+1) )
	970
	971	DO i = 0, nx + 1
	972	DO j = 0, ny + 1
	973	zu_s_inner(i,j) = zu(nzb_local(j,i))
	974	zw_w_inner(i,j) = zw(nzb_local(j,i))
	975	ENDDO
	976	ENDDO
	977
	978	ENDIF
[1221]	979	!
	980	!-- Set flag arrays to be used for masking of grid points
	981	DO i = nxlg, nxrg
	982	DO j = nysg, nyng
	983	DO k = nzb, nzt+1
[1353]	984	IF ( k <= nzb_s_inner(j,i) ) rflags_s_inner(k,j,i) = 0.0_wp
	985	IF ( k <= nzb_s_inner(j,i) ) rflags_invers(j,i,k) = 0.0_wp
[1221]	986	ENDDO
	987	ENDDO
	988	ENDDO
[1804]	989
[1]	990	ENDIF
	991
	992	!
	993	!-- Set the individual index arrays which define the k index from which on
	994	!-- the usual finite difference form (which does not use surface fluxes) is
	995	!-- applied
[1691]	996	IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN
[1]	997	nzb_diff_u = nzb_u_inner + 2
	998	nzb_diff_v = nzb_v_inner + 2
	999	nzb_diff_s_inner = nzb_s_inner + 2
	1000	nzb_diff_s_outer = nzb_s_outer + 2
	1001	ELSE
	1002	nzb_diff_u = nzb_u_inner + 1
	1003	nzb_diff_v = nzb_v_inner + 1
	1004	nzb_diff_s_inner = nzb_s_inner + 1
	1005	nzb_diff_s_outer = nzb_s_outer + 1
	1006	ENDIF
	1007
	1008	!
	1009	!-- Calculation of wall switches and factors required by diffusion_u/v.f90 and
	1010	!-- for limitation of near-wall mixing length l_wall further below
	1011	corner_nl = 0
	1012	corner_nr = 0
	1013	corner_sl = 0
	1014	corner_sr = 0
	1015	wall_l = 0
	1016	wall_n = 0
	1017	wall_r = 0
	1018	wall_s = 0
	1019
	1020	DO i = nxl, nxr
	1021	DO j = nys, nyn
	1022	!
	1023	!-- u-component
	1024	IF ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) THEN
[1353]	1025	wall_u(j,i) = 1.0_wp ! north wall (location of adjacent fluid)
	1026	fym(j,i) = 0.0_wp
	1027	fyp(j,i) = 1.0_wp
[1]	1028	ELSEIF ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) THEN
[1353]	1029	wall_u(j,i) = 1.0_wp ! south wall (location of adjacent fluid)
	1030	fym(j,i) = 1.0_wp
	1031	fyp(j,i) = 0.0_wp
[1]	1032	ENDIF
	1033	!
	1034	!-- v-component
	1035	IF ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) THEN
[1353]	1036	wall_v(j,i) = 1.0_wp ! rigth wall (location of adjacent fluid)
	1037	fxm(j,i) = 0.0_wp
	1038	fxp(j,i) = 1.0_wp
[1]	1039	ELSEIF ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) ) THEN
[1353]	1040	wall_v(j,i) = 1.0_wp ! left wall (location of adjacent fluid)
	1041	fxm(j,i) = 1.0_wp
	1042	fxp(j,i) = 0.0_wp
[1]	1043	ENDIF
	1044	!
	1045	!-- w-component, also used for scalars, separate arrays for shear
	1046	!-- production of tke
	1047	IF ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) THEN
[1353]	1048	wall_e_y(j,i) = 1.0_wp ! north wall (location of adjacent fluid)
	1049	wall_w_y(j,i) = 1.0_wp
	1050	fwym(j,i) = 0.0_wp
	1051	fwyp(j,i) = 1.0_wp
[1]	1052	ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) THEN
[1353]	1053	wall_e_y(j,i) = -1.0_wp ! south wall (location of adjacent fluid)
	1054	wall_w_y(j,i) = 1.0_wp
	1055	fwym(j,i) = 1.0_wp
	1056	fwyp(j,i) = 0.0_wp
[1]	1057	ENDIF
	1058	IF ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) THEN
[1353]	1059	wall_e_x(j,i) = 1.0_wp ! right wall (location of adjacent fluid)
	1060	wall_w_x(j,i) = 1.0_wp
	1061	fwxm(j,i) = 0.0_wp
	1062	fwxp(j,i) = 1.0_wp
[1]	1063	ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) THEN
[1353]	1064	wall_e_x(j,i) = -1.0_wp ! left wall (location of adjacent fluid)
	1065	wall_w_x(j,i) = 1.0_wp
	1066	fwxm(j,i) = 1.0_wp
	1067	fwxp(j,i) = 0.0_wp
[1]	1068	ENDIF
	1069	!
	1070	!-- Wall and corner locations inside buildings for limitation of
	1071	!-- near-wall mixing length l_wall
	1072	IF ( nzb_s_inner(j,i) > nzb_s_inner(j+1,i) ) THEN
	1073
	1074	wall_n(j,i) = nzb_s_inner(j+1,i) + 1 ! North wall
	1075
	1076	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
	1077	corner_nl(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northleft corner
	1078	nzb_s_inner(j,i-1) ) + 1
	1079	ENDIF
	1080
	1081	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
	1082	corner_nr(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northright corner
	1083	nzb_s_inner(j,i+1) ) + 1
	1084	ENDIF
	1085
	1086	ENDIF
	1087
	1088	IF ( nzb_s_inner(j,i) > nzb_s_inner(j-1,i) ) THEN
	1089
	1090	wall_s(j,i) = nzb_s_inner(j-1,i) + 1 ! South wall
	1091	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
	1092	corner_sl(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southleft corner
	1093	nzb_s_inner(j,i-1) ) + 1
	1094	ENDIF
	1095
	1096	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
	1097	corner_sr(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southright corner
	1098	nzb_s_inner(j,i+1) ) + 1
	1099	ENDIF
	1100
	1101	ENDIF
	1102
	1103	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
	1104	wall_l(j,i) = nzb_s_inner(j,i-1) + 1 ! Left wall
	1105	ENDIF
	1106
	1107	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
	1108	wall_r(j,i) = nzb_s_inner(j,i+1) + 1 ! Right wall
	1109	ENDIF
	1110
	1111	ENDDO
	1112	ENDDO
	1113
	1114	!
[114]	1115	!-- Calculate wall flag arrays for the multigrid method
[1575]	1116	IF ( psolver(1:9) == 'multigrid' ) THEN
[114]	1117	!
	1118	!-- Gridpoint increment of the current level
	1119	inc = 1
	1120
	1121	DO l = maximum_grid_level, 1 , -1
	1122
	1123	nxl_l = nxl_mg(l)
	1124	nxr_l = nxr_mg(l)
	1125	nys_l = nys_mg(l)
	1126	nyn_l = nyn_mg(l)
	1127	nzt_l = nzt_mg(l)
	1128
	1129	!
	1130	!-- Assign the flag level to be calculated
	1131	SELECT CASE ( l )
	1132	CASE ( 1 )
	1133	flags => wall_flags_1
	1134	CASE ( 2 )
	1135	flags => wall_flags_2
	1136	CASE ( 3 )
	1137	flags => wall_flags_3
	1138	CASE ( 4 )
	1139	flags => wall_flags_4
	1140	CASE ( 5 )
	1141	flags => wall_flags_5
	1142	CASE ( 6 )
	1143	flags => wall_flags_6
	1144	CASE ( 7 )
	1145	flags => wall_flags_7
	1146	CASE ( 8 )
	1147	flags => wall_flags_8
	1148	CASE ( 9 )
	1149	flags => wall_flags_9
	1150	CASE ( 10 )
	1151	flags => wall_flags_10
	1152	END SELECT
	1153
	1154	!
	1155	!-- Depending on the grid level, set the respective bits in case of
	1156	!-- neighbouring walls
	1157	!-- Bit 0: wall to the bottom
	1158	!-- Bit 1: wall to the top (not realized in remaining PALM code so far)
	1159	!-- Bit 2: wall to the south
	1160	!-- Bit 3: wall to the north
	1161	!-- Bit 4: wall to the left
	1162	!-- Bit 5: wall to the right
[116]	1163	!-- Bit 6: inside building
[114]	1164
	1165	flags = 0
	1166
[927]	1167	!
	1168	!-- In case of masking method, flags are not set and multigrid method
	1169	!-- works like FFT-solver
[1902]	1170	IF ( psolver == 'multigrid' .AND. .NOT. masking_method ) THEN
[927]	1171
	1172	DO i = nxl_l-1, nxr_l+1
	1173	DO j = nys_l-1, nyn_l+1
	1174	DO k = nzb, nzt_l+1
[114]	1175
	1176	!
[927]	1177	!-- Inside/outside building (inside building does not need
	1178	!-- further tests for walls)
	1179	IF ( kinc <= nzb_local(jinc,i*inc) ) THEN
[114]	1180
[927]	1181	flags(k,j,i) = IBSET( flags(k,j,i), 6 )
[114]	1182
[927]	1183	ELSE
[114]	1184	!
[927]	1185	!-- Bottom wall
	1186	IF ( (k-1)inc <= nzb_local(jinc,i*inc) ) THEN
	1187	flags(k,j,i) = IBSET( flags(k,j,i), 0 )
	1188	ENDIF
[114]	1189	!
[927]	1190	!-- South wall
	1191	IF ( kinc <= nzb_local((j-1)inc,i*inc) ) THEN
	1192	flags(k,j,i) = IBSET( flags(k,j,i), 2 )
	1193	ENDIF
[114]	1194	!
[927]	1195	!-- North wall
	1196	IF ( kinc <= nzb_local((j+1)inc,i*inc) ) THEN
	1197	flags(k,j,i) = IBSET( flags(k,j,i), 3 )
	1198	ENDIF
[114]	1199	!
[927]	1200	!-- Left wall
	1201	IF ( kinc <= nzb_local(jinc,(i-1)*inc) ) THEN
	1202	flags(k,j,i) = IBSET( flags(k,j,i), 4 )
	1203	ENDIF
[114]	1204	!
[927]	1205	!-- Right wall
	1206	IF ( kinc <= nzb_local(jinc,(i+1)*inc) ) THEN
	1207	flags(k,j,i) = IBSET( flags(k,j,i), 5 )
	1208	ENDIF
	1209
[114]	1210	ENDIF
	1211
[927]	1212	ENDDO
[114]	1213	ENDDO
	1214	ENDDO
	1215
[927]	1216	ENDIF
	1217
[114]	1218	!
	1219	!-- Test output of flag arrays
[145]	1220	! i = nxl_l
	1221	! WRITE (9,*) ' '
	1222	! WRITE (9,) ' mg level ', l, ' *', mg_switch_to_pe0_level
	1223	! WRITE (9,*) ' inc=', inc, ' i =', nxl_l
	1224	! WRITE (9,*) ' nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
	1225	! DO k = nzt_l+1, nzb, -1
	1226	! WRITE (9,'(194(1X,I2))') ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
	1227	! ENDDO
[114]	1228
	1229	inc = inc * 2
	1230
	1231	ENDDO
	1232
	1233	ENDIF
[861]	1234	!
	1235	!-- Allocate flags needed for masking walls.
[1677]	1236	ALLOCATE( wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
	1237	wall_flags_00(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1221]	1238	wall_flags_0 = 0
	1239	wall_flags_00 = 0
[114]	1240
[1886]	1241	IF ( scalar_advec == 'ws-scheme' .OR. &
[1557]	1242	scalar_advec == 'ws-scheme-mono' ) THEN
[114]	1243	!
[861]	1244	!-- Set flags to steer the degradation of the advection scheme in advec_ws
	1245	!-- near topography, inflow- and outflow boundaries as well as bottom and
	1246	!-- top of model domain. wall_flags_0 remains zero for all non-prognostic
	1247	!-- grid points.
	1248	DO i = nxl, nxr
	1249	DO j = nys, nyn
	1250	DO k = nzb_s_inner(j,i)+1, nzt
	1251	!
	1252	!-- scalar - x-direction
	1253	!-- WS1 (0), WS3 (1), WS5 (2)
[1886]	1254	IF ( ( k <= nzb_s_inner(j,i+1) .OR. k <= nzb_s_inner(j,i+2) &
	1255	.OR. k <= nzb_s_inner(j,i-1) ) &
	1256	.OR. ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
[1762]	1257	.AND. i == nxl ) .OR. &
	1258	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
	1259	.AND. i == nxr ) ) &
	1260	THEN
[861]	1261	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
[1886]	1262	ELSEIF ( ( k <= nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i+1)&
	1263	.AND. k > nzb_s_inner(j,i+2)&
	1264	.AND. k > nzb_s_inner(j,i-1)&
	1265	) .OR. &
	1266	( k <= nzb_s_inner(j,i-2) .AND. k > nzb_s_inner(j,i+1)&
	1267	.AND. k > nzb_s_inner(j,i+2)&
	1268	.AND. k > nzb_s_inner(j,i-1)&
	1269	) &
[1762]	1270	.OR. &
	1271	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
	1272	.AND. i == nxr-1 ) .OR. &
	1273	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
	1274	.AND. i == nxlu ) ) &
	1275	THEN
[861]	1276	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
[1886]	1277	ELSEIF ( k > nzb_s_inner(j,i+1) .AND. k > nzb_s_inner(j,i+2) &
	1278	.AND. k > nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i-1) &
	1279	.AND. k > nzb_s_inner(j,i-2) ) &
	1280	THEN
[861]	1281	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 )
	1282	ENDIF
	1283	!
	1284	!-- scalar - y-direction
	1285	!-- WS1 (3), WS3 (4), WS5 (5)
[1886]	1286	IF ( ( k <= nzb_s_inner(j+1,i) .OR. k <= nzb_s_inner(j+2,i) &
	1287	.OR. k <= nzb_s_inner(j-1,i) ) &
	1288	.OR. &
[1762]	1289	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
	1290	.AND. j == nys ) .OR. &
	1291	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
	1292	.AND. j == nyn ) ) &
	1293	THEN
[861]	1294	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 )
[1762]	1295	!
[861]	1296	!-- WS3
[1886]	1297	ELSEIF ( ( k <= nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j+1,i)&
	1298	.AND. k > nzb_s_inner(j+2,i)&
	1299	.AND. k > nzb_s_inner(j-1,i)&
	1300	) .OR. &
	1301	( k <= nzb_s_inner(j-2,i) .AND. k > nzb_s_inner(j+1,i)&
	1302	.AND. k > nzb_s_inner(j+2,i)&
	1303	.AND. k > nzb_s_inner(j-1,i)&
	1304	) &
	1305	.OR. &
[1762]	1306	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
	1307	.AND. j == nysv ) .OR. &
	1308	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
	1309	.AND. j == nyn-1 ) ) &
	1310	THEN
[861]	1311	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 )
[1762]	1312	!
[861]	1313	!-- WS5
[1886]	1314	ELSEIF ( k > nzb_s_inner(j+1,i) .AND. k > nzb_s_inner(j+2,i) &
	1315	.AND. k > nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j-1,i) &
	1316	.AND. k > nzb_s_inner(j-2,i) ) &
	1317	THEN
[861]	1318	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 )
	1319	ENDIF
	1320	!
	1321	!-- scalar - z-direction
	1322	!-- WS1 (6), WS3 (7), WS5 (8)
	1323	flag_set = .FALSE.
	1324	IF ( k == nzb_s_inner(j,i) + 1 .OR. k == nzt ) THEN
	1325	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 )
	1326	flag_set = .TRUE.
	1327	ELSEIF ( k == nzb_s_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN
	1328	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 7 )
	1329	flag_set = .TRUE.
	1330	ELSEIF ( k > nzb_s_inner(j,i) .AND. .NOT. flag_set ) THEN
	1331	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 )
	1332	ENDIF
[1886]	1333
[861]	1334	ENDDO
	1335	ENDDO
	1336	ENDDO
	1337	ENDIF
	1338
	1339	IF ( momentum_advec == 'ws-scheme' ) THEN
	1340	!
	1341	!-- Set wall_flags_0 to steer the degradation of the advection scheme in advec_ws
	1342	!-- near topography, inflow- and outflow boundaries as well as bottom and
	1343	!-- top of model domain. wall_flags_0 remains zero for all non-prognostic
	1344	!-- grid points.
	1345	DO i = nxl, nxr
	1346	DO j = nys, nyn
[1580]	1347	DO k = nzb+1, nzt
[861]	1348	!
[1580]	1349	!-- At first, set flags to WS1.
	1350	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
	1351	!-- in order to handle the left/south flux.
	1352	!-- near vertical walls.
	1353	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
	1354	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
	1355	!
[861]	1356	!-- u component - x-direction
	1357	!-- WS1 (9), WS3 (10), WS5 (11)
[1886]	1358	IF ( k <= nzb_u_inner(j,i+1) .OR. &
[1762]	1359	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
	1360	.AND. i <= nxlu ) .OR. &
	1361	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
	1362	.AND. i == nxr ) ) &
	1363	THEN
[1580]	1364	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
[1886]	1365	ELSEIF ( ( k <= nzb_u_inner(j,i+2) .AND. &
	1366	k > nzb_u_inner(j,i+1) ) .OR. &
	1367	k <= nzb_u_inner(j,i-1) &
	1368	.OR. &
[1762]	1369	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
	1370	.AND. i == nxr-1 ) .OR. &
	1371	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
	1372	.AND. i == nxlu+1) ) &
	1373	THEN
[861]	1374	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 )
[1580]	1375	!
	1376	!-- Clear flag for WS1
	1377	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 )
[1886]	1378	ELSEIF ( k > nzb_u_inner(j,i+1) .AND. k > nzb_u_inner(j,i+2) &
	1379	.AND. k > nzb_u_inner(j,i-1) ) &
	1380	THEN
[861]	1381	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 )
[1580]	1382	!
	1383	!-- Clear flag for WS1
	1384	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 )
[861]	1385	ENDIF
	1386	!
	1387	!-- u component - y-direction
	1388	!-- WS1 (12), WS3 (13), WS5 (14)
[1762]	1389	IF ( k <= nzb_u_inner(j+1,i) .OR. &
	1390	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
	1391	.AND. j == nys ) .OR. &
	1392	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
	1393	.AND. j == nyn ) ) &
	1394	THEN
	1395	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
[1886]	1396	ELSEIF ( ( k <= nzb_u_inner(j+2,i) .AND. &
	1397	k > nzb_u_inner(j+1,i) ) .OR. &
	1398	k <= nzb_u_inner(j-1,i) &
	1399	.OR. &
[1762]	1400	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
	1401	.AND. j == nysv ) .OR. &
	1402	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
	1403	.AND. j == nyn-1 ) ) &
	1404	THEN
[861]	1405	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
[1580]	1406	!
	1407	!-- Clear flag for WS1
	1408	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 )
[1886]	1409	ELSEIF ( k > nzb_u_inner(j+1,i) .AND. k > nzb_u_inner(j+2,i) &
	1410	.AND. k > nzb_u_inner(j-1,i) ) &
	1411	THEN
[861]	1412	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 )
[1580]	1413	!
	1414	!-- Clear flag for WS1
	1415	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 )
[861]	1416	ENDIF
	1417	!
	1418	!-- u component - z-direction
	1419	!-- WS1 (15), WS3 (16), WS5 (17)
	1420	flag_set = .FALSE.
	1421	IF ( k == nzb_u_inner(j,i) + 1 .OR. k == nzt ) THEN
	1422	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 )
	1423	flag_set = .TRUE.
	1424	ELSEIF ( k == nzb_u_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN
	1425	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 )
	1426	flag_set = .TRUE.
	1427	ELSEIF ( k > nzb_u_inner(j,i) .AND. .NOT. flag_set ) THEN
	1428	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 )
	1429	ENDIF
	1430
	1431	ENDDO
	1432	ENDDO
	1433	ENDDO
	1434
	1435	DO i = nxl, nxr
	1436	DO j = nys, nyn
[1580]	1437	DO k = nzb+1, nzt
[861]	1438	!
[1580]	1439	!-- At first, set flags to WS1.
	1440	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
	1441	!-- in order to handle the left/south flux.
	1442	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
	1443	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
	1444	!
[861]	1445	!-- v component - x-direction
	1446	!-- WS1 (18), WS3 (19), WS5 (20)
[1762]	1447	IF ( k <= nzb_v_inner(j,i+1) .OR. &
	1448	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
	1449	.AND. i == nxl ) .OR. &
	1450	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
	1451	.AND. i == nxr ) ) &
	1452	THEN
	1453	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
	1454	!
[861]	1455	!-- WS3
[1886]	1456	ELSEIF ( ( k <= nzb_v_inner(j,i+2) .AND. &
	1457	k > nzb_v_inner(j,i+1) ) .OR. &
	1458	k <= nzb_v_inner(j,i-1) &
[1762]	1459	.OR. &
	1460	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
	1461	.AND. i == nxr-1 ) .OR. &
	1462	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
	1463	.AND. i == nxlu ) ) &
	1464	THEN
[861]	1465	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 )
[1580]	1466	!
	1467	!-- Clear flag for WS1
	1468	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 )
[1886]	1469	ELSEIF ( k > nzb_v_inner(j,i+1) .AND. k > nzb_v_inner(j,i+2) &
	1470	.AND. k > nzb_v_inner(j,i-1) ) &
	1471	THEN
[861]	1472	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 )
[1580]	1473	!
	1474	!-- Clear flag for WS1
	1475	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 )
[861]	1476	ENDIF
	1477	!
	1478	!-- v component - y-direction
	1479	!-- WS1 (21), WS3 (22), WS5 (23)
[1762]	1480	IF ( k <= nzb_v_inner(j+1,i) .OR. &
	1481	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
	1482	.AND. j <= nysv ) .OR. &
	1483	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
	1484	.AND. j == nyn ) ) &
	1485	THEN
	1486	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
[1886]	1487	ELSEIF ( ( k <= nzb_v_inner(j+2,i) .AND. &
	1488	k > nzb_v_inner(j+1,i) ) .OR. &
	1489	k <= nzb_v_inner(j-1,i) &
	1490	.OR. &
[1762]	1491	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
	1492	.AND. j == nysv+1) .OR. &
	1493	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
	1494	.AND. j == nyn-1 ) ) &
	1495	THEN
[861]	1496	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
[1580]	1497	!
	1498	!-- Clear flag for WS1
	1499	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 )
[1886]	1500	ELSEIF ( k > nzb_v_inner(j+1,i) .AND. k > nzb_v_inner(j+2,i) &
	1501	.AND. k > nzb_v_inner(j-1,i) ) &
	1502	THEN
[861]	1503	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
[1580]	1504	!
	1505	!-- Clear flag for WS1
	1506	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 )
[861]	1507	ENDIF
	1508	!
	1509	!-- v component - z-direction
	1510	!-- WS1 (24), WS3 (25), WS5 (26)
	1511	flag_set = .FALSE.
	1512	IF ( k == nzb_v_inner(j,i) + 1 .OR. k == nzt ) THEN
	1513	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 )
	1514	flag_set = .TRUE.
	1515	ELSEIF ( k == nzb_v_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN
	1516	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 )
	1517	flag_set = .TRUE.
	1518	ELSEIF ( k > nzb_v_inner(j,i) .AND. .NOT. flag_set ) THEN
	1519	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 )
	1520	ENDIF
	1521
	1522	ENDDO
	1523	ENDDO
	1524	ENDDO
	1525	DO i = nxl, nxr
	1526	DO j = nys, nyn
[1580]	1527	DO k = nzb+1, nzt
[861]	1528	!
[1580]	1529	!-- At first, set flags to WS1.
	1530	!-- Since fluxes are swapped in advec_ws.f90, this is necessary to
	1531	!-- in order to handle the left/south flux.
	1532	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
	1533	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
	1534	!
[861]	1535	!-- w component - x-direction
	1536	!-- WS1 (27), WS3 (28), WS5 (29)
[1762]	1537	IF ( k <= nzb_w_inner(j,i+1) .OR. &
	1538	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
	1539	.AND. i == nxl ) .OR. &
	1540	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
	1541	.AND. i == nxr ) ) &
	1542	THEN
	1543	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
[1886]	1544	ELSEIF ( ( k <= nzb_w_inner(j,i+2) .AND. &
	1545	k > nzb_w_inner(j,i+1) ) .OR. &
	1546	k <= nzb_w_inner(j,i-1) &
	1547	.OR. &
[1762]	1548	( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) &
	1549	.AND. i == nxr-1 ) .OR. &
	1550	( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) &
	1551	.AND. i == nxlu ) ) &
	1552	THEN
[861]	1553	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 )
[1580]	1554	!
	1555	!-- Clear flag for WS1
	1556	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 )
[1886]	1557	ELSEIF ( k > nzb_w_inner(j,i+1) .AND. k > nzb_w_inner(j,i+2) &
	1558	.AND. k > nzb_w_inner(j,i-1) ) &
	1559	THEN
[861]	1560	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i),29 )
[1580]	1561	!
	1562	!-- Clear flag for WS1
	1563	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 )
[861]	1564	ENDIF
	1565	!
	1566	!-- w component - y-direction
	1567	!-- WS1 (30), WS3 (31), WS5 (32)
[1762]	1568	IF ( k <= nzb_w_inner(j+1,i) .OR. &
	1569	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
	1570	.AND. j == nys ) .OR. &
	1571	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
	1572	.AND. j == nyn ) ) &
	1573	THEN
	1574	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
[1886]	1575	ELSEIF ( ( k <= nzb_w_inner(j+2,i) .AND. &
	1576	k > nzb_w_inner(j+1,i) ) .OR. &
	1577	k <= nzb_w_inner(j-1,i) &
	1578	.OR. &
[1762]	1579	( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) &
	1580	.AND. j == nysv ) .OR. &
	1581	( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) &
	1582	.AND. j == nyn-1 ) ) &
	1583	THEN
[861]	1584	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 )
[1580]	1585	!
	1586	!-- Clear flag for WS1
	1587	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 )
[1886]	1588	ELSEIF ( k > nzb_w_inner(j+1,i) .AND. k > nzb_w_inner(j+2,i) &
	1589	.AND. k > nzb_w_inner(j-1,i) ) &
	1590	THEN
[1221]	1591	wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 0 )
[1580]	1592	!
	1593	!-- Clear flag for WS1
	1594	wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 )
[861]	1595	ENDIF
	1596	!
	1597	!-- w component - z-direction
[1679]	1598	!-- WS1 (33), WS3 (34), WS5 (35)
[861]	1599	flag_set = .FALSE.
	1600	IF ( k == nzb_w_inner(j,i) .OR. k == nzb_w_inner(j,i) + 1 &
	1601	.OR. k == nzt ) THEN
	1602	!
	1603	!-- Please note, at k == nzb_w_inner(j,i) a flag is explictely
	1604	!-- set, although this is not a prognostic level. However,
	1605	!-- contrary to the advection of u,v and s this is necessary
	1606	!-- because flux_t(nzb_w_inner(j,i)) is used for the tendency
	1607	!-- at k == nzb_w_inner(j,i)+1.
[1221]	1608	wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 1 )
[861]	1609	flag_set = .TRUE.
	1610	ELSEIF ( k == nzb_w_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN
[1221]	1611	wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 2 )
[861]	1612	flag_set = .TRUE.
	1613	ELSEIF ( k > nzb_w_inner(j,i) .AND. .NOT. flag_set ) THEN
[1221]	1614	wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 3 )
[861]	1615	ENDIF
	1616
	1617	ENDDO
	1618	ENDDO
	1619	ENDDO
	1620
	1621	ENDIF
	1622
	1623	!
[1677]	1624	!-- Exchange 3D integer wall_flags.
	1625	IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' &
	1626	.OR. scalar_advec == 'ws-scheme-mono' ) THEN
	1627	!
	1628	!-- Exchange ghost points for advection flags
	1629	CALL exchange_horiz_int( wall_flags_0, nbgp )
	1630	CALL exchange_horiz_int( wall_flags_00, nbgp )
	1631	!
	1632	!-- Set boundary flags at inflow and outflow boundary in case of
	1633	!-- non-cyclic boundary conditions.
[1762]	1634	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN
[1677]	1635	wall_flags_0(:,:,nxl-1) = wall_flags_0(:,:,nxl)
	1636	wall_flags_00(:,:,nxl-1) = wall_flags_00(:,:,nxl)
	1637	ENDIF
	1638
[1762]	1639	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN
[1677]	1640	wall_flags_0(:,:,nxr+1) = wall_flags_0(:,:,nxr)
	1641	wall_flags_00(:,:,nxr+1) = wall_flags_00(:,:,nxr)
	1642	ENDIF
	1643
[1762]	1644	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN
[1677]	1645	wall_flags_0(:,nyn+1,:) = wall_flags_0(:,nyn,:)
	1646	wall_flags_00(:,nyn+1,:) = wall_flags_00(:,nyn,:)
	1647	ENDIF
	1648
[1762]	1649	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN
[1677]	1650	wall_flags_0(:,nys-1,:) = wall_flags_0(:,nys,:)
	1651	wall_flags_00(:,nys-1,:) = wall_flags_00(:,nys,:)
	1652	ENDIF
	1653
	1654	ENDIF
	1655
	1656	!
[1]	1657	!-- In case of topography: limit near-wall mixing length l_wall further:
	1658	!-- Go through all points of the subdomain one by one and look for the closest
	1659	!-- surface
	1660	IF ( TRIM(topography) /= 'flat' ) THEN
	1661	DO i = nxl, nxr
	1662	DO j = nys, nyn
	1663
	1664	nzb_si = nzb_s_inner(j,i)
	1665	vi = vertical_influence(nzb_si)
	1666
	1667	IF ( wall_n(j,i) > 0 ) THEN
	1668	!
	1669	!-- North wall (y distance)
	1670	DO k = wall_n(j,i), nzb_si
[1353]	1671	l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5_wp * dy )
[1]	1672	ENDDO
	1673	!
	1674	!-- Above North wall (yz distance)
	1675	DO k = nzb_si + 1, nzb_si + vi
[1353]	1676	l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), &
	1677	SQRT( 0.25_wp * dy**2 + &
[1]	1678	( zu(k) - zw(nzb_si) )**2 ) )
	1679	ENDDO
	1680	!
	1681	!-- Northleft corner (xy distance)
	1682	IF ( corner_nl(j,i) > 0 ) THEN
	1683	DO k = corner_nl(j,i), nzb_si
	1684	l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
[1353]	1685	0.5_wp * SQRT( dx2 + dy2 ) )
[1]	1686	ENDDO
	1687	!
	1688	!-- Above Northleft corner (xyz distance)
	1689	DO k = nzb_si + 1, nzb_si + vi
[1353]	1690	l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
	1691	SQRT( 0.25_wp * (dx2 + dy2) + &
	1692	( zu(k) - zw(nzb_si) )**2 ) )
[1]	1693	ENDDO
	1694	ENDIF
	1695	!
	1696	!-- Northright corner (xy distance)
	1697	IF ( corner_nr(j,i) > 0 ) THEN
	1698	DO k = corner_nr(j,i), nzb_si
[1353]	1699	l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
	1700	0.5_wp * SQRT( dx2 + dy2 ) )
[1]	1701	ENDDO
	1702	!
	1703	!-- Above northright corner (xyz distance)
	1704	DO k = nzb_si + 1, nzb_si + vi
[1353]	1705	l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
	1706	SQRT( 0.25_wp * (dx2 + dy2) + &
	1707	( zu(k) - zw(nzb_si) )**2 ) )
[1]	1708	ENDDO
	1709	ENDIF
	1710	ENDIF
	1711
	1712	IF ( wall_s(j,i) > 0 ) THEN
	1713	!
	1714	!-- South wall (y distance)
	1715	DO k = wall_s(j,i), nzb_si
[1353]	1716	l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), 0.5_wp * dy )
[1]	1717	ENDDO
	1718	!
	1719	!-- Above south wall (yz distance)
[1353]	1720	DO k = nzb_si + 1, nzb_si + vi
	1721	l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), &
	1722	SQRT( 0.25_wp * dy**2 + &
[1]	1723	( zu(k) - zw(nzb_si) )**2 ) )
	1724	ENDDO
	1725	!
	1726	!-- Southleft corner (xy distance)
	1727	IF ( corner_sl(j,i) > 0 ) THEN
	1728	DO k = corner_sl(j,i), nzb_si
[1353]	1729	l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
	1730	0.5_wp * SQRT( dx2 + dy2 ) )
[1]	1731	ENDDO
	1732	!
	1733	!-- Above southleft corner (xyz distance)
	1734	DO k = nzb_si + 1, nzb_si + vi
[1353]	1735	l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
	1736	SQRT( 0.25_wp * (dx2 + dy2) + &
	1737	( zu(k) - zw(nzb_si) )**2 ) )
[1]	1738	ENDDO
	1739	ENDIF
	1740	!
	1741	!-- Southright corner (xy distance)
	1742	IF ( corner_sr(j,i) > 0 ) THEN
	1743	DO k = corner_sr(j,i), nzb_si
[1353]	1744	l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
	1745	0.5_wp * SQRT( dx2 + dy2 ) )
[1]	1746	ENDDO
	1747	!
	1748	!-- Above southright corner (xyz distance)
	1749	DO k = nzb_si + 1, nzb_si + vi
[1353]	1750	l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
	1751	SQRT( 0.25_wp * (dx2 + dy2) + &
	1752	( zu(k) - zw(nzb_si) )**2 ) )
[1]	1753	ENDDO
	1754	ENDIF
	1755
	1756	ENDIF
	1757
	1758	IF ( wall_l(j,i) > 0 ) THEN
	1759	!
	1760	!-- Left wall (x distance)
	1761	DO k = wall_l(j,i), nzb_si
[1353]	1762	l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), 0.5_wp * dx )
[1]	1763	ENDDO
	1764	!
	1765	!-- Above left wall (xz distance)
	1766	DO k = nzb_si + 1, nzb_si + vi
[1353]	1767	l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), &
	1768	SQRT( 0.25_wp * dx**2 + &
	1769	( zu(k) - zw(nzb_si) )**2 ) )
[1]	1770	ENDDO
	1771	ENDIF
	1772
	1773	IF ( wall_r(j,i) > 0 ) THEN
	1774	!
	1775	!-- Right wall (x distance)
	1776	DO k = wall_r(j,i), nzb_si
[1353]	1777	l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5_wp * dx )
[1]	1778	ENDDO
	1779	!
	1780	!-- Above right wall (xz distance)
	1781	DO k = nzb_si + 1, nzb_si + vi
[1353]	1782	l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), &
	1783	SQRT( 0.25_wp * dx**2 + &
[1]	1784	( zu(k) - zw(nzb_si) )**2 ) )
	1785	ENDDO
	1786
	1787	ENDIF
	1788
	1789	ENDDO
	1790	ENDDO
	1791
	1792	ENDIF
	1793
	1794	!
	1795	!-- Multiplication with wall_adjustment_factor
	1796	l_wall = wall_adjustment_factor * l_wall
	1797
	1798	!
[709]	1799	!-- Set lateral boundary conditions for l_wall
[667]	1800	CALL exchange_horiz( l_wall, nbgp )
	1801
[1]	1802	DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, &
	1803	nzb_tmp, vertical_influence, wall_l, wall_n, wall_r, wall_s )
	1804
	1805
	1806	END SUBROUTINE init_grid

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |