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

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

preliminary update for changes concerning non-cyclic boundary conditions

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