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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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