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

Last change on this file since 865 was 865, checked in by gryschka, 12 years ago

last commit documented

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