Home

Context Navigation

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

Last change on this file since 786 was 760, checked in by raasch, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 41.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |