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

Last change on this file since 810 was 810, checked in by maronga, 10 years ago

last commit documented

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