Changeset 667 for palm/trunk/SOURCE

Timestamp:

Dec 23, 2010 12:06:00 PM (12 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:

• . (modified) (1 prop)
• Makefile (modified) (7 diffs)
• advec_particles.f90 (modified) (3 diffs)
• advec_w_pw.f90 (modified) (2 diffs)
• advec_ws.f90 (copied) (copied from palm/branches/suehring/advec_ws.f90)
• asselin_filter.f90 (modified) (2 diffs)
• average_3d_data.f90 (modified) (24 diffs)
• boundary_conds.f90 (modified) (18 diffs)
• calc_liquid_water_content.f90 (modified) (2 diffs)
• calc_precipitation.f90 (modified) (1 diff)
• calc_spectra.f90 (modified) (2 diffs)
• check_for_restart.f90 (modified) (3 diffs)
• check_open.f90 (modified) (9 diffs)
• check_parameters.f90 (modified) (27 diffs)
• data_output_2d.f90 (modified) (51 diffs)
• data_output_3d.f90 (modified) (13 diffs)
• data_output_mask.f90 (modified) (6 diffs)
• diffusion_e.f90 (modified) (3 diffs)
• diffusion_s.f90 (modified) (3 diffs)
• diffusion_u.f90 (modified) (4 diffs)
• diffusion_v.f90 (modified) (3 diffs)
• diffusion_w.f90 (modified) (3 diffs)
• diffusivities.f90 (modified) (5 diffs)
• disturb_field.f90 (modified) (8 diffs)
• exchange_horiz.f90 (modified) (7 diffs)
• exchange_horiz_2d.f90 (modified) (10 diffs)
• flow_statistics.f90 (modified) (23 diffs)
• global_min_max.f90 (modified) (5 diffs)
• header.f90 (modified) (6 diffs)
• inflow_turbulence.f90 (modified) (6 diffs)
• init_1d_model.f90 (modified) (2 diffs)
• init_3d_model.f90 (modified) (45 diffs)
• init_advec.f90 (modified) (1 diff)
• init_coupling.f90 (modified) (4 diffs)
• init_grid.f90 (modified) (14 diffs)
• init_particles.f90 (modified) (3 diffs)
• init_pegrid.f90 (modified) (9 diffs)
• init_pt_anomaly.f90 (modified) (2 diffs)
• init_rankine.f90 (modified) (3 diffs)
• init_slope.f90 (modified) (3 diffs)
• local_stop.f90 (modified) (2 diffs)
• modules.f90 (modified) (14 diffs)
• palm.f90 (modified) (2 diffs)
• parin.f90 (modified) (4 diffs)
• poisfft.f90 (modified) (6 diffs)
• poisfft_hybrid.f90 (modified) (3 diffs)
• poismg.f90 (modified) (8 diffs)
• prandtl_fluxes.f90 (modified) (16 diffs)
• pres.f90 (modified) (27 diffs)
• production_e.f90 (modified) (2 diffs)
• prognostic_equations.f90 (modified) (30 diffs)
• read_3d_binary.f90 (modified) (18 diffs)
• sor.f90 (modified) (6 diffs)
• subsidence.f90 (modified) (3 diffs)
• sum_up_3d_data.f90 (modified) (50 diffs)
• surface_coupler.f90 (modified) (8 diffs)
• time_integration.f90 (modified) (6 diffs)
• timestep.f90 (modified) (4 diffs)
• user_3d_data_averaging.f90 (modified) (3 diffs)
• user_actions.f90 (modified) (2 diffs)
• user_data_output_2d.f90 (modified) (4 diffs)
• user_data_output_3d.f90 (modified) (4 diffs)
• user_init.f90 (modified) (2 diffs)
• user_init_plant_canopy.f90 (modified) (2 diffs)
• user_read_restart_data.f90 (modified) (2 diffs)
• wall_fluxes.f90 (modified) (2 diffs)
• write_compressed.f90 (modified) (3 diffs)

palm/trunk/SOURCE/Makefile

@@ -4 +4 @@
 # Current revisions:
 # ------------------
-# 
+# +advec_ws
 #
 # Former revisions:
@@ -57 +57 @@
 
 RCS = advec_particles.f90 advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 \
-        advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 advec_u_ups.f90 \
-        advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 advec_w_pw.f90 \
-        advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 average_3d_data.f90 \
-        boundary_conds.f90 buoyancy.f90 calc_liquid_water_content.f90 \
-        calc_precipitation.f90 calc_radiation.f90 calc_spectra.f90 \
-        check_for_restart.f90 check_open.f90 check_parameters.f90 \
-        close_file.f90 compute_vpt.f90 coriolis.f90 cpu_log.f90 \
-        cpu_statistics.f90 data_log.f90 data_output_dvrp.f90 \
-        data_output_mask.f90 data_output_profiles.f90 data_output_ptseries.f90 \
-        data_output_spectra.f90 data_output_tseries.f90 data_output_2d.f90 \
-        data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 diffusion_u.f90 \
-        diffusion_v.f90 diffusion_w.f90 diffusivities.f90 disturb_field.f90 \
-        disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 \
-        exchange_horiz_2d.f90 \
+        advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \
+        advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \
+        advec_w_pw.f90 advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 \
+        average_3d_data.f90 boundary_conds.f90 buoyancy.f90 \
+        calc_liquid_water_content.f90 calc_precipitation.f90 \
+        calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \
+        check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \
+        coriolis.f90 cpu_log.f90 cpu_statistics.f90 data_log.f90 \
+        data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \
+        data_output_ptseries.f90 data_output_spectra.f90 data_output_tseries.f90 \
+        data_output_2d.f90 data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 \
+        diffusion_u.f90 diffusion_v.f90 diffusion_w.f90 diffusivities.f90 \
+        disturb_field.f90 disturb_heatflux.f90 eqn_state_seawater.f90 \
+        exchange_horiz.f90 exchange_horiz_2d.f90 \
         fft_xy.f90 flow_statistics.f90 global_min_max.f90 \
         header.f90 impact_of_latent_heat.f90 inflow_turbulence.f90 \
@@ -77 +77 @@
         init_masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \
         init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
-        interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \
-        local_stop.f90 local_system.f90 local_tremain.f90 \
+        interaction_droplets_ptq.f90 local_flush.f90 \
+        local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
         local_tremain_ini.f90 message.f90 modules.f90 netcdf.f90 \
         package_parin.f90 palm.f90 parin.f90 particle_boundary_conds.f90 \
@@ -106 +106 @@
 OBJS =  advec_particles.o advec_s_bc.o advec_s_pw.o advec_s_up.o \
         advec_s_ups.o advec_u_pw.o advec_u_up.o advec_u_ups.o \
-        advec_v_pw.o advec_v_up.o advec_v_ups.o advec_w_pw.o \
+        advec_ws.o advec_v_pw.o advec_v_up.o advec_v_ups.o advec_w_pw.o \
         advec_w_up.o advec_w_ups.o asselin_filter.o average_3d_data.o \
         boundary_conds.o buoyancy.o calc_liquid_water_content.o \
@@ -184 +184 @@
 advec_v_up.o: modules.o
 advec_v_ups.o: modules.o
+advec_ws.o: modules.o
 advec_w_pw.o: modules.o
 advec_w_up.o: modules.o
@@ -230 +231 @@
 inflow_turbulence.o: modules.o
 init_1d_model.o: modules.o
-init_3d_model.o: modules.o random_function.o
+init_3d_model.o: modules.o random_function.o advec_ws.o
 init_advec.o: modules.o
 init_cloud_physics.o: modules.o
@@ -264 +265 @@
 production_e.o: modules.o wall_fluxes.o
 prognostic_equations.o: modules.o advec_s_pw.o advec_s_up.o advec_u_pw.o \
+        advec_ws.o \
         advec_u_up.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o  \
         buoyancy.o calc_precipitation.o calc_radiation.o coriolis.o \

palm/trunk/SOURCE/advec_particles.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
+! Declaration of de_dx, de_dy, de_dz adapted to additional 
+! ghost points. Furthermore the calls of exchange_horiz were modified.
 ! 
 ! TEST: PRINT statements on unit 9 (commented out)
@@ -153 +155 @@
     REAL    ::  location(1:30,1:3)
 
-    REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  de_dx, de_dy, de_dz
+    REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  de_dx, de_dy, de_dz
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::  trlpt, trnpt, trrpt, trspt
@@ -768 +770 @@
 !
 !--    Lateral boundary conditions
-       CALL exchange_horiz( de_dx )
-       CALL exchange_horiz( de_dy )
-       CALL exchange_horiz( de_dz )
-       CALL exchange_horiz( diss  )
+       CALL exchange_horiz( de_dx, nbgp )
+       CALL exchange_horiz( de_dy, nbgp )
+       CALL exchange_horiz( de_dz, nbgp )
+       CALL exchange_horiz( diss, nbgp  )
 
 !

palm/trunk/SOURCE/advec_w_pw.f90

@@ -89 +89 @@
        REAL    ::  gu, gv
 
-
        gu = 2.0 * u_gtrans
        gv = 2.0 * v_gtrans
@@ -103 +102 @@
                                                 )
        ENDDO
-
     END SUBROUTINE advec_w_pw_ij
 

palm/trunk/SOURCE/asselin_filter.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
 ! Former revisions:
@@ -49 +49 @@
 !$OMP PARALLEL PRIVATE (i,j,k)
 !$OMP DO
-    DO  i = nxl-1, nxr+1
-       DO  j = nys-1, nyn+1
+    DO  i = nxlg, nxrg
+       DO  j = nysg, nyng
 
           DO  k = nzb_2d(j,i), nzt+1

palm/trunk/SOURCE/average_3d_data.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
 ! Former revisions:
@@ -58 +58 @@
 
           CASE ( 'e' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       e_av(k,j,i) = e_av(k,j,i) / REAL( average_count_3d )
@@ -67 +67 @@
 
           CASE ( 'qsws*' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    qsws_av(j,i) = qsws_av(j,i) / REAL( average_count_3d )
                 ENDDO
@@ -74 +74 @@
 
           CASE ( 'lwp*' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    lwp_av(j,i) = lwp_av(j,i) / REAL( average_count_3d )
                 ENDDO
@@ -81 +81 @@
 
           CASE ( 'p' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       p_av(k,j,i) = p_av(k,j,i) / REAL( average_count_3d )
@@ -108 +108 @@
 
           CASE ( 'prr*' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    precipitation_rate_av(j,i) = precipitation_rate_av(j,i) / &
                                                 REAL( average_count_3d )
@@ -116 +116 @@
 
           CASE ( 'pt' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       pt_av(k,j,i) = pt_av(k,j,i) / REAL( average_count_3d )
@@ -125 +125 @@
 
           CASE ( 'q' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d )
@@ -134 +134 @@
 
           CASE ( 'ql' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d )
@@ -143 +143 @@
 
           CASE ( 'ql_c' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       ql_c_av(k,j,i) = ql_c_av(k,j,i) / REAL( average_count_3d )
@@ -152 +152 @@
 
           CASE ( 'ql_v' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       ql_v_av(k,j,i) = ql_v_av(k,j,i) / REAL( average_count_3d )
@@ -161 +161 @@
 
           CASE ( 'ql_vp' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       ql_vp_av(k,j,i) = ql_vp_av(k,j,i) / &
@@ -171 +171 @@
 
           CASE ( 'qv' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       qv_av(k,j,i) = qv_av(k,j,i) / REAL( average_count_3d )
@@ -180 +180 @@
 
           CASE ( 'rho' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       rho_av(k,j,i) = rho_av(k,j,i) / REAL( average_count_3d )
@@ -189 +189 @@
 
           CASE ( 's' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       s_av(k,j,i) = s_av(k,j,i) / REAL( average_count_3d )
@@ -198 +198 @@
 
           CASE ( 'sa' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d )
@@ -207 +207 @@
 
          CASE ( 'shf*' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    shf_av(j,i) = shf_av(j,i) / REAL( average_count_3d )
                 ENDDO
@@ -214 +214 @@
 
           CASE ( 't*' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    ts_av(j,i) = ts_av(j,i) / REAL( average_count_3d )
                 ENDDO
@@ -221 +221 @@
 
           CASE ( 'u' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d )
@@ -230 +230 @@
 
           CASE ( 'u*' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    us_av(j,i) = us_av(j,i) / REAL( average_count_3d )
                 ENDDO
@@ -237 +237 @@
 
           CASE ( 'v' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d )
@@ -246 +246 @@
 
           CASE ( 'vpt' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       vpt_av(k,j,i) = vpt_av(k,j,i) / REAL( average_count_3d )
@@ -255 +255 @@
 
           CASE ( 'w' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d )
@@ -264 +264 @@
 
           CASE ( 'z0*' )
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    z0_av(j,i) = z0_av(j,i) / REAL( average_count_3d )
                 ENDDO

palm/trunk/SOURCE/boundary_conds.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
+!
+! 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 
+!  
 ! Former revisions:
 ! -----------------
@@ -70 +76 @@
     IF ( range == 'main')  THEN
 !
-!--    Bottom boundary
-       IF ( ibc_uv_b == 0 )  THEN
-!
-!--       Satisfying the Dirichlet condition with an extra layer below the
-!--       surface where the u and v component change their sign
-          u_p(nzb,:,:) = -u_p(nzb+1,:,:)
-          v_p(nzb,:,:) = -v_p(nzb+1,:,:)
-       ELSE
+!--    Bottom boundary 
+       IF ( ibc_uv_b == 1 )  THEN
           u_p(nzb,:,:) = u_p(nzb+1,:,:)
           v_p(nzb,:,:) = v_p(nzb+1,:,:)
        ENDIF
-       DO  i = nxl-1, nxr+1
-          DO  j = nys-1, nyn+1
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              w_p(nzb_w_inner(j,i),j,i) = 0.0
           ENDDO
@@ -90 +90 @@
 !--    Top boundary
        IF ( ibc_uv_t == 0 )  THEN
-          u_p(nzt+1,:,:) = ug(nzt+1)
-          v_p(nzt+1,:,:) = vg(nzt+1)
+           u_p(nzt+1,:,:) = ug(nzt+1)
+           v_p(nzt+1,:,:) = vg(nzt+1)
        ELSE
-          u_p(nzt+1,:,:) = u_p(nzt,:,:)
-          v_p(nzt+1,:,:) = v_p(nzt,:,:)
+           u_p(nzt+1,:,:) = u_p(nzt,:,:)
+           v_p(nzt+1,:,:) = v_p(nzt,:,:)
        ENDIF
        w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
@@ -103 +103 @@
 !--    the sea surface temperature of the coupled ocean model.
        IF ( ibc_pt_b == 0 )  THEN
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
              ENDDO
           ENDDO
        ELSEIF ( ibc_pt_b == 1 )  THEN
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
              ENDDO
@@ -119 +119 @@
 !--    Temperature at top boundary
        IF ( ibc_pt_t == 0 )  THEN
-          pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
+           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
        ELSEIF ( ibc_pt_t == 1 )  THEN
-          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
+           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
        ELSEIF ( ibc_pt_t == 2 )  THEN
-          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
+           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
        ENDIF
 
@@ -130 +130 @@
 !--    Generally Neumann conditions with de/dz=0 are assumed
        IF ( .NOT. constant_diffusion )  THEN
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
              ENDDO
@@ -144 +144 @@
 !--       Bottom boundary: Neumann condition because salinity flux is always
 !--       given
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
              ENDDO
@@ -153 +153 @@
 !--       Top boundary: Dirichlet or Neumann
           IF ( ibc_sa_t == 0 )  THEN
-             sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
+              sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
           ELSEIF ( ibc_sa_t == 1 )  THEN
-             sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
+              sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
           ENDIF
 
@@ -167 +167 @@
 !--       Surface conditions for constant_humidity_flux
           IF ( ibc_q_b == 0 ) THEN
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
                 ENDDO
              ENDDO
           ELSE
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
                 ENDDO
@@ -182 +182 @@
 !--       Top boundary
           q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
+
+
        ENDIF
 
@@ -226 +228 @@
        c_max = dy / dt_3d
 
-       DO i = nxl-1, nxr+1
+       DO i = nxlg, nxrg
           DO k = nzb+1, nzt+1
 
@@ -299 +301 @@
 !--    Bottom boundary at the outflow
        IF ( ibc_uv_b == 0 )  THEN
-          u_p(nzb,-1,:) = -u_p(nzb+1,-1,:) 
-          v_p(nzb,0,:)  = -v_p(nzb+1,0,:)  
+          u_p(nzb,-1,:) = 0.0 
+          v_p(nzb,0,:)  = 0.0  
        ELSE                    
           u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
@@ -324 +326 @@
        c_max = dy / dt_3d
 
-       DO i = nxl-1, nxr+1
+       DO i = nxlg, nxrg
           DO k = nzb+1, nzt+1
 
@@ -397 +399 @@
 !--    Bottom boundary at the outflow
        IF ( ibc_uv_b == 0 )  THEN
-          u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:) 
-          v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:)  
+          u_p(nzb,ny+1,:) = 0.0 
+          v_p(nzb,ny+1,:) = 0.0   
        ELSE                    
           u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
@@ -422 +424 @@
        c_max = dx / dt_3d
 
-       DO j = nys-1, nyn+1
+       DO j = nysg, nyng
           DO k = nzb+1, nzt+1
 
@@ -495 +497 @@
 !--    Bottom boundary at the outflow
        IF ( ibc_uv_b == 0 )  THEN
-          u_p(nzb,:,-1) = -u_p(nzb+1,:,-1) 
-          v_p(nzb,:,-1) = -v_p(nzb+1,:,-1)  
+          u_p(nzb,:,0)  = 0.0 
+          v_p(nzb,:,-1) = 0.0
        ELSE                    
-          u_p(nzb,:,-1) =  u_p(nzb+1,:,-1)
+          u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
           v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
        ENDIF
@@ -520 +522 @@
        c_max = dx / dt_3d
 
-       DO j = nys-1, nyn+1
+       DO j = nysg, nyng
           DO k = nzb+1, nzt+1
 
@@ -593 +595 @@
 !--    Bottom boundary at the outflow
        IF ( ibc_uv_b == 0 )  THEN
-          u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1) 
-          v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1)  
+          u_p(nzb,:,nx+1) = 0.0
+          v_p(nzb,:,nx+1) = 0.0 
        ELSE                    
           u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)

palm/trunk/SOURCE/calc_liquid_water_content.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
 ! Former revisions:
@@ -47 +47 @@
     REAL :: alpha, e_s, q_s, t_l
 
-    DO  i = nxl-1, nxr+1
-       DO  j = nys-1, nyn+1
+    DO  i = nxlg, nxrg
+       DO  j = nysg, nyng
           DO  k = nzb_2d(j,i)+1, nzt
 

palm/trunk/SOURCE/calc_precipitation.f90

@@ -8 +8 @@
 ! Former revisions:
 ! -----------------
-! $Id$
+! $Id: calc_precipitation.f90 484 2010-02-05 07:36:54Z raasch 
 !
 ! 403 2009-10-22 13:57:16Z franke

palm/trunk/SOURCE/calc_spectra.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
+! of tend
 !
 ! Former revisions:
@@ -152 +153 @@
     IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
        DEALLOCATE( tend )
-       ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+       ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ENDIF
 

palm/trunk/SOURCE/check_for_restart.f90

@@ -5 +5 @@
 ! -----------------
 ! 
+! Exchange of terminate_coupled between ocean and atmosphere by PE0
 !
 ! Former revisions:
@@ -93 +94 @@
        terminate_coupled = 3
 #if defined( __parallel )
-       CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,          &
-                          target_id, 0,                                      &
-                          terminate_coupled_remote, 1, MPI_INTEGER,          &
-                          target_id, 0,                                      &
-                          comm_inter, status, ierr )
+       IF ( myid == 0 ) THEN
+          CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,          &
+                             target_id, 0,                                      &
+                             terminate_coupled_remote, 1, MPI_INTEGER,          &
+                             target_id, 0,                                      &
+                             comm_inter, status, ierr )
+       ENDIF
+       CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)
 #endif
     ENDIF
@@ -140 +144 @@
              ENDIF
 #if defined( __parallel )
-             CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,    &
-                                target_id,  0,                               &
-                                terminate_coupled_remote, 1, MPI_INTEGER,    &
-                                target_id,  0,                               &
-                                comm_inter, status, ierr )
+             IF ( myid == 0 ) THEN
+                CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,    &
+                                   target_id,  0,                               &
+                                   terminate_coupled_remote, 1, MPI_INTEGER,    &
+                                   target_id,  0,                               &
+                                   comm_inter, status, ierr )   
+             ENDIF
+             CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)  
+            
 #endif
           ENDIF

palm/trunk/SOURCE/check_open.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of total array size was adapted to nbgp.
 !
 ! Former revisions:
@@ -278 +278 @@
 !--          Output for combine_plot_fields
              IF ( data_output_2d_on_each_pe  .AND.  myid_char /= '' )  THEN
-                WRITE (21)  -1, nx+1, -1, ny+1    ! total array size
+                WRITE (21)  -nbgp, nx+nbgp, -nbgp, ny+nbgp    ! total array size
                 WRITE (21)   0, nx+1,  0, ny+1    ! output part
              ENDIF
@@ -319 +319 @@
 !--          Output for combine_plot_fields
              IF ( data_output_2d_on_each_pe  .AND.  myid_char /= '' )  THEN
-                WRITE (22)  -1, nx+1, 0, nz+1    ! total array size
+                WRITE (22)  -nbgp, nx+nbgp, 0, nz+1    ! total array size
                 WRITE (22)   0, nx+1, 0, nz+1    ! output part
              ENDIF
@@ -357 +357 @@
 !--          Output for combine_plot_fields
              IF ( data_output_2d_on_each_pe  .AND.  myid_char /= '' )  THEN
-                WRITE (23)  -1, ny+1, 0, nz+1    ! total array size
+                WRITE (23)  -nbgp, ny+nbgp, 0, nz+1    ! total array size
                 WRITE (23)   0, ny+1, 0, nz+1    ! output part
              ENDIF
@@ -392 +392 @@
 !--          Specifications for combine_plot_fields
              IF ( .NOT. do3d_compress )  THEN
-                WRITE ( 30 )  -1,nx+1,-1,ny+1,0,nz_do3d
+                WRITE ( 30 )  -nbgp,nx+nbgp,-nbgp,ny+nbgp, 0 ,nz_do3d
                 WRITE ( 30 )  0,nx+1,0,ny+1,0,nz_do3d
              ENDIF
@@ -503 +503 @@
                 openfile(33)%opened = .TRUE.
                 WRITE ( 33, 3300 )  TRIM( avs_coor_file ), &
-                                    TRIM( avs_coor_file ), (nx+2)*4, &
-                                    TRIM( avs_coor_file ), (nx+2)*4+(ny+2)*4
+                                    TRIM( avs_coor_file ), (nx+2*nbgp)*4, &
+                                    TRIM( avs_coor_file ), (nx+2*nbgp)*4+(ny+2*nbgp)*4
             
 
@@ -548 +548 @@
 !--       corresponding partial array of a PE only once at the top of the file
           IF ( avs_output  .AND.  do3d_compress )  THEN
-             WRITE ( 30 )  nxl-1, nxr+1, nys-1, nyn+1, nzb, nz_do3d
+             WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb, nz_do3d
           ENDIF
 
@@ -929 +929 @@
 #endif
              ENDIF
-
+             
              CALL handle_netcdf_error( 'check_open', 26 )
 !
@@ -1325 +1325 @@
 #endif
              ENDIF
-
+             
              CALL handle_netcdf_error( 'check_open', 43 )  
 

palm/trunk/SOURCE/check_parameters.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
+!
 ! 
+! 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
+!
+! Adapt error messages.
 !
 ! Former revisions:
@@ -180 +193 @@
 !
 !-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
-    IF ( coupling_mode /= 'uncoupled' )  THEN
+    IF ( coupling_mode /= 'uncoupled')  THEN
 
        IF ( dt_coupling == 9999999.9 )  THEN
@@ -189 +202 @@
 
 #if defined( __parallel )
-       CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
-                      ierr )
-       CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &
-                      status, ierr )
+       IF ( myid == 0 ) THEN
+          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
+                         ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &
+                         status, ierr )
+       ENDIF
+       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
+       
        IF ( dt_coupling /= remote )  THEN
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -200 +217 @@
        ENDIF
        IF ( dt_coupling <= 0.0 )  THEN
-          CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
-          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &
-                         status, ierr )
+          IF ( myid == 0  ) THEN
+             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
+             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &
+                            status, ierr )
+          ENDIF   
+          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
+          
           dt_coupling = MAX( dt_max, remote )
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -209 +230 @@
           CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
        ENDIF
-
-       CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
-                      ierr )
-       CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &
-                      status, ierr )
+       IF ( myid == 0 ) THEN
+          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
+                         ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &
+                         status, ierr )
+       ENDIF 
+       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
+      
        IF ( restart_time /= remote )  THEN
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -220 +244 @@
           CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
        ENDIF
-
-       CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
-                      ierr )
-       CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &
-                      status, ierr )
+       IF ( myid == 0 ) THEN
+          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
+                         ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &
+                         status, ierr )
+       ENDIF   
+       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
+      
        IF ( dt_restart /= remote )  THEN
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -233 +260 @@
 
        simulation_time_since_reference = end_time - coupling_start_time
-       CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
-                      14, comm_inter, ierr )
-       CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
-                      status, ierr )
+       IF  ( myid == 0 ) THEN
+          CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
+                         14, comm_inter, ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
+                         status, ierr )    
+       ENDIF 
+       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
+      
        IF ( simulation_time_since_reference /= remote )  THEN
           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
@@ -245 +276 @@
        ENDIF
 
-       CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
-       CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &
-                      status, ierr )
-       IF ( dx /= remote )  THEN
-          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
-                 '":  dx = ', dx, '& is not equal to dx_remote = ', remote
-          CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
-       CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &
-                      status, ierr )
-       IF ( dy /= remote )  THEN
-          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
-                 '":  dy = ', dy, '& is not equal to dy_remote = ', remote
-          CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       CALL MPI_SEND( nx, 1, MPI_INTEGER, target_id, 17, comm_inter, ierr )
-       CALL MPI_RECV( iremote, 1, MPI_INTEGER, target_id, 17, comm_inter, &
-                      status, ierr )
-       IF ( nx /= iremote )  THEN
-          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
-                 '": nx = ', nx, '& is not equal to nx_remote = ', iremote
-          CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
-       ENDIF
-
-       CALL MPI_SEND( ny, 1, MPI_INTEGER, target_id, 18, comm_inter, ierr )
-       CALL MPI_RECV( iremote, 1, MPI_INTEGER, target_id, 18, comm_inter, &
-                      status, ierr )
-       IF ( ny /= iremote )  THEN
-          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
-                 '": ny = ', ny, '& is not equal to ny_remote = ', iremote
-          CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
+  
+       IF ( myid == 0 ) THEN
+          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &
+                                                             status, ierr )
+       ENDIF
+       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
+
+
+       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
+
+          IF ( dx < remote ) THEN
+             WRITE( message_string, * ) 'coupling mode "', &
+                   TRIM( coupling_mode ),                  &
+           '": dx in Atmosphere is not equal to or not larger then dx in ocean'
+             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
+             WRITE( message_string, * ) 'coupling mode "', &
+                    TRIM( coupling_mode ), &
+             '": Domain size in x-direction is not equal in ocean and atmosphere'
+             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
+          ENDIF
+
+       ENDIF
+
+       IF ( myid == 0) THEN
+          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
+          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &
+                         status, ierr )
+       ENDIF
+       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
+
+       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
+
+          IF ( dy < remote )  THEN
+             WRITE( message_string, * ) 'coupling mode "', &
+                    TRIM( coupling_mode ), &
+                 '": dy in Atmosphere is not equal to or not larger then dy in ocean'
+             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
+             WRITE( message_string, * ) 'coupling mode "', &
+                   TRIM( coupling_mode ), &
+             '": Domain size in y-direction is not equal in ocean and atmosphere'
+             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
+             WRITE( message_string, * ) 'coupling mode "', &
+                   TRIM( coupling_mode ), &
+             '": nx+1 in ocean is not divisible without remainder with nx+1 in', & 
+             ' atmosphere'
+             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
+          ENDIF
+
+          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
+             WRITE( message_string, * ) 'coupling mode "', &
+                   TRIM( coupling_mode ), &
+             '": ny+1 in ocean is not divisible without remainder with ny+1 in', & 
+             ' atmosphere'
+             CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
+          ENDIF
+
        ENDIF
 #else
@@ -290 +353 @@
 !
 !-- Exchange via intercommunicator
-    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
+    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
        CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, &
                       ierr )
-    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
+    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
        CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, &
                       comm_inter, status, ierr )
     ENDIF
+    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
+    
 #endif
 
@@ -372 +437 @@
           CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
        ENDIF
+       IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' ) &
+       THEN
+          message_string = 'topography is still not allowed with ' // &
+                           'momentum_advec = "' // TRIM( momentum_advec ) //  &
+                           '"or scalar_advec = "' // TRIM( scalar_advec ) //'"'
+   ! message number still needs modification
+           CALL message( 'check_parameters', 'PA0341', 1, 2, 0, 6, 0 )
+       END IF
+          
 !
 !--    In case of non-flat topography, check whether the convention how to 
@@ -492 +566 @@
        CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
     ENDIF
-
+    
+    IF( momentum_advec == 'ws-scheme' .AND. & 
+        call_psolver_at_all_substeps == .FALSE. ) THEN
+        message_string = 'psolver must be called at each RK3 substep when "'//&
+                      TRIM(momentum_advec) // ' "is used for momentum_advec'
+        CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
+    END IF
 !
 !-- Advection schemes:
-    IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ups-scheme' ) &
-    THEN
+    IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' .AND. &
+         momentum_advec /= 'ups-scheme' ) THEN
        message_string = 'unknown advection scheme: momentum_advec = "' // &
                         TRIM( momentum_advec ) // '"'
        CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
     ENDIF
-    IF ( ( momentum_advec == 'ups-scheme'  .OR.  scalar_advec == 'ups-scheme' )&
-                                      .AND.  timestep_scheme /= 'euler' )  THEN
-       message_string = 'momentum_advec = "' // TRIM( momentum_advec ) // &
-                        '" is not allowed with timestep_scheme = "' //    &
-                        TRIM( timestep_scheme ) // '"'
+    IF ((( momentum_advec == 'ups-scheme'  .OR.  scalar_advec == 'ups-scheme' )&
+           .AND.  timestep_scheme /= 'euler' ) .OR. (( momentum_advec == 'ws-scheme'&
+           .OR.  scalar_advec == 'ws-scheme') .AND. (timestep_scheme == 'euler' .OR. &
+           timestep_scheme == 'leapfrog+euler' .OR. timestep_scheme == 'leapfrog'    &
+           .OR. timestep_scheme == 'runge-kutta-2'))) THEN
+       message_string = 'momentum_advec or scalar_advec = "' &
+         // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &
+         TRIM( timestep_scheme ) // '"'
        CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
     ENDIF
-
-    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'bc-scheme'  .AND.&
-         scalar_advec /= 'ups-scheme' )  THEN
+    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
+        scalar_advec /= 'bc-scheme'  .AND.  scalar_advec /= 'ups-scheme' )  THEN
        message_string = 'unknown advection scheme: scalar_advec = "' // &
                         TRIM( scalar_advec ) // '"'
@@ -563 +645 @@
     ENDIF
 
-    IF ( momentum_advec /= 'pw-scheme' .AND. timestep_scheme(1:5) == 'runge' ) &
-    THEN
+    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme') &
+         .AND. timestep_scheme(1:5) == 'runge' ) THEN
        message_string = 'momentum advection scheme "' // &
                         TRIM( momentum_advec ) // '" & does not work with ' // &
@@ -712 +794 @@
           ug_vertical_gradient_level_ind(1) = nzt+1
           ug(nzt+1) = ug_surface
-          DO  k = nzt, 0, -1
+          DO  k = nzt, nzb, -1
              IF ( i < 11 ) THEN
                 IF ( ug_vertical_gradient_level(i) > zu(k)  .AND. &
@@ -778 +860 @@
           vg_vertical_gradient_level_ind(1) = nzt+1
           vg(nzt+1) = vg_surface
-          DO  k = nzt, 0, -1
+          DO  k = nzt, nzb, -1
              IF ( i < 11 ) THEN
                 IF ( vg_vertical_gradient_level(i) > zu(k)  .AND. &
@@ -1020 +1102 @@
  
              
+
 !
 !-- Compute Coriolis parameter
@@ -1159 +1242 @@
 !
 !-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
-!-- Willimas advection scheme. Several schemes and tools do not work with
-!-- non-cyclic boundary conditions.
+!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
+!-- and tools do not work with non-cyclic boundary conditions.
     IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
        IF ( psolver /= 'multigrid' )  THEN
@@ -1167 +1250 @@
           CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
        ENDIF
-       IF ( momentum_advec /= 'pw-scheme' )  THEN
+       IF ( momentum_advec /= 'pw-scheme' .AND. &
+            momentum_advec /= 'ws-scheme')  THEN
           message_string = 'non-cyclic lateral boundaries do not allow ' // &
                            'momentum_advec = "' // TRIM( momentum_advec ) // '"'
           CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
        ENDIF
-       IF ( scalar_advec /= 'pw-scheme' )  THEN
+       IF ( scalar_advec /= 'pw-scheme' .AND. &
+            scalar_advec /= 'ws-scheme' )  THEN
           message_string = 'non-cyclic lateral boundaries do not allow ' // &
                            'scalar_advec = "' // TRIM( scalar_advec ) // '"'
           CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
        ENDIF
+       IF ( (scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' ) &
+          .AND. loop_optimization == 'vector' ) THEN
+          message_string = 'non-cyclic lateral boundaries do not allow ' // &
+                           'loop_optimization = vector and ' //  &
+                           'scalar_advec = "' // TRIM( scalar_advec ) // '"' 
+  ! The error message number still needs modification. 
+          CALL message( 'check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
+       END IF
        IF ( galilei_transformation )  THEN
           message_string = 'non-cyclic lateral boundaries do not allow ' // &
@@ -1407 +1500 @@
                         TRIM( bc_uv_b ) // '"'
        CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
+    ENDIF
+!
+!-- In case of coupled simulations u and v at the ground in atmosphere will be
+!-- assigned with the u and v values of the ocean surface
+    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
+       ibc_uv_b = 2
     ENDIF
 
@@ -2109 +2208 @@
              hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
 
+
           CASE ( 'u"pt"' )
              dopr_index(i) = 58
@@ -2244 +2344 @@
 
        END SELECT
+
 !
 !--    Check to which of the predefined coordinate systems the profile belongs
@@ -2584 +2685 @@
 !-- Upper plot limit (grid point value) for 1D profiles
     IF ( z_max_do1d == -1.0 )  THEN
+
        nz_do1d = nzt+1
+
     ELSE
        DO  k = nzb+1, nzt+1
@@ -2737 +2840 @@
 
 !
+
 !-- Check netcdf precison
     ldum = .FALSE.
@@ -3070 +3174 @@
     IF ( conserve_volume_flow )  THEN
        IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
-          IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
-             conserve_volume_flow_mode = 'inflow_profile'
-          ELSE
-             conserve_volume_flow_mode = 'initial_profiles'
-          ENDIF
+
+          conserve_volume_flow_mode = 'initial_profiles'
+
        ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
             TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' .AND.  &
@@ -3082 +3184 @@
           CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
        ENDIF
-       IF ( ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' ) .AND. &
-            TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' )  THEN
-          WRITE( message_string, * )  'noncyclic boundary conditions ', &
-               'require & conserve_volume_flow_mode = ''inflow_profile'''
+       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND. &
+          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
+          WRITE( message_string, * )  'non-cyclic boundary conditions ', &
+               'require  conserve_volume_flow_mode = ''initial_profiles'''
           CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
        ENDIF
@@ -3091 +3193 @@
             TRIM( conserve_volume_flow_mode ) == 'inflow_profile' )  THEN
           WRITE( message_string, * )  'cyclic boundary conditions ', &
-               'require & conserve_volume_flow_mode = ''initial_profiles''', &
+               'require conserve_volume_flow_mode = ''initial_profiles''', &
                ' or ''bulk_velocity'''
           CALL message( 'check_parameters', 'PA0156', 1, 2, 0, 6, 0 )
@@ -3100 +3202 @@
          TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
        WRITE( message_string, * )  'nonzero bulk velocity requires ', &
-            'conserve_volume_flow = .T. and & ', &
+            'conserve_volume_flow = .T. and ', &
             'conserve_volume_flow_mode = ''bulk_velocity'''
        CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
@@ -3139 +3241 @@
 
 
+
  END SUBROUTINE check_parameters

palm/trunk/SOURCE/data_output_2d.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! 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 modiefied.
 !
 ! Former revisions:
@@ -112 +114 @@
 
        CASE ( 'xy' )
-
           s = 1
-          ALLOCATE( level_z(0:nzt+1), local_2d(nxl-1:nxr+1,nys-1:nyn+1) )
+          ALLOCATE( level_z(nzb:nzt+1), local_2d(nxlg:nxrg,nysg:nyng) )
 
 !
@@ -130 +131 @@
                 IF ( iso2d_output )  CALL check_open( 21 )
 #if defined( __parallel )
-                ALLOCATE( total_2d(-1:nx+1,-1:ny+1) )
+                ALLOCATE( total_2d(-nbgp:nx+nbgp,-nbgp:ny+nbgp) )
 #endif
              ENDIF
@@ -136 +137 @@
 
        CASE ( 'xz' )
-
           s = 2
-          ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) )
+          ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
 
 !
@@ -154 +154 @@
                 IF ( iso2d_output )  CALL check_open( 22 )
 #if defined( __parallel )
-                ALLOCATE( total_2d(-1:nx+1,nzb:nzt+1) )
+                ALLOCATE( total_2d(-nbgp:nx+nbgp,nzb:nzt+1) )
 #endif
              ENDIF
@@ -162 +162 @@
 
           s = 3
-          ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) )
+          ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
 
 !
@@ -178 +178 @@
                 IF ( iso2d_output )  CALL check_open( 23 )
 #if defined( __parallel )
-                ALLOCATE( total_2d(-1:ny+1,nzb:nzt+1) )
+                ALLOCATE( total_2d(-nbgp:ny+nbgp,nzb:nzt+1) )
 #endif
              ENDIF
@@ -192 +192 @@
 !
 !-- Allocate a temporary array for resorting (kji -> ijk).
-    ALLOCATE( local_pf(nxl-1:nxr+1,nys-1:nyn+1,nzb:nzt+1) )
+    ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nzt+1) )
 
 !
@@ -219 +219 @@
              CASE ( 'lwp*_xy' )        ! 2d-array
                 IF ( av == 0 )  THEN
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) * &
                                                     dzw(1:nzt+1) )
@@ -226 +226 @@
                    ENDDO
                 ELSE
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) = lwp_av(j,i)
                       ENDDO
@@ -248 +248 @@
                    IF ( simulated_time >= particle_advection_start )  THEN
                       tend = prt_count
-                      CALL exchange_horiz( tend )
+                      CALL exchange_horiz( tend, nbgp )
                    ELSE
                       tend = 0.0
                    ENDIF
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          DO  k = nzb, nzt+1
                             local_pf(i,j,k) = tend(k,j,i)
@@ -261 +261 @@
                    resorted = .TRUE.
                 ELSE
-                   CALL exchange_horiz( pc_av )
+                   CALL exchange_horiz( pc_av, nbgp )
                    to_be_resorted => pc_av
                 ENDIF
@@ -287 +287 @@
                          ENDDO
                       ENDDO
-                      CALL exchange_horiz( tend )
+                      CALL exchange_horiz( tend, nbgp )
                    ELSE
                       tend = 0.0
-                   ENDIF
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   END IF
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          DO  k = nzb, nzt+1
                             local_pf(i,j,k) = tend(k,j,i)
@@ -300 +300 @@
                    resorted = .TRUE.
                 ELSE
-                   CALL exchange_horiz( pr_av )
+                   CALL exchange_horiz( pr_av, nbgp )
                    to_be_resorted => pr_av
                 ENDIF
@@ -306 +306 @@
              CASE ( 'pra*_xy' )        ! 2d-array / integral quantity => no av
                 CALL exchange_horiz_2d( precipitation_amount )
-                DO  i = nxl-1, nxr+1
-                   DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                       local_pf(i,j,nzb+1) =  precipitation_amount(j,i)
                    ENDDO
@@ -319 +319 @@
                 IF ( av == 0 )  THEN
                    CALL exchange_horiz_2d( precipitation_rate )
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1 
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  precipitation_rate(j,i)
                       ENDDO
@@ -326 +326 @@
                 ELSE
                    CALL exchange_horiz_2d( precipitation_rate_av )
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1 
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  precipitation_rate_av(j,i)
                       ENDDO
@@ -341 +341 @@
                       to_be_resorted => pt
                    ELSE
-                      DO  i = nxl-1, nxr+1
-                         DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                             DO  k = nzb, nzt+1
                                local_pf(i,j,k) = pt(k,j,i) + l_d_cp *    &
@@ -399 +399 @@
              CASE ( 'qsws*_xy' )        ! 2d-array
                 IF ( av == 0 ) THEN
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1 
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  qsws(j,i)
                       ENDDO
                    ENDDO
                 ELSE
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1 
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng 
                          local_pf(i,j,nzb+1) =  qsws_av(j,i)
                       ENDDO
@@ -417 +417 @@
              CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
                 IF ( av == 0 )  THEN
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          DO  k = nzb, nzt+1
                             local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
@@ -453 +453 @@
              CASE ( 'shf*_xy' )        ! 2d-array
                 IF ( av == 0 ) THEN
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1 
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  shf(j,i)
                       ENDDO
                    ENDDO
                 ELSE
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1 
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  shf_av(j,i)
                       ENDDO
@@ -471 +471 @@
              CASE ( 't*_xy' )        ! 2d-array
                 IF ( av == 0 )  THEN
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) = ts(j,i)
                       ENDDO
                    ENDDO
                 ELSE
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) = ts_av(j,i)
                       ENDDO
@@ -503 +503 @@
              CASE ( 'u*_xy' )        ! 2d-array
                 IF ( av == 0 )  THEN
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) = us(j,i)
                       ENDDO
                    ENDDO
                 ELSE
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) = us_av(j,i)
                       ENDDO
@@ -551 +551 @@
              CASE ( 'z0*_xy' )        ! 2d-array
                 IF ( av == 0 ) THEN
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1 
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  z0(j,i)
                       ENDDO
                    ENDDO
                 ELSE
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1 
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          local_pf(i,j,nzb+1) =  z0_av(j,i)
                       ENDDO
@@ -593 +593 @@
 !--       Resort the array to be output, if not done above
           IF ( .NOT. resorted )  THEN
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
                       local_pf(i,j,k) = to_be_resorted(k,j,i)
@@ -647 +647 @@
 !--                   Carry out the averaging (all data are on the PE)
                       DO  k = nzb, nzt+1
-                         DO  j = nys-1, nyn+1
-                            DO  i = nxl-1, nxr+1
+                         DO  j = nysg, nyng
+                            DO  i = nxlg, nxrg
                                local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
                             ENDDO
@@ -654 +654 @@
                       ENDDO
 
-                      local_2d = local_2d / ( nzt -nzb + 2.0 )
+                      local_2d = local_2d / ( nzt -nzb + 2.0)
 
                    ELSE
@@ -723 +723 @@
                          ENDIF
 #endif
-                         WRITE ( 21 )  nxl-1, nxr+1, nys-1, nyn+1
+                         WRITE ( 21 )  nxlg, nxrg, nysg, nyng
                          WRITE ( 21 )  local_2d
 
@@ -734 +734 @@
                          CALL MPI_BARRIER( comm2d, ierr )
 
-                         ngp = ( nxr-nxl+3 ) * ( nyn-nys+3 )
+                         ngp = ( nxrg-nxlg+1 ) * ( nyng-nysg+1 )
                          IF ( myid == 0 )  THEN
 !
 !--                         Local array can be relocated directly.
-                            total_2d(nxl-1:nxr+1,nys-1:nyn+1) = local_2d
+                            total_2d(nxlg:nxrg,nysg:nyng) = local_2d
 !
 !--                         Receive data from all other PEs.
@@ -760 +760 @@
 !--                         Output of the total cross-section.
                             IF ( iso2d_output )  THEN
-                               WRITE (21)  total_2d(0:nx+1,0:ny+1)
+                               WRITE (21)  total_2d(-nbgp:nx+nbgp,-nbgp:ny+nbgp)
                             ENDIF
 !
 !--                         Relocate the local array for the next loop increment
                             DEALLOCATE( local_2d )
-                            ALLOCATE( local_2d(nxl-1:nxr+1,nys-1:nyn+1) )
+                            ALLOCATE( local_2d(nxlg:nxrg,nysg:nyng) )
 
 #if defined( __netcdf )
@@ -789 +789 @@
 !
 !--                         First send the local index limits to PE0
-                            ind(1) = nxl-1; ind(2) = nxr+1
-                            ind(3) = nys-1; ind(4) = nyn+1
+                            ind(1) = nxlg; ind(2) = nxrg
+                            ind(3) = nysg; ind(4) = nyng
                             CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
                                            comm2d, ierr )
 !
 !--                         Send data to PE0
-                            CALL MPI_SEND( local_2d(nxl-1,nys-1), ngp, &
+                            CALL MPI_SEND( local_2d(nxlg,nysg), ngp, &
                                            MPI_REAL, 0, 1, comm2d, ierr )
                          ENDIF
@@ -882 +882 @@
 
                    ENDIF
+
 !
 !--                If required, carry out averaging along y
                    IF ( section(is,s) == -1 )  THEN
 
-                      ALLOCATE( local_2d_l(nxl-1:nxr+1,nzb:nzt+1) )
+                      ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
                       local_2d_l = 0.0
-                      ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
+                      ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
 !
 !--                   First local averaging on the PE
                       DO  k = nzb, nzt+1
                          DO  j = nys, nyn
-                            DO  i = nxl-1, nxr+1
+                            DO  i = nxlg, nxrg
                                local_2d_l(i,k) = local_2d_l(i,k) + &
                                                  local_pf(i,j,k)
@@ -903 +904 @@
 !--                   Now do the averaging over all PEs along y
                       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-                      CALL MPI_ALLREDUCE( local_2d_l(nxl-1,nzb),              &
-                                          local_2d(nxl-1,nzb), ngp, MPI_REAL, &
+                      CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb),              &
+                                          local_2d(nxlg,nzb), ngp, MPI_REAL, &
                                           MPI_SUM, comm1dy, ierr )
 #else
@@ -936 +937 @@
 !--                   BEGIN WORKAROUND---------------------------------------
                       IF ( npey /= 1  .AND.  section(is,s) /= -1)  THEN
-                         ALLOCATE( local_2d_l(nxl-1:nxr+1,nzb:nzt+1) )
+                         ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
                          local_2d_l = 0.0
                          IF ( section(is,s) >= nys .AND. section(is,s) <= nyn )&
@@ -945 +946 @@
 !
 !--                      Distribute data over all PEs along y
-                         ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
+                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
                          IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
-                         CALL MPI_ALLREDUCE( local_2d_l(nxl-1,nzb),            &
-                                             local_2d(nxl-1,nzb), ngp,         &
+                         CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb),            &
+                                             local_2d(nxlg,nzb), ngp,         &
                                              MPI_REAL, MPI_SUM, comm1dy, ierr )
 #else
@@ -1022 +1023 @@
                               ( section(is,s) == -1  .AND.  nys-1 == -1 ) )  &
                          THEN
-                            WRITE (22)  nxl-1, nxr+1, nzb, nzt+1
+                            WRITE (22)  nxlg, nxrg, nzb, nzt+1
                             WRITE (22)  local_2d
                          ELSE
@@ -1036 +1037 @@
                          CALL MPI_BARRIER( comm2d, ierr )
 
-                         ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
+                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
                          IF ( myid == 0 )  THEN
 !
@@ -1044 +1045 @@
                                  ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
                             THEN
-                               total_2d(nxl-1:nxr+1,nzb:nzt+1) = local_2d
+                               total_2d(nxlg:nxrg,nzb:nzt+1) = local_2d
                             ENDIF
 !
@@ -1073 +1074 @@
 !--                         Output of the total cross-section.
                             IF ( iso2d_output )  THEN
-                               WRITE (22)  total_2d(0:nx+1,nzb:nzt+1)
+                               WRITE (22)  total_2d(-nbgp:nx+nbgp,nzb:nzt+1)
                             ENDIF
 !
 !--                         Relocate the local array for the next loop increment
                             DEALLOCATE( local_2d )
-                            ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) )
+                            ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
 
 #if defined( __netcdf )
@@ -1099 +1100 @@
                                  ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
                             THEN
-                               ind(1) = nxl-1; ind(2) = nxr+1
+                               ind(1) = nxlg; ind(2) = nxrg
                                ind(3) = nzb;   ind(4) = nzt+1
                             ELSE
@@ -1110 +1111 @@
 !--                         If applicable, send data to PE0.
                             IF ( ind(1) /= -9999 )  THEN
-                               CALL MPI_SEND( local_2d(nxl-1,nzb), ngp, &
+                               CALL MPI_SEND( local_2d(nxlg,nzb), ngp, &
                                               MPI_REAL, 0, 1, comm2d, ierr )
                             ENDIF
@@ -1187 +1188 @@
                    IF ( section(is,s) == -1 )  THEN
 
-                      ALLOCATE( local_2d_l(nys-1:nyn+1,nzb:nzt+1) )
+                      ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
                       local_2d_l = 0.0
-                      ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
+                      ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
 !
 !--                   First local averaging on the PE
                       DO  k = nzb, nzt+1
-                         DO  j = nys-1, nyn+1
+                         DO  j = nysg, nyng
                             DO  i = nxl, nxr
                                local_2d_l(j,k) = local_2d_l(j,k) + &
@@ -1204 +1205 @@
 !--                   Now do the averaging over all PEs along x
                       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-                      CALL MPI_ALLREDUCE( local_2d_l(nys-1,nzb),              &
-                                          local_2d(nys-1,nzb), ngp, MPI_REAL, &
+                      CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb),              &
+                                          local_2d(nysg,nzb), ngp, MPI_REAL, &
                                           MPI_SUM, comm1dx, ierr )
 #else
@@ -1237 +1238 @@
 !--                   BEGIN WORKAROUND---------------------------------------
                       IF ( npex /= 1  .AND.  section(is,s) /= -1)  THEN
-                         ALLOCATE( local_2d_l(nys-1:nyn+1,nzb:nzt+1) )
+                         ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
                          local_2d_l = 0.0
                          IF ( section(is,s) >= nxl .AND. section(is,s) <= nxr )&
@@ -1246 +1247 @@
 !
 !--                      Distribute data over all PEs along x
-                         ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
+                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 )
                          IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
-                         CALL MPI_ALLREDUCE( local_2d_l(nys-1,nzb),            &
-                                             local_2d(nys-1,nzb), ngp,         &
+                         CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb),            &
+                                             local_2d(nysg,nzb), ngp,         &
                                              MPI_REAL, MPI_SUM, comm1dx, ierr )
 #else
@@ -1323 +1324 @@
                               ( section(is,s) == -1  .AND.  nxl-1 == -1 ) )  &
                          THEN
-                            WRITE (23)  nys-1, nyn+1, nzb, nzt+1
+                            WRITE (23)  nysg, nyng, nzb, nzt+1
                             WRITE (23)  local_2d
                          ELSE
@@ -1337 +1338 @@
                          CALL MPI_BARRIER( comm2d, ierr )
 
-                         ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
+                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
                          IF ( myid == 0 )  THEN
 !
@@ -1345 +1346 @@
                                  ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
                             THEN
-                               total_2d(nys-1:nyn+1,nzb:nzt+1) = local_2d
+                               total_2d(nysg:nyng,nzb:nzt+1) = local_2d
                             ENDIF
 !
@@ -1379 +1380 @@
 !--                         Relocate the local array for the next loop increment
                             DEALLOCATE( local_2d )
-                            ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) )
+                            ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
 
 #if defined( __netcdf )
@@ -1400 +1401 @@
                                  ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
                             THEN
-                               ind(1) = nys-1; ind(2) = nyn+1
+                               ind(1) = nysg; ind(2) = nyng
                                ind(3) = nzb;   ind(4) = nzt+1
                             ELSE
@@ -1411 +1412 @@
 !--                         If applicable, send data to PE0.
                             IF ( ind(1) /= -9999 )  THEN
-                               CALL MPI_SEND( local_2d(nys-1,nzb), ngp, &
+                               CALL MPI_SEND( local_2d(nysg,nzb), ngp, &
                                               MPI_REAL, 0, 1, comm2d, ierr )
                             ENDIF

palm/trunk/SOURCE/data_output_3d.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! 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.
 !
 ! Former revisions:
@@ -102 +104 @@
 !
 !-- Allocate a temporary array with the desired output dimensions.
-    ALLOCATE( local_pf(nxl-1:nxr+1,nys-1:nyn+1,nzb:nz_do3d) )
+    ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) )
 
 !
@@ -157 +159 @@
              IF ( av == 0 )  THEN
                 tend = prt_count
-                CALL exchange_horiz( tend )
-                DO  i = nxl-1, nxr+1
-                   DO  j = nys-1, nyn+1
+                CALL exchange_horiz( tend, nbgp )
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
                       DO  k = nzb, nz_do3d
                          local_pf(i,j,k) = tend(k,j,i)
@@ -167 +169 @@
                 resorted = .TRUE.
              ELSE
-                CALL exchange_horiz( pc_av )
+                CALL exchange_horiz( pc_av, nbgp )
                 to_be_resorted => pc_av
              ENDIF
@@ -192 +194 @@
                    ENDDO
                 ENDDO
-                CALL exchange_horiz( tend )
-                DO  i = nxl-1, nxr+1
-                   DO  j = nys-1, nyn+1
+                CALL exchange_horiz( tend, nbgp )
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
                          local_pf(i,j,k) = tend(k,j,i)
@@ -202 +204 @@
                 resorted = .TRUE.
              ELSE
-                CALL exchange_horiz( pr_av )
+                CALL exchange_horiz( pr_av, nbgp )
                 to_be_resorted => pr_av
              ENDIF
@@ -211 +213 @@
                    to_be_resorted => pt
                 ELSE
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          DO  k = nzb, nz_do3d
                             local_pf(i,j,k) = pt(k,j,i) + l_d_cp *    &
@@ -263 +265 @@
           CASE ( 'qv' )
              IF ( av == 0 )  THEN
-                DO  i = nxl-1, nxr+1
-                   DO  j = nys-1, nyn+1
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
                       DO  k = nzb, nz_do3d
                          local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
@@ -342 +344 @@
 !--    Resort the array to be output, if not done above
        IF ( .NOT. resorted )  THEN
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 DO  k = nzb, nz_do3d
                    local_pf(i,j,k) = to_be_resorted(k,j,i)
@@ -376 +378 @@
 !--       Determine the Skip-value for the next array. Record end and start 
 !--       require 4 byte each.
-          skip_do_avs = skip_do_avs + ( ((nx+2)*(ny+2)*(nz_do3d+1)) * 4 + 8 )
+          skip_do_avs = skip_do_avs + ( ((nx+2*nbgp)*(ny+2*nbgp)*(nz_do3d+1)) * 4 + 8 )
        ENDIF
 
@@ -386 +388 @@
 !--       of compressed data.
           CALL write_compressed( local_pf, 30, 33, myid, nxl, nxr, nyn, nys, &
-                                 nzb, nz_do3d, prec )
+                                 nzb, nz_do3d, prec, nbgp )
        ELSE
 !
@@ -400 +402 @@
                    WRITE ( 30 )  simulated_time, do3d_time_count(av), av
                 ENDIF
-                WRITE ( 30 )  nxl-1, nxr+1, nys-1, nyn+1, nzb, nz_do3d
+                WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb, nz_do3d
                 WRITE ( 30 )  local_pf
 
@@ -411 +413 @@
                 IF ( nxr == nx  .AND.  nyn /= ny )  THEN
                    nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
-                                  local_pf(nxl:nxr+1,nys:nyn,nzb:nz_do3d),    &
+                                  local_pf(nxl:nxrg,nys:nyn,nzb:nz_do3d),    &
                       start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
-                      count = (/ nxr-nxl+2, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
+                      count = (/ nxr-nxl+1+nbgp, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
                 ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
                    nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
-                                  local_pf(nxl:nxr,nys:nyn+1,nzb:nz_do3d),    &
+                                  local_pf(nxl:nxr,nys:nyng,nzb:nz_do3d),    &
                       start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
-                      count = (/ nxr-nxl+1, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
+                      count = (/ nxr-nxl+1, nyn-nys+1+nbgp, nz_do3d-nzb+1, 1 /) )
                 ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
                    nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
-                                  local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),  &
+                                  local_pf(nxl:nxrg,nys:nyng,nzb:nz_do3d),  &
                       start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
-                      count = (/ nxr-nxl+2, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
+                      count = (/ nxr-nxl+1+nbgp, nyn-nys+1+nbgp, nz_do3d-nzb+1, 1 /) )
                 ELSE
                    nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &

palm/trunk/SOURCE/data_output_mask.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! Calls of exchange_horiz are modified.
 !
 ! Former revisions:
@@ -123 +123 @@
              IF ( av == 0 )  THEN
                 tend = prt_count
-                CALL exchange_horiz( tend )
+                CALL exchange_horiz( tend, nbgp )
                 DO  i = 1, mask_size_l(mid,1)
                    DO  j = 1, mask_size_l(mid,2)
@@ -134 +134 @@
                 resorted = .TRUE.
              ELSE
-                CALL exchange_horiz( pc_av )
+                CALL exchange_horiz( pc_av, nbgp )
                 to_be_resorted => pc_av
              ENDIF
@@ -159 +159 @@
                    ENDDO
                 ENDDO
-                CALL exchange_horiz( tend )
+                CALL exchange_horiz( tend, nbgp )
                 DO  i = 1, mask_size_l(mid,1)
                    DO  j = 1, mask_size_l(mid,2)
@@ -170 +170 @@
                 resorted = .TRUE.
              ELSE
-                CALL exchange_horiz( pr_av )
+                CALL exchange_horiz( pr_av, nbgp )
                 to_be_resorted => pr_av
              ENDIF
@@ -439 +439 @@
 
        if = if + 1
+
     ENDDO
 

palm/trunk/SOURCE/diffusion_e.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
 ! Former revisions:
@@ -68 +68 @@
        REAL            ::  dvar_dz, l_stable, phi_m, var_reference
        REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
-                           l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1)
-       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
+                           l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
+       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
        REAL, DIMENSION(:,:), POINTER   ::  rif
        REAL, DIMENSION(:,:,:), POINTER ::  e, km, var
@@ -289 +289 @@
        REAL            ::  dvar_dz, l_stable, phi_m, var_reference
        REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
-                           l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1)
-       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  diss, tend
+                           l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
+       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
        REAL, DIMENSION(:,:), POINTER   ::  rif
        REAL, DIMENSION(:,:,:), POINTER ::  e, km, var

palm/trunk/SOURCE/diffusion_s.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
 ! Former revisions:
@@ -64 +64 @@
        REAL    ::  vertical_gridspace
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
-       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL    ::  wall_s_flux(0:4)
        REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t
@@ -176 +176 @@
        REAL    ::  vertical_gridspace
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
-       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL    ::  wall_s_flux(0:4)
        REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t

palm/trunk/SOURCE/diffusion_u.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
 ! Former revisions:
@@ -72 +72 @@
        INTEGER ::  i, j, k
        REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
-       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)
+       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -177 +177 @@
                       &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
                       &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
-                      &            )                                         &
+                      &            )                                          &
                       &   ) * ddzw(k)
              ENDDO
@@ -248 +248 @@
        INTEGER ::  i, j, k
        REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
-       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)
+       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(nzb:nzt+1)      ::  usvs
        REAL, DIMENSION(:,:),   POINTER ::  usws, uswst

palm/trunk/SOURCE/diffusion_v.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
 ! Former revisions:
@@ -70 +70 @@
        INTEGER ::  i, j, k
        REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
-       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)
+       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -246 +246 @@
        INTEGER ::  i, j, k
        REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
-       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)
+       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(nzb:nzt+1)      ::  vsus
        REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst

palm/trunk/SOURCE/diffusion_w.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
 ! Former revisions:
@@ -64 +64 @@
        REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
                    kmyp_z
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
-                   km_damp_y(nys-1:nyn+1)
-       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        & 
+                   km_damp_y(nysg:nyng)
+       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
        REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
@@ -211 +211 @@
        REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
                    kmyp_z
-       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
-                   km_damp_y(nys-1:nyn+1)
-       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        & 
+                   km_damp_y(nysg:nyng)
+       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
        REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w

palm/trunk/SOURCE/diffusivities.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 ! 
 ! Former revisions:
@@ -54 +54 @@
     REAL, SAVE ::  phi_m = 1.0
 
-    REAL    ::  var(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+    REAL    ::  var(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
 
     REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e
@@ -73 +73 @@
 
     !$OMP DO
-    DO  i = nxl-1, nxr+1
-       DO  j = nys-1, nyn+1
+    DO  i = nxlg, nxrg
+       DO  j = nysg, nyng
 
 !
@@ -157 +157 @@
 
     sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
-                                                ! data output
+                                                  ! data output
 
     !$OMP END PARALLEL
@@ -167 +167 @@
 !-- values of the diffusivities are not needed
     !$OMP PARALLEL DO
-    DO  i = nxl-1, nxr+1
-       DO  j = nys-1, nyn+1
+    DO  i = nxlg, nxrg
+       DO  j = nysg, nyng
           km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
           km(nzt+1,j,i)            = km(nzt,j,i)

palm/trunk/SOURCE/disturb_field.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
+! Calls of exchange_horiz are modified.
 !
 ! Former revisions:
@@ -44 +45 @@
 
     INTEGER ::  i, j, k
-    INTEGER ::  nzb_uv_inner(nys-1:nyn+1,nxl-1:nxr+1)
+    INTEGER ::  nzb_uv_inner(nysg:nyng,nxlg:nxrg)
 
     REAL    ::  randomnumber,                             &
-                dist1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                field(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+                dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+                field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::  dist2
 
@@ -57 +58 @@
 !-- Create an additional temporary array and initialize the arrays needed
 !-- to store the disturbance
-    ALLOCATE( dist2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     dist1 = 0.0
     dist2 = 0.0
@@ -102 +103 @@
 !-- Exchange of ghost points for the random perturbation
 
-    CALL exchange_horiz( dist1 )
-
+    CALL exchange_horiz( dist1, nbgp )
 !
 !-- Applying the Shuman filter in order to smooth the perturbations.
@@ -128 +128 @@
 !-- Exchange of ghost points for the filtered perturbation.
 !-- Afterwards, filter operation and exchange of ghost points are repeated.
-    CALL exchange_horiz( dist2 )
+    CALL exchange_horiz( dist2, nbgp )
 
     DO  i = nxl, nxr
@@ -141 +141 @@
     ENDDO
 
-    CALL exchange_horiz( dist1 )
+    CALL exchange_horiz( dist1, nbgp )
 
 !
@@ -148 +148 @@
 !-- (diffusion criterion))
     IF ( TRIM( topography ) /= 'flat' )  THEN
-       DO  i = nxl-1, nxr+1
-          DO  j = nys-1, nyn+1
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              dist1(nzb:nzb_uv_inner(j,i)+1,j,i) = 0.0
           ENDDO
@@ -157 +157 @@
 !
 !-- Random perturbation is added to the array to be disturbed.
-    DO  i = nxl-1, nxr+1
-       DO  j = nys-1, nyn+1
+    DO  i = nxlg, nxrg
+       DO  j = nysg, nyng
           DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
              field(k,j,i) = field(k,j,i) + dist1(k,j,i)

palm/trunk/SOURCE/exchange_horiz.f90

@@ -1 +Â @@
- SUBROUTINE exchange_horiz( ar )
+ SUBROUTINE exchange_horiz( ar, nbgp_local)
 
 !------------------------------------------------------------------------------!
 ! Current revisions:
 ! -----------------
-! 
+! 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. 
 !
 ! Former revisions:
@@ -41 +46 @@
     INTEGER, DIMENSION(MPI_STATUS_SIZE,4) ::  wait_stat
 #endif
-
-    REAL ::  ar(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
-
+    INTEGER :: i,nbgp_local
+    REAL, DIMENSION(nzb:nzt+1,nys-nbgp_local:nyn+nbgp_local, &
+                    nxl-nbgp_local:nxr+nbgp_local) ::  ar
 
     CALL cpu_log( log_point_s(2), 'exchange_horiz', 'start' )
 
+    IF ( exchange_mg == .TRUE. ) THEN
+      i = grid_level
+    ELSE
+      i = 0
+    END IF
 #if defined( __parallel )
 
@@ -56 +66 @@
 !--    within the PE memory
        IF ( bc_lr == 'cyclic' )  THEN
-          ar(:,nys:nyn,nxl-1) = ar(:,nys:nyn,nxr)
-          ar(:,nys:nyn,nxr+1) = ar(:,nys:nyn,nxl)
+          ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)
+          ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1)
        ENDIF
 
@@ -65 +75 @@
 !
 !--    Send left boundary, receive right one
-       CALL MPI_ISEND(                                                     &
-               ar(nzb,nys-1,nxl), ngp_yz(grid_level), MPI_REAL, pleft,  0, &
-                          comm2d, req(1), ierr )
-       CALL MPI_IRECV(                                                       &
-               ar(nzb,nys-1,nxr+1), ngp_yz(grid_level), MPI_REAL, pright, 0, &
-                          comm2d, req(2), ierr )
+       CALL MPI_ISEND(ar(nzb,nys-nbgp_local,nxl),1,type_yz(i),pleft,0,comm2d,&
+                      req(1),ierr)
+       CALL MPI_IRECV(ar(nzb,nys-nbgp_local,nxr+1),1,type_yz(i),pright,0,&
+                     comm2d,req(2),ierr)
 !
 !--    Send right boundary, receive left one
-       CALL MPI_ISEND(                                                     &
-               ar(nzb,nys-1,nxr), ngp_yz(grid_level), MPI_REAL, pright, 1, &
-                          comm2d, req(3), ierr )
-       CALL MPI_IRECV(                                                       &
-               ar(nzb,nys-1,nxl-1), ngp_yz(grid_level), MPI_REAL, pleft,  1, &
-                          comm2d, req(4), ierr )
+
+
+       CALL MPI_ISEND(ar(nzb,nys-nbgp_local,nxr+1-nbgp_local),1,type_yz(i),pright, 1,  &
+                      comm2d, req(3), ierr )
+       CALL MPI_IRECV(ar(nzb,nys-nbgp_local,nxl-nbgp_local),1,type_yz(i),pleft,1,&
+                      comm2d,req(4), ierr)
+
        CALL MPI_WAITALL( 4, req, wait_stat, ierr )
 
@@ -89 +98 @@
 !--    within the PE memory
        IF ( bc_ns == 'cyclic' )  THEN
-          ar(:,nys-1,:) = ar(:,nyn,:)
-          ar(:,nyn+1,:) = ar(:,nys,:)
+          ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)
+          ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:)
        ENDIF
 
@@ -98 +107 @@
 !
 !--    Send front boundary, receive rear one
-       CALL MPI_ISEND( ar(nzb,nys,nxl-1),   1, type_xz(grid_level), psouth, 0, &
+!--    MPI_ISEND initial send adress changed, type_xz() is sendet nbgp times
+
+       CALL MPI_ISEND( ar(nzb,nys,nxl-nbgp_local),1, type_xz(i), psouth, 0, &
                        comm2d, req(1), ierr )
-       CALL MPI_IRECV( ar(nzb,nyn+1,nxl-1), 1, type_xz(grid_level), pnorth, 0, &
+       CALL MPI_IRECV( ar(nzb,nyn+1,nxl-nbgp_local),1, type_xz(i), pnorth, 0, &
                        comm2d, req(2), ierr )
 !
 !--    Send rear boundary, receive front one
-       CALL MPI_ISEND( ar(nzb,nyn,nxl-1),   1, type_xz(grid_level), pnorth, 1, &
+       CALL MPI_ISEND( ar(nzb,nyn-nbgp_local+1,nxl-nbgp_local),1, type_xz(i), pnorth, 1, &
                        comm2d, req(3), ierr )
-       CALL MPI_IRECV( ar(nzb,nys-1,nxl-1), 1, type_xz(grid_level), psouth, 1, &
+       CALL MPI_IRECV( ar(nzb,nys-nbgp_local,nxl-nbgp_local),1, type_xz(i), psouth, 1, &
                        comm2d, req(4), ierr )
        call MPI_WAITALL( 4, req, wait_stat, ierr )
 
     ENDIF
-
 
 #else
@@ -118 +128 @@
 !-- Lateral boundary conditions in the non-parallel case
     IF ( bc_lr == 'cyclic' )  THEN
-       ar(:,nys:nyn,nxl-1) = ar(:,nys:nyn,nxr)
-       ar(:,nys:nyn,nxr+1) = ar(:,nys:nyn,nxl)
+        ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)
+        ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1)
     ENDIF
 
     IF ( bc_ns == 'cyclic' )  THEN
-       ar(:,nys-1,:) = ar(:,nyn,:)
-       ar(:,nyn+1,:) = ar(:,nys,:)
+        ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)
+        ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:)
     ENDIF
 
 #endif
-
     CALL cpu_log( log_point_s(2), 'exchange_horiz', 'stop' )
 
+
  END SUBROUTINE exchange_horiz

palm/trunk/SOURCE/exchange_horiz_2d.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! 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.
 !
 ! Former revisions:
@@ -37 +38 @@
     IMPLICIT NONE
 
-    REAL ::  ar(nys-1:nyn+1,nxl-1:nxr+1)
+    REAL ::  ar(nysg:nyng,nxlg:nxrg)
+    INTEGER :: i
 
 
@@ -51 +53 @@
 !--    One-dimensional decomposition along y, boundary values can be exchanged
 !--    within the PE memory
-       ar(nys:nyn,nxl-1) = ar(nys:nyn,nxr)
-       ar(nys:nyn,nxr+1) = ar(nys:nyn,nxl)
+       ar(:,nxl-nbgp:nxl-1) = ar(:,nxr-nbgp+1:nxr)
+       ar(:,nxr+1:nxr+nbgp) = ar(:,nxl:nxl+nbgp-1)
 
     ELSE
 !
 !--    Send left boundary, receive right one
-       CALL MPI_SENDRECV( ar(nys,nxl),   ngp_y, MPI_REAL, pleft,  0, &
-                          ar(nys,nxr+1), ngp_y, MPI_REAL, pright, 0, &
+
+       CALL MPI_SENDRECV( ar(nysg,nxl), 1, type_y, pleft,  0, &
+                          ar(nysg,nxr+1), 1, type_y, pright, 0, &
                           comm2d, status, ierr )
 !
 !--    Send right boundary, receive left one
-       CALL MPI_SENDRECV( ar(nys,nxr),   ngp_y, MPI_REAL, pright,  1, &
-                          ar(nys,nxl-1), ngp_y, MPI_REAL, pleft,   1, &
+       CALL MPI_SENDRECV( ar(nysg,nxr+1-nbgp), 1, type_y, pright,  1, &
+                          ar(nysg,nxlg), 1, type_y, pleft,   1, &
                           comm2d, status, ierr )
     ENDIF
@@ -71 +74 @@
 !--    One-dimensional decomposition along x, boundary values can be exchanged
 !--    within the PE memory
-       ar(nys-1,:) = ar(nyn,:)
-       ar(nyn+1,:) = ar(nys,:)
+       ar(nys-nbgp:nys-1,:) = ar(nyn-nbgp+1:nyn,:)
+       ar(nyn+1:nyn+nbgp,:) = ar(nys:nys+nbgp-1,:)
 
     ELSE
 !
 !--    Send front boundary, receive rear one
-       CALL MPI_SENDRECV( ar(nys,nxl-1),   1, type_x, psouth, 0, &
-                          ar(nyn+1,nxl-1), 1, type_x, pnorth, 0, &
+
+       CALL MPI_SENDRECV( ar(nys,nxlg), 1, type_x, psouth, 0, &        !replace number of sended elements from 
+                          ar(nyn+1,nxlg), 1, type_x, pnorth, 0, &      ! nbgp to 1
                           comm2d, status, ierr )
 !
 !--    Send rear boundary, receive front one
-       CALL MPI_SENDRECV( ar(nyn,nxl-1),   1, type_x, pnorth, 1, &
-                          ar(nys-1,nxl-1), 1, type_x, psouth, 1, &
-                          comm2d, status, ierr )
+       CALL MPI_SENDRECV( ar(nyn+1-nbgp,nxlg), 1, type_x, pnorth, 1, &
+                          ar(nysg,nxlg), 1, type_x, psouth, 1, &
+                          comm2d, status, ierr )
+
     ENDIF
 
@@ -92 +97 @@
 !-- Lateral boundary conditions in the non-parallel case
     IF ( bc_lr == 'cyclic' )  THEN
-       ar(nys:nyn,nxl-1) = ar(nys:nyn,nxr)
-       ar(nys:nyn,nxr+1) = ar(nys:nyn,nxl)
+       ar(:,nxl-nbgp:nxl-1) = ar(:,nxr-nbgp+1:nxr)
+       ar(:,nxr+1:nxr+nbgp) = ar(:,nxl:nxl+nbgp-1)
     ENDIF
 
     IF ( bc_ns == 'cyclic' )  THEN
-       ar(nys-1,:) = ar(nyn,:)
-       ar(nyn+1,:) = ar(nys,:)
-    ENDIF
+       ar(nys-nbgp:nys-1,:) = ar(nyn-nbgp+1:nyn,:)
+       ar(nyn+1:nyn+nbgp,:) = ar(nys:nys+nbgp-1,:)
+    ENDIF
+
 
 #endif
@@ -106 +112 @@
 !-- Neumann-conditions at inflow/outflow in case of non-cyclic boundary
 !-- conditions
-    IF ( inflow_l .OR. outflow_l )  ar(:,nxl-1) = ar(:,nxl)
-    IF ( inflow_r .OR. outflow_r )  ar(:,nxr+1) = ar(:,nxr)
-    IF ( inflow_s .OR. outflow_s )  ar(nys-1,:) = ar(nys,:)
-    IF ( inflow_n .OR. outflow_n )  ar(nyn+1,:) = ar(nyn,:)
-
+    IF ( inflow_l .OR. outflow_l )  THEN
+       DO i=nbgp, 1, -1
+         ar(:,nxl-i) = ar(:,nxl)
+       END DO
+    END IF
+    IF ( inflow_r .OR. outflow_r )  THEN 
+       DO i=1, nbgp
+          ar(:,nxr+i) = ar(:,nxr)
+       END DO
+    END IF
+    IF ( inflow_s .OR. outflow_s )  THEN
+       DO i=nbgp, 1, -1
+         ar(nys-i,:) = ar(nys,:)
+       END DO
+    END IF
+    IF ( inflow_n .OR. outflow_n )  THEN
+       DO i=1, nbgp
+         ar(nyn+i,:) = ar(nyn,:)
+       END DO
+    END IF
     CALL cpu_log( log_point_s(13), 'exchange_horiz_2d', 'stop' )
 
@@ -134 +155 @@
     IMPLICIT NONE
 
-    INTEGER ::  ar(nys-1:nyn+1,nxl-1:nxr+1)
+    REAL ::  ar(nysg:nyng,nxlg:nxrg)
+    INTEGER :: i
+
 
 
@@ -154 +177 @@
 !
 !--    Send left boundary, receive right one
-       CALL MPI_SENDRECV( ar(nys,nxl),   ngp_y, MPI_INTEGER, pleft,  0, &
-                          ar(nys,nxr+1), ngp_y, MPI_INTEGER, pright, 0, &
+       CALL MPI_SENDRECV( ar(nysg,nxl), 1, type_y_int, pleft,  0, &
+                          ar(nysg,nxr+1), 1, type_y_int, pright, 0, &
                           comm2d, status, ierr )
 !
 !--    Send right boundary, receive left one
-       CALL MPI_SENDRECV( ar(nys,nxr),   ngp_y, MPI_INTEGER, pright,  1, &
-                          ar(nys,nxl-1), ngp_y, MPI_INTEGER, pleft,   1, &
-                          comm2d, status, ierr )
+       CALL MPI_SENDRECV( ar(nysg,nxr+1-nbgp), 1, type_y_int, pright,  1, &
+                          ar(nysg,nxlg), 1, type_y_int, pleft,   1, &
+                          comm2d, status, ierr )
+
     ENDIF
 
@@ -168 +192 @@
 !--    One-dimensional decomposition along x, boundary values can be exchanged
 !--    within the PE memory
-       ar(nys-1,:) = ar(nyn,:)
-       ar(nyn+1,:) = ar(nys,:)
+       ar(nysg:nys-1,:) = ar(nyn+1-nbgp:nyn,:)
+       ar(nyn+1:nyng,:) = ar(nys:nys-1+nbgp,:)
+
 
     ELSE
 !
 !--    Send front boundary, receive rear one
-       CALL MPI_SENDRECV( ar(nys,nxl-1),   1, type_x_int, psouth, 0, &
-                          ar(nyn+1,nxl-1), 1, type_x_int, pnorth, 0, &
-                          comm2d, status, ierr )
+       CALL MPI_SENDRECV( ar(nys,nxlg),1, type_x_int, psouth, 0, &
+                          ar(nyn+1,nxlg),1, type_x_int, pnorth, 0, &
+                          comm2d, status, ierr )
+
 !
 !--    Send rear boundary, receive front one
-       CALL MPI_SENDRECV( ar(nyn,nxl-1),   1, type_x_int, pnorth, 1, &
-                          ar(nys-1,nxl-1), 1, type_x_int, psouth, 1, &
-                          comm2d, status, ierr )
+       CALL MPI_SENDRECV( ar(nyn+1-nbgp,nxlg), nbgp, type_x_int, pnorth, 1, &
+                          ar(nysg,nxlg), nbgp, type_x_int, psouth, 1, &
+                          comm2d, status, ierr )
+
     ENDIF
 
@@ -194 +221 @@
 
     IF ( bc_ns == 'cyclic' )  THEN
-       ar(nys-1,:) = ar(nyn,:)
-       ar(nyn+1,:) = ar(nys,:)
+       ar(nysg:nys-1,:) = ar(nyn+1-nbgp:nyn,:)
+       ar(nyn+1:nyng,:) = ar(nys:nys-1+nbgp,:)
     ENDIF
 

palm/trunk/SOURCE/flow_statistics.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! When advection is computed with ws-scheme, turbulent fluxes are already
+! computed in the respective advection routines and buffered in arrays
+! sums_xx_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.
+!
 !
 ! Former revisions:
@@ -102 +111 @@
     REAL    ::  sums_ll(nzb:nzt+1,2)
 
-
     CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
 
@@ -133 +141 @@
        sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
 
+!
+!--    Copy the turbulent quantities, evaluated in the advection routines to 
+!--    the local array sums_l() for further computations
+       IF ( ws_scheme_mom )  THEN
+!       
+!--       Boundary condition for u'u' and v'v', because below the surface no
+!--       computation for these quantities is done.
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                sums_us2_ws_l(nzb_u_inner(j,i),sr) =                          & 
+                    sums_us2_ws_l(nzb_u_inner(j,i)+1,sr)
+                sums_vs2_ws_l(nzb_v_inner(j,i),sr) =                          &  
+                    sums_vs2_ws_l(nzb_v_inner(j,i)+1,sr)
+             ENDDO
+          ENDDO
+!          
+!--       Swap the turbulent quantities evaluated in advec_ws.
+          sums_l(:,13,0) = sums_wsus_ws_l(:,sr)       ! w*u* 
+          sums_l(:,15,0) = sums_wsvs_ws_l(:,sr)       ! w*v* 
+          sums_l(:,30,0) = sums_us2_ws_l(:,sr)        ! u*2 
+          sums_l(:,31,0) = sums_vs2_ws_l(:,sr)        ! v*2 
+          sums_l(:,32,0) = sums_ws2_ws_l(:,sr)        ! w*2 
+          sums_l(:,34,0) = sums_l(:,34,0) + 0.5 *                             & 
+                (sums_us2_ws_l(:,sr) + sums_vs2_ws_l(:,sr)                    &
+                + sums_ws2_ws_l(:,sr))                      ! e*
+          DO  k = nzb, nzt
+             sums_l(nzb+5,pr_palm,0) = sums_l(nzb+5,pr_palm,0) + 0.5 * (      &
+                sums_us2_ws_l(k,sr) + sums_vs2_ws_l(k,sr) +                   &
+                sums_ws2_ws_l(k,sr))
+          ENDDO
+       ENDIF
+       IF ( ws_scheme_sca )  THEN
+          sums_l(:,17,0) = sums_wspts_ws_l(:,sr)      ! w*pts* from advec_s_ws
+          IF ( ocean ) sums_l(:,66,0) = sums_wssas_ws_l(:,sr) ! w*sa*
+          IF ( humidity  .OR.  passive_scalar ) sums_l(:,49,0) =              &
+                                                   sums_wsqs_ws_l(:,sr) !w*q* 
+       ENDIF
 ! 
 !--    Horizontally averaged profiles of horizontal velocities and temperature.
@@ -138 +183 @@
 !--    for other horizontal averages.
        tn = 0
+
        !$OMP PARALLEL PRIVATE( i, j, k, tn )
 #if defined( __intel_openmp_bug )
@@ -215 +261 @@
        ENDIF
        !$OMP END PARALLEL
-
 !
 !--    Summation of thread sums
@@ -304 +349 @@
        hom(:,1,2,sr) = sums(:,2)             ! v
        hom(:,1,4,sr) = sums(:,4)             ! pt
+
 
 !
@@ -354 +400 @@
           DO  j =  nys, nyn
              sums_l_etot = 0.0
-             sums_l_eper = 0.0
              DO  k = nzb_s_inner(j,i), nzt+1
-                u2   = u(k,j,i)**2
-                v2   = v(k,j,i)**2
-                w2   = w(k,j,i)**2
-                ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
-                vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
 !
 !--             Prognostic and diagnostic variables
@@ -369 +409 @@
                 sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
 
-!
-!--             Variances
-                sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
-                sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
-                sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
                 sums_l(k,33,tn) = sums_l(k,33,tn) + &
                                   ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
@@ -381 +416 @@
                                   ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
                 ENDIF
-!
-!--             Higher moments
-!--             (Computation of the skewness of w further below)
-                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
-                                                    rmask(j,i,sr)
-!
-!--             Perturbation energy
-                sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 * ( ust2 + vst2 + w2 ) &
-                                                    * rmask(j,i,sr)
+
                 sums_l_etot  = sums_l_etot + &
-                                        0.5 * ( u2 + v2 + w2 ) * rmask(j,i,sr)
-                sums_l_eper  = sums_l_eper + &
-                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
+                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
+                                        w(k,j,i)**2 ) * rmask(j,i,sr)
              ENDDO
 !
@@ -401 +427 @@
 !--          allow vectorization of that loop.
              sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
-             sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) + sums_l_eper
 !
 !--          2D-arrays (being collected in the last column of sums_l)
@@ -420 +445 @@
 
 !
+!--    Computation of statistics when ws-scheme is not used. Else these 
+!--    quantities are evaluated in the advection routines.
+       IF ( .NOT. ws_scheme_mom )  THEN
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                sums_l_eper = 0.0
+                DO  k = nzb_s_inner(j,i), nzt+1
+                   u2   = u(k,j,i)**2
+                   v2   = v(k,j,i)**2
+                   w2   = w(k,j,i)**2
+                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
+                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
+
+                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
+                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
+                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
+!
+!--   Higher moments
+!--  (Computation of the skewness of w further below)
+                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
+                                                    rmask(j,i,sr)
+!
+!--             Perturbation energy
+
+                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
+                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
+                   sums_l_eper  = sums_l_eper + &
+                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
+
+                ENDDO
+                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
+                     + sums_l_eper
+             ENDDO
+          ENDDO
+       ELSE
+          !$OMP DO
+          DO  i = nxl, nxr
+             DO  j =  nys, nyn
+                DO  k = nzb_s_inner(j,i), nzt + 1
+                   w2   = w(k,j,i)**2
+!
+!--                Higher moments
+!--                (Computation of the skewness of w further below)
+                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
+                                                    rmask(j,i,sr)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+
+!
 !--    Horizontally averaged profiles of the vertical fluxes
+
        !$OMP DO
        DO  i = nxl, nxr
@@ -587 +665 @@
                 pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
                               pt(k+1,j,i) - hom(k+1,1,4,sr) )
-!
-!--             Momentum flux w*u*
-                sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                     &
-                                                    ( w(k,j,i-1) + w(k,j,i) ) &
-                                                    * ust * rmask(j,i,sr)
-!
-!--             Momentum flux w*v*
-                sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                     &
-                                                    ( w(k,j-1,i) + w(k,j,i) ) &
-                                                    * vst * rmask(j,i,sr)
-!
-!--             Heat flux w*pt*
-!--             The following formula (comment line, not executed) does not
-!--             work if applied to subregions
-!                sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 *                     &
-!                                                    ( pt(k,j,i)+pt(k+1,j,i) ) &
-!                                                    * w(k,j,i) * rmask(j,i,sr)
-                sums_l(k,17,tn) = sums_l(k,17,tn) + pts * w(k,j,i) * &
-                                                    rmask(j,i,sr)
-!
+
 !--             Higher moments
                 sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
@@ -617 +676 @@
 !--             but so far there is no other suitable place to calculate)
                 IF ( ocean )  THEN
-                   pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
+                   IF( .NOT. ws_scheme_sca )  THEN
+                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
                                  sa(k+1,j,i) - hom(k+1,1,23,sr) )
-                   sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
+                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
                                                        rmask(j,i,sr)
+                   ENDIF
                    sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
                                                        rmask(j,i,sr)
@@ -635 +696 @@
                    sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
                                                        rmask(j,i,sr)
-                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
-                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
-                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
-                                                       rmask(j,i,sr)
+
                    IF ( cloud_physics  .OR.  cloud_droplets )  THEN
                       pts = 0.5 *                                            &
@@ -652 +710 @@
 !
 !--             Passive scalar flux
-                IF ( passive_scalar )  THEN
+                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca ))  THEN
                    pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
                                  q(k+1,j,i) - hom(k+1,1,41,sr) )
@@ -661 +719 @@
 !
 !--             Energy flux w*e*
-                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *           &
-                                              ( ust**2 + vst**2 + w(k,j,i)**2 )&
-                                              * rmask(j,i,sr)
-          
+!--             has to be adjusted
+                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
+                                             ( ust**2 + vst**2 + w(k,j,i)**2 )&
+                                             * rmask(j,i,sr)
              ENDDO
           ENDDO
        ENDDO
+!-     for reasons of speed optimization the loop is splitted, to avoid if-else
+!-     statements inside the loop
+!-     Fluxes which have been computed in part directly inside the advection routines
+!-     treated seperatly.
+!-     First treat the momentum fluxes
+       IF ( .NOT. ws_scheme_mom )  THEN
+         !$OMP DO
+         DO  i = nxl, nxr
+            DO  j = nys, nyn
+               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
+                  ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
+                              u(k+1,j,i) - hom(k+1,1,1,sr) )
+                  vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
+                              v(k+1,j,i) - hom(k+1,1,2,sr) )
+!                              
+!--               Momentum flux w*u*
+                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
+                                                    ( w(k,j,i-1) + w(k,j,i) ) &
+                                                    * ust * rmask(j,i,sr)
+!
+!--               Momentum flux w*v*
+                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
+                                                    ( w(k,j-1,i) + w(k,j,i) ) &
+                                                    * vst * rmask(j,i,sr)
+               ENDDO
+            ENDDO
+         ENDDO
+
+       ENDIF
+       IF ( .NOT. ws_scheme_sca )  THEN
+         !$OMP DO
+         DO  i = nxl, nxr
+            DO  j = nys, nyn
+               DO  k = nzb_diff_s_inner(j,i) - 1, nzt_diff
+!-                vertical heat flux
+                  sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
+                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
+                           pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
+                           * w(k,j,i) * rmask(j,i,sr)
+                  IF ( humidity )  THEN
+                     pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
+                                q(k+1,j,i) - hom(k+1,1,41,sr) )
+                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
+                                                      rmask(j,i,sr)
+                  ENDIF
+               ENDDO
+            ENDDO
+         ENDDO
+
+       ENDIF
+
 
 !
@@ -814 +923 @@
 
 #if defined( __parallel )
+
 !
 !--    Compute total sum from local sums
@@ -843 +953 @@
           sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
        ENDDO
+
 !--    Upstream-parts
        sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
@@ -868 +979 @@
           ENDDO
        ENDIF
-
 !
 !--    Collect horizontal average in hom.
@@ -934 +1044 @@
                                        ! upstream-parts u_x, u_y, u_z, v_x,
                                        ! v_y, usw. (in last but one profile)
-       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)  
+       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
                                        ! u*, w'u', w'v', t* (in last profile)
 
@@ -951 +1061 @@
        z_i(1) = 0.0
        first = .TRUE.
+
        IF ( ocean )  THEN
           DO  k = nzt, nzb+1, -1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
+             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
              IF ( hom(k,1,18,sr) < 0.0  .AND. &
+                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
                   hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
                 IF ( zw(k) < 1.5 * height )  THEN
@@ -969 +1082 @@
        ELSE
           DO  k = nzb, nzt-1
-             IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
+             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
+                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
                 first = .FALSE.
                 height = zw(k)
              ENDIF
-             IF ( hom(k,1,18,sr) < 0.0  .AND. &
+             IF ( hom(k,1,18,sr) < 0.0  .AND. & 
+                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
                   hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
                 IF ( zw(k) < 1.5 * height )  THEN
@@ -1026 +1141 @@
 !--    the characteristic convective boundary layer temperature.
 !--    The horizontal average at nzb+1 is input for the average temperature.
-       IF ( hom(nzb,1,18,sr) > 0.0  .AND.  z_i(1) /= 0.0 )  THEN
+       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
+           .AND.  z_i(1) /= 0.0 )  THEN
           hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
                                        hom(nzb,1,18,sr) *      &

palm/trunk/SOURCE/global_min_max.f90

@@ -5 +5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adapting of the index arrays, because MINLOC assumes lowerbound at 1 and not 
+! at nbgp.
 !
 ! Former revisions:
@@ -59 +60 @@
 !--    Determine the local minimum
        fmin_ijk_l = MINLOC( ar )
-       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - 1    ! MINLOC assumes lowerbound = 1
-       fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - 1
+       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - nbgp    ! MINLOC assumes lowerbound = 1
+       fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
        fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - 1
        fmin_l(1)  = ar(fmin_ijk_l(1),fmin_ijk_l(2),fmin_ijk_l(3))
@@ -100 +101 @@
 !--    Determine the local maximum
        fmax_ijk_l = MAXLOC( ar )
-       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - 1    ! MAXLOC assumes lowerbound = 1
-       fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - 1
+       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - nbgp    ! MAXLOC assumes lowerbound = 1
+       fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
        fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - 1
        fmax_l(1) = ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3))
@@ -221 +222 @@
           IF ( fmax_ijk(1) < 0 )  THEN
              value        = -value
-             value_ijk(1) = -value_ijk(1) - 10
+             value_ijk(1) = -value_ijk(1) - 10         !???
           ENDIF
 
@@ -228 +229 @@
 !
 !-- Limit index values to the range 0..nx, 0..ny
-    IF ( value_ijk(3) ==   -1 )  value_ijk(3) = nx
-    IF ( value_ijk(3) == nx+1 )  value_ijk(3) =  0
-    IF ( value_ijk(2) ==   -1 )  value_ijk(2) = ny
-    IF ( value_ijk(2) == ny+1 )  value_ijk(2) =  0
+    IF ( value_ijk(3) < 0  ) value_ijk(3) = nx +1 + value_ijk(3)
+    IF ( value_ijk(3) > nx ) value_ijk(3) = value_ijk(3) - (nx+1)
+    IF ( value_ijk(2) < 0  ) value_ijk(2) = ny +1 + value_ijk(2)
+    IF ( value_ijk(2) > ny ) value_ijk(2) = value_ijk(2) - (ny+1)
 
 

palm/trunk/SOURCE/header.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output of advection scheme.
+! Modified output of Prandtl-layer height.
 ! 
 ! Former revisions:
@@ -250 +251 @@
     IF ( momentum_advec == 'pw-scheme' )  THEN
        WRITE ( io, 113 )
-    ELSE
+    ELSEIF (momentum_advec == 'ws-scheme' ) THEN
+       WRITE ( io, 503 )
+    ELSEIF (momentum_advec == 'ups-scheme' ) THEN
        WRITE ( io, 114 )
        IF ( cut_spline_overshoot )  WRITE ( io, 124 )
@@ -267 +270 @@
     IF ( scalar_advec == 'pw-scheme' )  THEN
        WRITE ( io, 116 )
+    ELSEIF ( scalar_advec == 'ws-scheme' )  THEN
+       WRITE ( io, 504 )
     ELSEIF ( scalar_advec == 'ups-scheme' )  THEN
        WRITE ( io, 117 )
@@ -575 +580 @@
     ELSEIF( ibc_pt_t == 2 )  THEN
        roben  = TRIM( roben  ) // ' pt(nzt+1) = pt(nzt) + dpt/dz_ini'
+
     ENDIF
 
@@ -662 +668 @@
 
     IF ( prandtl_layer )  THEN
-       WRITE ( io, 305 )  0.5 * (zu(1)-zu(0)), roughness_length, kappa, &
+       WRITE ( io, 305 )  (zu(1)-zu(0)), roughness_length, kappa, &
                           rif_min, rif_max
        IF ( .NOT. constant_heatflux )  WRITE ( io, 308 )
@@ -1981 +1987 @@
             '    Dissipation calculation:           ',A/)
 502 FORMAT ('    Damping layer starts from ',F7.1,' m (GP ',I4,')'/)
+503 FORMAT (' --> Momentum advection via Wicker-Skamarock-Scheme 5th order')
+504 FORMAT (' --> Scalar advection via Wicker-Skamarock-Scheme 5th order')
 
 

palm/trunk/SOURCE/inflow_turbulence.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! Using nbgp recycling planes for a better resolution of the turbulent flow 
+! near the inflow.
 !
 ! Former revisions:
@@ -35 +36 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, imax, j, k, ngp_ifd, ngp_pr
+    INTEGER ::  i, imax, j, k, l, ngp_ifd, ngp_pr
 
     REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
-    REAL, DIMENSION(nzb:nzt+1,5) ::  avpr, avpr_l
-    REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,5) ::  inflow_dist
+    REAL, DIMENSION(nzb:nzt+1,5,nbgp) ::  avpr, avpr_l
+    REAL, DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) ::  inflow_dist
 
     CALL cpu_log( log_point(40), 'inflow_turbulence', 'start' )
 
 !
-!-- Carry out horizontal averaging in the recycling plane
+!-- Carry out spanwise averaging in the recycling plane
     avpr_l = 0.0
-    ngp_pr = ( nzt - nzb + 2 ) * 5
-    ngp_ifd = ngp_pr * ( nyn - nys + 3 )
+    ngp_pr = ( nzt - nzb + 2 ) * 5 * nbgp
+    ngp_ifd = ngp_pr * ( nyn - nys + 1 + 2 * nbgp )
 
 !
 !-- First, local averaging within the recycling domain
-    IF ( recycling_plane >= nxl )  THEN
-
-       imax = MIN( nxr, recycling_plane )
-
-       DO  i = nxl, imax
+
+    i = recycling_plane
+
+#if defined( __parallel )
+    IF ( myidx == id_recycling )  THEN
+       
+       DO  l = 1, nbgp
           DO  j = nys, nyn
-             DO  k = nzb, nzt+1
-
-                avpr_l(k,1) = avpr_l(k,1) + u(k,j,i)
-                avpr_l(k,2) = avpr_l(k,2) + v(k,j,i)
-                avpr_l(k,3) = avpr_l(k,3) + w(k,j,i)
-                avpr_l(k,4) = avpr_l(k,4) + pt(k,j,i)
-                avpr_l(k,5) = avpr_l(k,5) + e(k,j,i)
+             DO  k = nzb, nzt + 1
+
+                avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
+                avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
+                avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
+                avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
+                avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
 
              ENDDO
           ENDDO
-       ENDDO
-
-    ENDIF
-
-!    WRITE (9,*) '*** averaged profiles avpr_l'
-!    DO  k = nzb, nzt+1
-!       WRITE (9,'(F5.1,1X,F5.1,1X,F5.1,1X,F6.1,1X,F7.2)') avpr_l(k,1),avpr_l(k,2),avpr_l(k,3),avpr_l(k,4),avpr_l(k,5)
-!    ENDDO
-!    WRITE (9,*) ' '
-
-#if defined( __parallel )
+          i = i + 1
+       ENDDO
+
+    ENDIF
 !
 !-- Now, averaging over all PEs
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( avpr_l(nzb,1), avpr(nzb,1), ngp_pr, MPI_REAL, MPI_SUM, &
-                        comm2d, ierr )
+    CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr,  & 
+                    MPI_REAL, MPI_SUM, comm2d, ierr )
+
 #else
+    DO  l = 1, nbgp
+       DO  j = nys, nyn
+          DO  k = nzb, nzt + 1
+
+             avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
+             avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
+             avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
+             avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
+             avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
+
+          ENDDO
+       ENDDO
+       i = i + 1 
+    ENDDO
+    
     avpr = avpr_l
 #endif
 
-    avpr = avpr / ( ( ny + 1 ) * ( recycling_plane + 1 ) )
-
-!    WRITE (9,*) '*** averaged profiles'
-!    DO  k = nzb, nzt+1
-!       WRITE (9,'(F5.1,1X,F5.1,1X,F5.1,1X,F6.1,1X,F7.2)') avpr(k,1),avpr(k,2),avpr(k,3),avpr(k,4),avpr(k,5)
-!    ENDDO
-!    WRITE (9,*) ' '
-
+    avpr = avpr / ( ny + 1 )
 !
 !-- Calculate the disturbances at the recycling plane
@@ -101 +106 @@
 #if defined( __parallel )
     IF ( myidx == id_recycling )  THEN
-
-       DO  j = nys-1, nyn+1
+       DO  l = 1, nbgp
+          DO  j = nysg, nyng
+             DO  k = nzb, nzt + 1
+
+                inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
+                inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
+                inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
+                inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
+                inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
+              
+            ENDDO
+          ENDDO
+          i = i + 1
+       ENDDO
+
+    ENDIF
+#else
+    DO  l = 1, nbgp
+       DO  j = nysg, nyng
           DO  k = nzb, nzt+1
 
-              inflow_dist(k,j,1) = u(k,j,i+1) - avpr(k,1)
-              inflow_dist(k,j,2) = v(k,j,i)   - avpr(k,2)
-              inflow_dist(k,j,3) = w(k,j,i)   - avpr(k,3)
-              inflow_dist(k,j,4) = pt(k,j,i)  - avpr(k,4)
-              inflow_dist(k,j,5) = e(k,j,i)   - avpr(k,5)
-
-          ENDDO
-       ENDDO
-
-    ENDIF
-#else
-    DO  j = nys-1, nyn+1
-       DO  k = nzb, nzt+1
-
-          inflow_dist(k,j,1) = u(k,j,i+1) - avpr(k,1)
-          inflow_dist(k,j,2) = v(k,j,i)   - avpr(k,2)
-          inflow_dist(k,j,3) = w(k,j,i)   - avpr(k,3)
-          inflow_dist(k,j,4) = pt(k,j,i)  - avpr(k,4)
-          inflow_dist(k,j,5) = e(k,j,i)   - avpr(k,5)
-
-       ENDDO
+             inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
+             inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
+             inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
+             inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
+             inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
+              
+          ENDDO
+       ENDDO
+       i = i + 1
     ENDDO
 #endif
@@ -134 +144 @@
     IF ( myidx == id_recycling  .AND.  myidx /= id_inflow )  THEN
 
-       CALL MPI_SEND( inflow_dist(nzb,nys-1,1), ngp_ifd, MPI_REAL, &
+       CALL MPI_SEND( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
                       id_inflow, 1, comm1dx, ierr )
 
@@ -140 +150 @@
 
        inflow_dist = 0.0
-       CALL MPI_RECV( inflow_dist(nzb,nys-1,1), ngp_ifd, MPI_REAL, &
+       CALL MPI_RECV( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
                       id_recycling, 1, comm1dx, status, ierr )
 
@@ -150 +160 @@
     IF ( nxl == 0 )  THEN
 
-       DO  j = nys-1, nyn+1
-          DO  k = nzb, nzt+1
-
-!              WRITE (9,*) 'j=',j,' k=',k
-!              WRITE (9,*) 'mean_u = ', mean_inflow_profiles(k,1), ' dist_u = ',&
-!                          inflow_dist(k,j,1)
-!              WRITE (9,*) 'mean_v = ', mean_inflow_profiles(k,2), ' dist_v = ',&
-!                          inflow_dist(k,j,2)
-!              WRITE (9,*) 'mean_w = 0.0', ' dist_w = ',&
-!                          inflow_dist(k,j,3)
-!              WRITE (9,*) 'mean_pt = ', mean_inflow_profiles(k,4), ' dist_pt = ',&
-!                          inflow_dist(k,j,4)
-!              WRITE (9,*) 'mean_e = ', mean_inflow_profiles(k,5), ' dist_e = ',&
-!                          inflow_dist(k,j,5)
-              u(k,j,0)   = mean_inflow_profiles(k,1) + &
-                           inflow_dist(k,j,1) * inflow_damping_factor(k)
-              v(k,j,-1)  = mean_inflow_profiles(k,2) + &
-                           inflow_dist(k,j,2) * inflow_damping_factor(k)
-              w(k,j,-1)  = inflow_dist(k,j,3) * inflow_damping_factor(k)
-              pt(k,j,-1) = mean_inflow_profiles(k,4) + &
-                           inflow_dist(k,j,4) * inflow_damping_factor(k)
-              e(k,j,-1)  = mean_inflow_profiles(k,5) + &
-                           inflow_dist(k,j,5) * inflow_damping_factor(k)
-              e(k,j,-1)  = MAX( e(k,j,-1), 0.0 )
+       DO  j = nysg, nyng
+          DO  k = nzb, nzt + 1
+
+              u(k,j,-nbgp+1:0)   = mean_inflow_profiles(k,1) + &
+                           inflow_dist(k,j,1,1:nbgp) * inflow_damping_factor(k)
+              v(k,j,-nbgp:-1)  = mean_inflow_profiles(k,2) + &
+                           inflow_dist(k,j,2,1:nbgp) * inflow_damping_factor(k)
+              w(k,j,-nbgp:-1)  = inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)
+              pt(k,j,-nbgp:-1) = mean_inflow_profiles(k,4) + &
+                           inflow_dist(k,j,4,1:nbgp) * inflow_damping_factor(k)
+              e(k,j,-nbgp:-1)  = mean_inflow_profiles(k,5) + &
+                           inflow_dist(k,j,5,1:nbgp) * inflow_damping_factor(k)
+              e(k,j,-nbgp:-1)  = MAX( e(k,j,-nbgp:-1), 0.0 )
 
           ENDDO

palm/trunk/SOURCE/init_1d_model.f90

@@ -63 +63 @@
               l1d_m(nzb:nzt+1),  rif1d(nzb:nzt+1),   te_e(nzb:nzt+1),  &
               te_em(nzb:nzt+1),  te_u(nzb:nzt+1),    te_um(nzb:nzt+1), &
-              te_v(nzb:nzt+1),   te_vm(nzb:nzt+1),   u1d(nzb:nzt+1),   &
+              te_v(nzb:nzt+1),   te_vm(nzb:nzt+1),    u1d(nzb:nzt+1),   &
               u1d_m(nzb:nzt+1),  u1d_p(nzb:nzt+1),   v1d(nzb:nzt+1),   &
               v1d_m(nzb:nzt+1),  v1d_p(nzb:nzt+1) )
@@ -385 +385 @@
 !--       boundary condition applies to u and v.
 !--       The boundary condition for e is set further below ( (u*/cm)**2 ).
-          u1d_p(nzb) = -u1d_p(nzb+1)
-          v1d_p(nzb) = -v1d_p(nzb+1)
+         ! u1d_p(nzb) = -u1d_p(nzb+1)
+         ! v1d_p(nzb) = -v1d_p(nzb+1)
+
+          u1d_p(nzb) = 0.0
+          v1d_p(nzb) = 0.0
 
 !

palm/trunk/SOURCE/init_3d_model.f90

@@ -7 +7 @@
 ! Current revisions:
 ! -----------------
-! Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)
+! 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) 
 !
 ! Former revisions:
@@ -101 +115 @@
 !------------------------------------------------------------------------------!
 
+    USE advec_ws
     USE arrays_3d
     USE averaging
@@ -147 +162 @@
               ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions),               &
               ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions),             &
-              rmask(nys-1:nyn+1,nxl-1:nxr+1,0:statistic_regions),           &
+              rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions),           &
               sums(nzb:nzt+1,pr_palm+max_pr_user),                          &
               sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),   &
@@ -154 +169 @@
               sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions),                &
               ts_value(dots_max,0:statistic_regions) )
-    ALLOCATE( km_damp_x(nxl-1:nxr+1), km_damp_y(nys-1:nyn+1) )
-
-    ALLOCATE( rif_1(nys-1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &
-              ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1),  &
-              us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1),   &
-              uswst_1(nys-1:nyn+1,nxl-1:nxr+1),                               &
-              vsws_1(nys-1:nyn+1,nxl-1:nxr+1),                                &
-              vswst_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
+    ALLOCATE( km_damp_x(nxlg:nxrg), km_damp_y(nysg:nyng) )
+
+    ALLOCATE( rif_1(nysg:nyng,nxlg:nxrg), shf_1(nysg:nyng,nxlg:nxrg), &
+              ts(nysg:nyng,nxlg:nxrg), tswst_1(nysg:nyng,nxlg:nxrg),  &
+              us(nysg:nyng,nxlg:nxrg), usws_1(nysg:nyng,nxlg:nxrg),   &
+              uswst_1(nysg:nyng,nxlg:nxrg),                               &
+              vsws_1(nysg:nyng,nxlg:nxrg),                                &
+              vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg) )
 
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
 !
 !--    Leapfrog scheme needs two timelevels of diffusion quantities
-       ALLOCATE( rif_2(nys-1:nyn+1,nxl-1:nxr+1),   &
-                 shf_2(nys-1:nyn+1,nxl-1:nxr+1),   &
-                 tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
-                 usws_2(nys-1:nyn+1,nxl-1:nxr+1),  &
-                 uswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
-                 vswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
-                 vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
+       ALLOCATE( rif_2(nysg:nyng,nxlg:nxrg),   &
+                 shf_2(nysg:nyng,nxlg:nxrg),   &
+                 tswst_2(nysg:nyng,nxlg:nxrg), &
+                 usws_2(nysg:nyng,nxlg:nxrg),  &
+                 uswst_2(nysg:nyng,nxlg:nxrg), &
+                 vswst_2(nysg:nyng,nxlg:nxrg), &
+                 vsws_2(nysg:nyng,nxlg:nxrg) )
     ENDIF
 
     ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra),         &
-              e_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-              e_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-              e_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-              kh_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-              km_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-              p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),    &
-              pt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-              pt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-              pt_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-              tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-              u_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-              u_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-              u_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-              v_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-              v_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-              v_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-              w_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-              w_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-              w_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+              e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+              e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+              e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+              kh_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+              km_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+              p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
+              pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+              pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+              pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+              tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+              u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+              u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+              u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+              v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+              v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+              v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+              w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-       ALLOCATE( kh_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                 km_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+       ALLOCATE( kh_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+                 km_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ENDIF
 
@@ -204 +219 @@
 !
 !--    2D-humidity/scalar arrays
-       ALLOCATE ( qs(nys-1:nyn+1,nxl-1:nxr+1),     &
-                  qsws_1(nys-1:nyn+1,nxl-1:nxr+1), &
-                  qswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
+       ALLOCATE ( qs(nysg:nyng,nxlg:nxrg),     &
+                  qsws_1(nysg:nyng,nxlg:nxrg), &
+                  qswst_1(nysg:nyng,nxlg:nxrg) )
 
        IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-          ALLOCATE( qsws_2(nys-1:nyn+1,nxl-1:nxr+1), &
-                    qswst_2(nys-1:nyn+1,nxl-1:nxr+1) )
+          ALLOCATE( qsws_2(nysg:nyng,nxlg:nxrg), &
+                    qswst_2(nysg:nyng,nxlg:nxrg) )
        ENDIF
 !
 !--    3D-humidity/scalar arrays
-       ALLOCATE( q_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                 q_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                 q_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+       ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+                 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+                 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
 !
 !--    3D-arrays needed for humidity only
        IF ( humidity )  THEN
-          ALLOCATE( vpt_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+          ALLOCATE( vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
           IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-             ALLOCATE( vpt_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+             ALLOCATE( vpt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
           ENDIF
 
@@ -230 +245 @@
 !
 !--          Liquid water content
-             ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 !
 !--          Precipitation amount and rate (only needed if output is switched)
-             ALLOCATE( precipitation_amount(nys-1:nyn+1,nxl-1:nxr+1), &
-                       precipitation_rate(nys-1:nyn+1,nxl-1:nxr+1) )
+             ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &
+                       precipitation_rate(nysg:nyng,nxlg:nxrg) )
           ENDIF
 
@@ -241 +256 @@
 !--          Liquid water content, change in liquid water content,
 !--          real volume of particles (with weighting), volume of particles
-             ALLOCATE ( ql_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                        ql_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                        ql_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                        ql_vp(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+                        ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+                        ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+                        ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
           ENDIF
 
@@ -252 +267 @@
 
     IF ( ocean )  THEN
-       ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
-                 saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
-       ALLOCATE( prho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                 rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-                 sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
-                 sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
-                 sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+       ALLOCATE( saswsb_1(nysg:nyng,nxlg:nxrg), &
+                 saswst_1(nysg:nyng,nxlg:nxrg) )
+       ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
+                 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+                 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
+                 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
+                 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
        prho => prho_1
        rho  => rho_1  ! routines calc_mean_profile and diffusion_e require
                       ! density to be apointer
        IF ( humidity_remote )  THEN
-          ALLOCATE( qswst_remote(nys-1:nyn+1,nxl-1:nxr+1) )
+          ALLOCATE( qswst_remote(nysg:nyng,nxlg:nxrg))
           qswst_remote = 0.0
        ENDIF
@@ -272 +287 @@
 !-- particle velocities
     IF ( use_sgs_for_particles )  THEN
-       ALLOCATE ( diss(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+       ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     ELSE
        ALLOCATE ( diss(2,2,2) )  ! required because diss is used as a
@@ -288 +303 @@
 !-- 3D-arrays for the leaf area density and the canopy drag coefficient
     IF ( plant_canopy ) THEN
-       ALLOCATE ( lad_s(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-                  lad_u(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-                  lad_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-                  lad_w(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-                  cdc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+       ALLOCATE ( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+                  lad_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+                  lad_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+                  lad_w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
+                  cdc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
        IF ( passive_scalar ) THEN
-          ALLOCATE ( sls(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
-                     sec(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
+          ALLOCATE ( sls(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
+                     sec(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 
        ENDIF
 
        IF ( cthf /= 0.0 ) THEN
-          ALLOCATE ( lai(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
-                     canopy_heat_flux(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+          ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg),   &
+                     canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
        ENDIF
 
@@ -309 +324 @@
 !-- 4D-array for storing the Rif-values at vertical walls
     IF ( topography /= 'flat' )  THEN
-       ALLOCATE( rif_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1,1:4) )
+       ALLOCATE( rif_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg,1:4) )
        rif_wall = 0.0
-    ENDIF
-
-!
-!-- Velocities at nzb+1 needed for volume flow control
-    IF ( conserve_volume_flow )  THEN
-       ALLOCATE( u_nzb_p1_for_vfc(nys:nyn), v_nzb_p1_for_vfc(nxl:nxr) )
-       u_nzb_p1_for_vfc = 0.0
-       v_nzb_p1_for_vfc = 0.0
     ENDIF
 
@@ -325 +332 @@
 !-- are needed for radiation boundary conditions
     IF ( outflow_l )  THEN
-       ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,1:2), &
-                 v_m_l(nzb:nzt+1,nys-1:nyn+1,0:1), &
-                 w_m_l(nzb:nzt+1,nys-1:nyn+1,0:1) )
+       ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &
+                 v_m_l(nzb:nzt+1,nysg:nyng,0:1), &
+                 w_m_l(nzb:nzt+1,nysg:nyng,0:1) )
     ENDIF
     IF ( outflow_r )  THEN
-       ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
-                 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
-                 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx) )
+       ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
+                 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &
+                 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) )
     ENDIF
     IF ( outflow_l  .OR.  outflow_r )  THEN
-       ALLOCATE( c_u(nzb:nzt+1,nys-1:nyn+1), c_v(nzb:nzt+1,nys-1:nyn+1), &
-                 c_w(nzb:nzt+1,nys-1:nyn+1) )
+       ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &
+                 c_w(nzb:nzt+1,nysg:nyng) )
     ENDIF
     IF ( outflow_s )  THEN
-       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1), &
-                 v_m_s(nzb:nzt+1,1:2,nxl-1:nxr+1), &
-                 w_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1) )
+       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &
+                 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &
+                 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) )
     ENDIF
     IF ( outflow_n )  THEN
-       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
-                 v_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
-                 w_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1) )
+       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
+                 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &
+                 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) )
     ENDIF
     IF ( outflow_s  .OR.  outflow_n )  THEN
-       ALLOCATE( c_u(nzb:nzt+1,nxl-1:nxr+1), c_v(nzb:nzt+1,nxl-1:nxr+1), &
-                 c_w(nzb:nzt+1,nxl-1:nxr+1) )
+       ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &
+                 c_w(nzb:nzt+1,nxlg:nxrg) )
     ENDIF
 
@@ -435 +442 @@
 !
 !--       Transfer initial profiles to the arrays of the 3D model
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 e(:,j,i)  = e1d
                 kh(:,j,i) = kh1d
@@ -447 +454 @@
 
           IF ( humidity  .OR.  passive_scalar )  THEN
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    q(:,j,i) = q_init
                 ENDDO
@@ -455 +462 @@
 
           IF ( .NOT. constant_diffusion )  THEN
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    e(:,j,i)  = e1d
                 ENDDO
@@ -505 +512 @@
                 ENDDO
              ENDDO
-             IF ( conserve_volume_flow )  THEN
-                IF ( nxr == nx )  THEN
-                   DO  j = nys, nyn
-                      DO  k = nzb + 1, nzb_u_inner(j,nx)
-                         u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &
-                                               u1d(k) * dzu(k)
-                      ENDDO
-                   ENDDO
-                ENDIF
-                IF ( nyn == ny )  THEN
-                   DO  i = nxl, nxr
-                      DO  k = nzb + 1, nzb_v_inner(ny,i)
-                         v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &
-                                               v1d(k) * dzu(k)
-                      ENDDO
-                   ENDDO
-                ENDIF
-             ENDIF
+             
 !
 !--          WARNING: The extra boundary conditions set after running the
@@ -530 +520 @@
 !--          ---------  advection scheme: keep u and v zero one layer below
 !--                     the topography.
-             IF ( ibc_uv_b == 0 )  THEN
-!
-!--             Satisfying the Dirichlet condition with an extra layer below
-!--             the surface where the u and v component change their sign.
-                DO  i = nxl-1, nxr+1
-                   DO  j = nys-1, nyn+1
-                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
-                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
-                   ENDDO
-                ENDDO
-
-             ELSE
+!
+!--           Following was removed, because mirror boundary condition are
+!--           replaced by dirichlet boundary conditions 
+!
+!             IF ( ibc_uv_b == 0 )  THEN
+!!
+!!--             Satisfying the Dirichlet condition with an extra layer below
+!!--             the surface where the u and v component change their sign.
+!                DO  i = nxl-1, nxr+1
+!                   DO  j = nys-1, nyn+1
+!                      IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)
+!                      IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)
+!                   ENDDO
+!                ENDDO
+!
+!             ELSE
+             IF ( ibc_uv_b == 1 )  THEN
 !
 !--             Neumann condition
@@ -560 +555 @@
 !--       Use constructed initial profiles (velocity constant with height,
 !--       temperature profile with constant gradient)
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 pt(:,j,i) = pt_init
                 u(:,j,i)  = u_init
@@ -574 +569 @@
 !--       in the limiting formula!). The original values are stored to be later
 !--       used for volume flow control.
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0
                 v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0
              ENDDO
           ENDDO
-          IF ( conserve_volume_flow )  THEN
-             IF ( nxr == nx )  THEN
-                DO  j = nys, nyn
-                   DO  k = nzb + 1, nzb_u_inner(j,nx) + 1
-                      u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &
-                                            u_init(k) * dzu(k)
-                   ENDDO
-                ENDDO
-             ENDIF
-             IF ( nyn == ny )  THEN
-                DO  i = nxl, nxr
-                   DO  k = nzb + 1, nzb_v_inner(ny,i) + 1
-                      v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &
-                                            v_init(k) * dzu(k)
-                   ENDDO
-                ENDDO
-             ENDIF
-          ENDIF
 
           IF ( humidity  .OR.  passive_scalar )  THEN
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    q(:,j,i) = q_init
                 ENDDO
@@ -608 +585 @@
 
           IF ( ocean )  THEN
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    sa(:,j,i) = sa_init
                 ENDDO
@@ -660 +637 @@
 
        ENDIF
+!
+!--    Bottom boundary
+       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2  )  THEN
+          u(nzb,:,:) = 0.0
+          v(nzb,:,:) = 0.0
+       ENDIF
 
 !
@@ -683 +666 @@
        hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 )
        hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 )
-       IF ( ibc_uv_b == 0 )  THEN
-          hom(nzb,1,5,:) = -hom(nzb+1,1,5,:)  ! due to satisfying the Dirichlet 
-          hom(nzb,1,6,:) = -hom(nzb+1,1,6,:)  ! condition with an extra layer 
-              ! below the surface where the u and v component change their sign
+       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2)  THEN
+          hom(nzb,1,5,:) = 0.0   
+          hom(nzb,1,6,:) = 0.0  
        ENDIF
        hom(:,1,7,:)  = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 )
@@ -733 +715 @@
 !--             Over topography surface_heatflux is replaced by wall_heatflux(0)
                 IF ( TRIM( topography ) /= 'flat' )  THEN
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          IF ( nzb_s_inner(j,i) /= 0 )  THEN
                             shf(j,i) = wall_heatflux(0)
@@ -755 +737 @@
                 IF ( TRIM( topography ) /= 'flat' )  THEN
                    wall_qflux = wall_humidityflux
-                   DO  i = nxl-1, nxr+1
-                      DO  j = nys-1, nyn+1
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
                          IF ( nzb_s_inner(j,i) /= 0 )  THEN
                             qsws(j,i) = wall_qflux(0)
@@ -827 +809 @@
        ENDIF
 
-!
-!--    Calculate the initial volume flow at the right and north boundary
-       IF ( conserve_volume_flow )  THEN
-
-          volume_flow_initial_l = 0.0
-          volume_flow_area_l    = 0.0
-  
-          IF ( nxr == nx )  THEN
-             DO  j = nys, nyn
-                DO  k = nzb_2d(j,nx) + 1, nzt
-                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
-                                              u(k,j,nx) * dzu(k)
-                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzu(k)
-                ENDDO
-!
-!--             Correction if velocity at nzb+1 has been set zero further above
-                volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
-                                           u_nzb_p1_for_vfc(j)
-             ENDDO
-          ENDIF
-
-          IF ( nyn == ny )  THEN
-             DO  i = nxl, nxr
-                DO  k = nzb_2d(ny,i) + 1, nzt  
-                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
-                                              v(k,ny,i) * dzu(k)
-                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzu(k)
-                ENDDO
-!
-!--             Correction if velocity at nzb+1 has been set zero further above
-                volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
-                                           v_nzb_p1_for_vfc(i)
-             ENDDO
-          ENDIF
-
-#if defined( __parallel )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
-                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
-                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
-#else
-          volume_flow_initial = volume_flow_initial_l
-          volume_flow_area    = volume_flow_area_l
-#endif
-!
-!--       In case of 'bulk_velocity' mode, volume_flow_initial is overridden
-!--       and calculated from u|v_bulk instead.
-          IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
-             volume_flow_initial(1) = u_bulk * volume_flow_area(1)
-             volume_flow_initial(2) = v_bulk * volume_flow_area(2)
-          ENDIF
-
-       ENDIF
 
 !
@@ -968 +895 @@
           sa_p  = sa
        ENDIF
-
+       
 
     ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.    &
-             TRIM( initializing_actions ) == 'cyclic_fill' )  &
+         TRIM( initializing_actions ) == 'cyclic_fill' )  &
     THEN
 !
@@ -978 +905 @@
        IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
 
-          WRITE (9,*) 'before read_parts_of_var_list'
-          CALL local_flush( 9 )
           CALL read_parts_of_var_list
-          WRITE (9,*) 'after read_parts_of_var_list'
-          CALL local_flush( 9 )
           CALL close_file( 13 )
 
@@ -1002 +925 @@
 !--          conditions are used)
              IF ( inflow_l )  THEN
-                DO  j = nys-1, nyn+1
+                DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      u(k,j,-1)  = mean_inflow_profiles(k,1)
-                      v(k,j,-1)  = mean_inflow_profiles(k,2)
-                      w(k,j,-1)  = 0.0
-                      pt(k,j,-1) = mean_inflow_profiles(k,4)
-                      e(k,j,-1)  = mean_inflow_profiles(k,5)
+                      u(k,j,nxlg:-1)  = mean_inflow_profiles(k,1)
+                      v(k,j,nxlg:-1)  = mean_inflow_profiles(k,2)
+                      w(k,j,nxlg:-1)  = 0.0
+                      pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)
+                      e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
                    ENDDO
                 ENDDO
@@ -1064 +987 @@
 !
 !--    Read binary data from restart file
-          WRITE (9,*) 'before read_3d_binary'
-          CALL local_flush( 9 )
+
        CALL read_3d_binary
-          WRITE (9,*) 'after read_3d_binary'
-          CALL local_flush( 9 )
 
 !
@@ -1074 +994 @@
        IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.  &
             topography /= 'flat' )  THEN
-!
-!--       Correction of initial volume flow
-          IF ( conserve_volume_flow )  THEN
-             IF ( nxr == nx )  THEN
-                DO  j = nys, nyn
-                   DO  k = nzb + 1, nzb_u_inner(j,nx)
-                      u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &
-                                            u(k,j,nx) * dzu(k)
-                   ENDDO
-                ENDDO
-             ENDIF
-             IF ( nyn == ny )  THEN
-                DO  i = nxl, nxr
-                   DO  k = nzb + 1, nzb_v_inner(ny,i)
-                      v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &
-                                            v(k,ny,i) * dzu(k)
-                   ENDDO
-                ENDDO
-             ENDIF
-          ENDIF
-
 !
 !--       Inside buildings set velocities and TKE back to zero.
@@ -1100 +999 @@
 !--       maybe revise later.
           IF ( timestep_scheme(1:5) == 'runge' )  THEN
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    u  (nzb:nzb_u_inner(j,i),j,i) = 0.0
                    v  (nzb:nzb_v_inner(j,i),j,i) = 0.0
@@ -1118 +1017 @@
              ENDDO
           ELSE
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    u  (nzb:nzb_u_inner(j,i),j,i) = 0.0
                    v  (nzb:nzb_v_inner(j,i),j,i) = 0.0
@@ -1139 +1038 @@
 
 !
-!--    Calculate the initial volume flow at the right and north boundary
-       IF ( conserve_volume_flow  .AND.  &
-            TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
-
-          volume_flow_initial_l = 0.0
-          volume_flow_area_l    = 0.0
-  
-          IF ( nxr == nx )  THEN
-             DO  j = nys, nyn
-                DO  k = nzb_2d(j,nx) + 1, nzt
-                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
-                                              u(k,j,nx) * dzu(k)
-                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzu(k)
-                ENDDO
-!
-!--             Correction if velocity inside buildings has been set to zero
-!--             further above
-                volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
-                                           u_nzb_p1_for_vfc(j)
-             ENDDO
-          ENDIF
-
-          IF ( nyn == ny )  THEN
-             DO  i = nxl, nxr
-                DO  k = nzb_2d(ny,i) + 1, nzt  
-                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
-                                              v(k,ny,i) * dzu(k)
-                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzu(k)
-                ENDDO
-!
-!--             Correction if velocity inside buildings has been set to zero
-!--             further above
-                volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
-                                           v_nzb_p1_for_vfc(i)
-             ENDDO
-          ENDIF
-
-#if defined( __parallel )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
-                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
-                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
-#else
-          volume_flow_initial = volume_flow_initial_l
-          volume_flow_area    = volume_flow_area_l
-#endif  
-       ENDIF
-
-
-!
 !--    Calculate initial temperature field and other constants used in case
 !--    of a sloping surface
@@ -1243 +1090 @@
           w_m_n(:,:,:) = w(:,ny-1:ny,:)
        ENDIF
-
-    ENDIF
-
+       
+    ENDIF
+!
+!-- Calculate the initial volume flow at the right and north boundary
+    IF ( conserve_volume_flow ) THEN
+
+       volume_flow_initial_l = 0.0
+       volume_flow_area_l    = 0.0
+  
+       IF  ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
+
+          IF ( nxr == nx )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb_2d(j,nx) + 1, nzt
+                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
+                                              hom_sum(k,1,0) * dzw(k)
+                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
+                ENDDO
+             ENDDO
+          ENDIF
+          
+          IF ( nyn == ny )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb_2d(ny,i) + 1, nzt  
+                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
+                                               hom_sum(k,2,0) * dzw(k)
+                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
+                ENDDO
+             ENDDO
+          ENDIF
+
+       ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+
+          IF ( nxr == nx )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb_2d(j,nx) + 1, nzt
+                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
+                                               u(k,j,nx) * dzw(k)
+                   volume_flow_area_l(1)    = volume_flow_area_l(1) + dzw(k)
+                ENDDO
+             ENDDO
+          ENDIF
+          
+          IF ( nyn == ny )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb_2d(ny,i) + 1, nzt  
+                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
+                                              v(k,ny,i) * dzw(k)
+                   volume_flow_area_l(2)    = volume_flow_area_l(2) + dzw(k)
+                ENDDO
+             ENDDO
+          ENDIF
+
+       ENDIF
+
+#if defined( __parallel )
+          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
+                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
+          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
+                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
+
+          CALL MPI_ALLREDUCE( volume_flow_initial_l(2), volume_flow_initial(2),&
+                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
+          CALL MPI_ALLREDUCE( volume_flow_area_l(2), volume_flow_area(2),      &
+                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
+
+#else
+          volume_flow_initial = volume_flow_initial_l
+          volume_flow_area    = volume_flow_area_l
+#endif  
+
+!
+!--       In case of 'bulk_velocity' mode, volume_flow_initial is overridden
+!--       and calculated from u|v_bulk instead.
+          IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
+             volume_flow_initial(1) = u_bulk * volume_flow_area(1)
+             volume_flow_initial(2) = v_bulk * volume_flow_area(2)
+          ENDIF
+
+       ENDIF
 !
 !-- Initialization of the leaf area density
@@ -1254 +1178 @@
           CASE( 'block' )
 
-             DO  i = nxl-1, nxr+1
-                DO  j = nys-1, nyn+1
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
                    lad_s(:,j,i) = lad(:)
                    cdc(:,j,i)   = drag_coefficient
@@ -1277 +1201 @@
           END SELECT
 
-       CALL exchange_horiz( lad_s )
-       CALL exchange_horiz( cdc )
+       CALL exchange_horiz( lad_s, nbgp )
+       CALL exchange_horiz( cdc, nbgp )
 
        IF ( passive_scalar ) THEN
-          CALL exchange_horiz( sls )
-          CALL exchange_horiz( sec )
+          CALL exchange_horiz( sls, nbgp )
+          CALL exchange_horiz( sec, nbgp )
        ENDIF
 
@@ -1311 +1235 @@
        lad_w(nzt+1,:,:)     = lad_w(nzt,:,:)
 
-       CALL exchange_horiz( lad_u )
-       CALL exchange_horiz( lad_v )
-       CALL exchange_horiz( lad_w )
+       CALL exchange_horiz( lad_u, nbgp )
+       CALL exchange_horiz( lad_v, nbgp )
+       CALL exchange_horiz( lad_w, nbgp )
 
 !
@@ -1322 +1246 @@
 !--       integration of the leaf area density
           lai(:,:,:) = 0.0
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 DO  k = pch_index-1, 0, -1
                    lai(k,j,i) = lai(k+1,j,i) +                   &
@@ -1339 +1263 @@
 !--       Evaluation of the upward kinematic vertical heat flux within the 
 !--       canopy
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 DO  k = 0, pch_index
                    canopy_heat_flux(k,j,i) = cthf *                    &
@@ -1384 +1308 @@
 !-- Initialize quantities for special advections schemes
     CALL init_advec
+    IF ( momentum_advec == 'ws-scheme' .OR.  &
+         scalar_advec == 'ws-scheme' ) CALL ws_init
 
 !
@@ -1439 +1365 @@
     IF ( bc_lr == 'dirichlet/radiation' )  THEN
 
-       DO  i = nxl-1, nxr+1
+       DO  i = nxlg, nxrg
           IF ( i >= nx - outflow_damping_width )  THEN
              km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
@@ -1452 +1378 @@
     ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
 
-       DO  i = nxl-1, nxr+1
+       DO  i = nxlg, nxrg
           IF ( i <= outflow_damping_width )  THEN
              km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
@@ -1467 +1393 @@
     IF ( bc_ns == 'dirichlet/radiation' )  THEN
 
-       DO  j = nys-1, nyn+1
+       DO  j = nysg, nyng
           IF ( j >= ny - outflow_damping_width )  THEN
              km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
@@ -1480 +1406 @@
     ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
 
-       DO  j = nys-1, nyn+1
+       DO  j = nysg, nyng
           IF ( j <= outflow_damping_width )  THEN
              km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
@@ -1594 +1520 @@
 !-- buoyancy, etc. A zero value will occur for cases where all grid points of
 !-- the respective subdomain lie below the surface topography
-    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   )
+    ngp_2dh_outer   = MAX( 1, ngp_2dh_outer(:,:)   ) 
     ngp_3d_inner    = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )),            &
                            ngp_3d_inner(:) )
-    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) )
+    ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 
 
     DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l, ngp_3d_inner_tmp )

palm/trunk/SOURCE/init_advec.f90

@@ -8 +8 @@
 ! Former revisions:
 ! -----------------
-! $Id$
+! $Id: init_advec.f90 484 2010-02-05 07:36:54Z raasch 
 ! RCS Log replace by Id keyword, revision history cleaned up
 !

palm/trunk/SOURCE/init_coupling.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+!
+! determination of target_id's moved to init_pegrid
+!
 !
 ! Former revisions:
@@ -24 +26 @@
     USE pegrid
     USE control_parameters
+    USE indices
 
     IMPLICIT NONE
@@ -48 +51 @@
 !-- ATTENTION: numprocs will be reset according to the new communicators
 #if defined ( __parallel )
+
+!myid_absolut = myid
     IF ( myid == 0 )  THEN
        READ (*,*,ERR=10,END=10)  coupling_mode, bc_data(1), bc_data(2)
@@ -94 +99 @@
 
        IF ( myid < bc_data(1) ) THEN
-          inter_color = 0
-          numprocs = bc_data(1)
+          inter_color     = 0
+          numprocs        = bc_data(1)
+          coupling_mode   = 'atmosphere_to_ocean'
        ELSE
-          inter_color = 1
-          numprocs = bc_data(2)
+          inter_color     = 1
+          numprocs        = bc_data(2)
+          coupling_mode   = 'ocean_to_atmosphere'
        ENDIF
-!
-!--    Calculate the target PE for coupling and set up the new communicators.
-!--    Currently only 1:1 topologies are supported.
-       target_id = myid - ISIGN( numprocs, inter_color - 1 )
-       IF ( inter_color == 0 ) THEN
-          coupling_mode = 'atmosphere_to_ocean'
-       ELSE
-          coupling_mode = 'ocean_to_atmosphere'
-       ENDIF
+
        CALL MPI_COMM_SPLIT( MPI_COMM_WORLD, inter_color, 0, comm_palm, ierr )
        comm2d = comm_palm

palm/trunk/SOURCE/init_grid.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! 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.
 !
 ! Former revisions:
@@ -76 +81 @@
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::  distance
-
+    
+!
+!   Computation of the array bounds.
+    nxlg = nxl - nbgp
+    nxrg = nxr + nbgp
+    nysg = nys - nbgp
+    nyng = nyn + nbgp
 !
 !-- Allocate grid arrays
     ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &
-              dzw(1:nzt+1), l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1) )
+              dzw(1:nzt+1), l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) )
 
 !
@@ -97 +108 @@
 !
 !--    Grid for atmosphere with surface at z=0 (k=0, w-grid).
-!--    Since the w-level lies on the surface, the first u-level (staggered!)
-!--    lies below the surface (used for "mirror" boundary condition).
 !--    The first u-level above the surface corresponds to the top of the 
 !--    Prandtl-layer.
-       zu(0) = - dz * 0.5
+
+       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
+          zu(0) = 0.0
+      !    zu(0) = - dz * 0.5
+       ELSE
+          zu(0) = - dz * 0.5
+       ENDIF
        zu(1) =   dz * 0.5
 
@@ -173 +188 @@
        dd2zu(k) = 1.0 / ( dzu(k) + dzu(k+1) )
     ENDDO
+    
+!    
+!-- In case of FFT method or SOR swap inverse grid lenght ddzu to ddzu_fft and 
+!-- modify the lowest entry. This is necessary for atmosphere runs, because
+!-- zu(0) and so ddzu(1) changed. Accompanied with this modified arrays 
+!-- pressure solver uses wrong values and this causes kinks in the profiles
+!-- of turbulent quantities.  
+    IF ( psolver /= 'multigrid' ) THEN
+       ALLOCATE( ddzu_pres(1:nzt+1) )
+       ddzu_pres = ddzu
+       IF( .NOT. ocean ) ddzu_pres(1) = ddzu_pres(2)
+    ENDIF   
 
 !
@@ -187 +214 @@
 
        dzu_mg(:,maximum_grid_level) = dzu
+!       
+!--    To ensure a equally spaced grid. For ocean runs this is not necessary, 
+!--    because zu(0) does not changed so far. Also this would cause errors 
+!--    if a vertical stretching for ocean runs is used.
+       IF ( .NOT. ocean ) dzu_mg(1,maximum_grid_level) = dzu(2)
        dzw_mg(:,maximum_grid_level) = dzw
        nzt_l = nzt
@@ -239 +271 @@
 !-- the flag arrays needed for the multigrid method
     gls = 2**( maximum_grid_level )
+
     ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr),       &
               corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr),       &
-              nzb_local(-gls:ny+gls,-gls:nx+gls), nzb_tmp(-1:ny+1,-1:nx+1), &
+              nzb_local(-gls:ny+gls,-gls:nx+gls),                                   &
+              nzb_tmp(-nbgp:ny+nbgp,-nbgp:nx+nbgp),                         &
               wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr),             &
               wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
-    ALLOCATE( fwxm(nys-1:nyn+1,nxl-1:nxr+1), fwxp(nys-1:nyn+1,nxl-1:nxr+1), &
-              fwym(nys-1:nyn+1,nxl-1:nxr+1), fwyp(nys-1:nyn+1,nxl-1:nxr+1), &
-              fxm(nys-1:nyn+1,nxl-1:nxr+1), fxp(nys-1:nyn+1,nxl-1:nxr+1),   &
-              fym(nys-1:nyn+1,nxl-1:nxr+1), fyp(nys-1:nyn+1,nxl-1:nxr+1),   &
-              nzb_s_inner(nys-1:nyn+1,nxl-1:nxr+1),                         &
-              nzb_s_outer(nys-1:nyn+1,nxl-1:nxr+1),                         &
-              nzb_u_inner(nys-1:nyn+1,nxl-1:nxr+1),                         &
-              nzb_u_outer(nys-1:nyn+1,nxl-1:nxr+1),                         &
-              nzb_v_inner(nys-1:nyn+1,nxl-1:nxr+1),                         &
-              nzb_v_outer(nys-1:nyn+1,nxl-1:nxr+1),                         &
-              nzb_w_inner(nys-1:nyn+1,nxl-1:nxr+1),                         &
-              nzb_w_outer(nys-1:nyn+1,nxl-1:nxr+1),                         &
-              nzb_diff_s_inner(nys-1:nyn+1,nxl-1:nxr+1),                    &
-              nzb_diff_s_outer(nys-1:nyn+1,nxl-1:nxr+1),                    &
-              nzb_diff_u(nys-1:nyn+1,nxl-1:nxr+1),                          &
-              nzb_diff_v(nys-1:nyn+1,nxl-1:nxr+1),                          &
-              nzb_2d(nys-1:nyn+1,nxl-1:nxr+1),                              &
-              wall_e_x(nys-1:nyn+1,nxl-1:nxr+1),                            &
-              wall_e_y(nys-1:nyn+1,nxl-1:nxr+1),                            &
-              wall_u(nys-1:nyn+1,nxl-1:nxr+1),                              &
-              wall_v(nys-1:nyn+1,nxl-1:nxr+1),                              &
-              wall_w_x(nys-1:nyn+1,nxl-1:nxr+1),                            &
-              wall_w_y(nys-1:nyn+1,nxl-1:nxr+1) )
-
-    ALLOCATE( l_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+    ALLOCATE( fwxm(nysg:nyng,nxlg:nxrg), fwxp(nysg:nyng,nxlg:nxrg),         &
+              fwym(nysg:nyng,nxlg:nxrg), fwyp(nysg:nyng,nxlg:nxrg),         &
+              fxm(nysg:nyng,nxlg:nxrg), fxp(nysg:nyng,nxlg:nxrg),           &
+              fym(nysg:nyng,nxlg:nxrg), fyp(nysg:nyng,nxlg:nxrg),           &
+              nzb_s_inner(nysg:nyng,nxlg:nxrg),                             &
+              nzb_s_outer(nysg:nyng,nxlg:nxrg),                             &
+              nzb_u_inner(nysg:nyng,nxlg:nxrg),                             &
+              nzb_u_outer(nysg:nyng,nxlg:nxrg),                             &
+              nzb_v_inner(nysg:nyng,nxlg:nxrg),                             &
+              nzb_v_outer(nysg:nyng,nxlg:nxrg),                             &
+              nzb_w_inner(nysg:nyng,nxlg:nxrg),                             &
+              nzb_w_outer(nysg:nyng,nxlg:nxrg),                             &
+              nzb_diff_s_inner(nysg:nyng,nxlg:nxrg),                        &
+              nzb_diff_s_outer(nysg:nyng,nxlg:nxrg),                        &
+              nzb_diff_u(nysg:nyng,nxlg:nxrg),                              &
+              nzb_diff_v(nysg:nyng,nxlg:nxrg),                              &
+              nzb_2d(nysg:nyng,nxlg:nxrg),                                  &
+              wall_e_x(nysg:nyng,nxlg:nxrg),                                &
+              wall_e_y(nysg:nyng,nxlg:nxrg),                                &
+              wall_u(nysg:nyng,nxlg:nxrg),                                  &
+              wall_v(nysg:nyng,nxlg:nxrg),                                  &
+              wall_w_x(nysg:nyng,nxlg:nxrg),                                &
+              wall_w_y(nysg:nyng,nxlg:nxrg) )
+
+
+
+    ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 
     nzb_s_inner = nzb;  nzb_s_outer = nzb
@@ -327 +363 @@
     vertical_influence(0) = vertical_influence(1)
 
-    DO  i = nxl-1, nxr+1
-       DO  j = nys-1, nyn+1
+    DO  i = nxlg, nxrg
+       DO  j = nysg, nyng
           DO  k = nzb_s_inner(j,i) + 1, &
                   nzb_s_inner(j,i) + vertical_influence(nzb_s_inner(j,i))
@@ -477 +513 @@
           nzb_local(:,-gls:-1)        = nzb_local(:,nx-gls+1:nx)
           nzb_local(:,nx+1:nx+gls)    = nzb_local(:,0:gls-1)
+
+
       
           GOTO 12
@@ -588 +626 @@
 !--    nzb_s_outer: 
 !--    extend nzb_local east-/westwards first, then north-/southwards
-       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
+       nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
        DO  j = -1, ny + 1
           DO  i = 0, nx
@@ -620 +658 @@
 !--    nzb_u_inner: 
 !--    extend nzb_local rightwards only
-       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
+       nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
        DO  j = -1, ny + 1
           DO  i = 0, nx + 1
@@ -626 +664 @@
           ENDDO
        ENDDO
-       nzb_u_inner = nzb_tmp(nys-1:nyn+1,nxl-1:nxr+1)
+       nzb_u_inner = nzb_tmp(nysg:nyng,nxlg:nxrg)
 
 !
@@ -652 +690 @@
 !--    nzb_v_inner:
 !--    extend nzb_local northwards only
-       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
+       nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
        DO  i = -1, nx + 1
           DO  j = 0, ny + 1
@@ -658 +696 @@
           ENDDO
        ENDDO
-       nzb_v_inner = nzb_tmp(nys-1:nyn+1,nxl-1:nxr+1)
+       nzb_v_inner = nzb_tmp(nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp)
 
 !
@@ -1096 +1134 @@
 !
 !-- Need to set lateral boundary conditions for l_wall
-    CALL exchange_horiz( l_wall )
+
+    CALL exchange_horiz( l_wall, nbgp )
 
     DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, &

palm/trunk/SOURCE/init_particles.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation 
+! of arrays.
 !
 ! Former revisions:
@@ -185 +186 @@
        ENDIF
 
-       ALLOCATE( prt_count(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),       &
-                 prt_start_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
+       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                 prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
                  particle_mask(maximum_number_of_particles),         &
                  particles(maximum_number_of_particles) )
@@ -209 +210 @@
 !--    particles, which can be also periodically released at later times.
 !--    Also allocate array for particle tail coordinates, if needed.
-       ALLOCATE( prt_count(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),       &
-                 prt_start_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
+       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                 prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
                  particle_mask(maximum_number_of_particles),         &
                  particles(maximum_number_of_particles) )

palm/trunk/SOURCE/init_pegrid.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
+!
+! 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!!!!!!!!
@@ -79 +94 @@
 
 
+
     IMPLICIT NONE
 
@@ -88 +104 @@
 
     INTEGER, DIMENSION(:), ALLOCATABLE ::  ind_all, nxlf, nxrf, nynf, nysf
+
+    INTEGER, DIMENSION(2) :: pdims_remote
 
     LOGICAL ::  found
@@ -103 +121 @@
 
 #if defined( __parallel )
+
 !
 !-- Determine the processor topology or check it, if prescribed by the user
@@ -624 +643 @@
 #endif
 
+! 
+!-- Determine the number of ghost points
+    IF (scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme') THEN
+       nbgp = 3
+    ELSE
+       nbgp = 1
+    END IF 
+
 !
 !-- In case of coupled runs, create a new MPI derived datatype for the
 !-- exchange of surface (xy) data .
 !-- Gridpoint number for the exchange of ghost points (xy-plane)
-    ngp_xy  = ( nxr - nxl + 3 ) * ( nyn - nys + 3 )
+
+    ngp_xy  = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
 
 !
@@ -635 +663 @@
     CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
     CALL MPI_TYPE_COMMIT( type_xy, ierr )
+
+
+    IF ( TRIM( coupling_mode ) .NE. 'uncoupled' ) THEN
+   
+!
+!--    Pass the number of grid points of the atmosphere model to
+!--    the ocean model and vice versa
+       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
+
+          nx_a = nx
+          ny_a = ny
+
+          IF ( myid == 0 ) THEN
+             CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, &
+                            comm_inter, ierr )
+             CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, &
+                            comm_inter, ierr )
+             CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, &
+                            comm_inter, ierr )
+             CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, &
+                            comm_inter, status, ierr )
+             CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, &
+                            comm_inter, status, ierr )
+             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &
+                            comm_inter, status, ierr )
+          ENDIF
+
+          CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr)
+          CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr) 
+          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
+       
+       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
+
+          nx_o = nx
+          ny_o = ny 
+
+          IF ( myid == 0 ) THEN
+             CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, &
+                            comm_inter, status, ierr )
+             CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, &
+                            comm_inter, status, ierr )
+             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, &
+                            comm_inter, status, ierr )
+             CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, &
+                            comm_inter, ierr )
+             CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, &
+                            comm_inter, ierr )
+             CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, &
+                            comm_inter, ierr )
+          ENDIF
+
+          CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
+          CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr) 
+          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr) 
+
+       ENDIF
+  
+       ngp_a = (nx_a+1+2*nbgp)*(ny_a+1+2*nbgp)
+       ngp_o = (nx_o+1+2*nbgp)*(ny_o+1+2*nbgp)
+
+!
+!--    determine if the horizontal grid and the number of PEs 
+!--    in ocean and atmosphere is same or not
+!--    (different number of PEs still not implemented)
+       IF ( nx_o == nx_a .AND. ny_o == ny_a .AND.  &
+            pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
+       THEN
+          coupling_topology = 0
+       ELSE
+          coupling_topology = 1
+       ENDIF 
+
+!
+!--    Determine the target PEs for the exchange between ocean and
+!--    atmosphere (comm2d)
+       IF ( coupling_topology == 0) THEN
+          IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
+             target_id = myid + numprocs
+          ELSE
+             target_id = myid 
+          ENDIF
+
+       ELSE
+
+!
+!--       In case of nonequivalent topology in ocean and atmosphere only for
+!--       PE0 in ocean and PE0 in atmosphere a target_id is needed, since
+!--       data echxchange between ocean and atmosphere will be done only by
+!--       those PEs.    
+          IF ( myid == 0 ) THEN
+             IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
+                target_id = numprocs 
+             ELSE
+                target_id = 0
+             ENDIF
+ print*, coupling_mode, myid, " -> ", target_id, "numprocs: ", numprocs 
+          ENDIF
+       ENDIF
+
+    ENDIF
+
+
 #endif
 
@@ -854 +984 @@
     ELSE
 
-       maximum_grid_level = 1
+       maximum_grid_level = 0
 
     ENDIF
@@ -863 +993 @@
 !
 !-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
-    ngp_y  = nyn - nys + 1
+    ngp_y  = nyn - nys + 1 + 2 * nbgp
 
 !
 !-- Define a new MPI derived datatype for the exchange of ghost points in
 !-- y-direction for 2D-arrays (line)