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

Last change on this file since 667 was 667, checked in by suehring, 11 years ago

summary:


Gryschka:

  • Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
  • Exchange of u and v from ocean surface to atmosphere surface
  • Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
  • Inflow turbulence is now defined by flucuations around spanwise mean
  • Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

  • New advection added ( Wicker and Skamarock 5th order ), therefore:
    • New module advec_ws.f90
    • Modified exchange of ghost boundaries.
    • Modified evaluation of turbulent fluxes
    • New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90


Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90


Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90


Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90


When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90


Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90


Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90


Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90


Steering parameter dissipation_control added in inipar.

Makefile


Module advec_ws added.

Modules


Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90


New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90


Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90


Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90


Output of total array size was adapted to nbgp.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90


Calls of exchange_horiz are modified.

diffusion_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90


Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90


Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90


determination of target_id's moved to init_pegrid

init_pt_anomaly.f90


Call of exchange_horiz are modified.

init_rankine.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90


Output of advection scheme.

poismg.f90


Calls of exchange_horiz are modified.

prandtl_fluxes.f90


Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90


+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90


Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90


Allocation with nbgp.

wall_fluxes.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90


Array bounds and nx, ny adapted with nbgp.

sor.f90


bug removed in declaration of ddzw(), nz replaced by nzt+1

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