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

Last change on this file since 77 was 77, checked in by raasch, 14 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

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