Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Dec 23, 2010 12:06:00 PM (13 years ago)
Author:
suehring
Message:

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

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE

    • Property svn:mergeinfo set to (toggle deleted branches)
      /palm/branches/suehring423-666
      /palm/branches/letzel/masked_output/SOURCE296-409

  • palm/trunk/SOURCE/init_3d_model.f90

    r631 r667  
    77! Current revisions:
    88! -----------------
    9 ! Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)
     9! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
     10! allocation of arrays. Calls of exchange_horiz are modified.
     11! Call ws_init to initialize arrays needed for statistical evaluation and
     12! optimization when ws-scheme is used.
     13!
     14!
     15! Initial volume flow is now calculated by using the variable hom_sum.
     16! Therefore the correction of initial volume flow for non-flat topography
     17! removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
     18!
     19! Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
     20! mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
     21! representative for the height z0
     22!
     23! Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)
    1024!
    1125! Former revisions:
    … …  
    101115!------------------------------------------------------------------------------!
    102116
     117    USE advec_ws
    103118    USE arrays_3d
    104119    USE averaging
    … …  
    147162              ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),               &
    148163              ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),             &
    149               rmask(nys-1:nyn+1,nxl-1:nxr+1,0:statistic_regions),           &
     164              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),           &
    150165              sums(nzb:nzt+1,pr_palm+max_pr_user),                          &
    151166              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),   &
    … …  
    154169              sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
    155170              ts_value(dots_max,0:statistic_regions) )
    156     ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
    157 
    158     ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &
    159               ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1),  &
    160               us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1),   &
    161               uswst_1(nys-1:nyn+1,nxl-1:nxr+1),                               &
    162               vsws_1(nys-1:nyn+1,nxl-1:nxr+1),                                &
    163               vswst_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
     171    ALLOCATE( km_damp_x(nxlg:nxrg), km_damp_y(nysg:nyng) )
     172
     173    ALLOCATE( rif_1(nysg:nyng,nxlg:nxrg), shf_1(nysg:nyng,nxlg:nxrg), &
     174              ts(nysg:nyng,nxlg:nxrg), tswst_1(nysg:nyng,nxlg:nxrg),  &
     175              us(nysg:nyng,nxlg:nxrg), usws_1(nysg:nyng,nxlg:nxrg),   &
     176              uswst_1(nysg:nyng,nxlg:nxrg),                               &
     177              vsws_1(nysg:nyng,nxlg:nxrg),                                &
     178              vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg) )
    164179
    165180    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    166181!
    167182!--    Leapfrog scheme needs two timelevels of diffusion quantities
    168        ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1),   &
    169                  shf_2(nys-1:nyn+1,nxl-1:nxr+1),   &
    170                  tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
    171                  usws_2(nys-1:nyn+1,nxl-1:nxr+1),  &
    172                  uswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
    173                  vswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
    174                  vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
     183       ALLOCATE( rif_2(nysg:nyng,nxlg:nxrg),   &
     184                 shf_2(nysg:nyng,nxlg:nxrg),   &
     185                 tswst_2(nysg:nyng,nxlg:nxrg), &
     186                 usws_2(nysg:nyng,nxlg:nxrg),  &
     187                 uswst_2(nysg:nyng,nxlg:nxrg), &
     188                 vswst_2(nysg:nyng,nxlg:nxrg), &
     189                 vsws_2(nysg:nyng,nxlg:nxrg) )
    175190    ENDIF
    176191
    177192    ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra),         &
    178               e_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    179               e_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    180               e_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    181               kh_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    182               km_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    183               p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),    &
    184               pt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    185               pt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    186               pt_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    187               tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    188               u_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    189               u_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    190               u_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    191               v_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    192               v_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    193               v_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    194               w_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    195               w_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    196               w_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     193              e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     194              e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     195              e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     196              kh_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     197              km_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     198              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
     199              pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     200              pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     201              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     202              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     203              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     204              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     205              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     206              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     207              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     208              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     209              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     210              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     211              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    197212
    198213    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    199        ALLOCATE( kh_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    200                  km_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     214       ALLOCATE( kh_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     215                 km_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    201216    ENDIF
    202217
    … …  
    204219!
    205220!--    2D-humidity/scalar arrays
    206        ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1),     &
    207                   qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
    208                   qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
     221       ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),     &
     222                  qsws_1(nysg:nyng,nxlg:nxrg), &
     223                  qswst_1(nysg:nyng,nxlg:nxrg) )
    209224
    210225       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    211           ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
    212                     qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
     226          ALLOCATE( qsws_2(nysg:nyng,nxlg:nxrg), &
     227                    qswst_2(nysg:nyng,nxlg:nxrg) )
    213228       ENDIF
    214229!
    215230!--    3D-humidity/scalar arrays
    216        ALLOCATE( q_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    217                  q_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    218                  q_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     231       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     232                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     233                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    219234
    220235!
    221236!--    3D-arrays needed for humidity only
    222237       IF ( humidity )  THEN
    223           ALLOCATE( vpt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     238          ALLOCATE( vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    224239
    225240          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    226              ALLOCATE( vpt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     241             ALLOCATE( vpt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    227242          ENDIF
    228243
    … …  
    230245!
    231246!--          Liquid water content
    232              ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     247             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    233248!
    234249!--          Precipitation amount and rate (only needed if output is switched)
    235              ALLOCATE( precipitation_amount(nys-1:nyn+1,nxl-1:nxr+1), &
    236                        precipitation_rate(nys-1:nyn+1,nxl-1:nxr+1) )
     250             ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &
     251                       precipitation_rate(nysg:nyng,nxlg:nxrg) )
    237252          ENDIF
    238253
    … …  
    241256!--          Liquid water content, change in liquid water content,
    242257!--          real volume of particles (with weighting), volume of particles
    243              ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    244                         ql_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    245                         ql_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    246                         ql_vp(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     258             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     259                        ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     260                        ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     261                        ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    247262          ENDIF
    248263
    … …  
    252267
    253268    IF ( ocean )  THEN
    254        ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
    255                  saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
    256        ALLOCATE( prho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    257                  rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    258                  sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
    259                  sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
    260                  sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     269       ALLOCATE( saswsb_1(nysg:nyng,nxlg:nxrg), &
     270                 saswst_1(nysg:nyng,nxlg:nxrg) )
     271       ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     272                 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     273                 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
     274                 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
     275                 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    261276       prho => prho_1
    262277       rho  => rho_1  ! routines calc_mean_profile and diffusion_e require
    263278                      ! density to be apointer
    264279       IF ( humidity_remote )  THEN
    265           ALLOCATE( qswst_remote(nys-1:nyn+1,nxl-1:nxr+1) )
     280          ALLOCATE( qswst_remote(nysg:nyng,nxlg:nxrg))
    266281          qswst_remote = 0.0
    267282       ENDIF
    … …  
    272287!-- particle velocities
    273288    IF ( use_sgs_for_particles )  THEN
    274        ALLOCATE ( diss(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     289       ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    275290    ELSE
    276291       ALLOCATE ( diss(2,2,2) )  ! required because diss is used as a
    … …  
    288303!-- 3D-arrays for the leaf area density and the canopy drag coefficient
    289304    IF ( plant_canopy ) THEN
    290        ALLOCATE ( lad_s(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    291                   lad_u(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    292                   lad_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    293                   lad_w(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    294                   cdc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     305       ALLOCATE ( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     306                  lad_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     307                  lad_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     308                  lad_w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
     309                  cdc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    295310
    296311       IF ( passive_scalar ) THEN
    297           ALLOCATE ( sls(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
    298                      sec(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     312          ALLOCATE ( sls(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
     313                     sec(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    299314       ENDIF
    300315
    301316       IF ( cthf /= 0.0 ) THEN
    302           ALLOCATE ( lai(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
    303                      canopy_heat_flux(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     317          ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
     318                     canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    304319       ENDIF
    305320
    … …  
    309324!-- 4D-array for storing the Rif-values at vertical walls
    310325    IF ( topography /= 'flat' )  THEN
    311        ALLOCATE( rif_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1,1:4) )
     326       ALLOCATE( rif_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg,1:4) )
    312327       rif_wall = 0.0
    313     ENDIF
    314 
    315 !
    316 !-- Velocities at nzb+1 needed for volume flow control
    317     IF ( conserve_volume_flow )  THEN
    318        ALLOCATE( u_nzb_p1_for_vfc(nys:nyn), v_nzb_p1_for_vfc(nxl:nxr) )
    319        u_nzb_p1_for_vfc = 0.0
    320        v_nzb_p1_for_vfc = 0.0
    321328    ENDIF
    322329
    … …  
    325332!-- are needed for radiation boundary conditions
    326333    IF ( outflow_l )  THEN
    327        ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,1:2), &
    328                  v_m_l(nzb:nzt+1,nys-1:nyn+1,0:1), &
    329                  w_m_l(nzb:nzt+1,nys-1:nyn+1,0:1) )
     334       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &
     335                 v_m_l(nzb:nzt+1,nysg:nyng,0:1), &
     336                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
    330337    ENDIF
    331338    IF ( outflow_r )  THEN
    332        ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
    333                  v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
    334                  w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx) )
     339       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
     340                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
     341                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
    335342    ENDIF
    336343    IF ( outflow_l  .OR.  outflow_r )  THEN
    337        ALLOCATE( c_u(nzb:nzt+1,nys-1:nyn+1), c_v(nzb:nzt+1,nys-1:nyn+1), &
    338                  c_w(nzb:nzt+1,nys-1:nyn+1) )
     344       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &
     345                 c_w(nzb:nzt+1,nysg:nyng) )
    339346    ENDIF
    340347    IF ( outflow_s )  THEN
    341        ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1), &
    342                  v_m_s(nzb:nzt+1,1:2,nxl-1:nxr+1), &
    343                  w_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1) )
     348       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &
     349                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &
     350                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
    344351    ENDIF
    345352    IF ( outflow_n )  THEN
    346        ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
    347                  v_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
    348                  w_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1) )
     353       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
     354                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
     355                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
    349356    ENDIF
    350357    IF ( outflow_s  .OR.  outflow_n )  THEN
    351        ALLOCATE( c_u(nzb:nzt+1,nxl-1:nxr+1), c_v(nzb:nzt+1,nxl-1:nxr+1), &
    352                  c_w(nzb:nzt+1,nxl-1:nxr+1) )
     358       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &
     359                 c_w(nzb:nzt+1,nxlg:nxrg) )
    353360    ENDIF
    354361
    … …  
    435442!
    436443!--       Transfer initial profiles to the arrays of the 3D model
    437           DO  i = nxl-1, nxr+1
    438              DO  j = nys-1, nyn+1
     444          DO  i = nxlg, nxrg
     445             DO  j = nysg, nyng
    439446                e(:,j,i)  = e1d
    440447                kh(:,j,i) = kh1d
    … …  
    447454
    448455          IF ( humidity  .OR.  passive_scalar )  THEN
    449              DO  i = nxl-1, nxr+1
    450                 DO  j = nys-1, nyn+1
     456             DO  i = nxlg, nxrg
     457                DO  j = nysg, nyng
    451458                   q(:,j,i) = q_init
    452459                ENDDO
    … …  
    455462
    456463          IF ( .NOT. constant_diffusion )  THEN
    457              DO  i = nxl-1, nxr+1
    458                 DO  j = nys-1, nyn+1
     464             DO  i = nxlg, nxrg
     465                DO  j = nysg, nyng
    459466                   e(:,j,i)  = e1d
    460467                ENDDO
    … …  
    505512                ENDDO
    506513             ENDDO
    507              IF ( conserve_volume_flow )  THEN
    508                 IF ( nxr == nx )  THEN
    509                    DO  j = nys, nyn
    510                       DO  k = nzb + 1, nzb_u_inner(j,nx)
    511                          u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &
    512                                                u1d(k) * dzu(k)
    513                       ENDDO
    514                    ENDDO
    515                 ENDIF
    516                 IF ( nyn == ny )  THEN
    517                    DO  i = nxl, nxr
    518                       DO  k = nzb + 1, nzb_v_inner(ny,i)
    519                          v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &
    520                                                v1d(k) * dzu(k)
    521                       ENDDO
    522                    ENDDO
    523                 ENDIF
    524              ENDIF
     514             
    525515!
    526516!--          WARNING: The extra boundary conditions set after running the
    … …  
    530520!--          ---------  advection scheme: keep u and v zero one layer below
    531521!--                     the topography.
    532              IF ( ibc_uv_b == 0 )  THEN
    533 !
    534 !--             Satisfying the Dirichlet condition with an extra layer below
    535 !--             the surface where the u and v component change their sign.
    536                 DO  i = nxl-1, nxr+1
    537                    DO  j = nys-1, nyn+1
    538                       IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
    539                       IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
    540                    ENDDO
    541                 ENDDO
    542 
    543              ELSE
     522!
     523!--           Following was removed, because mirror boundary condition are
     524!--           replaced by dirichlet boundary conditions
     525!
     526!             IF ( ibc_uv_b == 0 )  THEN
     527!!
     528!!--             Satisfying the Dirichlet condition with an extra layer below
     529!!--             the surface where the u and v component change their sign.
     530!                DO  i = nxl-1, nxr+1
     531!                   DO  j = nys-1, nyn+1
     532!                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
     533!                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
     534!                   ENDDO
     535!                ENDDO
     536!
     537!             ELSE
     538             IF ( ibc_uv_b == 1 )  THEN
    544539!
    545540!--             Neumann condition
    … …  
    560555!--       Use constructed initial profiles (velocity constant with height,
    561556!--       temperature profile with constant gradient)
    562           DO  i = nxl-1, nxr+1
    563              DO  j = nys-1, nyn+1
     557          DO  i = nxlg, nxrg
     558             DO  j = nysg, nyng
    564559                pt(:,j,i) = pt_init
    565560                u(:,j,i)  = u_init
    … …  
    574569!--       in the limiting formula!). The original values are stored to be later
    575570!--       used for volume flow control.
    576           DO  i = nxl-1, nxr+1
    577              DO  j = nys-1, nyn+1
     571          DO  i = nxlg, nxrg
     572             DO  j = nysg, nyng
    578573                u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
    579574                v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
    580575             ENDDO
    581576          ENDDO
    582           IF ( conserve_volume_flow )  THEN
    583              IF ( nxr == nx )  THEN
    584                 DO  j = nys, nyn
    585                    DO  k = nzb + 1, nzb_u_inner(j,nx) + 1
    586                       u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &
    587                                             u_init(k) * dzu(k)
    588                    ENDDO
    589                 ENDDO
    590              ENDIF
    591              IF ( nyn == ny )  THEN
    592                 DO  i = nxl, nxr
    593                    DO  k = nzb + 1, nzb_v_inner(ny,i) + 1
    594                       v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &
    595                                             v_init(k) * dzu(k)
    596                    ENDDO
    597                 ENDDO
    598              ENDIF
    599           ENDIF
    600577
    601578          IF ( humidity  .OR.  passive_scalar )  THEN
    602              DO  i = nxl-1, nxr+1
    603                 DO  j = nys-1, nyn+1
     579             DO  i = nxlg, nxrg
     580                DO  j = nysg, nyng
    604581                   q(:,j,i) = q_init
    605582                ENDDO
    … …  
    608585
    609586          IF ( ocean )  THEN
    610              DO  i = nxl-1, nxr+1
    611                 DO  j = nys-1, nyn+1
     587             DO  i = nxlg, nxrg
     588                DO  j = nysg, nyng
    612589                   sa(:,j,i) = sa_init
    613590                ENDDO
    … …  
    660637
    661638       ENDIF
     639!
     640!--    Bottom boundary
     641       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
     642          u(nzb,:,:) = 0.0
     643          v(nzb,:,:) = 0.0
     644       ENDIF
    662645
    663646!
    … …  
    683666       hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
    684667       hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
    685        IF ( ibc_uv_b == 0 )  THEN
    686           hom(nzb,1,5,:) = -hom(nzb+1,1,5,:)  ! due to satisfying the Dirichlet
    687           hom(nzb,1,6,:) = -hom(nzb+1,1,6,:)  ! condition with an extra layer
    688               ! below the surface where the u and v component change their sign
     668       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
     669          hom(nzb,1,5,:) = 0.0   
     670          hom(nzb,1,6,:) = 0.0 
    689671       ENDIF
    690672       hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
    … …  
    733715!--             Over topography surface_heatflux is replaced by wall_heatflux(0)
    734716                IF ( TRIM( topography ) /= 'flat' )  THEN
    735                    DO  i = nxl-1, nxr+1
    736                       DO  j = nys-1, nyn+1
     717                   DO  i = nxlg, nxrg
     718                      DO  j = nysg, nyng
    737719                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
    738720                            shf(j,i) = wall_heatflux(0)
    … …  
    755737                IF ( TRIM( topography ) /= 'flat' )  THEN
    756738                   wall_qflux = wall_humidityflux
    757                    DO  i = nxl-1, nxr+1
    758                       DO  j = nys-1, nyn+1
     739                   DO  i = nxlg, nxrg
     740                      DO  j = nysg, nyng
    759741                         IF ( nzb_s_inner(j,i) /= 0 )  THEN
    760742                            qsws(j,i) = wall_qflux(0)
    … …  
    827809       ENDIF
    828810
    829 !
    830 !--    Calculate the initial volume flow at the right and north boundary
    831        IF ( conserve_volume_flow )  THEN
    832 
    833           volume_flow_initial_l = 0.0
    834           volume_flow_area_l    = 0.0
    835  
    836           IF ( nxr == nx )  THEN
    837              DO  j = nys, nyn
    838                 DO  k = nzb_2d(j,nx) + 1, nzt
    839                    volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
    840                                               u(k,j,nx) * dzu(k)
    841                    volume_flow_area_l(1)    = volume_flow_area_l(1) + dzu(k)
    842                 ENDDO
    843 !
    844 !--             Correction if velocity at nzb+1 has been set zero further above
    845                 volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
    846                                            u_nzb_p1_for_vfc(j)
    847              ENDDO
    848           ENDIF
    849 
    850           IF ( nyn == ny )  THEN
    851              DO  i = nxl, nxr
    852                 DO  k = nzb_2d(ny,i) + 1, nzt 
    853                    volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
    854                                               v(k,ny,i) * dzu(k)
    855                    volume_flow_area_l(2)    = volume_flow_area_l(2) + dzu(k)
    856                 ENDDO
    857 !
    858 !--             Correction if velocity at nzb+1 has been set zero further above
    859                 volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
    860                                            v_nzb_p1_for_vfc(i)
    861              ENDDO
    862           ENDIF
    863 
    864 #if defined( __parallel )
    865           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    866           CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
    867                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
    868           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    869           CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
    870                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
    871 #else
    872           volume_flow_initial = volume_flow_initial_l
    873           volume_flow_area    = volume_flow_area_l
    874 #endif
    875 !
    876 !--       In case of 'bulk_velocity' mode, volume_flow_initial is overridden
    877 !--       and calculated from u|v_bulk instead.
    878           IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
    879              volume_flow_initial(1) = u_bulk * volume_flow_area(1)
    880              volume_flow_initial(2) = v_bulk * volume_flow_area(2)
    881           ENDIF
    882 
    883        ENDIF
    884811
    885812!
    … …  
    968895          sa_p  = sa
    969896       ENDIF
    970 
     897       
    971898
    972899    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.    &
    973              TRIM( initializing_actions ) == 'cyclic_fill' )  &
     900         TRIM( initializing_actions ) == 'cyclic_fill' )  &
    974901    THEN
    975902!
    … …  
    978905       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
    979906
    980           WRITE (9,*) 'before read_parts_of_var_list'
    981           CALL local_flush( 9 )
    982907          CALL read_parts_of_var_list
    983           WRITE (9,*) 'after read_parts_of_var_list'
    984           CALL local_flush( 9 )
    985908          CALL close_file( 13 )
    986909
    … …  
    1002925!--          conditions are used)
    1003926             IF ( inflow_l )  THEN
    1004                 DO  j = nys-1, nyn+1
     927                DO  j = nysg, nyng
    1005928                   DO  k = nzb, nzt+1
    1006                       u(k,j,-1)  = mean_inflow_profiles(k,1)
    1007                       v(k,j,-1)  = mean_inflow_profiles(k,2)
    1008                       w(k,j,-1)  = 0.0
    1009                       pt(k,j,-1) = mean_inflow_profiles(k,4)
    1010                       e(k,j,-1)  = mean_inflow_profiles(k,5)
     929                      u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
     930                      v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
     931                      w(k,j,nxlg:-1)  = 0.0
     932                      pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
     933                      e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
    1011934                   ENDDO
    1012935                ENDDO
    … …  
    1064987!
    1065988!--    Read binary data from restart file
    1066           WRITE (9,*) 'before read_3d_binary'
    1067           CALL local_flush( 9 )
     989
    1068990       CALL read_3d_binary
    1069           WRITE (9,*) 'after read_3d_binary'
    1070           CALL local_flush( 9 )
    1071991
    1072992!
    … …  
    1074994       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.  &
    1075995            topography /= 'flat' )  THEN
    1076 !
    1077 !--       Correction of initial volume flow
    1078           IF ( conserve_volume_flow )  THEN
    1079              IF ( nxr == nx )  THEN
    1080                 DO  j = nys, nyn
    1081                    DO  k = nzb + 1, nzb_u_inner(j,nx)
    1082                       u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &
    1083                                             u(k,j,nx) * dzu(k)
    1084                    ENDDO
    1085                 ENDDO
    1086              ENDIF
    1087              IF ( nyn == ny )  THEN
    1088                 DO  i = nxl, nxr
    1089                    DO  k = nzb + 1, nzb_v_inner(ny,i)
    1090                       v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &
    1091                                             v(k,ny,i) * dzu(k)
    1092                    ENDDO
    1093                 ENDDO
    1094              ENDIF
    1095           ENDIF
    1096 
    1097996!
    1098997!--       Inside buildings set velocities and TKE back to zero.
    … …  
    1100999!--       maybe revise later.
    11011000          IF ( timestep_scheme(1:5) == 'runge' )  THEN
    1102              DO  i = nxl-1, nxr+1
    1103                 DO  j = nys-1, nyn+1
     1001             DO  i = nxlg, nxrg
     1002                DO  j = nysg, nyng
    11041003                   u  (nzb:nzb_u_inner(j,i),j,i) = 0.0
    11051004                   v  (nzb:nzb_v_inner(j,i),j,i) = 0.0
    … …  
    11181017             ENDDO
    11191018          ELSE
    1120              DO  i = nxl-1, nxr+1
    1121                 DO  j = nys-1, nyn+1
     1019             DO  i = nxlg, nxrg
     1020                DO  j = nysg, nyng
    11221021                   u  (nzb:nzb_u_inner(j,i),j,i) = 0.0
    11231022                   v  (nzb:nzb_v_inner(j,i),j,i) = 0.0
    … …  
    11391038
    11401039!
    1141 !--    Calculate the initial volume flow at the right and north boundary
    1142        IF ( conserve_volume_flow  .AND.  &
    1143             TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
    1144 
    1145           volume_flow_initial_l = 0.0
    1146           volume_flow_area_l    = 0.0
    1147  
    1148           IF ( nxr == nx )  THEN
    1149              DO  j = nys, nyn
    1150                 DO  k = nzb_2d(j,nx) + 1, nzt
    1151                    volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
    1152                                               u(k,j,nx) * dzu(k)
    1153                    volume_flow_area_l(1)    = volume_flow_area_l(1) + dzu(k)
    1154                 ENDDO
    1155 !
    1156 !--             Correction if velocity inside buildings has been set to zero
    1157 !--             further above
    1158                 volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
    1159                                            u_nzb_p1_for_vfc(j)
    1160              ENDDO
    1161           ENDIF
    1162 
    1163           IF ( nyn == ny )  THEN
    1164              DO  i = nxl, nxr
    1165                 DO  k = nzb_2d(ny,i) + 1, nzt 
    1166                    volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
    1167                                               v(k,ny,i) * dzu(k)
    1168                    volume_flow_area_l(2)    = volume_flow_area_l(2) + dzu(k)
    1169                 ENDDO
    1170 !
    1171 !--             Correction if velocity inside buildings has been set to zero
    1172 !--             further above
    1173                 volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
    1174                                            v_nzb_p1_for_vfc(i)
    1175              ENDDO
    1176           ENDIF
    1177 
    1178 #if defined( __parallel )
    1179           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1180           CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
    1181                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
    1182           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1183           CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
    1184                               2, MPI_REAL, MPI_SUM, comm2d, ierr )
    1185 #else
    1186           volume_flow_initial = volume_flow_initial_l
    1187           volume_flow_area    = volume_flow_area_l
    1188 #endif 
    1189        ENDIF
    1190 
    1191 
    1192 !
    11931040!--    Calculate initial temperature field and other constants used in case
    11941041!--    of a sloping surface
    … …  
    12431090          w_m_n(:,:,:) = w(:,ny-1:ny,:)
    12441091       ENDIF
    1245 
    1246     ENDIF
    1247 
     1092       
     1093    ENDIF
     1094!
     1095!-- Calculate the initial volume flow at the right and north boundary
     1096    IF ( conserve_volume_flow ) THEN
     1097
     1098       volume_flow_initial_l = 0.0
     1099       volume_flow_area_l    = 0.0
     1100 
     1101       IF  ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
     1102
     1103          IF ( nxr == nx )  THEN
     1104             DO  j = nys, nyn
     1105                DO  k = nzb_2d(j,nx) + 1, nzt
     1106                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
     1107                                              hom_sum(k,1,0) * dzw(k)
     1108                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
     1109                ENDDO
     1110             ENDDO
     1111          ENDIF
     1112         
     1113          IF ( nyn == ny )  THEN
     1114             DO  i = nxl, nxr
     1115                DO  k = nzb_2d(ny,i) + 1, nzt 
     1116                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
     1117                                               hom_sum(k,2,0) * dzw(k)
     1118                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
     1119                ENDDO
     1120             ENDDO
     1121          ENDIF
     1122
     1123       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     1124
     1125          IF ( nxr == nx )  THEN
     1126             DO  j = nys, nyn
     1127                DO  k = nzb_2d(j,nx) + 1, nzt
     1128                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
     1129                                               u(k,j,nx) * dzw(k)
     1130                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
     1131                ENDDO
     1132             ENDDO
     1133          ENDIF
     1134         
     1135          IF ( nyn == ny )  THEN
     1136             DO  i = nxl, nxr
     1137                DO  k = nzb_2d(ny,i) + 1, nzt 
     1138                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
     1139                                              v(k,ny,i) * dzw(k)
     1140                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
     1141                ENDDO
     1142             ENDDO
     1143          ENDIF
     1144
     1145       ENDIF
     1146
     1147#if defined( __parallel )
     1148          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
     1149                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
     1150          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
     1151                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
     1152
     1153          CALL MPI_ALLREDUCE( volume_flow_initial_l(2), volume_flow_initial(2),&
     1154                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
     1155          CALL MPI_ALLREDUCE( volume_flow_area_l(2), volume_flow_area(2),      &
     1156                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
     1157
     1158#else
     1159          volume_flow_initial = volume_flow_initial_l
     1160          volume_flow_area    = volume_flow_area_l
     1161#endif 
     1162
     1163!
     1164!--       In case of 'bulk_velocity' mode, volume_flow_initial is overridden
     1165!--       and calculated from u|v_bulk instead.
     1166          IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
     1167             volume_flow_initial(1) = u_bulk * volume_flow_area(1)
     1168             volume_flow_initial(2) = v_bulk * volume_flow_area(2)
     1169          ENDIF
     1170
     1171       ENDIF
    12481172!
    12491173!-- Initialization of the leaf area density
    … …  
    12541178          CASE( 'block' )
    12551179
    1256              DO  i = nxl-1, nxr+1
    1257                 DO  j = nys-1, nyn+1
     1180             DO  i = nxlg, nxrg
     1181                DO  j = nysg, nyng
    12581182                   lad_s(:,j,i) = lad(:)
    12591183                   cdc(:,j,i)   = drag_coefficient
    … …  
    12771201          END SELECT
    12781202
    1279        CALL exchange_horiz( lad_s )
    1280        CALL exchange_horiz( cdc )
     1203       CALL exchange_horiz( lad_s, nbgp )
     1204       CALL exchange_horiz( cdc, nbgp )
    12811205
    12821206       IF ( passive_scalar ) THEN
    1283           CALL exchange_horiz( sls )
    1284           CALL exchange_horiz( sec )
     1207          CALL exchange_horiz( sls, nbgp )
     1208          CALL exchange_horiz( sec, nbgp )
    12851209       ENDIF
    12861210
    … …  
    13111235       lad_w(nzt+1,:,:)     = lad_w(nzt,:,:)
    13121236
    1313        CALL exchange_horiz( lad_u )
    1314        CALL exchange_horiz( lad_v )
    1315        CALL exchange_horiz( lad_w )
     1237       CALL exchange_horiz( lad_u, nbgp )
     1238       CALL exchange_horiz( lad_v, nbgp )
     1239       CALL exchange_horiz( lad_w, nbgp )
    13161240
    13171241!
    … …  
    13221246!--       integration of the leaf area density
    13231247          lai(:,:,:) = 0.0
    1324           DO  i = nxl-1, nxr+1
    1325              DO  j = nys-1, nyn+1
     1248          DO  i = nxlg, nxrg
     1249             DO  j = nysg, nyng
    13261250                DO  k = pch_index-1, 0, -1
    13271251                   lai(k,j,i) = lai(k+1,j,i) +                   &
    … …  
    13391263!--       Evaluation of the upward kinematic vertical heat flux within the
    13401264!--       canopy
    1341           DO  i = nxl-1, nxr+1
    1342              DO  j = nys-1, nyn+1
     1265          DO  i = nxlg, nxrg
     1266             DO  j = nysg, nyng
    13431267                DO  k = 0, pch_index
    13441268                   canopy_heat_flux(k,j,i) = cthf *                    &
    … …  
    13841308!-- Initialize quantities for special advections schemes
    13851309    CALL init_advec
     1310    IF ( momentum_advec == 'ws-scheme' .OR.  &
     1311         scalar_advec == 'ws-scheme' ) CALL ws_init
    13861312
    13871313!
    … …  
    14391365    IF ( bc_lr == 'dirichlet/radiation' )  THEN
    14401366
    1441        DO  i = nxl-1, nxr+1
     1367       DO  i = nxlg, nxrg
    14421368          IF ( i >= nx - outflow_damping_width )  THEN
    14431369             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
    … …  
    14521378    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
    14531379
    1454        DO  i = nxl-1, nxr+1
     1380       DO  i = nxlg, nxrg
    14551381          IF ( i <= outflow_damping_width )  THEN
    14561382             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
    … …  
    14671393    IF ( bc_ns == 'dirichlet/radiation' )  THEN
    14681394
    1469        DO  j = nys-1, nyn+1
     1395       DO  j = nysg, nyng
    14701396          IF ( j >= ny - outflow_damping_width )  THEN
    14711397             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
    … …  
    14801406    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
    14811407
    1482        DO  j = nys-1, nyn+1
     1408       DO  j = nysg, nyng
    14831409          IF ( j <= outflow_damping_width )  THEN
    14841410             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
    … …  
    15941520!-- buoyancy, etc. A zero value will occur for cases where all grid points of
    15951521!-- the respective subdomain lie below the surface topography
    1596     ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   )
     1522    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   ) 
    15971523    ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),            &
    15981524                           ngp_3d_inner(:) )
    1599     ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) )
     1525    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
    16001526
    16011527    DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l, ngp_3d_inner_tmp )
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |