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

Last change on this file since 979 was 979, checked in by fricke, 9 years ago

last commit documented

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