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

Last change on this file since 1069 was 1069, checked in by maronga, 11 years ago

allow usage of topography in combination with cloud physics, allow usage of topography in coupled ocean model, minor changes in mbuild and mrun

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