Home

Context Navigation

init_grid.f90 @ 707

Last change on this file since 707 was 707, checked in by raasch, 13 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:

Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:

localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

Property svn:keywords set to Id

File size: 41.3 KB

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

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

Context Navigation

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

Download in other formats: