Home

Context Navigation

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

Last change on this file since 251 was 240, checked in by letzel, 16 years ago
External pressure gradient (check_parameters, init_3d_model, header, modules, parin, prognostic_equations) New topography case 'single_street_canyon'
Property svn:keywords set to `Id`
File size: 39.6 KB

Rev	Line
[1]	1	SUBROUTINE init_grid
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
[240]	6	! new topography case 'single_street_canyon'
[139]	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Id: init_grid.f90 240 2009-02-18 17:50:38Z raasch $
	11	!
[226]	12	! 217 2008-12-09 18:00:48Z letzel
	13	! +topography_grid_convention
	14	!
[139]	15	! 134 2007-11-21 07:28:38Z letzel
[134]	16	! Redefine initial nzb_local as the actual total size of topography (later the
	17	! extent of topography in nzb_local is reduced by 1dx at the E topography walls
	18	! and by 1dy at the N topography walls to form the basis for nzb_s_inner);
	19	! for consistency redefine 'single_building' case.
[114]	20	! Calculation of wall flag arrays
[1]	21	!
[98]	22	! 94 2007-06-01 15:25:22Z raasch
	23	! Grid definition for ocean version
	24	!
[77]	25	! 75 2007-03-22 09:54:05Z raasch
	26	! storage of topography height arrays zu_s_inner and zw_s_inner,
	27	! 2nd+3rd argument removed from exchange horiz
	28	!
[39]	29	! 19 2007-02-23 04:53:48Z raasch
	30	! Setting of nzt_diff
	31	!
[3]	32	! RCS Log replace by Id keyword, revision history cleaned up
	33	!
[1]	34	! Revision 1.17 2006/08/22 14:00:05 raasch
	35	! +dz_max to limit vertical stretching,
	36	! bugfix in index array initialization for line- or point-like topography
	37	! structures
	38	!
	39	! Revision 1.1 1997/08/11 06:17:45 raasch
	40	! Initial revision (Testversion)
	41	!
	42	!
	43	! Description:
	44	! ------------
	45	! Creating grid depending constants
	46	!------------------------------------------------------------------------------!
	47
	48	USE arrays_3d
	49	USE control_parameters
	50	USE grid_variables
	51	USE indices
	52	USE pegrid
	53
	54	IMPLICIT NONE
	55
[240]	56	INTEGER :: bh, blx, bly, bxl, bxr, byn, bys, ch, cwx, cwy, cxl, cxr, cyn, &
	57	cys, gls, i, inc, i_center, j, j_center, k, l, nxl_l, nxr_l, &
	58	nyn_l, nys_l, nzb_si, nzt_l, vi
[1]	59
	60	INTEGER, DIMENSION(:), ALLOCATABLE :: vertical_influence
	61
	62	INTEGER, DIMENSION(:,:), ALLOCATABLE :: corner_nl, corner_nr, corner_sl, &
	63	corner_sr, wall_l, wall_n, wall_r,&
	64	wall_s, nzb_local, nzb_tmp
	65
	66	REAL :: dx_l, dy_l, dz_stretched
	67
	68	REAL, DIMENSION(0:ny,0:nx) :: topo_height
	69
	70	REAL, DIMENSION(:,:,:), ALLOCATABLE :: distance
	71
	72	!
	73	!-- Allocate grid arrays
	74	ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &
	75	dzw(1:nzt+1), l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1) )
	76
	77	!
	78	!-- Compute height of u-levels from constant grid length and dz stretch factors
	79	IF ( dz == -1.0 ) THEN
	80	IF ( myid == 0 ) PRINT*,'+++ init_grid: missing dz'
	81	CALL local_stop
	82	ELSEIF ( dz <= 0.0 ) THEN
	83	IF ( myid == 0 ) PRINT*,'+++ init_grid: dz=',dz,' <= 0.0'
	84	CALL local_stop
	85	ENDIF
[94]	86
[1]	87	!
[94]	88	!-- Define the vertical grid levels
	89	IF ( .NOT. ocean ) THEN
	90	!
	91	!-- Grid for atmosphere with surface at z=0 (k=0, w-grid).
	92	!-- Since the w-level lies on the surface, the first u-level (staggered!)
	93	!-- lies below the surface (used for "mirror" boundary condition).
	94	!-- The first u-level above the surface corresponds to the top of the
	95	!-- Prandtl-layer.
	96	zu(0) = - dz * 0.5
	97	zu(1) = dz * 0.5
[1]	98
[94]	99	dz_stretch_level_index = nzt+1
	100	dz_stretched = dz
	101	DO k = 2, nzt+1
	102	IF ( dz_stretch_level <= zu(k-1) .AND. dz_stretched < dz_max ) THEN
	103	dz_stretched = dz_stretched * dz_stretch_factor
	104	dz_stretched = MIN( dz_stretched, dz_max )
	105	IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
	106	ENDIF
	107	zu(k) = zu(k-1) + dz_stretched
	108	ENDDO
[1]	109
	110	!
[94]	111	!-- Compute the w-levels. They are always staggered half-way between the
	112	!-- corresponding u-levels. The top w-level is extrapolated linearly.
	113	zw(0) = 0.0
	114	DO k = 1, nzt
	115	zw(k) = ( zu(k) + zu(k+1) ) * 0.5
	116	ENDDO
	117	zw(nzt+1) = zw(nzt) + 2.0 * ( zu(nzt+1) - zw(nzt) )
[1]	118
[94]	119	ELSE
[1]	120	!
[94]	121	!-- Grid for ocean with solid surface at z=0 (k=0, w-grid). The free water
	122	!-- surface is at k=nzt (w-grid).
	123	!-- Since the w-level lies always on the surface, the first/last u-level
	124	!-- (staggered!) lies below the bottom surface / above the free surface.
	125	!-- It is used for "mirror" boundary condition.
	126	!-- The first u-level above the bottom surface corresponds to the top of the
	127	!-- Prandtl-layer.
	128	zu(nzt+1) = dz * 0.5
	129	zu(nzt) = - dz * 0.5
	130
	131	dz_stretch_level_index = 0
	132	dz_stretched = dz
	133	DO k = nzt-1, 0, -1
	134	IF ( dz_stretch_level <= ABS( zu(k+1) ) .AND. &
	135	dz_stretched < dz_max ) THEN
	136	dz_stretched = dz_stretched * dz_stretch_factor
	137	dz_stretched = MIN( dz_stretched, dz_max )
	138	IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
	139	ENDIF
	140	zu(k) = zu(k+1) - dz_stretched
	141	ENDDO
	142
	143	!
	144	!-- Compute the w-levels. They are always staggered half-way between the
	145	!-- corresponding u-levels.
	146	!-- The top w-level (nzt+1) is not used but set for consistency, since
	147	!-- w and all scalar variables are defined up tp nzt+1.
	148	zw(nzt+1) = dz
	149	zw(nzt) = 0.0
	150	DO k = 0, nzt
	151	zw(k) = ( zu(k) + zu(k+1) ) * 0.5
	152	ENDDO
	153
	154	ENDIF
	155
	156	!
[1]	157	!-- Compute grid lengths.
	158	DO k = 1, nzt+1
	159	dzu(k) = zu(k) - zu(k-1)
	160	ddzu(k) = 1.0 / dzu(k)
	161	dzw(k) = zw(k) - zw(k-1)
	162	ddzw(k) = 1.0 / dzw(k)
	163	ENDDO
	164
	165	DO k = 1, nzt
	166	dd2zu(k) = 1.0 / ( dzu(k) + dzu(k+1) )
	167	ENDDO
	168
	169	!
	170	!-- In case of multigrid method, compute grid lengths and grid factors for the
	171	!-- grid levels
	172	IF ( psolver == 'multigrid' ) THEN
	173
	174	ALLOCATE( ddx2_mg(maximum_grid_level), ddy2_mg(maximum_grid_level), &
	175	dzu_mg(nzb+1:nzt+1,maximum_grid_level), &
	176	dzw_mg(nzb+1:nzt+1,maximum_grid_level), &
	177	f1_mg(nzb+1:nzt,maximum_grid_level), &
	178	f2_mg(nzb+1:nzt,maximum_grid_level), &
	179	f3_mg(nzb+1:nzt,maximum_grid_level) )
	180
	181	dzu_mg(:,maximum_grid_level) = dzu
	182	dzw_mg(:,maximum_grid_level) = dzw
	183	nzt_l = nzt
	184	DO l = maximum_grid_level-1, 1, -1
	185	dzu_mg(nzb+1,l) = 2.0 * dzu_mg(nzb+1,l+1)
	186	dzw_mg(nzb+1,l) = 2.0 * dzw_mg(nzb+1,l+1)
	187	nzt_l = nzt_l / 2
	188	DO k = 2, nzt_l+1
	189	dzu_mg(k,l) = dzu_mg(2k-2,l+1) + dzu_mg(2k-1,l+1)
	190	dzw_mg(k,l) = dzw_mg(2k-2,l+1) + dzw_mg(2k-1,l+1)
	191	ENDDO
	192	ENDDO
	193
	194	nzt_l = nzt
	195	dx_l = dx
	196	dy_l = dy
	197	DO l = maximum_grid_level, 1, -1
	198	ddx2_mg(l) = 1.0 / dx_l**2
	199	ddy2_mg(l) = 1.0 / dy_l**2
	200	DO k = nzb+1, nzt_l
	201	f2_mg(k,l) = 1.0 / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
	202	f3_mg(k,l) = 1.0 / ( dzu_mg(k,l) * dzw_mg(k,l) )
	203	f1_mg(k,l) = 2.0 * ( ddx2_mg(l) + ddy2_mg(l) ) + &
	204	f2_mg(k,l) + f3_mg(k,l)
	205	ENDDO
	206	nzt_l = nzt_l / 2
	207	dx_l = dx_l * 2.0
	208	dy_l = dy_l * 2.0
	209	ENDDO
	210
	211	ENDIF
	212
	213	!
	214	!-- Compute the reciprocal values of the horizontal grid lengths.
	215	ddx = 1.0 / dx
	216	ddy = 1.0 / dy
	217	dx2 = dx * dx
	218	dy2 = dy * dy
	219	ddx2 = 1.0 / dx2
	220	ddy2 = 1.0 / dy2
	221
	222	!
	223	!-- Compute the grid-dependent mixing length.
	224	DO k = 1, nzt
	225	l_grid(k) = ( dx * dy * dzw(k) )**0.33333333333333
	226	ENDDO
	227
	228	!
	229	!-- Allocate outer and inner index arrays for topography and set
[114]	230	!-- defaults.
	231	!-- nzb_local has to contain additional layers of ghost points for calculating
	232	!-- the flag arrays needed for the multigrid method
	233	gls = 2**( maximum_grid_level )
	234	ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), &
	235	corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), &
	236	nzb_local(-gls:ny+gls,-gls:nx+gls), nzb_tmp(-1:ny+1,-1:nx+1), &
	237	wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr), &
[1]	238	wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
	239	ALLOCATE( fwxm(nys-1:nyn+1,nxl-1:nxr+1), fwxp(nys-1:nyn+1,nxl-1:nxr+1), &
	240	fwym(nys-1:nyn+1,nxl-1:nxr+1), fwyp(nys-1:nyn+1,nxl-1:nxr+1), &
	241	fxm(nys-1:nyn+1,nxl-1:nxr+1), fxp(nys-1:nyn+1,nxl-1:nxr+1), &
	242	fym(nys-1:nyn+1,nxl-1:nxr+1), fyp(nys-1:nyn+1,nxl-1:nxr+1), &
	243	nzb_s_inner(nys-1:nyn+1,nxl-1:nxr+1), &
	244	nzb_s_outer(nys-1:nyn+1,nxl-1:nxr+1), &
	245	nzb_u_inner(nys-1:nyn+1,nxl-1:nxr+1), &
	246	nzb_u_outer(nys-1:nyn+1,nxl-1:nxr+1), &
	247	nzb_v_inner(nys-1:nyn+1,nxl-1:nxr+1), &
	248	nzb_v_outer(nys-1:nyn+1,nxl-1:nxr+1), &
	249	nzb_w_inner(nys-1:nyn+1,nxl-1:nxr+1), &
	250	nzb_w_outer(nys-1:nyn+1,nxl-1:nxr+1), &
	251	nzb_diff_s_inner(nys-1:nyn+1,nxl-1:nxr+1), &
	252	nzb_diff_s_outer(nys-1:nyn+1,nxl-1:nxr+1), &
	253	nzb_diff_u(nys-1:nyn+1,nxl-1:nxr+1), &
	254	nzb_diff_v(nys-1:nyn+1,nxl-1:nxr+1), &
	255	nzb_2d(nys-1:nyn+1,nxl-1:nxr+1), &
	256	wall_e_x(nys-1:nyn+1,nxl-1:nxr+1), &
	257	wall_e_y(nys-1:nyn+1,nxl-1:nxr+1), &
	258	wall_u(nys-1:nyn+1,nxl-1:nxr+1), &
	259	wall_v(nys-1:nyn+1,nxl-1:nxr+1), &
	260	wall_w_x(nys-1:nyn+1,nxl-1:nxr+1), &
	261	wall_w_y(nys-1:nyn+1,nxl-1:nxr+1) )
	262
	263	ALLOCATE( l_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
	264
	265	nzb_s_inner = nzb; nzb_s_outer = nzb
	266	nzb_u_inner = nzb; nzb_u_outer = nzb
	267	nzb_v_inner = nzb; nzb_v_outer = nzb
	268	nzb_w_inner = nzb; nzb_w_outer = nzb
	269
	270	!
[19]	271	!-- Define vertical gridpoint from (or to) which on the usual finite difference
[1]	272	!-- form (which does not use surface fluxes) is applied
	273	IF ( prandtl_layer .OR. use_surface_fluxes ) THEN
	274	nzb_diff = nzb + 2
	275	ELSE
	276	nzb_diff = nzb + 1
	277	ENDIF
[19]	278	IF ( use_top_fluxes ) THEN
	279	nzt_diff = nzt - 1
	280	ELSE
	281	nzt_diff = nzt
	282	ENDIF
[1]	283
	284	nzb_diff_s_inner = nzb_diff; nzb_diff_s_outer = nzb_diff
	285	nzb_diff_u = nzb_diff; nzb_diff_v = nzb_diff
	286
	287	wall_e_x = 0.0; wall_e_y = 0.0; wall_u = 0.0; wall_v = 0.0
	288	wall_w_x = 0.0; wall_w_y = 0.0
	289	fwxp = 1.0; fwxm = 1.0; fwyp = 1.0; fwym = 1.0
	290	fxp = 1.0; fxm = 1.0; fyp = 1.0; fym = 1.0
	291
	292	!
	293	!-- Initialize near-wall mixing length l_wall only in the vertical direction
	294	!-- for the moment,
	295	!-- multiplication with wall_adjustment_factor near the end of this routine
	296	l_wall(nzb,:,:) = l_grid(1)
	297	DO k = nzb+1, nzt
	298	l_wall(k,:,:) = l_grid(k)
	299	ENDDO
	300	l_wall(nzt+1,:,:) = l_grid(nzt)
	301
	302	ALLOCATE ( vertical_influence(nzb:nzt) )
	303	DO k = 1, nzt
	304	vertical_influence(k) = MIN ( INT( l_grid(k) / &
	305	( wall_adjustment_factor * dzw(k) ) + 0.5 ), nzt - k )
	306	ENDDO
	307
	308	DO k = 1, MAXVAL( nzb_s_inner )
	309	IF ( l_grid(k) > 1.5 * dx * wall_adjustment_factor .OR. &
	310	l_grid(k) > 1.5 * dy * wall_adjustment_factor ) THEN
	311	IF ( myid == 0 ) THEN
	312	PRINT*, '+++ WARNING: grid anisotropy exceeds '// &
	313	'threshold given by only local'
	314	PRINT*, ' horizontal reduction of near_wall '// &
	315	'mixing length l_wall'
	316	PRINT*, ' starting from height level k = ', k, '.'
	317	ENDIF
	318	EXIT
	319	ENDIF
	320	ENDDO
	321	vertical_influence(0) = vertical_influence(1)
	322
	323	DO i = nxl-1, nxr+1
	324	DO j = nys-1, nyn+1
	325	DO k = nzb_s_inner(j,i) + 1, &
	326	nzb_s_inner(j,i) + vertical_influence(nzb_s_inner(j,i))
	327	l_wall(k,j,i) = zu(k) - zw(nzb_s_inner(j,i))
	328	ENDDO
	329	ENDDO
	330	ENDDO
	331
	332	!
	333	!-- Set outer and inner index arrays for non-flat topography.
	334	!-- Here consistency checks concerning domain size and periodicity are
	335	!-- necessary.
	336	!-- Within this SELECT CASE structure only nzb_local is initialized
	337	!-- individually depending on the chosen topography type, all other index
	338	!-- arrays are initialized further below.
	339	SELECT CASE ( TRIM( topography ) )
	340
	341	CASE ( 'flat' )
	342	!
	343	!-- No actions necessary
	344
	345	CASE ( 'single_building' )
	346	!
	347	!-- Single rectangular building, by default centered in the middle of the
	348	!-- total domain
	349	blx = NINT( building_length_x / dx )
	350	bly = NINT( building_length_y / dy )
	351	bh = NINT( building_height / dz )
	352
	353	IF ( building_wall_left == 9999999.9 ) THEN
	354	building_wall_left = ( nx + 1 - blx ) / 2 * dx
	355	ENDIF
	356	bxl = NINT( building_wall_left / dx )
	357	bxr = bxl + blx
	358
	359	IF ( building_wall_south == 9999999.9 ) THEN
	360	building_wall_south = ( ny + 1 - bly ) / 2 * dy
	361	ENDIF
	362	bys = NINT( building_wall_south / dy )
	363	byn = bys + bly
	364
	365	!
	366	!-- Building size has to meet some requirements
	367	IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR. &
	368	( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) ) THEN
	369	IF ( myid == 0 ) THEN
	370	PRINT*, '+++ init_grid: inconsistent building parameters:'
	371	PRINT*, ' bxl=', bxl, 'bxr=', bxr, 'bys=', bys, &
	372	'byn=', byn, 'nx=', nx, 'ny=', ny
	373	ENDIF
	374	CALL local_stop
	375	ENDIF
	376
	377	!
[217]	378	!-- Define the building.
[1]	379	nzb_local = 0
[134]	380	nzb_local(bys:byn,bxl:bxr) = bh
[1]	381
[240]	382	CASE ( 'single_street_canyon' )
	383	!
	384	!-- Single quasi-2D street canyon of infinite length in x or y direction.
	385	!-- The canyon is centered in the other direction by default.
	386	IF ( canyon_width_x /= 9999999.9 ) THEN
	387	!
	388	!-- Street canyon in y direction
	389	cwx = NINT( canyon_width_x / dx )
	390	IF ( canyon_wall_left == 9999999.9 ) THEN
	391	canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
	392	ENDIF
	393	cxl = NINT( canyon_wall_left / dx )
	394	cxr = cxl + cwx
	395
	396	ELSEIF ( canyon_width_y /= 9999999.9 ) THEN
	397	!
	398	!-- Street canyon in x direction
	399	cwy = NINT( canyon_width_y / dy )
	400	IF ( canyon_wall_south == 9999999.9 ) THEN
	401	canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
	402	ENDIF
	403	cys = NINT( canyon_wall_south / dy )
	404	cyn = cys + cwy
	405
	406	ELSE
	407	IF ( myid == 0 ) THEN
	408	PRINT*, '+++ user_init_grid: no street canyon width given'
	409	ENDIF
	410	CALL local_stop
	411	ENDIF
	412
	413	ch = NINT( canyon_height / dz )
	414	dp_level_ind_b = ch
	415	!
	416	!-- Street canyon size has to meet some requirements
	417	IF ( canyon_width_x /= 9999999.9 ) THEN
	418	IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR. &
	419	( ch < 3 ) ) THEN
	420	IF ( myid == 0 ) THEN
	421	PRINT*, '+++ user_init_grid: inconsistent canyon parameters:'
	422	PRINT*, ' cxl=', cxl, 'cxr=', cxr, &
	423	'cwx=', cwx, &
	424	'ch=', ch, 'nx=', nx, 'ny=', ny
	425	ENDIF
	426	CALL local_stop
	427	ENDIF
	428	ELSEIF ( canyon_width_y /= 9999999.9 ) THEN
	429	IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR. &
	430	( ch < 3 ) ) THEN
	431	IF ( myid == 0 ) THEN
	432	PRINT*, '+++ user_init_grid: inconsistent canyon parameters:'
	433	PRINT*, ' cys=', cys, 'cyn=', cyn, &
	434	'cwy=', cwy, &
	435	'ch=', ch, 'nx=', nx, 'ny=', ny
	436	ENDIF
	437	CALL local_stop
	438	ENDIF
	439	ENDIF
	440	IF ( canyon_width_x /= 9999999.9 .AND. canyon_width_y /= 9999999.9 ) &
	441	THEN
	442	IF ( myid == 0 ) THEN
	443	PRINT*, '+++ user_init_grid: inconsistent canyon parameters:'
	444	PRINT*, ' street canyon can only be oriented'
	445	PRINT*, ' either in x- or in y-direction'
	446	ENDIF
	447	CALL local_stop
	448	ENDIF
	449
	450	nzb_local = ch
	451	IF ( canyon_width_x /= 9999999.9 ) THEN
	452	nzb_local(:,cxl+1:cxr-1) = 0
	453	ELSEIF ( canyon_width_y /= 9999999.9 ) THEN
	454	nzb_local(cys+1:cyn-1,:) = 0
	455	ENDIF
	456
[1]	457	CASE ( 'read_from_file' )
	458	!
	459	!-- Arbitrary irregular topography data in PALM format (exactly matching
	460	!-- the grid size and total domain size)
	461	OPEN( 90, FILE='TOPOGRAPHY_DATA', STATUS='OLD', FORM='FORMATTED', &
	462	ERR=10 )
	463	DO j = ny, 0, -1
	464	READ( 90, *, ERR=11, END=11 ) ( topo_height(j,i), i = 0, nx )
	465	ENDDO
	466	!
	467	!-- Calculate the index height of the topography
	468	DO i = 0, nx
	469	DO j = 0, ny
	470	nzb_local(j,i) = NINT( topo_height(j,i) / dz )
	471	ENDDO
	472	ENDDO
[114]	473	!
	474	!-- Add cyclic boundaries (additional layers are for calculating flag
	475	!-- arrays needed for the multigrid sover)
	476	nzb_local(-gls:-1,0:nx) = nzb_local(ny-gls+1:ny,0:nx)
	477	nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,0:nx)
	478	nzb_local(:,-gls:-1) = nzb_local(:,nx-gls+1:nx)
	479	nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
[1]	480
	481	GOTO 12
	482
	483	10 IF ( myid == 0 ) THEN
	484	PRINT*, '+++ init_grid: file TOPOGRAPHY_DATA does not exist'
	485	ENDIF
	486	CALL local_stop
	487
	488	11 IF ( myid == 0 ) THEN
	489	PRINT*, '+++ init_grid: errors in file TOPOGRAPHY_DATA'
	490	ENDIF
	491	CALL local_stop
	492
	493	12 CLOSE( 90 )
	494
	495	CASE DEFAULT
	496	!
	497	!-- The DEFAULT case is reached either if the parameter topography
[217]	498	!-- contains a wrong character string or if the user has defined a special
[1]	499	!-- case in the user interface. There, the subroutine user_init_grid
	500	!-- checks which of these two conditions applies.
[114]	501	CALL user_init_grid( gls, nzb_local )
[1]	502
	503	END SELECT
	504
	505	!
[114]	506	!-- Test output of nzb_local -1:ny+1,-1:nx+1
[145]	507	! WRITE (9,) ' nzb_local *'
	508	! DO j = ny+1, -1, -1
	509	! WRITE (9,'(194(1X,I2))') ( nzb_local(j,i), i = -1, nx+1 )
	510	! ENDDO
[114]	511
	512	!
[1]	513	!-- Consistency checks and index array initialization are only required for
[217]	514	!-- non-flat topography, also the initialization of topography height arrays
[49]	515	!-- zu_s_inner and zw_w_inner
[1]	516	IF ( TRIM( topography ) /= 'flat' ) THEN
	517
	518	!
	519	!-- Consistency checks
	520	IF ( MINVAL( nzb_local ) < 0 .OR. MAXVAL( nzb_local ) > nz + 1 ) THEN
	521	IF ( myid == 0 ) THEN
	522	PRINT*, '+++ init_grid: nzb_local values are outside the', &
	523	'model domain'
	524	PRINT*, ' MINVAL( nzb_local ) = ', MINVAL(nzb_local)
	525	PRINT*, ' MAXVAL( nzb_local ) = ', MAXVAL(nzb_local)
	526	ENDIF
	527	CALL local_stop
	528	ENDIF
	529
	530	IF ( bc_lr == 'cyclic' ) THEN
	531	IF ( ANY( nzb_local(:,-1) /= nzb_local(:,nx) ) .OR. &
	532	ANY( nzb_local(:,0) /= nzb_local(:,nx+1) ) ) THEN
	533	IF ( myid == 0 ) THEN
	534	PRINT*, '+++ init_grid: nzb_local does not fulfill cyclic', &
	535	' boundary condition in x-direction'
	536	ENDIF
	537	CALL local_stop
	538	ENDIF
	539	ENDIF
	540	IF ( bc_ns == 'cyclic' ) THEN
	541	IF ( ANY( nzb_local(-1,:) /= nzb_local(ny,:) ) .OR. &
	542	ANY( nzb_local(0,:) /= nzb_local(ny+1,:) ) ) THEN
	543	IF ( myid == 0 ) THEN
	544	PRINT*, '+++ init_grid: nzb_local does not fulfill cyclic', &
	545	' boundary condition in y-direction'
	546	ENDIF
	547	CALL local_stop
	548	ENDIF
	549	ENDIF
	550
[217]	551	IF ( topography_grid_convention == 'cell_edge' ) THEN
[134]	552	!
[217]	553	!-- The array nzb_local as defined using the 'cell_edge' convention
	554	!-- describes the actual total size of topography which is defined at the
	555	!-- cell edges where u=0 on the topography walls in x-direction and v=0
	556	!-- on the topography walls in y-direction. However, PALM uses individual
	557	!-- arrays nzb_u\|v\|w\|s_inner\|outer that are based on nzb_s_inner.
	558	!-- Therefore, the extent of topography in nzb_local is now reduced by
	559	!-- 1dx at the E topography walls and by 1dy at the N topography walls
	560	!-- to form the basis for nzb_s_inner.
	561	DO j = -gls, ny + gls
	562	DO i = -gls, nx
	563	nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j,i+1) )
	564	ENDDO
[134]	565	ENDDO
[217]	566	!-- apply cyclic boundary conditions in x-direction
	567	!(ist das erforderlich? Ursache von Seung Bus Fehler?)
	568	nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
	569	DO i = -gls, nx + gls
	570	DO j = -gls, ny
	571	nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j+1,i) )
	572	ENDDO
[134]	573	ENDDO
[217]	574	!-- apply cyclic boundary conditions in y-direction
	575	!(ist das erforderlich? Ursache von Seung Bus Fehler?)
	576	nzb_local(ny+1:ny+gls,:) = nzb_local(0:gls-1,:)
	577	ENDIF
[134]	578
[1]	579	!
	580	!-- Initialize index arrays nzb_s_inner and nzb_w_inner
	581	nzb_s_inner = nzb_local(nys-1:nyn+1,nxl-1:nxr+1)
	582	nzb_w_inner = nzb_local(nys-1:nyn+1,nxl-1:nxr+1)
	583
	584	!
	585	!-- Initialize remaining index arrays:
	586	!-- first pre-initialize them with nzb_s_inner...
	587	nzb_u_inner = nzb_s_inner
	588	nzb_u_outer = nzb_s_inner
	589	nzb_v_inner = nzb_s_inner
	590	nzb_v_outer = nzb_s_inner
	591	nzb_w_outer = nzb_s_inner
	592	nzb_s_outer = nzb_s_inner
	593
	594	!
	595	!-- ...then extend pre-initialized arrays in their according directions
	596	!-- based on nzb_local using nzb_tmp as a temporary global index array
	597
	598	!
	599	!-- nzb_s_outer:
	600	!-- extend nzb_local east-/westwards first, then north-/southwards
[114]	601	nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
[1]	602	DO j = -1, ny + 1
	603	DO i = 0, nx
	604	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i), &
	605	nzb_local(j,i+1) )
	606	ENDDO
	607	ENDDO
	608	DO i = nxl, nxr
	609	DO j = nys, nyn
	610	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
	611	nzb_tmp(j+1,i) )
	612	ENDDO
	613	!
	614	!-- non-cyclic boundary conditions (overwritten by call of
	615	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
	616	IF ( nys == 0 ) THEN
	617	j = -1
	618	nzb_s_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
	619	ENDIF
	620	IF ( nys == ny ) THEN
	621	j = ny + 1
	622	nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
	623	ENDIF
	624	ENDDO
	625	!
	626	!-- nzb_w_outer:
	627	!-- identical to nzb_s_outer
	628	nzb_w_outer = nzb_s_outer
	629
	630	!
	631	!-- nzb_u_inner:
	632	!-- extend nzb_local rightwards only
[114]	633	nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
[1]	634	DO j = -1, ny + 1
	635	DO i = 0, nx + 1
	636	nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )
	637	ENDDO
	638	ENDDO
	639	nzb_u_inner = nzb_tmp(nys-1:nyn+1,nxl-1:nxr+1)
	640
	641	!
	642	!-- nzb_u_outer:
	643	!-- extend current nzb_tmp (nzb_u_inner) north-/southwards
	644	DO i = nxl, nxr
	645	DO j = nys, nyn
	646	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &
	647	nzb_tmp(j+1,i) )
	648	ENDDO
	649	!
	650	!-- non-cyclic boundary conditions (overwritten by call of
	651	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
	652	IF ( nys == 0 ) THEN
	653	j = -1
	654	nzb_u_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )
	655	ENDIF
	656	IF ( nys == ny ) THEN
	657	j = ny + 1
	658	nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )
	659	ENDIF
	660	ENDDO
	661
	662	!
	663	!-- nzb_v_inner:
	664	!-- extend nzb_local northwards only
[114]	665	nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
[1]	666	DO i = -1, nx + 1
	667	DO j = 0, ny + 1
	668	nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )
	669	ENDDO
	670	ENDDO
	671	nzb_v_inner = nzb_tmp(nys-1:nyn+1,nxl-1:nxr+1)
	672
	673	!
	674	!-- nzb_v_outer:
	675	!-- extend current nzb_tmp (nzb_v_inner) right-/leftwards
	676	DO j = nys, nyn
	677	DO i = nxl, nxr
	678	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i), &
	679	nzb_tmp(j,i+1) )
	680	ENDDO
	681	!
	682	!-- non-cyclic boundary conditions (overwritten by call of
	683	!-- exchange_horiz_2d_int below in case of cyclic boundary conditions)
	684	IF ( nxl == 0 ) THEN
	685	i = -1
	686	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i+1), nzb_tmp(j,i) )
	687	ENDIF
	688	IF ( nxr == nx ) THEN
	689	i = nx + 1
	690	nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i) )
	691	ENDIF
	692	ENDDO
	693
	694	!
	695	!-- Exchange of lateral boundary values (parallel computers) and cyclic
	696	!-- boundary conditions, if applicable.
	697	!-- Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local
	698	!-- they do not require exchange and are not included here.
	699	CALL exchange_horiz_2d_int( nzb_u_inner )
	700	CALL exchange_horiz_2d_int( nzb_u_outer )
	701	CALL exchange_horiz_2d_int( nzb_v_inner )
	702	CALL exchange_horiz_2d_int( nzb_v_outer )
	703	CALL exchange_horiz_2d_int( nzb_w_outer )
	704	CALL exchange_horiz_2d_int( nzb_s_outer )
	705
[49]	706	!
	707	!-- Allocate and set the arrays containing the topography height
	708	IF ( myid == 0 ) THEN
	709
	710	ALLOCATE( zu_s_inner(0:nx+1,0:ny+1), zw_w_inner(0:nx+1,0:ny+1) )
	711
	712	DO i = 0, nx + 1
	713	DO j = 0, ny + 1
	714	zu_s_inner(i,j) = zu(nzb_local(j,i))
	715	zw_w_inner(i,j) = zw(nzb_local(j,i))
	716	ENDDO
	717	ENDDO
	718
	719	ENDIF
	720
[1]	721	ENDIF
	722
	723	!
	724	!-- Preliminary: to be removed after completion of the topography code!
	725	!-- Set the former default k index arrays nzb_2d
	726	nzb_2d = nzb
	727
	728	!
	729	!-- Set the individual index arrays which define the k index from which on
	730	!-- the usual finite difference form (which does not use surface fluxes) is
	731	!-- applied
	732	IF ( prandtl_layer .OR. use_surface_fluxes ) THEN
	733	nzb_diff_u = nzb_u_inner + 2
	734	nzb_diff_v = nzb_v_inner + 2
	735	nzb_diff_s_inner = nzb_s_inner + 2
	736	nzb_diff_s_outer = nzb_s_outer + 2
	737	ELSE
	738	nzb_diff_u = nzb_u_inner + 1
	739	nzb_diff_v = nzb_v_inner + 1
	740	nzb_diff_s_inner = nzb_s_inner + 1
	741	nzb_diff_s_outer = nzb_s_outer + 1
	742	ENDIF
	743
	744	!
	745	!-- Calculation of wall switches and factors required by diffusion_u/v.f90 and
	746	!-- for limitation of near-wall mixing length l_wall further below
	747	corner_nl = 0
	748	corner_nr = 0
	749	corner_sl = 0
	750	corner_sr = 0
	751	wall_l = 0
	752	wall_n = 0
	753	wall_r = 0
	754	wall_s = 0
	755
	756	DO i = nxl, nxr
	757	DO j = nys, nyn
	758	!
	759	!-- u-component
	760	IF ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) THEN
	761	wall_u(j,i) = 1.0 ! north wall (location of adjacent fluid)
	762	fym(j,i) = 0.0
	763	fyp(j,i) = 1.0
	764	ELSEIF ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) THEN
	765	wall_u(j,i) = 1.0 ! south wall (location of adjacent fluid)
	766	fym(j,i) = 1.0
	767	fyp(j,i) = 0.0
	768	ENDIF
	769	!
	770	!-- v-component
	771	IF ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) THEN
	772	wall_v(j,i) = 1.0 ! rigth wall (location of adjacent fluid)
	773	fxm(j,i) = 0.0
	774	fxp(j,i) = 1.0
	775	ELSEIF ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) ) THEN
	776	wall_v(j,i) = 1.0 ! left wall (location of adjacent fluid)
	777	fxm(j,i) = 1.0
	778	fxp(j,i) = 0.0
	779	ENDIF
	780	!
	781	!-- w-component, also used for scalars, separate arrays for shear
	782	!-- production of tke
	783	IF ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) THEN
	784	wall_e_y(j,i) = 1.0 ! north wall (location of adjacent fluid)
	785	wall_w_y(j,i) = 1.0
	786	fwym(j,i) = 0.0
	787	fwyp(j,i) = 1.0
	788	ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) THEN
	789	wall_e_y(j,i) = -1.0 ! south wall (location of adjacent fluid)
	790	wall_w_y(j,i) = 1.0
	791	fwym(j,i) = 1.0
	792	fwyp(j,i) = 0.0
	793	ENDIF
	794	IF ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) THEN
	795	wall_e_x(j,i) = 1.0 ! right wall (location of adjacent fluid)
	796	wall_w_x(j,i) = 1.0
	797	fwxm(j,i) = 0.0
	798	fwxp(j,i) = 1.0
	799	ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) THEN
	800	wall_e_x(j,i) = -1.0 ! left wall (location of adjacent fluid)
	801	wall_w_x(j,i) = 1.0
	802	fwxm(j,i) = 1.0
	803	fwxp(j,i) = 0.0
	804	ENDIF
	805	!
	806	!-- Wall and corner locations inside buildings for limitation of
	807	!-- near-wall mixing length l_wall
	808	IF ( nzb_s_inner(j,i) > nzb_s_inner(j+1,i) ) THEN
	809
	810	wall_n(j,i) = nzb_s_inner(j+1,i) + 1 ! North wall
	811
	812	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
	813	corner_nl(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northleft corner
	814	nzb_s_inner(j,i-1) ) + 1
	815	ENDIF
	816
	817	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
	818	corner_nr(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northright corner
	819	nzb_s_inner(j,i+1) ) + 1
	820	ENDIF
	821
	822	ENDIF
	823
	824	IF ( nzb_s_inner(j,i) > nzb_s_inner(j-1,i) ) THEN
	825
	826	wall_s(j,i) = nzb_s_inner(j-1,i) + 1 ! South wall
	827	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
	828	corner_sl(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southleft corner
	829	nzb_s_inner(j,i-1) ) + 1
	830	ENDIF
	831
	832	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
	833	corner_sr(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southright corner
	834	nzb_s_inner(j,i+1) ) + 1
	835	ENDIF
	836
	837	ENDIF
	838
	839	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN
	840	wall_l(j,i) = nzb_s_inner(j,i-1) + 1 ! Left wall
	841	ENDIF
	842
	843	IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN
	844	wall_r(j,i) = nzb_s_inner(j,i+1) + 1 ! Right wall
	845	ENDIF
	846
	847	ENDDO
	848	ENDDO
	849
	850	!
[114]	851	!-- Calculate wall flag arrays for the multigrid method
	852	IF ( psolver == 'multigrid' ) THEN
	853	!
	854	!-- Gridpoint increment of the current level
	855	inc = 1
	856
	857	DO l = maximum_grid_level, 1 , -1
	858
	859	nxl_l = nxl_mg(l)
	860	nxr_l = nxr_mg(l)
	861	nys_l = nys_mg(l)
	862	nyn_l = nyn_mg(l)
	863	nzt_l = nzt_mg(l)
	864
	865	!
	866	!-- Assign the flag level to be calculated
	867	SELECT CASE ( l )
	868	CASE ( 1 )
	869	flags => wall_flags_1
	870	CASE ( 2 )
	871	flags => wall_flags_2
	872	CASE ( 3 )
	873	flags => wall_flags_3
	874	CASE ( 4 )
	875	flags => wall_flags_4
	876	CASE ( 5 )
	877	flags => wall_flags_5
	878	CASE ( 6 )
	879	flags => wall_flags_6
	880	CASE ( 7 )
	881	flags => wall_flags_7
	882	CASE ( 8 )
	883	flags => wall_flags_8
	884	CASE ( 9 )
	885	flags => wall_flags_9
	886	CASE ( 10 )
	887	flags => wall_flags_10
	888	END SELECT
	889
	890	!
	891	!-- Depending on the grid level, set the respective bits in case of
	892	!-- neighbouring walls
	893	!-- Bit 0: wall to the bottom
	894	!-- Bit 1: wall to the top (not realized in remaining PALM code so far)
	895	!-- Bit 2: wall to the south
	896	!-- Bit 3: wall to the north
	897	!-- Bit 4: wall to the left
	898	!-- Bit 5: wall to the right
[116]	899	!-- Bit 6: inside building
[114]	900
	901	flags = 0
	902
	903	DO i = nxl_l-1, nxr_l+1
	904	DO j = nys_l-1, nyn_l+1
	905	DO k = nzb, nzt_l+1
	906
	907	!
	908	!-- Inside/outside building (inside building does not need
	909	!-- further tests for walls)
	910	IF ( kinc <= nzb_local(jinc,i*inc) ) THEN
	911
	912	flags(k,j,i) = IBSET( flags(k,j,i), 6 )
	913
	914	ELSE
	915	!
	916	!-- Bottom wall
	917	IF ( (k-1)inc <= nzb_local(jinc,i*inc) ) THEN
	918	flags(k,j,i) = IBSET( flags(k,j,i), 0 )
	919	ENDIF
	920	!
	921	!-- South wall
	922	IF ( kinc <= nzb_local((j-1)inc,i*inc) ) THEN
	923	flags(k,j,i) = IBSET( flags(k,j,i), 2 )
	924	ENDIF
	925	!
	926	!-- North wall
	927	IF ( kinc <= nzb_local((j+1)inc,i*inc) ) THEN
	928	flags(k,j,i) = IBSET( flags(k,j,i), 3 )
	929	ENDIF
	930	!
	931	!-- Left wall
	932	IF ( kinc <= nzb_local(jinc,(i-1)*inc) ) THEN
	933	flags(k,j,i) = IBSET( flags(k,j,i), 4 )
	934	ENDIF
	935	!
	936	!-- Right wall
	937	IF ( kinc <= nzb_local(jinc,(i+1)*inc) ) THEN
	938	flags(k,j,i) = IBSET( flags(k,j,i), 5 )
	939	ENDIF
	940
	941	ENDIF
	942
	943	ENDDO
	944	ENDDO
	945	ENDDO
	946
	947	!
	948	!-- Test output of flag arrays
[145]	949	! i = nxl_l
	950	! WRITE (9,*) ' '
	951	! WRITE (9,) ' mg level ', l, ' *', mg_switch_to_pe0_level
	952	! WRITE (9,*) ' inc=', inc, ' i =', nxl_l
	953	! WRITE (9,*) ' nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
	954	! DO k = nzt_l+1, nzb, -1
	955	! WRITE (9,'(194(1X,I2))') ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
	956	! ENDDO
[114]	957
	958	inc = inc * 2
	959
	960	ENDDO
	961
	962	ENDIF
	963
	964	!
[1]	965	!-- In case of topography: limit near-wall mixing length l_wall further:
	966	!-- Go through all points of the subdomain one by one and look for the closest
	967	!-- surface
	968	IF ( TRIM(topography) /= 'flat' ) THEN
	969	DO i = nxl, nxr
	970	DO j = nys, nyn
	971
	972	nzb_si = nzb_s_inner(j,i)
	973	vi = vertical_influence(nzb_si)
	974
	975	IF ( wall_n(j,i) > 0 ) THEN
	976	!
	977	!-- North wall (y distance)
	978	DO k = wall_n(j,i), nzb_si
	979	l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5 * dy )
	980	ENDDO
	981	!
	982	!-- Above North wall (yz distance)
	983	DO k = nzb_si + 1, nzb_si + vi
	984	l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), &
	985	SQRT( 0.25 * dy**2 + &
	986	( zu(k) - zw(nzb_si) )**2 ) )
	987	ENDDO
	988	!
	989	!-- Northleft corner (xy distance)
	990	IF ( corner_nl(j,i) > 0 ) THEN
	991	DO k = corner_nl(j,i), nzb_si
	992	l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
	993	0.5 * SQRT( dx2 + dy2 ) )
	994	ENDDO
	995	!
	996	!-- Above Northleft corner (xyz distance)
	997	DO k = nzb_si + 1, nzb_si + vi
	998	l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
	999	SQRT( 0.25 * (dx2 + dy2) + &
	1000	( zu(k) - zw(nzb_si) )**2 ) )
	1001	ENDDO
	1002	ENDIF
	1003	!
	1004	!-- Northright corner (xy distance)
	1005	IF ( corner_nr(j,i) > 0 ) THEN
	1006	DO k = corner_nr(j,i), nzb_si
	1007	l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
	1008	0.5 * SQRT( dx2 + dy2 ) )
	1009	ENDDO
	1010	!
	1011	!-- Above northright corner (xyz distance)
	1012	DO k = nzb_si + 1, nzb_si + vi
	1013	l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
	1014	SQRT( 0.25 * (dx2 + dy2) + &
	1015	( zu(k) - zw(nzb_si) )**2 ) )
	1016	ENDDO
	1017	ENDIF
	1018	ENDIF
	1019
	1020	IF ( wall_s(j,i) > 0 ) THEN
	1021	!
	1022	!-- South wall (y distance)
	1023	DO k = wall_s(j,i), nzb_si
	1024	l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), 0.5 * dy )
	1025	ENDDO
	1026	!
	1027	!-- Above south wall (yz distance)
	1028	DO k = nzb_si + 1, &
	1029	nzb_si + vi
	1030	l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), &
	1031	SQRT( 0.25 * dy**2 + &
	1032	( zu(k) - zw(nzb_si) )**2 ) )
	1033	ENDDO
	1034	!
	1035	!-- Southleft corner (xy distance)
	1036	IF ( corner_sl(j,i) > 0 ) THEN
	1037	DO k = corner_sl(j,i), nzb_si
	1038	l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
	1039	0.5 * SQRT( dx2 + dy2 ) )
	1040	ENDDO
	1041	!
	1042	!-- Above southleft corner (xyz distance)
	1043	DO k = nzb_si + 1, nzb_si + vi
	1044	l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
	1045	SQRT( 0.25 * (dx2 + dy2) + &
	1046	( zu(k) - zw(nzb_si) )**2 ) )
	1047	ENDDO
	1048	ENDIF
	1049	!
	1050	!-- Southright corner (xy distance)
	1051	IF ( corner_sr(j,i) > 0 ) THEN
	1052	DO k = corner_sr(j,i), nzb_si
	1053	l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
	1054	0.5 * SQRT( dx2 + dy2 ) )
	1055	ENDDO
	1056	!
	1057	!-- Above southright corner (xyz distance)
	1058	DO k = nzb_si + 1, nzb_si + vi
	1059	l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
	1060	SQRT( 0.25 * (dx2 + dy2) + &
	1061	( zu(k) - zw(nzb_si) )**2 ) )
	1062	ENDDO
	1063	ENDIF
	1064
	1065	ENDIF
	1066
	1067	IF ( wall_l(j,i) > 0 ) THEN
	1068	!
	1069	!-- Left wall (x distance)
	1070	DO k = wall_l(j,i), nzb_si
	1071	l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), 0.5 * dx )
	1072	ENDDO
	1073	!
	1074	!-- Above left wall (xz distance)
	1075	DO k = nzb_si + 1, nzb_si + vi
	1076	l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), &
	1077	SQRT( 0.25 * dx**2 + &
	1078	( zu(k) - zw(nzb_si) )**2 ) )
	1079	ENDDO
	1080	ENDIF
	1081
	1082	IF ( wall_r(j,i) > 0 ) THEN
	1083	!
	1084	!-- Right wall (x distance)
	1085	DO k = wall_r(j,i), nzb_si
	1086	l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5 * dx )
	1087	ENDDO
	1088	!
	1089	!-- Above right wall (xz distance)
	1090	DO k = nzb_si + 1, nzb_si + vi
	1091	l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), &
	1092	SQRT( 0.25 * dx**2 + &
	1093	( zu(k) - zw(nzb_si) )**2 ) )
	1094	ENDDO
	1095
	1096	ENDIF
	1097
	1098	ENDDO
	1099	ENDDO
	1100
	1101	ENDIF
	1102
	1103	!
	1104	!-- Multiplication with wall_adjustment_factor
	1105	l_wall = wall_adjustment_factor * l_wall
	1106
	1107	!
	1108	!-- Need to set lateral boundary conditions for l_wall
[75]	1109	CALL exchange_horiz( l_wall )
[1]	1110
	1111	DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, &
	1112	nzb_tmp, vertical_influence, wall_l, wall_n, wall_r, wall_s )
	1113
	1114
	1115	END SUBROUTINE init_grid

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |