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

Last change on this file since 1017 was 1017, checked in by raasch, 11 years ago

last commit documented

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