Home

Context Navigation

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

Last change on this file since 927 was 927, checked in by raasch, 12 years ago
masking method can be switched on for mg-solver using inipar parameter masking_method
Property svn:keywords set to `Id`
File size: 53.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |