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

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

last commit documented

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