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

Last change on this file since 927 was 927, checked in by raasch, 9 years ago

masking method can be switched on for mg-solver using inipar parameter masking_method

  • Property svn:keywords set to Id
File size: 53.7 KB
Line 
1 SUBROUTINE init_grid
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! Wall flags are not set for multigrid method in case of masking method
7!
8! Former revisions:
9! -----------------
10! $Id: init_grid.f90 927 2012-06-06 19:15:04Z raasch $
11!
12! 864 2012-03-27 15:10:33Z gryschka
13! In case of ocean and Dirichlet bottom bc for u and v dzu_mg and ddzu_pres
14! were not correctly defined 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!
995!--       In case of masking method, flags are not set and multigrid method
996!--       works like FFT-solver
997          IF ( .NOT. masking_method )  THEN
998
999             DO  i = nxl_l-1, nxr_l+1
1000                DO  j = nys_l-1, nyn_l+1
1001                   DO  k = nzb, nzt_l+1
1002                         
1003!
1004!--                   Inside/outside building (inside building does not need
1005!--                   further tests for walls)
1006                      IF ( k*inc <= nzb_local(j*inc,i*inc) )  THEN
1007
1008                         flags(k,j,i) = IBSET( flags(k,j,i), 6 )
1009
1010                      ELSE
1011!
1012!--                      Bottom wall
1013                         IF ( (k-1)*inc <= nzb_local(j*inc,i*inc) )  THEN
1014                            flags(k,j,i) = IBSET( flags(k,j,i), 0 )
1015                         ENDIF
1016!
1017!--                      South wall
1018                         IF ( k*inc <= nzb_local((j-1)*inc,i*inc) )  THEN
1019                            flags(k,j,i) = IBSET( flags(k,j,i), 2 )
1020                         ENDIF
1021!
1022!--                      North wall
1023                         IF ( k*inc <= nzb_local((j+1)*inc,i*inc) )  THEN
1024                            flags(k,j,i) = IBSET( flags(k,j,i), 3 )
1025                         ENDIF
1026!
1027!--                      Left wall
1028                         IF ( k*inc <= nzb_local(j*inc,(i-1)*inc) )  THEN
1029                            flags(k,j,i) = IBSET( flags(k,j,i), 4 )
1030                         ENDIF
1031!
1032!--                      Right wall
1033                         IF ( k*inc <= nzb_local(j*inc,(i+1)*inc) )  THEN
1034                            flags(k,j,i) = IBSET( flags(k,j,i), 5 )
1035                         ENDIF
1036
1037                      ENDIF
1038                           
1039                   ENDDO
1040                ENDDO
1041             ENDDO
1042
1043          ENDIF
1044
1045!
1046!--       Test output of flag arrays
1047!          i = nxl_l
1048!          WRITE (9,*)  ' '
1049!          WRITE (9,*)  '*** mg level ', l, ' ***', mg_switch_to_pe0_level
1050!          WRITE (9,*)  '    inc=', inc, '  i =', nxl_l
1051!          WRITE (9,*)  '    nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
1052!          DO  k = nzt_l+1, nzb, -1
1053!             WRITE (9,'(194(1X,I2))')  ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
1054!          ENDDO
1055
1056          inc = inc * 2
1057
1058       ENDDO
1059
1060    ENDIF
1061!
1062!-- Allocate flags needed for masking walls.
1063    ALLOCATE( wall_flags_0(nzb:nzt,nys:nyn,nxl:nxr) )
1064    wall_flags_0 = 0
1065
1066    IF ( scalar_advec == 'ws-scheme' )  THEN
1067!
1068!--    Set flags to steer the degradation of the advection scheme in advec_ws
1069!--    near topography, inflow- and outflow boundaries as well as bottom and
1070!--    top of model domain. wall_flags_0 remains zero for all non-prognostic
1071!--    grid points.
1072       DO  i = nxl, nxr
1073          DO  j = nys, nyn
1074             DO  k = nzb_s_inner(j,i)+1, nzt
1075!
1076!--             scalar - x-direction
1077!--             WS1 (0), WS3 (1), WS5 (2)
1078                IF ( k <= nzb_s_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) &
1079                     .OR. ( outflow_r .AND. i == nxr ) )  THEN
1080                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
1081                ELSEIF ( k <= nzb_s_inner(j,i+2) .OR. k <= nzb_s_inner(j,i-1)  &
1082                         .OR. ( outflow_r .AND. i == nxr-1 )                   &
1083                         .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
1084                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
1085                ELSE
1086                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 )
1087                ENDIF
1088!
1089!--             scalar - y-direction
1090!--             WS1 (3), WS3 (4), WS5 (5)
1091                IF ( k <= nzb_s_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) &
1092                     .OR. ( outflow_n .AND. j == nyn ) )  THEN
1093                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 )
1094!--             WS3
1095                ELSEIF ( k <= nzb_s_inner(j+2,i) .OR. k <= nzb_s_inner(j-1,i)  &
1096                        .OR. ( outflow_s .AND. j == nysv )                     &
1097                        .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
1098                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 )
1099!--             WS5
1100                ELSE
1101                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 )
1102                ENDIF
1103!
1104!--             scalar - z-direction
1105!--             WS1 (6), WS3 (7), WS5 (8)
1106                flag_set = .FALSE.
1107                IF ( k == nzb_s_inner(j,i) + 1 .OR. k == nzt )  THEN
1108                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 )
1109                   flag_set = .TRUE.
1110                ELSEIF ( k == nzb_s_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
1111                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 7 )
1112                   flag_set = .TRUE.
1113                ELSEIF ( k > nzb_s_inner(j,i) .AND. .NOT. flag_set )  THEN
1114                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 )
1115                ENDIF
1116             ENDDO
1117          ENDDO
1118       ENDDO
1119    ENDIF
1120
1121    IF ( momentum_advec == 'ws-scheme' )  THEN
1122!
1123!--    Set wall_flags_0 to steer the degradation of the advection scheme in advec_ws
1124!--    near topography, inflow- and outflow boundaries as well as bottom and
1125!--    top of model domain. wall_flags_0 remains zero for all non-prognostic
1126!--    grid points.
1127       DO  i = nxl, nxr
1128          DO  j = nys, nyn
1129             DO  k = nzb_u_inner(j,i)+1, nzt
1130!
1131!--             u component - x-direction
1132!--             WS1 (9), WS3 (10), WS5 (11)
1133                IF ( k <= nzb_u_inner(j,i+1)                                  &
1134                     .OR. ( outflow_r .AND. i == nxr ) )  THEN
1135                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
1136                ELSEIF ( k <= nzb_u_inner(j,i+2) .OR. k <= nzb_u_inner(j,i-1) &
1137                         .OR. ( outflow_r .AND. i == nxr-1 )                  &
1138                         .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
1139                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 )
1140                ELSE
1141                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 )
1142                ENDIF
1143!
1144!--             u component - y-direction
1145!--             WS1 (12), WS3 (13), WS5 (14)
1146                IF ( k <= nzb_u_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) &
1147                    .OR. ( outflow_n .AND. j == nyn ) )  THEN
1148                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
1149                ELSEIF ( k <= nzb_u_inner(j+2,i) .OR. k <= nzb_u_inner(j-1,i)  &
1150                    .OR. ( outflow_s .AND. j == nysv )                         &
1151                    .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
1152                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
1153                ELSE
1154                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 )
1155                ENDIF
1156!
1157!--             u component - z-direction
1158!--             WS1 (15), WS3 (16), WS5 (17)
1159                flag_set = .FALSE.
1160                IF ( k == nzb_u_inner(j,i) + 1 .OR. k == nzt )  THEN
1161                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 )
1162                   flag_set = .TRUE.
1163                ELSEIF ( k == nzb_u_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
1164                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 )
1165                   flag_set = .TRUE.
1166                ELSEIF ( k > nzb_u_inner(j,i) .AND. .NOT. flag_set )  THEN
1167                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 )
1168                ENDIF
1169
1170             ENDDO
1171          ENDDO
1172       ENDDO
1173
1174       DO  i = nxl, nxr
1175          DO  j = nys, nyn
1176             DO  k = nzb_v_inner(j,i)+1, nzt
1177!
1178!--             v component - x-direction
1179!--             WS1 (18), WS3 (19), WS5 (20)
1180                IF ( k <= nzb_v_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) &
1181                    .OR. ( outflow_r .AND. i == nxr ) )  THEN
1182                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
1183!--             WS3
1184                ELSEIF ( k <= nzb_v_inner(j,i+2) .OR. k <= nzb_v_inner(j,i-1)  &
1185                    .OR. ( outflow_r .AND. i == nxr-1 )                        &
1186                    .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
1187                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 )
1188                ELSE
1189                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 )
1190                ENDIF
1191!
1192!--             v component - y-direction
1193!--             WS1 (21), WS3 (22), WS5 (23)
1194                IF ( k <= nzb_v_inner(j+1,i)                                   &
1195                     .OR. ( outflow_n .AND. j == nyn ) )  THEN
1196                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
1197                ELSEIF ( k <= nzb_v_inner(j+2,i) .OR. k <= nzb_v_inner(j-1,i)  &
1198                         .OR. ( outflow_s .AND. j == nysv )                    &
1199                         .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
1200                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
1201                ELSE
1202                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
1203                ENDIF
1204!
1205!--             v component - z-direction
1206!--             WS1 (24), WS3 (25), WS5 (26)
1207                flag_set = .FALSE.
1208                IF ( k == nzb_v_inner(j,i) + 1 .OR. k == nzt )  THEN
1209                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 )
1210                   flag_set = .TRUE.
1211                ELSEIF ( k == nzb_v_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
1212                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 )
1213                   flag_set = .TRUE.
1214                ELSEIF ( k > nzb_v_inner(j,i) .AND. .NOT. flag_set )  THEN
1215                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 )
1216                ENDIF
1217
1218             ENDDO
1219          ENDDO
1220       ENDDO
1221       DO  i = nxl, nxr
1222          DO  j = nys, nyn
1223             DO  k = nzb_w_inner(j,i)+1, nzt
1224!
1225!--             w component - x-direction
1226!--             WS1 (27), WS3 (28), WS5 (29)
1227                IF ( k <= nzb_w_inner(j,i+1) .OR. ( outflow_l .AND. i == nxl ) &
1228                    .OR. ( outflow_r .AND. i == nxr ) )  THEN
1229                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
1230                ELSEIF ( k <= nzb_w_inner(j,i+2) .OR. k <= nzb_w_inner(j,i-1)  &
1231                         .OR. ( outflow_r .AND. i == nxr-1 )                   &
1232                         .OR. ( outflow_l .AND. i == nxlu  ) )  THEN
1233                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 )
1234                ELSE
1235                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i),29 )
1236                ENDIF
1237!
1238!--             w component - y-direction
1239!--             WS1 (30), WS3 (31), WS5 (32)
1240                IF ( k <= nzb_w_inner(j+1,i) .OR. ( outflow_s .AND. j == nys ) &
1241                     .OR. ( outflow_n .AND. j == nyn ) )  THEN
1242                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
1243                ELSEIF ( k <= nzb_w_inner(j+2,i) .OR. k <= nzb_w_inner(j-1,i)  &
1244                         .OR. ( outflow_s .AND. j == nysv )                    &
1245                         .OR. ( outflow_n .AND. j == nyn-1  ) )  THEN
1246                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 )
1247                ELSE
1248                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 32 )
1249                ENDIF
1250!
1251!--             w component - z-direction
1252!--             WS1 (33), WS3 (34), WS5 (35)
1253                flag_set = .FALSE.
1254                IF ( k == nzb_w_inner(j,i) .OR. k == nzb_w_inner(j,i) + 1      &
1255                                           .OR. k == nzt )  THEN
1256!
1257!--                Please note, at k == nzb_w_inner(j,i) a flag is explictely
1258!--                set, although this is not a prognostic level. However,
1259!--                contrary to the advection of u,v and s this is necessary
1260!--                because flux_t(nzb_w_inner(j,i)) is used for the tendency
1261!--                at k == nzb_w_inner(j,i)+1.
1262                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 33 )
1263                   flag_set = .TRUE.
1264                ELSEIF ( k == nzb_w_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
1265                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 34 )
1266                   flag_set = .TRUE.
1267                ELSEIF ( k > nzb_w_inner(j,i) .AND. .NOT. flag_set )  THEN
1268                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 35 )
1269                ENDIF
1270
1271             ENDDO
1272          ENDDO
1273       ENDDO
1274
1275    ENDIF
1276
1277!
1278!-- In case of topography: limit near-wall mixing length l_wall further:
1279!-- Go through all points of the subdomain one by one and look for the closest
1280!-- surface
1281    IF ( TRIM(topography) /= 'flat' )  THEN
1282       DO  i = nxl, nxr
1283          DO  j = nys, nyn
1284
1285             nzb_si = nzb_s_inner(j,i)
1286             vi     = vertical_influence(nzb_si)
1287
1288             IF ( wall_n(j,i) > 0 )  THEN
1289!
1290!--             North wall (y distance)
1291                DO  k = wall_n(j,i), nzb_si
1292                   l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5 * dy )
1293                ENDDO
1294!
1295!--             Above North wall (yz distance)
1296                DO  k = nzb_si + 1, nzb_si + vi
1297                   l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i),     &
1298                                          SQRT( 0.25 * dy**2 + &
1299                                          ( zu(k) - zw(nzb_si) )**2 ) )
1300                ENDDO
1301!
1302!--             Northleft corner (xy distance)
1303                IF ( corner_nl(j,i) > 0 )  THEN
1304                   DO  k = corner_nl(j,i), nzb_si
1305                      l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), &
1306                                               0.5 * SQRT( dx**2 + dy**2 ) )
1307                   ENDDO
1308!
1309!--                Above Northleft corner (xyz distance)
1310                   DO  k = nzb_si + 1, nzb_si + vi
1311                      l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1),             &
1312                                               SQRT( 0.25 * (dx**2 + dy**2) + &
1313                                               ( zu(k) - zw(nzb_si) )**2 ) )
1314                   ENDDO
1315                ENDIF
1316!
1317!--             Northright corner (xy distance)
1318                IF ( corner_nr(j,i) > 0 )  THEN
1319                   DO  k = corner_nr(j,i), nzb_si
1320                       l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
1321                                                0.5 * SQRT( dx**2 + dy**2 ) )
1322                   ENDDO
1323!
1324!--                Above northright corner (xyz distance)
1325                   DO  k = nzb_si + 1, nzb_si + vi
1326                      l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), &
1327                                               SQRT( 0.25 * (dx**2 + dy**2) + &
1328                                               ( zu(k) - zw(nzb_si) )**2 ) )
1329                   ENDDO
1330                ENDIF
1331             ENDIF
1332
1333             IF ( wall_s(j,i) > 0 )  THEN
1334!
1335!--             South wall (y distance)
1336                DO  k = wall_s(j,i), nzb_si
1337                   l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), 0.5 * dy )
1338                ENDDO
1339!
1340!--             Above south wall (yz distance)
1341                DO  k = nzb_si + 1, &
1342                        nzb_si + vi
1343                   l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i),     &
1344                                          SQRT( 0.25 * dy**2 + &
1345                                          ( zu(k) - zw(nzb_si) )**2 ) )
1346                ENDDO
1347!
1348!--             Southleft corner (xy distance)
1349                IF ( corner_sl(j,i) > 0 )  THEN
1350                   DO  k = corner_sl(j,i), nzb_si
1351                      l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &
1352                                               0.5 * SQRT( dx**2 + dy**2 ) )
1353                   ENDDO
1354!
1355!--                Above southleft corner (xyz distance)
1356                   DO  k = nzb_si + 1, nzb_si + vi
1357                      l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1),             &
1358                                               SQRT( 0.25 * (dx**2 + dy**2) + &
1359                                               ( zu(k) - zw(nzb_si) )**2 ) )
1360                   ENDDO
1361                ENDIF
1362!
1363!--             Southright corner (xy distance)
1364                IF ( corner_sr(j,i) > 0 )  THEN
1365                   DO  k = corner_sr(j,i), nzb_si
1366                      l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &
1367                                               0.5 * SQRT( dx**2 + dy**2 ) )
1368                   ENDDO
1369!
1370!--                Above southright corner (xyz distance)
1371                   DO  k = nzb_si + 1, nzb_si + vi
1372                      l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1),             &
1373                                               SQRT( 0.25 * (dx**2 + dy**2) + &
1374                                               ( zu(k) - zw(nzb_si) )**2 ) )
1375                   ENDDO
1376                ENDIF
1377
1378             ENDIF
1379
1380             IF ( wall_l(j,i) > 0 )  THEN
1381!
1382!--             Left wall (x distance)
1383                DO  k = wall_l(j,i), nzb_si
1384                   l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), 0.5 * dx )
1385                ENDDO
1386!
1387!--             Above left wall (xz distance)
1388                DO  k = nzb_si + 1, nzb_si + vi
1389                   l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1),     &
1390                                          SQRT( 0.25 * dx**2 + &
1391                                          ( zu(k) - zw(nzb_si) )**2 ) )
1392                ENDDO
1393             ENDIF
1394
1395             IF ( wall_r(j,i) > 0 )  THEN
1396!
1397!--             Right wall (x distance)
1398                DO  k = wall_r(j,i), nzb_si
1399                   l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5 * dx )
1400                ENDDO
1401!
1402!--             Above right wall (xz distance)
1403                DO  k = nzb_si + 1, nzb_si + vi
1404                   l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1),     &
1405                                          SQRT( 0.25 * dx**2 + &
1406                                          ( zu(k) - zw(nzb_si) )**2 ) )
1407                ENDDO
1408
1409             ENDIF
1410
1411          ENDDO
1412       ENDDO
1413
1414    ENDIF
1415
1416!
1417!-- Multiplication with wall_adjustment_factor
1418    l_wall = wall_adjustment_factor * l_wall
1419
1420!
1421!-- Set lateral boundary conditions for l_wall
1422    CALL exchange_horiz( l_wall, nbgp )
1423
1424    DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, &
1425                nzb_tmp, vertical_influence, wall_l, wall_n, wall_r, wall_s )
1426
1427#endif
1428
1429 END SUBROUTINE init_grid
Note: See TracBrowser for help on using the repository browser.