source: palm/tags/release-3.8/SOURCE/init_grid.f90 @ 4343

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

bugfixes of r722 merged to release 3.8

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