Home

Context Navigation

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

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 55.3 KB

Rev	Line
[1]	1	SUBROUTINE init_grid
	2
	3	!------------------------------------------------------------------------------!
[254]	4	! Current revisions:
[1]	5	! -----------------
[1015]	6	! lower index for calculating wall_flags_0 set to nzb_w_inner instead of
	7	! nzb_w_inner+1
[863]	8	!
	9	! Former revisions:
	10	! -----------------
	11	! $Id: init_grid.f90 1015 2012-09-27 09:23:24Z raasch $
	12	!
[997]	13	! 996 2012-09-07 10:41:47Z raasch
	14	! little reformatting
	15	!
[979]	16	! 978 2012-08-09 08:28:32Z fricke
	17	! Bugfix: nzb_max is set to nzt at non-cyclic lateral boundaries
	18	! Bugfix: Set wall_flags_0 for inflow boundary
	19	!
[928]	20	! 927 2012-06-06 19:15:04Z raasch
	21	! Wall flags are not set for multigrid method in case of masking method
	22	!
[865]	23	! 864 2012-03-27 15:10:33Z gryschka
[927]	24	! In case of ocean and Dirichlet bottom bc for u and v dzu_mg and ddzu_pres
	25	! were not correctly defined for k=1.
[865]	26	!
[863]	27	! 861 2012-03-26 14:18:34Z suehring
[861]	28	! Set wall_flags_0. The array is needed for degradation in ws-scheme near walls,
	29	! inflow and outflow boundaries as well as near the bottom and the top of the
[863]	30	! model domain.!
[861]	31	! Initialization of nzb_s_inner and nzb_w_inner.
	32	! gls has to be at least nbgp to do not exceed the array bounds of nzb_local
	33	! while setting wall_flags_0
	34	!
[844]	35	! 843 2012-02-29 15:16:21Z gryschka
	36	! In case of ocean and dirichlet bc for u and v at the bottom
	37	! the first u-level ist defined at same height as the first w-level
	38	!
[819]	39	! 818 2012-02-08 16:11:23Z maronga
	40	! Bugfix: topo_height is only required if topography is used. It is thus now
	41	! allocated in the topography branch
	42	!
[810]	43	! 809 2012-01-30 13:32:58Z maronga
	44	! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
	45	!
[808]	46	! 807 2012-01-25 11:53:51Z maronga
	47	! New cpp directive "__check" implemented which is used by check_namelist_files
	48	!
[760]	49	! 759 2011-09-15 13:58:31Z raasch
	50	! Splitting of parallel I/O in blocks of PEs
	51	!
[723]	52	! 722 2011-04-11 06:21:09Z raasch
	53	! Bugfix: bc_lr/ns_cyc replaced by bc_lr/ns, because variables are not yet set
	54	! here
	55	!
[710]	56	! 709 2011-03-30 09:31:40Z raasch
	57	! formatting adjustments
	58	!
[708]	59	! 707 2011-03-29 11:39:40Z raasch
	60	! bc_lr/ns replaced by bc_lr/ns_cyc
	61	!
[668]	62	! 667 2010-12-23 12:06:00Z suehring/gryschka
[667]	63	! Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
	64	! Furthermore the allocation of arrays and steering of loops is done with these
	65	! parameters. Call of exchange_horiz are modified.
	66	! In case of dirichlet bounday condition at the bottom zu(0)=0.0
	67	! dzu_mg has to be set explicitly for a equally spaced grid near bottom.
	68	! ddzu_pres added to use a equally spaced grid near bottom.
[139]	69	!
[556]	70	! 555 2010-09-07 07:32:53Z raasch
[818]	71	! Bugfix: default setting of nzb_local for flat topography
[556]	72	!
[392]	73	! 274 2009-03-26 15:11:21Z heinze
	74	! Output of messages replaced by message handling routine.
	75	! new topography case 'single_street_canyon'
	76	!
[226]	77	! 217 2008-12-09 18:00:48Z letzel
	78	! +topography_grid_convention
	79	!
[139]	80	! 134 2007-11-21 07:28:38Z letzel
[134]	81	! Redefine initial nzb_local as the actual total size of topography (later the
	82	! extent of topography in nzb_local is reduced by 1dx at the E topography walls
	83	! and by 1dy at the N topography walls to form the basis for nzb_s_inner);
	84	! for consistency redefine 'single_building' case.
[114]	85	! Calculation of wall flag arrays
[1]	86	!
[98]	87	! 94 2007-06-01 15:25:22Z raasch
	88	! Grid definition for ocean version
	89	!
[77]	90	! 75 2007-03-22 09:54:05Z raasch
	91	! storage of topography height arrays zu_s_inner and zw_s_inner,
	92	! 2nd+3rd argument removed from exchange horiz
	93	!
[39]	94	! 19 2007-02-23 04:53:48Z raasch
	95	! Setting of nzt_diff
	96	!
[3]	97	! RCS Log replace by Id keyword, revision history cleaned up
	98	!
[1]	99	! Revision 1.17 2006/08/22 14:00:05 raasch
	100	! +dz_max to limit vertical stretching,
	101	! bugfix in index array initialization for line- or point-like topography
	102	! structures
	103	!
	104	! Revision 1.1 1997/08/11 06:17:45 raasch
	105	! Initial revision (Testversion)
	106	!
	107	!
	108	! Description:
	109	! ------------
	110	! Creating grid depending constants
	111	!------------------------------------------------------------------------------!
	112
	113	USE arrays_3d
	114	USE control_parameters
	115	USE grid_variables
	116	USE indices
	117	USE pegrid
	118
	119	IMPLICIT NONE
	120
[240]	121	INTEGER :: bh, blx, bly, bxl, bxr, byn, bys, ch, cwx, cwy, cxl, cxr, cyn, &
[759]	122	cys, gls, i, ii, inc, i_center, j, j_center, k, l, nxl_l, &
	123	nxr_l, nyn_l, nys_l, nzb_si, nzt_l, vi
[1]	124
	125	INTEGER, DIMENSION(:), ALLOCATABLE :: vertical_influence
	126
	127	INTEGER, DIMENSION(:,:), ALLOCATABLE :: corner_nl, corner_nr, corner_sl, &
	128	corner_sr, wall_l, wall_n, wall_r,&
	129	wall_s, nzb_local, nzb_tmp
	130
[861]	131	LOGICAL :: flag_set = .FALSE.
	132
[1]	133	REAL :: dx_l, dy_l, dz_stretched
	134
[818]	135	REAL, DIMENSION(:,:), ALLOCATABLE :: topo_height
[1]	136
	137	REAL, DIMENSION(:,:,:), ALLOCATABLE :: distance
[667]	138
[1]	139	!
[709]	140	!-- Calculation of horizontal array bounds including ghost layers
[667]	141	nxlg = nxl - nbgp
	142	nxrg = nxr + nbgp
	143	nysg = nys - nbgp
	144	nyng = nyn + nbgp
[709]	145
[667]	146	!
[1]	147	!-- Allocate grid arrays
	148	ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &
[667]	149	dzw(1:nzt+1), l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) )
[1]	150
	151	!
	152	!-- Compute height of u-levels from constant grid length and dz stretch factors
	153	IF ( dz == -1.0 ) THEN
[254]	154	message_string = 'missing dz'
	155	CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
[1]	156	ELSEIF ( dz <= 0.0 ) THEN
[254]	157	WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
	158	CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
[1]	159	ENDIF
[94]	160
[1]	161	!
[94]	162	!-- Define the vertical grid levels
	163	IF ( .NOT. ocean ) THEN
	164	!
	165	!-- Grid for atmosphere with surface at z=0 (k=0, w-grid).
[843]	166	!-- The second u-level (k=1) corresponds to the top of the
[94]	167	!-- Prandtl-layer.
[667]	168
	169	IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
	170	zu(0) = 0.0
	171	! zu(0) = - dz * 0.5
	172	ELSE
	173	zu(0) = - dz * 0.5
	174	ENDIF
[94]	175	zu(1) = dz * 0.5
[1]	176
[94]	177	dz_stretch_level_index = nzt+1
	178	dz_stretched = dz
	179	DO k = 2, nzt+1
	180	IF ( dz_stretch_level <= zu(k-1) .AND. dz_stretched < dz_max ) THEN
	181	dz_stretched = dz_stretched * dz_stretch_factor
	182	dz_stretched = MIN( dz_stretched, dz_max )
	183	IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
	184	ENDIF
	185	zu(k) = zu(k-1) + dz_stretched
	186	ENDDO
[1]	187
	188	!
[94]	189	!-- Compute the w-levels. They are always staggered half-way between the
[843]	190	!-- corresponding u-levels. In case of dirichlet bc for u and v at the
	191	!-- ground the first u- and w-level (k=0) are defined at same height (z=0).
	192	!-- The top w-level is extrapolated linearly.
[94]	193	zw(0) = 0.0
	194	DO k = 1, nzt
	195	zw(k) = ( zu(k) + zu(k+1) ) * 0.5
	196	ENDDO
	197	zw(nzt+1) = zw(nzt) + 2.0 * ( zu(nzt+1) - zw(nzt) )
[1]	198
[94]	199	ELSE
[1]	200	!
[843]	201	!-- Grid for ocean with free water surface is at k=nzt (w-grid).
	202	!-- In case of neumann bc at the ground the first first u-level (k=0) lies
	203	!-- below the first w-level (k=0). In case of dirichlet bc the first u- and
	204	!-- w-level are defined at same height, but staggered from the second level.
	205	!-- The second u-level (k=1) corresponds to the top of the Prandtl-layer.
[94]	206	zu(nzt+1) = dz * 0.5
	207	zu(nzt) = - dz * 0.5
	208
	209	dz_stretch_level_index = 0
	210	dz_stretched = dz
	211	DO k = nzt-1, 0, -1
	212	IF ( dz_stretch_level <= ABS( zu(k+1) ) .AND. &
	213	dz_stretched < dz_max ) THEN
	214	dz_stretched = dz_stretched * dz_stretch_factor
	215	dz_stretched = MIN( dz_stretched, dz_max )
	216	IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
	217	ENDIF
	218	zu(k) = zu(k+1) - dz_stretched
	219	ENDDO
	220
	221	!
	222	!-- Compute the w-levels. They are always staggered half-way between the
[843]	223	!-- corresponding u-levels, except in case of dirichlet bc for u and v
	224	!-- at the ground. In this case the first u- and w-level are defined at
	225	!-- same height. The top w-level (nzt+1) is not used but set for
	226	!-- consistency, since w and all scalar variables are defined up tp nzt+1.
[94]	227	zw(nzt+1) = dz
	228	zw(nzt) = 0.0
	229	DO k = 0, nzt
	230	zw(k) = ( zu(k) + zu(k+1) ) * 0.5
	231	ENDDO
	232
[843]	233	!
	234	!-- In case of dirichlet bc for u and v the first u- and w-level are defined
	235	!-- at same height.
	236	IF ( ibc_uv_b == 0 ) THEN
	237	zu(0) = zw(0)
	238	ENDIF
	239
[94]	240	ENDIF
	241
	242	!
[1]	243	!-- Compute grid lengths.
	244	DO k = 1, nzt+1
	245	dzu(k) = zu(k) - zu(k-1)
	246	ddzu(k) = 1.0 / dzu(k)
	247	dzw(k) = zw(k) - zw(k-1)
	248	ddzw(k) = 1.0 / dzw(k)
	249	ENDDO
	250
	251	DO k = 1, nzt
	252	dd2zu(k) = 1.0 / ( dzu(k) + dzu(k+1) )
	253	ENDDO
[667]	254
	255	!
[709]	256	!-- The FFT- SOR-pressure solvers assume grid spacings of a staggered grid
	257	!-- everywhere. For the actual grid, the grid spacing at the lowest level
	258	!-- is only dz/2, but should be dz. Therefore, an additional array
	259	!-- containing with appropriate grid information is created for these
	260	!-- solvers.
	261	IF ( psolver /= 'multigrid' ) THEN
[667]	262	ALLOCATE( ddzu_pres(1:nzt+1) )
	263	ddzu_pres = ddzu
[864]	264	ddzu_pres(1) = ddzu_pres(2) ! change for lowest level
[667]	265	ENDIF
[1]	266
	267	!
	268	!-- In case of multigrid method, compute grid lengths and grid factors for the
	269	!-- grid levels
	270	IF ( psolver == 'multigrid' ) THEN
	271
	272	ALLOCATE( ddx2_mg(maximum_grid_level), ddy2_mg(maximum_grid_level), &
	273	dzu_mg(nzb+1:nzt+1,maximum_grid_level), &
	274	dzw_mg(nzb+1:nzt+1,maximum_grid_level), &
	275	f1_mg(nzb+1:nzt,maximum_grid_level), &
	276	f2_mg(nzb+1:nzt,maximum_grid_level), &
	277	f3_mg(nzb+1:nzt,maximum_grid_level) )
	278
	279	dzu_mg(:,maximum_grid_level) = dzu
[667]	280	!
[864]	281	!-- Next line to ensure an equally spaced grid.
	282	dzu_mg(1,maximum_grid_level) = dzu(2)
[709]	283
[1]	284	dzw_mg(:,maximum_grid_level) = dzw
	285	nzt_l = nzt
	286	DO l = maximum_grid_level-1, 1, -1
	287	dzu_mg(nzb+1,l) = 2.0 * dzu_mg(nzb+1,l+1)
	288	dzw_mg(nzb+1,l) = 2.0 * dzw_mg(nzb+1,l+1)
	289	nzt_l = nzt_l / 2
	290	DO k = 2, nzt_l+1
	291	dzu_mg(k,l) = dzu_mg(2k-2,l+1) + dzu_mg(2k-1,l+1)
	292	dzw_mg(k,l) = dzw_mg(2k-2,l+1) + dzw_mg(2k-1,l+1)
	293	ENDDO
	294	ENDDO
	295
	296	nzt_l = nzt
	297	dx_l = dx
	298	dy_l = dy
	299	DO l = maximum_grid_level, 1, -1
	300	ddx2_mg(l) = 1.0 / dx_l**2
	301	ddy2_mg(l) = 1.0 / dy_l**2
	302	DO k = nzb+1, nzt_l
	303	f2_mg(k,l) = 1.0 / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
	304	f3_mg(k,l) = 1.0 / ( dzu_mg(k,l) * dzw_mg(k,l) )
	305	f1_mg(k,l) = 2.0 * ( ddx2_mg(l) + ddy2_mg(l) ) + &
	306	f2_mg(k,l) + f3_mg(k,l)
	307	ENDDO
	308	nzt_l = nzt_l / 2
	309	dx_l = dx_l * 2.0
	310	dy_l = dy_l * 2.0
	311	ENDDO
	312
	313	ENDIF
	314
	315	!
	316	!-- Compute the reciprocal values of the horizontal grid lengths.
	317	ddx = 1.0 / dx
	318	ddy = 1.0 / dy
	319	dx2 = dx * dx
	320	dy2 = dy * dy
	321	ddx2 = 1.0 / dx2
	322	ddy2 = 1.0 / dy2
	323
	324	!
	325	!-- Compute the grid-dependent mixing length.
	326	DO k = 1, nzt
	327	l_grid(k) = ( dx * dy * dzw(k) )**0.33333333333333
	328	ENDDO
	329
	330	!
	331	!-- Allocate outer and inner index arrays for topography and set
[114]	332	!-- defaults.
	333	!-- nzb_local has to contain additional layers of ghost points for calculating
	334	!-- the flag arrays needed for the multigrid method
	335	gls = 2**( maximum_grid_level )
[861]	336	IF ( gls < nbgp ) gls = nbgp
[667]	337
[114]	338	ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), &
	339	corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), &
[667]	340	nzb_local(-gls:ny+gls,-gls:nx+gls), &
	341	nzb_tmp(-nbgp:ny+nbgp,-nbgp:nx+nbgp), &
[114]	342	wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr), &
[1]	343	wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
[667]	344	ALLOCATE( fwxm(nysg:nyng,nxlg:nxrg), fwxp(nysg:nyng,nxlg:nxrg), &
	345	fwym(nysg:nyng,nxlg:nxrg), fwyp(nysg:nyng,nxlg:nxrg), &
	346	fxm(nysg:nyng,nxlg:nxrg), fxp(nysg:nyng,nxlg:nxrg), &
	347	fym(nysg:nyng,nxlg:nxrg), fyp(nysg:nyng,nxlg:nxrg), &
	348	nzb_s_inner(nysg:nyng,nxlg:nxrg), &
	349	nzb_s_outer(nysg:nyng,nxlg:nxrg), &
	350	nzb_u_inner(nysg:nyng,nxlg:nxrg), &
	351	nzb_u_outer(nysg:nyng,nxlg:nxrg), &
	352	nzb_v_inner(nysg:nyng,nxlg:nxrg), &
	353	nzb_v_outer(nysg:nyng,nxlg:nxrg), &
	354	nzb_w_inner(nysg:nyng,nxlg:nxrg), &
	355	nzb_w_outer(nysg:nyng,nxlg:nxrg), &
	356	nzb_diff_s_inner(nysg:nyng,nxlg:nxrg), &
	357	nzb_diff_s_outer(nysg:nyng,nxlg:nxrg), &
	358	nzb_diff_u(nysg:nyng,nxlg:nxrg), &
	359	nzb_diff_v(nysg:nyng,nxlg:nxrg), &
	360	nzb_2d(nysg:nyng,nxlg:nxrg), &
	361	wall_e_x(nysg:nyng,nxlg:nxrg), &
	362	wall_e_y(nysg:nyng,nxlg:nxrg), &
	363	wall_u(nysg:nyng,nxlg:nxrg), &
	364	wall_v(nysg:nyng,nxlg:nxrg), &
	365	wall_w_x(nysg:nyng,nxlg:nxrg), &
	366	wall_w_y(nysg:nyng,nxlg:nxrg) )
[1]	367
	368
[667]	369
	370	ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	371
[818]	372
[1]	373	nzb_s_inner = nzb; nzb_s_outer = nzb
	374	nzb_u_inner = nzb; nzb_u_outer = nzb
	375	nzb_v_inner = nzb; nzb_v_outer = nzb
	376	nzb_w_inner = nzb; nzb_w_outer = nzb
	377
	378	!
[19]	379	!-- Define vertical gridpoint from (or to) which on the usual finite difference
[1]	380	!-- form (which does not use surface fluxes) is applied
	381	IF ( prandtl_layer .OR. use_surface_fluxes ) THEN
	382	nzb_diff = nzb + 2
	383	ELSE
	384	nzb_diff = nzb + 1
	385	ENDIF
[19]	386	IF ( use_top_fluxes ) THEN
	387	nzt_diff = nzt - 1
	388	ELSE
	389	nzt_diff = nzt
	390	ENDIF
[1]	391
	392	nzb_diff_s_inner = nzb_diff; nzb_diff_s_outer = nzb_diff
	393	nzb_diff_u = nzb_diff; nzb_diff_v = nzb_diff
	394
	395	wall_e_x = 0.0; wall_e_y = 0.0; wall_u = 0.0; wall_v = 0.0
	396	wall_w_x = 0.0; wall_w_y = 0.0
	397	fwxp = 1.0; fwxm = 1.0; fwyp = 1.0; fwym = 1.0
	398	fxp = 1.0; fxm = 1.0; fyp = 1.0; fym = 1.0
	399
	400	!
	401	!-- Initialize near-wall mixing length l_wall only in the vertical direction
	402	!-- for the moment,
	403	!-- multiplication with wall_adjustment_factor near the end of this routine
	404	l_wall(nzb,:,:) = l_grid(1)
	405	DO k = nzb+1, nzt
	406	l_wall(k,:,:) = l_grid(k)
	407	ENDDO
	408	l_wall(nzt+1,:,:) = l_grid(nzt)
	409
	410	ALLOCATE ( vertical_influence(nzb:nzt) )
	411	DO k = 1, nzt
	412	vertical_influence(k) = MIN ( INT( l_grid(k) / &
	413	( wall_adjustment_factor * dzw(k) ) + 0.5 ), nzt - k )
	414	ENDDO
	415
	416	DO k = 1, MAXVAL( nzb_s_inner )
	417	IF ( l_grid(k) > 1.5 * dx * wall_adjustment_factor .OR. &
	418	l_grid(k) > 1.5 * dy * wall_adjustment_factor ) THEN
[254]	419	WRITE( message_string, * ) 'grid anisotropy exceeds ', &
	420	'threshold given by only local', &
	421	' &horizontal reduction of near_wall ', &
	422	'mixing length l_wall', &
	423	' &starting from height level k = ', k, '.'
	424	CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
[1]	425	EXIT
	426	ENDIF
	427	ENDDO
	428	vertical_influence(0) = vertical_influence(1)
	429
[667]	430	DO i = nxlg, nxrg
	431	DO j = nysg, nyng
[1]	432	DO k = nzb_s_inner(j,i) + 1, &
	433	nzb_s_inner(j,i) + vertical_influence(nzb_s_inner(j,i))
	434	l_wall(k,j,i) = zu(k) - zw(nzb_s_inner(j,i))
	435	ENDDO
	436	ENDDO
	437	ENDDO
	438
	439	!
	440	!-- Set outer and inner index arrays for non-flat topography.
	441	!-- Here consistency checks concerning domain size and periodicity are
	442	!-- necessary.
	443	!-- Within this SELECT CASE structure only nzb_local is initialized
	444	!-- individually depending on the chosen topography type, all other index
	445	!-- arrays are initialized further below.
	446	SELECT CASE ( TRIM( topography ) )
	447
	448	CASE ( 'flat' )
	449	!
[555]	450	!-- nzb_local is required for the multigrid solver
	451	nzb_local = 0
[1]	452
	453	CASE ( 'single_building' )
	454	!
	455	!-- Single rectangular building, by default centered in the middle of the
	456	!-- total domain
	457	blx = NINT( building_length_x / dx )
	458	bly = NINT( building_length_y / dy )
	459	bh = NINT( building_height / dz )
	460
	461	IF ( building_wall_left == 9999999.9 ) THEN
	462	building_wall_left = ( nx + 1 - blx ) / 2 * dx
	463	ENDIF
	464	bxl = NINT( building_wall_left / dx )
	465	bxr = bxl + blx
	466
	467	IF ( building_wall_south == 9999999.9 ) THEN
	468	building_wall_south = ( ny + 1 - bly ) / 2 * dy
	469	ENDIF
	470	bys = NINT( building_wall_south / dy )
	471	byn = bys + bly
	472
	473	!
	474	!-- Building size has to meet some requirements
	475	IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR. &
	476	( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) ) THEN
[274]	477	WRITE( message_string, * ) 'inconsistent building parameters:', &
	478	'& bxl=', bxl, 'bxr=', bxr, 'bys=', bys, &
	479	'byn=', byn, 'nx=', nx, 'ny=', ny
[254]	480	CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
[1]	481	ENDIF
	482
	483	!
[217]	484	!-- Define the building.
[1]	485	nzb_local = 0
[134]	486	nzb_local(bys:byn,bxl:bxr) = bh
[1]	487
[240]	488	CASE ( 'single_street_canyon' )
	489	!
	490	!-- Single quasi-2D street canyon of infinite length in x or y direction.
	491	!-- The canyon is centered in the other direction by default.
	492	IF ( canyon_width_x /= 9999999.9 ) THEN
	493	!
	494	!-- Street canyon in y direction
	495	cwx = NINT( canyon_width_x / dx )
	496	IF ( canyon_wall_left == 9999999.9 ) THEN
	497	canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
	498	ENDIF
	499	cxl = NINT( canyon_wall_left / dx )
	500	cxr = cxl + cwx
	501
	502	ELSEIF ( canyon_width_y /= 9999999.9 ) THEN
	503	!
	504	!-- Street canyon in x direction
	505	cwy = NINT( canyon_width_y / dy )
	506	IF ( canyon_wall_south == 9999999.9 ) THEN
	507	canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
	508	ENDIF
	509	cys = NINT( canyon_wall_south / dy )
	510	cyn = cys + cwy
	511
	512	ELSE
[254]	513
	514	message_string = 'no street canyon width given'
	515	CALL message( 'init_grid', 'PA0204', 1, 2, 0, 6, 0 )
	516
[240]	517	ENDIF
	518
	519	ch = NINT( canyon_height / dz )
	520	dp_level_ind_b = ch
	521	!
	522	!-- Street canyon size has to meet some requirements
	523	IF ( canyon_width_x /= 9999999.9 ) THEN
	524	IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR. &
	525	( ch < 3 ) ) THEN
[254]	526	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
[274]	527	'&cxl=', cxl, 'cxr=', cxr, &
	528	'cwx=', cwx, &
[254]	529	'ch=', ch, 'nx=', nx, 'ny=', ny
	530	CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 )
[240]	531	ENDIF
	532	ELSEIF ( canyon_width_y /= 9999999.9 ) THEN
	533	IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR. &
	534	( ch < 3 ) ) THEN
[254]	535	WRITE( message_string, * ) 'inconsistent canyon parameters:', &
[274]	536	'&cys=', cys, 'cyn=', cyn, &
	537	'cwy=', cwy, &
[254]	538	'ch=', ch, 'nx=', nx, 'ny=', ny
	539	CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 )
[240]	540	ENDIF
	541	ENDIF
[274]	542	IF ( canyon_width_x /= 9999999.9 .AND. canyon_width_y /= 9999999.9 ) &
[240]	543	THEN
[274]	544	message_string = 'inconsistent canyon parameters:' // &
[254]	545	'&street canyon can only be oriented' // &
	546	'&either in x- or in y-direction'
	547	CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 )
[240]	548	ENDIF
	549
	550	nzb_local = ch
	551	IF ( canyon_width_x /= 9999999.9 ) THEN
	552	nzb_local(:,cxl+1:cxr-1) = 0
	553	ELSEIF ( canyon_width_y /= 9999999.9 ) THEN
	554	nzb_local(cys+1:cyn-1,:) = 0
	555	ENDIF
	556
[1]	557	CASE ( 'read_from_file' )
[759]	558
[818]	559	ALLOCATE ( topo_height(0:ny,0:nx) )
	560
[759]	561	DO ii = 0, io_blocks-1
	562	IF ( ii == io_group ) THEN
	563
[1]	564	!
[759]	565	!-- Arbitrary irregular topography data in PALM format (exactly
	566	!-- matching the grid size and total domain size)
	567	OPEN( 90, FILE='TOPOGRAPHY_DATA', STATUS='OLD', &
	568	FORM='FORMATTED', ERR=10 )
	569	DO j = ny, 0, -1
	570	READ( 90, *, ERR=11, END=11 ) ( topo_height(j,i), i = 0,nx )
	571	ENDDO
	572
	573	GOTO 12
	574
	575	10 message_string = 'file TOPOGRAPHY_DATA does not exist'
	576	CALL message( 'init_grid', 'PA0208', 1, 2, 0, 6, 0 )
	577
	578	11 message_string = 'errors in file TOPOGRAPHY_DATA'
	579	CALL message( 'init_grid', 'PA0209', 1, 2, 0, 6, 0 )
	580
	581	12 CLOSE( 90 )
	582
	583	ENDIF
[809]	584	#if defined( __parallel ) && ! defined ( __check )
[759]	585	CALL MPI_BARRIER( comm2d, ierr )
	586	#endif
[559]	587	ENDDO
[759]	588
[1]	589	!
	590	!-- Calculate the index height of the topography
	591	DO i = 0, nx
	592	DO j = 0, ny
	593	nzb_local(j,i) = NINT( topo_height(j,i) / dz )
	594	ENDDO
	595	ENDDO
[818]	596
	597	DEALLOCATE ( topo_height )
[114]	598	!
[759]	599	!-- Add cyclic boundaries (additional layers are for calculating
	600	!-- flag arrays needed for the multigrid sover)
[114]	601	nzb_local(-gls:-1,0:nx) = nzb_local(ny-gls+1:ny,0:nx)
	602	nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,0:nx)
	603	nzb_local(:,-gls:-1) = nzb_local(:,nx-gls+1:nx)
	604	nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
[667]	605
[1]	606	CASE DEFAULT
	607	!
	608	!-- The DEFAULT case is reached either if the parameter topography
[217]	609	!-- contains a wrong character string or if the user has defined a special
[1]	610	!-- case in the user interface. There, the subroutine user_init_grid
	611	!-- checks which of these two conditions applies.
[114]	612	CALL user_init_grid( gls, nzb_local )
[1]	613
	614	END SELECT
	615	!
[861]	616	!-- Determine the maximum level of topography. Furthermore it is used for
	617	!-- steering the degradation of order of the applied advection scheme.
[978]	618	!-- In case of non-cyclic lateral boundaries, the order of the advection
[996]	619	!-- scheme have to be reduced up to nzt (required at the lateral boundaries).
[861]	620	nzb_max = MAXVAL( nzb_local )
[978]	621	IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR. &
	622	inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s ) THEN
	623	nzb_max = nzt
	624	ENDIF
	625
[861]	626	!
[1]	627	!-- Consistency checks and index array initialization are only required for
[217]	628	!-- non-flat topography, also the initialization of topography height arrays
[49]	629	!-- zu_s_inner and zw_w_inner
[1]	630	IF ( TRIM( topography ) /= 'flat' ) THEN
	631
	632	!
	633	!-- Consistency checks
	634	IF ( MINVAL( nzb_local ) < 0 .OR. MAXVAL( nzb_local ) > nz + 1 ) THEN
[274]	635	WRITE( message_string, * ) 'nzb_local values are outside the', &
	636	'model domain', &
	637	'&MINVAL( nzb_local ) = ', MINVAL(nzb_local), &
	638	'&MAXVAL( nzb_local ) = ', MAXVAL(nzb_local)
[254]	639	CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
[1]	640	ENDIF
	641
[722]	642	IF ( bc_lr == 'cyclic' ) THEN
[1]	643	IF ( ANY( nzb_local(:,-1) /= nzb_local(:,nx) ) .OR. &
	644	ANY( nzb_local(:,0) /= nzb_local(:,nx+1) ) ) THEN
[254]	645	message_string = 'nzb_local does not fulfill cyclic' // &
	646	' boundary condition in x-direction'
	647	CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
[1]	648	ENDIF
	649	ENDIF
[722]	650	IF ( bc_ns == 'cyclic' ) THEN
[1]	651	IF ( ANY( nzb_local(-1,:) /= nzb_local(ny,:) ) .OR. &
	652	ANY( nzb_local(0,:) /= nzb_local(ny+1,:) ) ) THEN
[254]	653	message_string = 'nzb_local does not fulfill cyclic' // &
	654	' boundary condition in y-direction'
	655	CALL message( 'init_grid', 'PA0212', 1, 2, 0, 6, 0 )
[1]	656	ENDIF
	657	ENDIF
	658
[217]	659	IF ( topography_grid_convention == 'cell_edge' ) THEN
[134]	660	!
[217]	661	!-- The array nzb_local as defined using the 'cell_edge' convention
	662	!-- describes the actual total size of topography which is defined at the
	663	!-- cell edges where u=0 on the topography walls in x-direction and v=0
	664	!-- on the topography walls in y-direction. However, PALM uses individual
	665	!-- arrays nzb_u\|v\|w\|s_inner\|outer that are based on nzb_s_inner.
	666	!-- Therefore, the extent of topography in nzb_local is now reduced by
	667	!-- 1dx at the E topography walls and by 1dy at the N topography walls
	668	!-- to form the basis for nzb_s_inner.
	669	DO j = -gls, ny + gls
	670	DO i = -gls, nx
	671	nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j,i+1) )
	672	ENDDO
[134]	673	ENDDO
[217]	674	!-- apply cyclic boundary conditions in x-direction
	675	!(ist das erforderlich? Ursache von Seung Bus Fehler?)
	676	nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
	677	DO i = -gls, nx + gls
	678	DO j = -gls, ny
	679	nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j+1,i) )
	680	ENDDO
[134]	681	ENDDO
[217]	682	!-- apply cyclic boundary conditions in y-direction
	683	!(ist das erforderlich? Ursache von Seung Bus Fehler?)
	684	nzb_local(ny+1:ny+gls,:) = nzb_local(0:gls-1,:)
	685	ENDIF
[134]	686
[1]	687	!
	688	!-- Initialize index arrays nzb_s_inner and nzb_w_inner
[861]	689	nzb_s_inner = nzb_local(nysg:nyng,nxlg:nxrg)
	690	nzb_w_inner = nzb_local(nysg:nyng,nxlg:nxrg)
[1]	691
	692	!
	693	!-- Initialize remaining index arrays:
	694	!-- first pre-initialize them with nzb_s_inner...
	695	nzb_u_inner = nzb_s_inner
	696	nzb_u_outer = nzb_s_inner
	697	nzb_v_inner = nzb_s_inner
	698	nzb_v_outer = nzb_s_inner
	699	nzb_w_outer = nzb_s_inner
	700	nzb_s_outer = nzb_s_inner
	701
	702	!
	703	!-- ...then extend pre-initialized arrays in their according directions
	704	!-- based on nzb_local using nzb_tmp as a temporary global index array
	705
	706	!
	707	!-- nzb_s_outer:
	708	!-- extend nzb_local east-/westwards first, then north-/southwards
[667]	709	nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
[1]	710	DO j = -1, ny + 1
	711	DO i = 0, nx
	712	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i), &
	713	nzb_local(j,i+1) )
	714	ENDDO
	715	ENDDO
	716	DO i = nxl, nxr
	717	DO j = nys, nyn
	718	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
	719	nzb_tmp(j+1,i) )
	720	ENDDO
	721	!
	722	!-- non-cyclic boundary conditions (overwritten by call of
	723	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
	724	IF ( nys == 0 ) THEN
	725	j = -1
	726	nzb_s_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
	727	ENDIF
	728	IF ( nys == ny ) THEN
	729	j = ny + 1
	730	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
	731	ENDIF
	732	ENDDO
	733	!
	734	!-- nzb_w_outer:
	735	!-- identical to nzb_s_outer
	736	nzb_w_outer = nzb_s_outer
	737
	738	!
	739	!-- nzb_u_inner:
	740	!-- extend nzb_local rightwards only
[667]	741	nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
[1]	742	DO j = -1, ny + 1
	743	DO i = 0, nx + 1
	744	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )
	745	ENDDO
	746	ENDDO
[667]	747	nzb_u_inner = nzb_tmp(nysg:nyng,nxlg:nxrg)
[1]	748
	749	!
	750	!-- nzb_u_outer:
	751	!-- extend current nzb_tmp (nzb_u_inner) north-/southwards
	752	DO i = nxl, nxr
	753	DO j = nys, nyn
	754	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
	755	nzb_tmp(j+1,i) )
	756	ENDDO
	757	!
	758	!-- non-cyclic boundary conditions (overwritten by call of
	759	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
	760	IF ( nys == 0 ) THEN
	761	j = -1
	762	nzb_u_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
	763	ENDIF
	764	IF ( nys == ny ) THEN
	765	j = ny + 1
	766	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
	767	ENDIF
	768	ENDDO
	769
	770	!
	771	!-- nzb_v_inner:
	772	!-- extend nzb_local northwards only
[667]	773	nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
[1]	774	DO i = -1, nx + 1
	775	DO j = 0, ny + 1
	776	nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )
	777	ENDDO
	778	ENDDO
[667]	779	nzb_v_inner = nzb_tmp(nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp)
[1]	780
	781	!
	782	!-- nzb_v_outer:
	783	!-- extend current nzb_tmp (nzb_v_inner) right-/leftwards
	784	DO j = nys, nyn
	785	DO i = nxl, nxr
	786	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i), &
	787	nzb_tmp(j,i+1) )
	788	ENDDO
	789	!
	790	!-- non-cyclic boundary conditions (overwritten by call of
	791	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
	792	IF ( nxl == 0 ) THEN
	793	i = -1
	794	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i+1), nzb_tmp(j,i) )
	795	ENDIF
	796	IF ( nxr == nx ) THEN
	797	i = nx + 1
	798	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i) )
	799	ENDIF
	800	ENDDO
[809]	801	#if ! defined ( __check )
[1]	802	!
	803	!-- Exchange of lateral boundary values (parallel computers) and cyclic
	804	!-- boundary conditions, if applicable.
	805	!-- Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local
	806	!-- they do not require exchange and are not included here.
	807	CALL exchange_horiz_2d_int( nzb_u_inner )
	808	CALL exchange_horiz_2d_int( nzb_u_outer )
	809	CALL exchange_horiz_2d_int( nzb_v_inner )
	810	CALL exchange_horiz_2d_int( nzb_v_outer )
	811	CALL exchange_horiz_2d_int( nzb_w_outer )
	812	CALL exchange_horiz_2d_int( nzb_s_outer )
	813
[49]	814	!
	815	!-- Allocate and set the arrays containing the topography height
	816	IF ( myid == 0 ) THEN
	817
	818	ALLOCATE( zu_s_inner(0:nx+1,0:ny+1), zw_w_inner(0:nx+1,0:ny+1) )
	819
	820	DO i = 0, nx + 1
	821	DO j = 0, ny + 1
	822	zu_s_inner(i,j) = zu(nzb_local(j,i))
	823	zw_w_inner(i,j) = zw(nzb_local(j,i))
	824	ENDDO
	825	ENDDO
	826
	827	ENDIF
[807]	828	#endif
[1]	829	ENDIF
	830
[809]	831	#if ! defined ( __check )
[1]	832	!
	833	!-- Preliminary: to be removed after completion of the topography code!
	834	!-- Set the former default k index arrays nzb_2d
	835	nzb_2d = nzb
	836
	837	!
	838	!-- Set the individual index arrays which define the k index from which on
	839	!-- the usual finite difference form (which does not use surface fluxes) is
	840	!-- applied
	841	IF ( prandtl_layer .OR. use_surface_fluxes ) THEN
	842	nzb_diff_u = nzb_u_inner + 2
	843	nzb_diff_v = nzb_v_inner + 2
	844	nzb_diff_s_inner = nzb_s_inner + 2
	845	nzb_diff_s_outer = nzb_s_outer + 2
	846	ELSE
	847	nzb_diff_u = nzb_u_inner + 1
	848	nzb_diff_v = nzb_v_inner + 1
	849	nzb_diff_s_inner = nzb_s_inner + 1
	850	nzb_diff_s_outer = nzb_s_outer + 1
	851	ENDIF
	852
	853	!
	854	!-- Calculation of wall switches and factors required by diffusion_u/v.f90 and
	855	!-- for limitation of near-wall mixing length l_wall further below
	856	corner_nl = 0
	857	corner_nr = 0
	858	corner_sl = 0
	859	corner_sr = 0
	860	wall_l = 0
	861	wall_n = 0
	862	wall_r = 0
	863	wall_s = 0
	864
	865	DO i = nxl, nxr
	866	DO j = nys, nyn
	867	!
	868	!-- u-component
	869	IF ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) THEN
	870	wall_u(j,i) = 1.0 ! north wall (location of adjacent fluid)
	871	fym(j,i) = 0.0
	872	fyp(j,i) = 1.0
	873	ELSEIF ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) THEN
	874	wall_u(j,i) = 1.0 ! south wall (location of adjacent fluid)
	875	fym(j,i) = 1.0
	876	fyp(j,i) = 0.0
	877	ENDIF
	878	!
	879	!-- v-component
	880	IF ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) THEN
	881	wall_v(j,i) = 1.0 ! rigth wall (location of adjacent fluid)
	882	fxm(j,i) = 0.0
	883	fxp(j,i) = 1.0
	884	ELSEIF ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) ) THEN
	885	wall_v(j,i) = 1.0 ! left wall (location of adjacent fluid)
	886	fxm(j,i) = 1.0
	887	fxp(j,i) = 0.0
	888	ENDIF
	889	!
	890	!-- w-component, also used for scalars, separate arrays for shear
	891	!-- production of tke
	892	IF ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) THEN
	893	wall_e_y(j,i) = 1.0 ! north wall (location of adjacent fluid)
	894	wall_w_y(j,i) = 1.0
	895	fwym(j,i) = 0.0
	896	fwyp(j,i) = 1.0
	897	ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) THEN
	898	wall_e_y(j,i) = -1.0 ! south wall (location of adjacent fluid)
	899	wall_w_y(j,i) = 1.0
	900	fwym(j,i) = 1.0
	901	fwyp(j,i) = 0.0
	902	ENDIF
	903	IF ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) THEN
	904	wall_e_x(j,i) = 1.0 ! right wall (location of adjacent fluid)
	905	wall_w_x(j,i) = 1.0
	906	fwxm(j,i) = 0.0
	907	fwxp(j,i) = 1.0
	908	ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) THEN
	909	wall_e_x(j,i) = -1.0 ! left wall (location of adjacent fluid)
	910	wall_w_x(j,i) = 1.0
	911	fwxm(j,i) = 1.0
	912	fwxp(j,i) = 0.0
	913	ENDIF
	914	!
	915	!-- Wall and corner locations inside buildings for limitation of
	916	!-- near-wall mixing length l_wall
	917	IF ( nzb_s_inner(j,i) > nzb_s_inner(j+1,i) ) THEN
	918
	919	wall_n(j,i) = nzb_s_inner(j+1,i) + 1 ! North wall
	920
	921	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
	922	corner_nl(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northleft corner
	923	nzb_s_inner(j,i-1) ) + 1
	924	ENDIF
	925
	926	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
	927	corner_nr(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northright corner
	928	nzb_s_inner(j,i+1) ) + 1
	929	ENDIF
	930
	931	ENDIF
	932
	933	IF ( nzb_s_inner(j,i) > nzb_s_inner(j-1,i) ) THEN
	934
	935	wall_s(j,i) = nzb_s_inner(j-1,i) + 1 ! South wall
	936	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
	937	corner_sl(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southleft corner
	938	nzb_s_inner(j,i-1) ) + 1
	939	ENDIF
	940
	941	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
	942	corner_sr(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southright corner
	943	nzb_s_inner(j,i+1) ) + 1
	944	ENDIF
	945
	946	ENDIF
	947
	948	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
	949	wall_l(j,i) = nzb_s_inner(j,i-1) + 1 ! Left wall
	950	ENDIF
	951
	952	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
	953	wall_r(j,i) = nzb_s_inner(j,i+1) + 1 ! Right wall
	954	ENDIF
	955
	956	ENDDO
	957	ENDDO
	958
	959	!
[114]	960	!-- Calculate wall flag arrays for the multigrid method
	961	IF ( psolver == 'multigrid' ) THEN
	962	!
	963	!-- Gridpoint increment of the current level
	964	inc = 1
	965
	966	DO l = maximum_grid_level, 1 , -1
	967
	968	nxl_l = nxl_mg(l)
	969	nxr_l = nxr_mg(l)
	970	nys_l = nys_mg(l)
	971	nyn_l = nyn_mg(l)
	972	nzt_l = nzt_mg(l)
	973
	974	!
	975	!-- Assign the flag level to be calculated
	976	SELECT CASE ( l )
	977	CASE ( 1 )
	978	flags => wall_flags_1
	979	CASE ( 2 )
	980	flags => wall_flags_2
	981	CASE ( 3 )
	982	flags => wall_flags_3
	983	CASE ( 4 )
	984	flags => wall_flags_4
	985	CASE ( 5 )
	986	flags => wall_flags_5
	987	CASE ( 6 )
	988	flags => wall_flags_6
	989	CASE ( 7 )
	990	flags => wall_flags_7
	991	CASE ( 8 )
	992	flags => wall_flags_8
	993	CASE ( 9 )
	994	flags => wall_flags_9
	995	CASE ( 10 )
	996	flags => wall_flags_10
	997	END SELECT
	998
	999	!
	1000	!-- Depending on the grid level, set the respective bits in case of
	1001	!-- neighbouring walls
	1002	!-- Bit 0: wall to the bottom
	1003	!-- Bit 1: wall to the top (not realized in remaining PALM code so far)
	1004	!-- Bit 2: wall to the south
	1005	!-- Bit 3: wall to the north
	1006	!-- Bit 4: wall to the left
	1007	!-- Bit 5: wall to the right
[116]	1008	!-- Bit 6: inside building
[114]	1009
	1010	flags = 0
	1011
[927]	1012	!
	1013	!-- In case of masking method, flags are not set and multigrid method
	1014	!-- works like FFT-solver
	1015	IF ( .NOT. masking_method ) THEN
	1016
	1017	DO i = nxl_l-1, nxr_l+1
	1018	DO j = nys_l-1, nyn_l+1
	1019	DO k = nzb, nzt_l+1
[114]	1020
	1021	!
[927]	1022	!-- Inside/outside building (inside building does not need
	1023	!-- further tests for walls)
	1024	IF ( kinc <= nzb_local(jinc,i*inc) ) THEN
[114]	1025
[927]	1026	flags(k,j,i) = IBSET( flags(k,j,i), 6 )
[114]	1027
[927]	1028	ELSE
[114]	1029	!
[927]	1030	!-- Bottom wall
	1031	IF ( (k-1)inc <= nzb_local(jinc,i*inc) ) THEN
	1032	flags(k,j,i) = IBSET( flags(k,j,i), 0 )
	1033	ENDIF
[114]	1034	!
[927]	1035	!-- South wall
	1036	IF ( kinc <= nzb_local((j-1)inc,i*inc) ) THEN
	1037	flags(k,j,i) = IBSET( flags(k,j,i), 2 )
	1038	ENDIF
[114]	1039	!
[927]	1040	!-- North wall
	1041	IF ( kinc <= nzb_local((j+1)inc,i*inc) ) THEN
	1042	flags(k,j,i) = IBSET( flags(k,j,i), 3 )
	1043	ENDIF
[114]	1044	!
[927]	1045	!-- Left wall
	1046	IF ( kinc <= nzb_local(jinc,(i-1)*inc) ) THEN
	1047	flags(k,j,i) = IBSET( flags(k,j,i), 4 )
	1048	ENDIF
[114]	1049	!
[927]	1050	!-- Right wall
	1051	IF ( kinc <= nzb_local(jinc,(i+1)*inc) ) THEN
	1052	flags(k,j,i) = IBSET( flags(k,j,i), 5 )
	1053	ENDIF
	1054
[114]	1055	ENDIF
	1056
[927]	1057	ENDDO
[114]	1058	ENDDO
	1059	ENDDO
	1060
[927]	1061	ENDIF
	1062
[114]	1063	!
	1064	!-- Test output of flag arrays
[145]	1065	! i = nxl_l
	1066	! WRITE (9,*) ' '
	1067	! WRITE (9,) ' mg level ', l, ' *', mg_switch_to_pe0_level
	1068	! WRITE (9,*) ' inc=', inc, ' i =', nxl_l
	1069	! WRITE (9,*) ' nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
	1070	! DO k = nzt_l+1, nzb, -1
	1071	! WRITE (9,'(194(1X,I2))') ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
	1072	! ENDDO
[114]	1073
	1074	inc = inc * 2
	1075
	1076	ENDDO
	1077
	1078	ENDIF
[861]	1079	!
	1080	!-- Allocate flags needed for masking walls.
	1081	ALLOCATE( wall_flags_0(nzb:nzt,nys:nyn,nxl:nxr) )
	1082	wall_flags_0 = 0
[114]	1083
[861]	1084	IF ( scalar_advec == 'ws-scheme' ) THEN
[114]	1085	!
[861]	1086	!-- Set flags to steer the degradation of the advection scheme in advec_ws
	1087	!-- near topography, inflow- and outflow boundaries as well as bottom and
	1088	!-- top of model domain. wall_flags_0 remains zero for all non-prognostic
	1089	!-- grid points.
	1090	DO i = nxl, nxr
	1091	DO j = nys, nyn
	1092	DO k = nzb_s_inner(j,i)+1, nzt
	1093	!
	1094	!-- scalar - x-direction
	1095	!-- WS1 (0), WS3 (1), WS5 (2)
[978]	1096	IF ( k <= nzb_s_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )&
	1097	.AND. i == nxl ) .OR. ( ( inflow_r .OR. outflow_r ) &
	1098	.AND. i == nxr ) ) THEN
[861]	1099	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
	1100	ELSEIF ( k <= nzb_s_inner(j,i+2) .OR. k <= nzb_s_inner(j,i-1) &
[978]	1101	.OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) &
	1102	.OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu ) &
	1103	) THEN
[861]	1104	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
	1105	ELSE
	1106	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 )
	1107	ENDIF
	1108	!
	1109	!-- scalar - y-direction
	1110	!-- WS1 (3), WS3 (4), WS5 (5)
[978]	1111	IF ( k <= nzb_s_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )&
	1112	.AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n ) &
	1113	.AND. j == nyn ) ) THEN
[861]	1114	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 )
	1115	!-- WS3
	1116	ELSEIF ( k <= nzb_s_inner(j+2,i) .OR. k <= nzb_s_inner(j-1,i) &
[978]	1117	.OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv ) &
	1118	.OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) &
	1119	) THEN
[861]	1120	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 )
	1121	!-- WS5
	1122	ELSE
	1123	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 )
	1124	ENDIF
	1125	!
	1126	!-- scalar - z-direction
	1127	!-- WS1 (6), WS3 (7), WS5 (8)
	1128	flag_set = .FALSE.
	1129	IF ( k == nzb_s_inner(j,i) + 1 .OR. k == nzt ) THEN
	1130	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 )
	1131	flag_set = .TRUE.
	1132	ELSEIF ( k == nzb_s_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN
	1133	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 7 )
	1134	flag_set = .TRUE.
	1135	ELSEIF ( k > nzb_s_inner(j,i) .AND. .NOT. flag_set ) THEN
	1136	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 )
	1137	ENDIF
	1138	ENDDO
	1139	ENDDO
	1140	ENDDO
	1141	ENDIF
	1142
	1143	IF ( momentum_advec == 'ws-scheme' ) THEN
	1144	!
	1145	!-- Set wall_flags_0 to steer the degradation of the advection scheme in advec_ws
	1146	!-- near topography, inflow- and outflow boundaries as well as bottom and
	1147	!-- top of model domain. wall_flags_0 remains zero for all non-prognostic
	1148	!-- grid points.
	1149	DO i = nxl, nxr
	1150	DO j = nys, nyn
	1151	DO k = nzb_u_inner(j,i)+1, nzt
	1152	!
	1153	!-- u component - x-direction
	1154	!-- WS1 (9), WS3 (10), WS5 (11)
	1155	IF ( k <= nzb_u_inner(j,i+1) &
[978]	1156	.OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu ) &
	1157	.OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr ) &
	1158	) THEN
[861]	1159	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
	1160	ELSEIF ( k <= nzb_u_inner(j,i+2) .OR. k <= nzb_u_inner(j,i-1) &
[978]	1161	.OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 )&
	1162	.OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu+1)&
	1163	) THEN
[861]	1164	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 )
	1165	ELSE
	1166	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 )
	1167	ENDIF
[978]	1168
[861]	1169	!
	1170	!-- u component - y-direction
	1171	!-- WS1 (12), WS3 (13), WS5 (14)
[978]	1172	IF ( k <= nzb_u_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )&
	1173	.AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n ) &
	1174	.AND. j == nyn ) ) THEN
[861]	1175	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
	1176	ELSEIF ( k <= nzb_u_inner(j+2,i) .OR. k <= nzb_u_inner(j-1,i) &
[978]	1177	.OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv ) &
	1178	.OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) &
	1179	) THEN
[861]	1180	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
	1181	ELSE
	1182	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 )
	1183	ENDIF
	1184	!
	1185	!-- u component - z-direction
	1186	!-- WS1 (15), WS3 (16), WS5 (17)
	1187	flag_set = .FALSE.
	1188	IF ( k == nzb_u_inner(j,i) + 1 .OR. k == nzt ) THEN
	1189	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 )
	1190	flag_set = .TRUE.
	1191	ELSEIF ( k == nzb_u_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN
	1192	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 )
	1193	flag_set = .TRUE.
	1194	ELSEIF ( k > nzb_u_inner(j,i) .AND. .NOT. flag_set ) THEN
	1195	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 )
	1196	ENDIF
	1197
	1198	ENDDO
	1199	ENDDO
	1200	ENDDO
	1201
	1202	DO i = nxl, nxr
	1203	DO j = nys, nyn
	1204	DO k = nzb_v_inner(j,i)+1, nzt
	1205	!
	1206	!-- v component - x-direction
	1207	!-- WS1 (18), WS3 (19), WS5 (20)
[978]	1208	IF ( k <= nzb_v_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )&
	1209	.AND. i == nxl ) .OR. (( inflow_r .OR. outflow_r ) &
	1210	.AND. i == nxr ) ) THEN
[861]	1211	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
	1212	!-- WS3
	1213	ELSEIF ( k <= nzb_v_inner(j,i+2) .OR. k <= nzb_v_inner(j,i-1) &
[978]	1214	.OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) &
	1215	.OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu ) &
	1216	) THEN
[861]	1217	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 )
	1218	ELSE
	1219	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 )
	1220	ENDIF
	1221	!
	1222	!-- v component - y-direction
	1223	!-- WS1 (21), WS3 (22), WS5 (23)
	1224	IF ( k <= nzb_v_inner(j+1,i) &
[978]	1225	.OR. ( ( inflow_s .OR. outflow_s ) .AND. i == nysv ) &
	1226	.OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn ) &
	1227	) THEN
[861]	1228	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
	1229	ELSEIF ( k <= nzb_v_inner(j+2,i) .OR. k <= nzb_v_inner(j-1,i) &
[978]	1230	.OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv+1 )&
	1231	.OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 )&
	1232	) THEN
[861]	1233	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
	1234	ELSE
	1235	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
	1236	ENDIF
	1237	!
	1238	!-- v component - z-direction
	1239	!-- WS1 (24), WS3 (25), WS5 (26)
	1240	flag_set = .FALSE.
	1241	IF ( k == nzb_v_inner(j,i) + 1 .OR. k == nzt ) THEN
	1242	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 )
	1243	flag_set = .TRUE.
	1244	ELSEIF ( k == nzb_v_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN
	1245	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 )
	1246	flag_set = .TRUE.
	1247	ELSEIF ( k > nzb_v_inner(j,i) .AND. .NOT. flag_set ) THEN
	1248	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 )
	1249	ENDIF
	1250
	1251	ENDDO
	1252	ENDDO
	1253	ENDDO
	1254	DO i = nxl, nxr
	1255	DO j = nys, nyn
[1015]	1256	DO k = nzb_w_inner(j,i), nzt
[861]	1257	!
	1258	!-- w component - x-direction
	1259	!-- WS1 (27), WS3 (28), WS5 (29)
[978]	1260	IF ( k <= nzb_w_inner(j,i+1) .OR. ( ( inflow_l .OR. outflow_l )&
	1261	.AND. i == nxl ) .OR. ( ( inflow_r .OR. outflow_r ) &
	1262	.AND. i == nxr ) ) THEN
[861]	1263	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
	1264	ELSEIF ( k <= nzb_w_inner(j,i+2) .OR. k <= nzb_w_inner(j,i-1) &
[978]	1265	.OR. ( ( inflow_r .OR. outflow_r ) .AND. i == nxr-1 ) &
	1266	.OR. ( ( inflow_l .OR. outflow_l ) .AND. i == nxlu ) &
	1267	) THEN
[861]	1268	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 )
	1269	ELSE
	1270	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i),29 )
	1271	ENDIF
	1272	!
	1273	!-- w component - y-direction
	1274	!-- WS1 (30), WS3 (31), WS5 (32)
[978]	1275	IF ( k <= nzb_w_inner(j+1,i) .OR. ( ( inflow_s .OR. outflow_s )&
	1276	.AND. j == nys ) .OR. ( ( inflow_n .OR. outflow_n ) &
	1277	.AND. j == nyn ) ) THEN
[861]	1278	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
	1279	ELSEIF ( k <= nzb_w_inner(j+2,i) .OR. k <= nzb_w_inner(j-1,i) &
[978]	1280	.OR. ( ( inflow_s .OR. outflow_s ) .AND. j == nysv ) &
	1281	.OR. ( ( inflow_n .OR. outflow_n ) .AND. j == nyn-1 ) &
	1282	) THEN
[861]	1283	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 )
	1284	ELSE
	1285	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 32 )
	1286	ENDIF
	1287	!
	1288	!-- w component - z-direction
	1289	!-- WS1 (33), WS3 (34), WS5 (35)
	1290	flag_set = .FALSE.
	1291	IF ( k == nzb_w_inner(j,i) .OR. k == nzb_w_inner(j,i) + 1 &
	1292	.OR. k == nzt ) THEN
	1293	!
	1294	!-- Please note, at k == nzb_w_inner(j,i) a flag is explictely
	1295	!-- set, although this is not a prognostic level. However,
	1296	!-- contrary to the advection of u,v and s this is necessary
	1297	!-- because flux_t(nzb_w_inner(j,i)) is used for the tendency
	1298	!-- at k == nzb_w_inner(j,i)+1.
	1299	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 33 )
	1300	flag_set = .TRUE.
	1301	ELSEIF ( k == nzb_w_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN
	1302	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 34 )
	1303	flag_set = .TRUE.
	1304	ELSEIF ( k > nzb_w_inner(j,i) .AND. .NOT. flag_set ) THEN
	1305	wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 35 )
	1306	ENDIF
	1307
	1308	ENDDO
	1309	ENDDO
	1310	ENDDO
	1311
	1312	ENDIF
	1313
	1314	!
[1]	1315	!-- In case of topography: limit near-wall mixing length l_wall further:
	1316	!-- Go through all points of the subdomain one by one and look for the closest
	1317	!-- surface
	1318	IF ( TRIM(topography) /= 'flat' ) THEN
	1319	DO i = nxl, nxr
	1320	DO j = nys, nyn
	1321
	1322	nzb_si = nzb_s_inner(j,i)
	1323	vi = vertical_influence(nzb_si)
	1324
	1325	IF ( wall_n(j,i) > 0 ) THEN
	1326	!
	1327	!-- North wall (y distance)
	1328	DO k = wall_n(j,i), nzb_si
	1329	l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5 * dy )
	1330	ENDDO
	1331	!
	1332	!-- Above North wall (yz distance)
	1333	DO k = nzb_si + 1, nzb_si + vi
	1334	l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), &
	1335	SQRT( 0.25 * dy**2 + &
	1336	( zu(k) - zw(nzb_si) )**2 ) )
	1337	ENDDO
	1338	!
	1339	!-- Northleft corner (xy distance)
	1340	IF ( corner_nl(j,i) > 0 ) THEN
	1341	DO k = corner_nl(j,i), nzb_si
	1342	l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
	1343	0.5 * SQRT( dx2 + dy2 ) )
	1344	ENDDO
	1345	!
	1346	!-- Above Northleft corner (xyz distance)
	1347	DO k = nzb_si + 1, nzb_si + vi
	1348	l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
	1349	SQRT( 0.25 * (dx2 + dy2) + &
	1350	( zu(k) - zw(nzb_si) )**2 ) )
	1351	ENDDO
	1352	ENDIF
	1353	!
	1354	!-- Northright corner (xy distance)
	1355	IF ( corner_nr(j,i) > 0 ) THEN
	1356	DO k = corner_nr(j,i), nzb_si
	1357	l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
	1358	0.5 * SQRT( dx2 + dy2 ) )
	1359	ENDDO
	1360	!
	1361	!-- Above northright corner (xyz distance)
	1362	DO k = nzb_si + 1, nzb_si + vi
	1363	l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
	1364	SQRT( 0.25 * (dx2 + dy2) + &
	1365	( zu(k) - zw(nzb_si) )**2 ) )
	1366	ENDDO
	1367	ENDIF
	1368	ENDIF
	1369
	1370	IF ( wall_s(j,i) > 0 ) THEN
	1371	!
	1372	!-- South wall (y distance)
	1373	DO k = wall_s(j,i), nzb_si
	1374	l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), 0.5 * dy )
	1375	ENDDO
	1376	!
	1377	!-- Above south wall (yz distance)
	1378	DO k = nzb_si + 1, &
	1379	nzb_si + vi
	1380	l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), &
	1381	SQRT( 0.25 * dy**2 + &
	1382	( zu(k) - zw(nzb_si) )**2 ) )
	1383	ENDDO
	1384	!
	1385	!-- Southleft corner (xy distance)
	1386	IF ( corner_sl(j,i) > 0 ) THEN
	1387	DO k = corner_sl(j,i), nzb_si
	1388	l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
	1389	0.5 * SQRT( dx2 + dy2 ) )
	1390	ENDDO
	1391	!
	1392	!-- Above southleft corner (xyz distance)
	1393	DO k = nzb_si + 1, nzb_si + vi
	1394	l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
	1395	SQRT( 0.25 * (dx2 + dy2) + &
	1396	( zu(k) - zw(nzb_si) )**2 ) )
	1397	ENDDO
	1398	ENDIF
	1399	!
	1400	!-- Southright corner (xy distance)
	1401	IF ( corner_sr(j,i) > 0 ) THEN
	1402	DO k = corner_sr(j,i), nzb_si
	1403	l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
	1404	0.5 * SQRT( dx2 + dy2 ) )
	1405	ENDDO
	1406	!
	1407	!-- Above southright corner (xyz distance)
	1408	DO k = nzb_si + 1, nzb_si + vi
	1409	l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
	1410	SQRT( 0.25 * (dx2 + dy2) + &
	1411	( zu(k) - zw(nzb_si) )**2 ) )
	1412	ENDDO
	1413	ENDIF
	1414
	1415	ENDIF
	1416
	1417	IF ( wall_l(j,i) > 0 ) THEN
	1418	!
	1419	!-- Left wall (x distance)
	1420	DO k = wall_l(j,i), nzb_si
	1421	l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), 0.5 * dx )
	1422	ENDDO
	1423	!
	1424	!-- Above left wall (xz distance)
	1425	DO k = nzb_si + 1, nzb_si + vi
	1426	l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), &
	1427	SQRT( 0.25 * dx**2 + &
	1428	( zu(k) - zw(nzb_si) )**2 ) )
	1429	ENDDO
	1430	ENDIF
	1431
	1432	IF ( wall_r(j,i) > 0 ) THEN
	1433	!
	1434	!-- Right wall (x distance)
	1435	DO k = wall_r(j,i), nzb_si
	1436	l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5 * dx )
	1437	ENDDO
	1438	!
	1439	!-- Above right wall (xz distance)
	1440	DO k = nzb_si + 1, nzb_si + vi
	1441	l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), &
	1442	SQRT( 0.25 * dx**2 + &
	1443	( zu(k) - zw(nzb_si) )**2 ) )
	1444	ENDDO
	1445
	1446	ENDIF
	1447
	1448	ENDDO
	1449	ENDDO
	1450
	1451	ENDIF
	1452
	1453	!
	1454	!-- Multiplication with wall_adjustment_factor
	1455	l_wall = wall_adjustment_factor * l_wall
	1456
	1457	!
[709]	1458	!-- Set lateral boundary conditions for l_wall
[667]	1459	CALL exchange_horiz( l_wall, nbgp )
	1460
[1]	1461	DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, &
	1462	nzb_tmp, vertical_influence, wall_l, wall_n, wall_r, wall_s )
	1463
[807]	1464	#endif
[1]	1465
	1466	END SUBROUTINE init_grid

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |