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

Last change on this file since 121 was 116, checked in by raasch, 17 years ago

further preliminary updates concerning particle sorting and documentation

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