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

Last change on this file since 928 was 928, checked in by raasch, 12 years ago

last commit documented

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