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

Last change on this file since 808 was 808, checked in by maronga, 10 years ago

last commit documented

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