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

Last change on this file since 861 was 861, checked in by suehring, 12 years ago

WS5 is available in combination with topography. Version number changed from 3.8 to 3.8a. Modification in init_pt_anomaly.

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