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

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

formatting adjustments

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