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

Last change on this file since 686 was 668, checked in by suehring, 13 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 41.2 KB
Line 
1 SUBROUTINE init_grid
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7! Former revisions:
8! -----------------
9! $Id: init_grid.f90 668 2010-12-23 13:22:58Z gryschka $
10!
11! 667 2010-12-23 12:06:00Z suehring/gryschka
12! Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
13! Furthermore the allocation of arrays and steering of loops is done with these
14! parameters. Call of exchange_horiz are modified.
15! In case of dirichlet bounday condition at the bottom zu(0)=0.0
16! dzu_mg has to be set explicitly for a equally spaced grid near bottom.
17! ddzu_pres added to use a equally spaced grid near bottom.
18!
19! 555 2010-09-07 07:32:53Z raasch
20! Bugfix: default setting of nzb_local for flat topographie
21!
22! 274 2009-03-26 15:11:21Z heinze
23! Output of messages replaced by message handling routine.
24! new topography case 'single_street_canyon'
25!
26! 217 2008-12-09 18:00:48Z letzel
27! +topography_grid_convention
28!
29! 134 2007-11-21 07:28:38Z letzel
30! Redefine initial nzb_local as the actual total size of topography (later the
31! extent of topography in nzb_local is reduced by 1dx at the E topography walls
32! and by 1dy at the N topography walls to form the basis for nzb_s_inner);
33! for consistency redefine 'single_building' case.
34! Calculation of wall flag arrays
35!
36! 94 2007-06-01 15:25:22Z raasch
37! Grid definition for ocean version
38!
39! 75 2007-03-22 09:54:05Z raasch
40! storage of topography height arrays zu_s_inner and zw_s_inner,
41! 2nd+3rd argument removed from exchange horiz
42!
43! 19 2007-02-23 04:53:48Z raasch
44! Setting of nzt_diff
45!
46! RCS Log replace by Id keyword, revision history cleaned up
47!
48! Revision 1.17  2006/08/22 14:00:05  raasch
49! +dz_max to limit vertical stretching,
50! bugfix in index array initialization for line- or point-like topography
51! structures
52!
53! Revision 1.1  1997/08/11 06:17:45  raasch
54! Initial revision (Testversion)
55!
56!
57! Description:
58! ------------
59! Creating grid depending constants
60!------------------------------------------------------------------------------!
61
62    USE arrays_3d
63    USE control_parameters
64    USE grid_variables
65    USE indices
66    USE pegrid
67
68    IMPLICIT NONE
69
70    INTEGER ::  bh, blx, bly, bxl, bxr, byn, bys, ch, cwx, cwy, cxl, cxr, cyn, &
71                cys, gls, i, inc, i_center, j, j_center, k, l, nxl_l, nxr_l, &
72                nyn_l, nys_l, nzb_si, nzt_l, vi
73
74    INTEGER, DIMENSION(:), ALLOCATABLE   ::  vertical_influence
75
76    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  corner_nl, corner_nr, corner_sl,  &
77                                             corner_sr, wall_l, wall_n, wall_r,&
78                                             wall_s, nzb_local, nzb_tmp
79
80    REAL    ::  dx_l, dy_l, dz_stretched
81
82    REAL, DIMENSION(0:ny,0:nx)          ::  topo_height
83
84    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  distance
85   
86!
87!   Computation of the array bounds.
88    nxlg = nxl - nbgp
89    nxrg = nxr + nbgp
90    nysg = nys - nbgp
91    nyng = nyn + nbgp
92!
93!-- Allocate grid arrays
94    ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &
95              dzw(1:nzt+1), l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) )
96
97!
98!-- Compute height of u-levels from constant grid length and dz stretch factors
99    IF ( dz == -1.0 )  THEN
100       message_string = 'missing dz'
101       CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 ) 
102    ELSEIF ( dz <= 0.0 )  THEN
103       WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
104       CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
105    ENDIF
106
107!
108!-- Define the vertical grid levels
109    IF ( .NOT. ocean )  THEN
110!
111!--    Grid for atmosphere with surface at z=0 (k=0, w-grid).
112!--    The first u-level above the surface corresponds to the top of the
113!--    Prandtl-layer.
114
115       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
116          zu(0) = 0.0
117      !    zu(0) = - dz * 0.5
118       ELSE
119          zu(0) = - dz * 0.5
120       ENDIF
121       zu(1) =   dz * 0.5
122
123       dz_stretch_level_index = nzt+1
124       dz_stretched = dz
125       DO  k = 2, nzt+1
126          IF ( dz_stretch_level <= zu(k-1)  .AND.  dz_stretched < dz_max )  THEN
127             dz_stretched = dz_stretched * dz_stretch_factor
128             dz_stretched = MIN( dz_stretched, dz_max )
129             IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1
130          ENDIF
131          zu(k) = zu(k-1) + dz_stretched
132       ENDDO
133
134!
135!--    Compute the w-levels. They are always staggered half-way between the
136!--    corresponding u-levels. The top w-level is extrapolated linearly.
137       zw(0) = 0.0
138       DO  k = 1, nzt
139          zw(k) = ( zu(k) + zu(k+1) ) * 0.5
140       ENDDO
141       zw(nzt+1) = zw(nzt) + 2.0 * ( zu(nzt+1) - zw(nzt) )
142
143    ELSE
144!
145!--    Grid for ocean with solid surface at z=0 (k=0, w-grid). The free water
146!--    surface is at k=nzt (w-grid).
147!--    Since the w-level lies always on the surface, the first/last u-level
148!--    (staggered!) lies below the bottom surface / above the free surface.
149!--    It is used for "mirror" boundary condition.
150!--    The first u-level above the bottom surface corresponds to the top of the
151!--    Prandtl-layer.
152       zu(nzt+1) =   dz * 0.5
153       zu(nzt)   = - dz * 0.5
154
155       dz_stretch_level_index = 0
156       dz_stretched = dz
157       DO  k = nzt-1, 0, -1
158          IF ( dz_stretch_level <= ABS( zu(k+1) )  .AND.  &
159               dz_stretched < dz_max )  THEN
160             dz_stretched = dz_stretched * dz_stretch_factor
161             dz_stretched = MIN( dz_stretched, dz_max )
162             IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1
163          ENDIF
164          zu(k) = zu(k+1) - dz_stretched
165       ENDDO
166
167!
168!--    Compute the w-levels. They are always staggered half-way between the
169!--    corresponding u-levels.
170!--    The top w-level (nzt+1) is not used but set for consistency, since
171!--    w and all scalar variables are defined up tp nzt+1.
172       zw(nzt+1) = dz
173       zw(nzt)   = 0.0
174       DO  k = 0, nzt
175          zw(k) = ( zu(k) + zu(k+1) ) * 0.5
176       ENDDO
177
178    ENDIF
179
180!
181!-- Compute grid lengths.
182    DO  k = 1, nzt+1
183       dzu(k)  = zu(k) - zu(k-1)
184       ddzu(k) = 1.0 / dzu(k)
185       dzw(k)  = zw(k) - zw(k-1)
186       ddzw(k) = 1.0 / dzw(k)
187    ENDDO
188
189    DO  k = 1, nzt
190       dd2zu(k) = 1.0 / ( dzu(k) + dzu(k+1) )
191    ENDDO
192   
193!   
194!-- In case of FFT method or SOR swap inverse grid lenght ddzu to ddzu_fft and
195!-- modify the lowest entry. This is necessary for atmosphere runs, because
196!-- zu(0) and so ddzu(1) changed. Accompanied with this modified arrays
197!-- pressure solver uses wrong values and this causes kinks in the profiles
198!-- of turbulent quantities. 
199    IF ( psolver /= 'multigrid' ) THEN
200       ALLOCATE( ddzu_pres(1:nzt+1) )
201       ddzu_pres = ddzu
202       IF( .NOT. ocean ) ddzu_pres(1) = ddzu_pres(2)
203    ENDIF   
204
205!
206!-- In case of multigrid method, compute grid lengths and grid factors for the
207!-- grid levels
208    IF ( psolver == 'multigrid' )  THEN
209
210       ALLOCATE( ddx2_mg(maximum_grid_level), ddy2_mg(maximum_grid_level), &
211                 dzu_mg(nzb+1:nzt+1,maximum_grid_level),                   &
212                 dzw_mg(nzb+1:nzt+1,maximum_grid_level),                   &
213                 f1_mg(nzb+1:nzt,maximum_grid_level),                      &
214                 f2_mg(nzb+1:nzt,maximum_grid_level),                      &
215                 f3_mg(nzb+1:nzt,maximum_grid_level) )
216
217       dzu_mg(:,maximum_grid_level) = dzu
218!       
219!--    To ensure a equally spaced grid. For ocean runs this is not necessary,
220!--    because zu(0) does not changed so far. Also this would cause errors
221!--    if a vertical stretching for ocean runs is used.
222       IF ( .NOT. ocean ) dzu_mg(1,maximum_grid_level) = dzu(2)
223       dzw_mg(:,maximum_grid_level) = dzw
224       nzt_l = nzt
225       DO  l = maximum_grid_level-1, 1, -1
226           dzu_mg(nzb+1,l) = 2.0 * dzu_mg(nzb+1,l+1)
227           dzw_mg(nzb+1,l) = 2.0 * dzw_mg(nzb+1,l+1)
228           nzt_l = nzt_l / 2
229           DO  k = 2, nzt_l+1
230              dzu_mg(k,l) = dzu_mg(2*k-2,l+1) + dzu_mg(2*k-1,l+1)
231              dzw_mg(k,l) = dzw_mg(2*k-2,l+1) + dzw_mg(2*k-1,l+1)
232           ENDDO
233       ENDDO
234
235       nzt_l = nzt
236       dx_l  = dx
237       dy_l  = dy
238       DO  l = maximum_grid_level, 1, -1
239          ddx2_mg(l) = 1.0 / dx_l**2
240          ddy2_mg(l) = 1.0 / dy_l**2
241          DO  k = nzb+1, nzt_l
242             f2_mg(k,l) = 1.0 / ( dzu_mg(k+1,l) * dzw_mg(k,l) )
243             f3_mg(k,l) = 1.0 / ( dzu_mg(k,l)   * dzw_mg(k,l) )
244             f1_mg(k,l) = 2.0 * ( ddx2_mg(l) + ddy2_mg(l) ) + &
245                          f2_mg(k,l) + f3_mg(k,l)
246          ENDDO
247          nzt_l = nzt_l / 2
248          dx_l  = dx_l * 2.0
249          dy_l  = dy_l * 2.0
250       ENDDO
251
252    ENDIF
253
254!
255!-- Compute the reciprocal values of the horizontal grid lengths.
256    ddx = 1.0 / dx
257    ddy = 1.0 / dy
258    dx2 = dx * dx
259    dy2 = dy * dy
260    ddx2 = 1.0 / dx2
261    ddy2 = 1.0 / dy2
262
263!
264!-- Compute the grid-dependent mixing length.
265    DO  k = 1, nzt
266       l_grid(k)  = ( dx * dy * dzw(k) )**0.33333333333333
267    ENDDO
268
269!
270!-- Allocate outer and inner index arrays for topography and set
271!-- defaults.
272!-- nzb_local has to contain additional layers of ghost points for calculating
273!-- the flag arrays needed for the multigrid method
274    gls = 2**( maximum_grid_level )
275
276    ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr),       &
277              corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr),       &
278              nzb_local(-gls:ny+gls,-gls:nx+gls),                                   &
279              nzb_tmp(-nbgp:ny+nbgp,-nbgp:nx+nbgp),                         &
280              wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr),             &
281              wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
282    ALLOCATE( fwxm(nysg:nyng,nxlg:nxrg), fwxp(nysg:nyng,nxlg:nxrg),         &
283              fwym(nysg:nyng,nxlg:nxrg), fwyp(nysg:nyng,nxlg:nxrg),         &
284              fxm(nysg:nyng,nxlg:nxrg), fxp(nysg:nyng,nxlg:nxrg),           &
285              fym(nysg:nyng,nxlg:nxrg), fyp(nysg:nyng,nxlg:nxrg),           &
286              nzb_s_inner(nysg:nyng,nxlg:nxrg),                             &
287              nzb_s_outer(nysg:nyng,nxlg:nxrg),                             &
288              nzb_u_inner(nysg:nyng,nxlg:nxrg),                             &
289              nzb_u_outer(nysg:nyng,nxlg:nxrg),                             &
290              nzb_v_inner(nysg:nyng,nxlg:nxrg),                             &
291              nzb_v_outer(nysg:nyng,nxlg:nxrg),                             &
292              nzb_w_inner(nysg:nyng,nxlg:nxrg),                             &
293              nzb_w_outer(nysg:nyng,nxlg:nxrg),                             &
294              nzb_diff_s_inner(nysg:nyng,nxlg:nxrg),                        &
295              nzb_diff_s_outer(nysg:nyng,nxlg:nxrg),                        &
296              nzb_diff_u(nysg:nyng,nxlg:nxrg),                              &
297              nzb_diff_v(nysg:nyng,nxlg:nxrg),                              &
298              nzb_2d(nysg:nyng,nxlg:nxrg),                                  &
299              wall_e_x(nysg:nyng,nxlg:nxrg),                                &
300              wall_e_y(nysg:nyng,nxlg:nxrg),                                &
301              wall_u(nysg:nyng,nxlg:nxrg),                                  &
302              wall_v(nysg:nyng,nxlg:nxrg),                                  &
303              wall_w_x(nysg:nyng,nxlg:nxrg),                                &
304              wall_w_y(nysg:nyng,nxlg:nxrg) )
305
306
307
308    ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
309
310    nzb_s_inner = nzb;  nzb_s_outer = nzb
311    nzb_u_inner = nzb;  nzb_u_outer = nzb
312    nzb_v_inner = nzb;  nzb_v_outer = nzb
313    nzb_w_inner = nzb;  nzb_w_outer = nzb
314
315!
316!-- Define vertical gridpoint from (or to) which on the usual finite difference
317!-- form (which does not use surface fluxes) is applied
318    IF ( prandtl_layer  .OR.  use_surface_fluxes )  THEN
319       nzb_diff = nzb + 2
320    ELSE
321       nzb_diff = nzb + 1
322    ENDIF
323    IF ( use_top_fluxes )  THEN
324       nzt_diff = nzt - 1
325    ELSE
326       nzt_diff = nzt
327    ENDIF
328
329    nzb_diff_s_inner = nzb_diff;  nzb_diff_s_outer = nzb_diff
330    nzb_diff_u = nzb_diff;  nzb_diff_v = nzb_diff
331
332    wall_e_x = 0.0;  wall_e_y = 0.0;  wall_u = 0.0;  wall_v = 0.0
333    wall_w_x = 0.0;  wall_w_y = 0.0
334    fwxp = 1.0;  fwxm = 1.0;  fwyp = 1.0;  fwym = 1.0
335    fxp  = 1.0;  fxm  = 1.0;  fyp  = 1.0;  fym  = 1.0
336
337!
338!-- Initialize near-wall mixing length l_wall only in the vertical direction
339!-- for the moment,
340!-- multiplication with wall_adjustment_factor near the end of this routine
341    l_wall(nzb,:,:)   = l_grid(1)
342    DO  k = nzb+1, nzt
343       l_wall(k,:,:)  = l_grid(k)
344    ENDDO
345    l_wall(nzt+1,:,:) = l_grid(nzt)
346
347    ALLOCATE ( vertical_influence(nzb:nzt) )
348    DO  k = 1, nzt
349       vertical_influence(k) = MIN ( INT( l_grid(k) / &
350                     ( wall_adjustment_factor * dzw(k) ) + 0.5 ), nzt - k )
351    ENDDO
352
353    DO  k = 1, MAXVAL( nzb_s_inner )
354       IF ( l_grid(k) > 1.5 * dx * wall_adjustment_factor .OR.  &
355            l_grid(k) > 1.5 * dy * wall_adjustment_factor )  THEN
356          WRITE( message_string, * ) 'grid anisotropy exceeds ', &
357                                     'threshold given by only local', &
358                                     ' &horizontal reduction of near_wall ', &
359                                     'mixing length l_wall', &
360                                     ' &starting from height level k = ', k, '.'
361          CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
362          EXIT
363       ENDIF
364    ENDDO
365    vertical_influence(0) = vertical_influence(1)
366
367    DO  i = nxlg, nxrg
368       DO  j = nysg, nyng
369          DO  k = nzb_s_inner(j,i) + 1, &
370                  nzb_s_inner(j,i) + vertical_influence(nzb_s_inner(j,i))
371             l_wall(k,j,i) = zu(k) - zw(nzb_s_inner(j,i))
372          ENDDO
373       ENDDO
374    ENDDO
375
376!
377!-- Set outer and inner index arrays for non-flat topography.
378!-- Here consistency checks concerning domain size and periodicity are
379!-- necessary.
380!-- Within this SELECT CASE structure only nzb_local is initialized
381!-- individually depending on the chosen topography type, all other index
382!-- arrays are initialized further below.
383    SELECT CASE ( TRIM( topography ) )
384
385       CASE ( 'flat' )
386!
387!--       nzb_local is required for the multigrid solver
388          nzb_local = 0
389
390       CASE ( 'single_building' )
391!
392!--       Single rectangular building, by default centered in the middle of the
393!--       total domain
394          blx = NINT( building_length_x / dx )
395          bly = NINT( building_length_y / dy )
396          bh  = NINT( building_height / dz )
397
398          IF ( building_wall_left == 9999999.9 )  THEN
399             building_wall_left = ( nx + 1 - blx ) / 2 * dx
400          ENDIF
401          bxl = NINT( building_wall_left / dx )
402          bxr = bxl + blx
403
404          IF ( building_wall_south == 9999999.9 )  THEN
405             building_wall_south = ( ny + 1 - bly ) / 2 * dy
406          ENDIF
407          bys = NINT( building_wall_south / dy )
408          byn = bys + bly
409
410!
411!--       Building size has to meet some requirements
412          IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR.  &
413               ( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) )  THEN
414             WRITE( message_string, * ) 'inconsistent building parameters:',   &
415                                      '& bxl=', bxl, 'bxr=', bxr, 'bys=', bys, &
416                                      'byn=', byn, 'nx=', nx, 'ny=', ny
417             CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
418          ENDIF
419
420!
421!--       Define the building.
422          nzb_local = 0
423          nzb_local(bys:byn,bxl:bxr) = bh
424
425       CASE ( 'single_street_canyon' )
426!
427!--       Single quasi-2D street canyon of infinite length in x or y direction.
428!--       The canyon is centered in the other direction by default.
429          IF ( canyon_width_x /= 9999999.9 )  THEN
430!
431!--          Street canyon in y direction
432             cwx = NINT( canyon_width_x / dx )
433             IF ( canyon_wall_left == 9999999.9 )  THEN
434                canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
435             ENDIF
436             cxl = NINT( canyon_wall_left / dx )
437             cxr = cxl + cwx
438
439          ELSEIF ( canyon_width_y /= 9999999.9 )  THEN
440!
441!--          Street canyon in x direction
442             cwy = NINT( canyon_width_y / dy )
443             IF ( canyon_wall_south == 9999999.9 )  THEN
444                canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
445             ENDIF
446             cys = NINT( canyon_wall_south / dy )
447             cyn = cys + cwy
448
449          ELSE
450             
451             message_string = 'no street canyon width given'
452             CALL message( 'init_grid', 'PA0204', 1, 2, 0, 6, 0 )
453 
454          ENDIF
455
456          ch             = NINT( canyon_height / dz )
457          dp_level_ind_b = ch
458!
459!--       Street canyon size has to meet some requirements
460          IF ( canyon_width_x /= 9999999.9 )  THEN
461             IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR.  &
462               ( ch < 3 ) )  THEN
463                WRITE( message_string, * ) 'inconsistent canyon parameters:', &
464                                           '&cxl=', cxl, 'cxr=', cxr,         &
465                                           'cwx=', cwx,                       &
466                                           'ch=', ch, 'nx=', nx, 'ny=', ny
467                CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 ) 
468             ENDIF
469          ELSEIF ( canyon_width_y /= 9999999.9 )  THEN
470             IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR.  &
471               ( ch < 3 ) )  THEN
472                WRITE( message_string, * ) 'inconsistent canyon parameters:', &
473                                           '&cys=', cys, 'cyn=', cyn,         &
474                                           'cwy=', cwy,                       &
475                                           'ch=', ch, 'nx=', nx, 'ny=', ny
476                CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 ) 
477             ENDIF
478          ENDIF
479          IF ( canyon_width_x /= 9999999.9 .AND. canyon_width_y /= 9999999.9 ) &
480               THEN
481             message_string = 'inconsistent canyon parameters:' //     & 
482                              '&street canyon can only be oriented' // &
483                              '&either in x- or in y-direction'
484             CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 )
485          ENDIF
486
487          nzb_local = ch
488          IF ( canyon_width_x /= 9999999.9 )  THEN
489             nzb_local(:,cxl+1:cxr-1) = 0
490          ELSEIF ( canyon_width_y /= 9999999.9 )  THEN
491             nzb_local(cys+1:cyn-1,:) = 0
492          ENDIF
493
494       CASE ( 'read_from_file' )
495!
496!--       Arbitrary irregular topography data in PALM format (exactly matching
497!--       the grid size and total domain size)
498          OPEN( 90, FILE='TOPOGRAPHY_DATA', STATUS='OLD', FORM='FORMATTED',  &
499               ERR=10 )
500          DO  j = ny, 0, -1
501             READ( 90, *, ERR=11, END=11 )  ( topo_height(j,i), i = 0, nx )
502          ENDDO
503!
504!--       Calculate the index height of the topography
505          DO  i = 0, nx
506             DO  j = 0, ny
507                nzb_local(j,i) = NINT( topo_height(j,i) / dz )
508             ENDDO
509          ENDDO
510!
511!--       Add cyclic boundaries (additional layers are for calculating flag
512!--       arrays needed for the multigrid sover)
513          nzb_local(-gls:-1,0:nx)     = nzb_local(ny-gls+1:ny,0:nx)
514          nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,0:nx)
515          nzb_local(:,-gls:-1)        = nzb_local(:,nx-gls+1:nx)
516          nzb_local(:,nx+1:nx+gls)    = nzb_local(:,0:gls-1)
517
518
519     
520          GOTO 12
521         
522 10       message_string = 'file TOPOGRAPHY_DATA does not exist'
523          CALL message( 'init_grid', 'PA0208', 1, 2, 0, 6, 0 )
524
525 11       message_string = 'errors in file TOPOGRAPHY_DATA'
526          CALL message( 'init_grid', 'PA0209', 1, 2, 0, 6, 0 )
527
528 12       CLOSE( 90 )
529
530       CASE DEFAULT
531!
532!--       The DEFAULT case is reached either if the parameter topography
533!--       contains a wrong character string or if the user has defined a special
534!--       case in the user interface. There, the subroutine user_init_grid
535!--       checks which of these two conditions applies.
536          CALL user_init_grid( gls, nzb_local )
537
538    END SELECT
539
540!
541!-- Test output of nzb_local -1:ny+1,-1:nx+1
542!    WRITE (9,*)  '*** nzb_local ***'
543!    DO  j = ny+1, -1, -1
544!       WRITE (9,'(194(1X,I2))')  ( nzb_local(j,i), i = -1, nx+1 )
545!    ENDDO
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 == 'cyclic' )  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 == 'cyclic' )  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!-- Need to set lateral boundary conditions for l_wall
1138
1139    CALL exchange_horiz( l_wall, nbgp )
1140
1141    DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, &
1142                nzb_tmp, vertical_influence, wall_l, wall_n, wall_r, wall_s )
1143
1144
1145 END SUBROUTINE init_grid
Note: See TracBrowser for help on using the repository browser.