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

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

last commit documented

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