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

Last change on this file since 1418 was 1418, checked in by fricke, 10 years ago

Bugfixes concerning grid stretching for the ocean and calculation of the salinity flux in routine surface_coupler

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