Home

Context Navigation

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

Last change on this file since 1418 was 1418, checked in by fricke, 10 years ago
Bugfixes concerning grid stretching for the ocean and calculation of the salinity flux in routine surface_coupler
Property svn:keywords set to `Id`
File size: 60.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |