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

Last change on this file since 843 was 843, checked in by gryschka, 12 years ago

Redefiend first u-level (k=0) in ocean in case of bc_uv_b='dirichlet'

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