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

Last change on this file since 759 was 759, checked in by raasch, 10 years ago

New:
---

The number of parallel I/O operations can be limited with new mrun-option -w.
(advec_particles, data_output_2d, data_output_3d, header, init_grid, init_pegrid, init_3d_model, modules, palm, parin, write_3d_binary)

Changed:


mrun option -T is obligatory

Errors:


Bugfix: No zero assignments to volume_flow_initial and volume_flow_area in
case of normal restart runs. (init_3d_model)

initialization of u_0, v_0. This is just to avoid access of uninitialized
memory in exchange_horiz_2d, which causes respective error messages
when the Intel thread checker (inspector) is used. (production_e)

Bugfix for ts limitation (prandtl_fluxes)

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