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

Last change on this file since 1070 was 1070, checked in by maronga, 9 years ago

last commit documented

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