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

Last change on this file since 1678 was 1678, checked in by boeske, 6 years ago

last commit documented

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