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

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

bugfixes for r707 concerning multigrid method for non-cyclic boundary conditions

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