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

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

unused variables remove from several routines

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