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

Last change on this file since 710 was 710, checked in by raasch, 10 years ago

last commit documented

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