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

Last change on this file since 863 was 863, checked in by suehring, 9 years ago

last commit documented

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