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

Last change on this file since 1310 was 1310, checked in by raasch, 7 years ago

update of GPL copyright

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