Home

Context Navigation

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

Last change on this file since 1001 was 997, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 55.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |