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

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

new utility check_namelist_files implemented

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