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

Last change on this file since 818 was 818, checked in by maronga, 9 years ago

bugfix: namelist file check now possible for topography and re-enabled. mrungui update (-z option)

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