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

Last change on this file since 2012 was 2001, checked in by knoop, 8 years ago

last commit documented

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