Home

Context Navigation

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

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 55.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |