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

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

last commit documented

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