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

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

New:


openACC porting of reduction operations
additional 3D-flag arrays for replacing the 2D-index arrays nzb_s_inner and nzb_diff_s_inner
(flow_statistics, init_grid, init_3d_model, modules, palm, pres, time_integration)

Changed:


for PGI/openACC performance reasons set default compile options for openACC to "-ta=nocache",
and set environment variable PGI_ACC_SYNCHRONOUS=1
(MAKE.inc.pgi.openacc, palm_simple_run)

wall_flags_0 changed to 32bit INTEGER, additional array wall_flags_00 introduced to hold
bits 32-63
(advec_ws, init_grid, modules, palm)

Errors:


dummy argument tri in 1d-routines replaced by tri_for_1d because of name
conflict with arry tri in module arrays_3d
(tridia_solver)

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