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

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

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:


Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:


localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

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