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

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

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

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