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

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

last commit documented

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