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

Last change on this file since 821 was 819, checked in by maronga, 13 years ago

last commit documented

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