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

Last change on this file since 1001 was 997, checked in by raasch, 12 years ago

last commit documented

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