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

Last change on this file since 2021 was 2021, checked in by suehring, 5 years ago

Bugfix, restore flags for nest boundaries in multigrid solver; bugfix: setting Neumann boundary conditions for topography arrays in case of non-cyclic boundary conditions

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