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

Last change on this file since 1037 was 1037, checked in by raasch, 11 years ago

last commit documented

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