Home


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

Changeset 667


Ignore:
Timestamp:
Dec 23, 2010 12:06:00 PM (14 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:
66 edited
1 copied

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE

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

  • palm/trunk/SOURCE/Makefile

    r482 r667  
    44# Current revisions:
    55# ------------------
    6 #
     6# +advec_ws
    77#
    88# Former revisions:
    … …  
    5757
    5858RCS = advec_particles.f90 advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 \
    59         advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 advec_u_ups.f90 \
    60         advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 advec_w_pw.f90 \
    61         advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 average_3d_data.f90 \
    62         boundary_conds.f90 buoyancy.f90 calc_liquid_water_content.f90 \
    63         calc_precipitation.f90 calc_radiation.f90 calc_spectra.f90 \
    64         check_for_restart.f90 check_open.f90 check_parameters.f90 \
    65         close_file.f90 compute_vpt.f90 coriolis.f90 cpu_log.f90 \
    66         cpu_statistics.f90 data_log.f90 data_output_dvrp.f90 \
    67         data_output_mask.f90 data_output_profiles.f90 data_output_ptseries.f90 \
    68         data_output_spectra.f90 data_output_tseries.f90 data_output_2d.f90 \
    69         data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 diffusion_u.f90 \
    70         diffusion_v.f90 diffusion_w.f90 diffusivities.f90 disturb_field.f90 \
    71         disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 \
    72         exchange_horiz_2d.f90 \
     59        advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \
     60        advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \
     61        advec_w_pw.f90 advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 \
     62        average_3d_data.f90 boundary_conds.f90 buoyancy.f90 \
     63        calc_liquid_water_content.f90 calc_precipitation.f90 \
     64        calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \
     65        check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \
     66        coriolis.f90 cpu_log.f90 cpu_statistics.f90 data_log.f90 \
     67        data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \
     68        data_output_ptseries.f90 data_output_spectra.f90 data_output_tseries.f90 \
     69        data_output_2d.f90 data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 \
     70        diffusion_u.f90 diffusion_v.f90 diffusion_w.f90 diffusivities.f90 \
     71        disturb_field.f90 disturb_heatflux.f90 eqn_state_seawater.f90 \
     72        exchange_horiz.f90 exchange_horiz_2d.f90 \
    7373        fft_xy.f90 flow_statistics.f90 global_min_max.f90 \
    7474        header.f90 impact_of_latent_heat.f90 inflow_turbulence.f90 \
    … …  
    7777        init_masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \
    7878        init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
    79         interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90 \
    80         local_stop.f90 local_system.f90 local_tremain.f90 \
     79        interaction_droplets_ptq.f90 local_flush.f90 \
     80        local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
    8181        local_tremain_ini.f90 message.f90 modules.f90 netcdf.f90 \
    8282        package_parin.f90 palm.f90 parin.f90 particle_boundary_conds.f90 \
    … …  
    106106OBJS =  advec_particles.o advec_s_bc.o advec_s_pw.o advec_s_up.o \
    107107        advec_s_ups.o advec_u_pw.o advec_u_up.o advec_u_ups.o \
    108         advec_v_pw.o advec_v_up.o advec_v_ups.o advec_w_pw.o \
     108        advec_ws.o advec_v_pw.o advec_v_up.o advec_v_ups.o advec_w_pw.o \
    109109        advec_w_up.o advec_w_ups.o asselin_filter.o average_3d_data.o \
    110110        boundary_conds.o buoyancy.o calc_liquid_water_content.o \
    … …  
    184184advec_v_up.o: modules.o
    185185advec_v_ups.o: modules.o
     186advec_ws.o: modules.o
    186187advec_w_pw.o: modules.o
    187188advec_w_up.o: modules.o
    … …  
    230231inflow_turbulence.o: modules.o
    231232init_1d_model.o: modules.o
    232 init_3d_model.o: modules.o random_function.o
     233init_3d_model.o: modules.o random_function.o advec_ws.o
    233234init_advec.o: modules.o
    234235init_cloud_physics.o: modules.o
    … …  
    264265production_e.o: modules.o wall_fluxes.o
    265266prognostic_equations.o: modules.o advec_s_pw.o advec_s_up.o advec_u_pw.o \
     267        advec_ws.o \
    266268        advec_u_up.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o  \
    267269        buoyancy.o calc_precipitation.o calc_radiation.o coriolis.o \

  • palm/trunk/SOURCE/advec_particles.f90

    r623 r667  
    44! Current revisions:
    55! -----------------
     6! Declaration of de_dx, de_dy, de_dz adapted to additional
     7! ghost points. Furthermore the calls of exchange_horiz were modified.
    68!
    79! TEST: PRINT statements on unit 9 (commented out)
    … …  
    153155    REAL    ::  location(1:30,1:3)
    154156
    155     REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  de_dx, de_dy, de_dz
     157    REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  de_dx, de_dy, de_dz
    156158
    157159    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  trlpt, trnpt, trrpt, trspt
    … …  
    768770!
    769771!--    Lateral boundary conditions
    770        CALL exchange_horiz( de_dx )
    771        CALL exchange_horiz( de_dy )
    772        CALL exchange_horiz( de_dz )
    773        CALL exchange_horiz( diss  )
     772       CALL exchange_horiz( de_dx, nbgp )
     773       CALL exchange_horiz( de_dy, nbgp )
     774       CALL exchange_horiz( de_dz, nbgp )
     775       CALL exchange_horiz( diss, nbgp  )
    774776
    775777!

  • palm/trunk/SOURCE/advec_w_pw.f90

    r484 r667  
    8989       REAL    ::  gu, gv
    9090
    91 
    9291       gu = 2.0 * u_gtrans
    9392       gv = 2.0 * v_gtrans
    … …  
    103102                                                )
    104103       ENDDO
    105 
    106104    END SUBROUTINE advec_w_pw_ij
    107105

  • palm/trunk/SOURCE/asselin_filter.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    77!
    88! Former revisions:
    … …  
    4949!$OMP PARALLEL PRIVATE (i,j,k)
    5050!$OMP DO
    51     DO  i = nxl-1, nxr+1
    52        DO  j = nys-1, nyn+1
     51    DO  i = nxlg, nxrg
     52       DO  j = nysg, nyng
    5353
    5454          DO  k = nzb_2d(j,i), nzt+1

  • palm/trunk/SOURCE/average_3d_data.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    77!
    88! Former revisions:
    … …  
    5858
    5959          CASE ( 'e' )
    60              DO  i = nxl-1, nxr+1
    61                 DO  j = nys-1, nyn+1
     60             DO  i = nxlg, nxrg
     61                DO  j = nysg, nyng
    6262                   DO  k = nzb, nzt+1
    6363                      e_av(k,j,i) = e_av(k,j,i) / REAL( average_count_3d )
    … …  
    6767
    6868          CASE ( 'qsws*' )
    69              DO  i = nxl-1, nxr+1
    70                 DO  j = nys-1, nyn+1
     69             DO  i = nxlg, nxrg
     70                DO  j = nysg, nyng
    7171                   qsws_av(j,i) = qsws_av(j,i) / REAL( average_count_3d )
    7272                ENDDO
    … …  
    7474
    7575          CASE ( 'lwp*' )
    76              DO  i = nxl-1, nxr+1
    77                 DO  j = nys-1, nyn+1
     76             DO  i = nxlg, nxrg
     77                DO  j = nysg, nyng
    7878                   lwp_av(j,i) = lwp_av(j,i) / REAL( average_count_3d )
    7979                ENDDO
    … …  
    8181
    8282          CASE ( 'p' )
    83              DO  i = nxl-1, nxr+1
    84                 DO  j = nys-1, nyn+1
     83             DO  i = nxlg, nxrg
     84                DO  j = nysg, nyng
    8585                   DO  k = nzb, nzt+1
    8686                      p_av(k,j,i) = p_av(k,j,i) / REAL( average_count_3d )
    … …  
    108108
    109109          CASE ( 'prr*' )
    110              DO  i = nxl-1, nxr+1
    111                 DO  j = nys-1, nyn+1
     110             DO  i = nxlg, nxrg
     111                DO  j = nysg, nyng
    112112                   precipitation_rate_av(j,i) = precipitation_rate_av(j,i) / &
    113113                                                REAL( average_count_3d )
    … …  
    116116
    117117          CASE ( 'pt' )
    118              DO  i = nxl-1, nxr+1
    119                 DO  j = nys-1, nyn+1
     118             DO  i = nxlg, nxrg
     119                DO  j = nysg, nyng
    120120                   DO  k = nzb, nzt+1
    121121                      pt_av(k,j,i) = pt_av(k,j,i) / REAL( average_count_3d )
    … …  
    125125
    126126          CASE ( 'q' )
    127              DO  i = nxl-1, nxr+1
    128                 DO  j = nys-1, nyn+1
     127             DO  i = nxlg, nxrg
     128                DO  j = nysg, nyng
    129129                   DO  k = nzb, nzt+1
    130130                      q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d )
    … …  
    134134
    135135          CASE ( 'ql' )
    136              DO  i = nxl-1, nxr+1
    137                 DO  j = nys-1, nyn+1
     136             DO  i = nxlg, nxrg
     137                DO  j = nysg, nyng
    138138                   DO  k = nzb, nzt+1
    139139                      ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d )
    … …  
    143143
    144144          CASE ( 'ql_c' )
    145              DO  i = nxl-1, nxr+1
    146                 DO  j = nys-1, nyn+1
     145             DO  i = nxlg, nxrg
     146                DO  j = nysg, nyng
    147147                   DO  k = nzb, nzt+1
    148148                      ql_c_av(k,j,i) = ql_c_av(k,j,i) / REAL( average_count_3d )
    … …  
    152152
    153153          CASE ( 'ql_v' )
    154              DO  i = nxl-1, nxr+1
    155                 DO  j = nys-1, nyn+1
     154             DO  i = nxlg, nxrg
     155                DO  j = nysg, nyng
    156156                   DO  k = nzb, nzt+1
    157157                      ql_v_av(k,j,i) = ql_v_av(k,j,i) / REAL( average_count_3d )
    … …  
    161161
    162162          CASE ( 'ql_vp' )
    163              DO  i = nxl-1, nxr+1
    164                 DO  j = nys-1, nyn+1
     163             DO  i = nxlg, nxrg
     164                DO  j = nysg, nyng
    165165                   DO  k = nzb, nzt+1
    166166                      ql_vp_av(k,j,i) = ql_vp_av(k,j,i) / &
    … …  
    171171
    172172          CASE ( 'qv' )
    173              DO  i = nxl-1, nxr+1
    174                 DO  j = nys-1, nyn+1
     173             DO  i = nxlg, nxrg
     174                DO  j = nysg, nyng
    175175                   DO  k = nzb, nzt+1
    176176                      qv_av(k,j,i) = qv_av(k,j,i) / REAL( average_count_3d )
    … …  
    180180
    181181          CASE ( 'rho' )
    182              DO  i = nxl-1, nxr+1
    183                 DO  j = nys-1, nyn+1
     182             DO  i = nxlg, nxrg
     183                DO  j = nysg, nyng
    184184                   DO  k = nzb, nzt+1
    185185                      rho_av(k,j,i) = rho_av(k,j,i) / REAL( average_count_3d )
    … …  
    189189
    190190          CASE ( 's' )
    191              DO  i = nxl-1, nxr+1
    192                 DO  j = nys-1, nyn+1
     191             DO  i = nxlg, nxrg
     192                DO  j = nysg, nyng
    193193                   DO  k = nzb, nzt+1
    194194                      s_av(k,j,i) = s_av(k,j,i) / REAL( average_count_3d )
    … …  
    198198
    199199          CASE ( 'sa' )
    200              DO  i = nxl-1, nxr+1
    201                 DO  j = nys-1, nyn+1
     200             DO  i = nxlg, nxrg
     201                DO  j = nysg, nyng
    202202                   DO  k = nzb, nzt+1
    203203                      sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d )
    … …  
    207207
    208208         CASE ( 'shf*' )
    209              DO  i = nxl-1, nxr+1
    210                 DO  j = nys-1, nyn+1
     209             DO  i = nxlg, nxrg
     210                DO  j = nysg, nyng
    211211                   shf_av(j,i) = shf_av(j,i) / REAL( average_count_3d )
    212212                ENDDO
    … …  
    214214
    215215          CASE ( 't*' )
    216              DO  i = nxl-1, nxr+1
    217                 DO  j = nys-1, nyn+1
     216             DO  i = nxlg, nxrg
     217                DO  j = nysg, nyng
    218218                   ts_av(j,i) = ts_av(j,i) / REAL( average_count_3d )
    219219                ENDDO
    … …  
    221221
    222222          CASE ( 'u' )
    223              DO  i = nxl-1, nxr+1
    224                 DO  j = nys-1, nyn+1
     223             DO  i = nxlg, nxrg
     224                DO  j = nysg, nyng
    225225                   DO  k = nzb, nzt+1
    226226                      u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d )
    … …  
    230230
    231231          CASE ( 'u*' )
    232              DO  i = nxl-1, nxr+1
    233                 DO  j = nys-1, nyn+1
     232             DO  i = nxlg, nxrg
     233                DO  j = nysg, nyng
    234234                   us_av(j,i) = us_av(j,i) / REAL( average_count_3d )
    235235                ENDDO
    … …  
    237237
    238238          CASE ( 'v' )
    239              DO  i = nxl-1, nxr+1
    240                 DO  j = nys-1, nyn+1
     239             DO  i = nxlg, nxrg
     240                DO  j = nysg, nyng
    241241                   DO  k = nzb, nzt+1
    242242                      v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d )
    … …  
    246246
    247247          CASE ( 'vpt' )
    248              DO  i = nxl-1, nxr+1
    249                 DO  j = nys-1, nyn+1
     248             DO  i = nxlg, nxrg
     249                DO  j = nysg, nyng
    250250                   DO  k = nzb, nzt+1
    251251                      vpt_av(k,j,i) = vpt_av(k,j,i) / REAL( average_count_3d )
    … …  
    255255
    256256          CASE ( 'w' )
    257              DO  i = nxl-1, nxr+1
    258                 DO  j = nys-1, nyn+1
     257             DO  i = nxlg, nxrg
     258                DO  j = nysg, nyng
    259259                   DO  k = nzb, nzt+1
    260260                      w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d )
    … …  
    264264
    265265          CASE ( 'z0*' )
    266              DO  i = nxl-1, nxr+1
    267                 DO  j = nys-1, nyn+1
     266             DO  i = nxlg, nxrg
     267                DO  j = nysg, nyng
    268268                   z0_av(j,i) = z0_av(j,i) / REAL( average_count_3d )
    269269                ENDDO

  • palm/trunk/SOURCE/boundary_conds.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
     6!
     7! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
     8!
    69!
    7 !
     10! Removed mirror boundary conditions for u and v at the bottom in case of
     11! ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
     12! in init_3d_model
     13! 
    814! Former revisions:
    915! -----------------
    … …  
    7076    IF ( range == 'main')  THEN
    7177!
    72 !--    Bottom boundary
    73        IF ( ibc_uv_b == 0 )  THEN
    74 !
    75 !--       Satisfying the Dirichlet condition with an extra layer below the
    76 !--       surface where the u and v component change their sign
    77           u_p(nzb,:,:) = -u_p(nzb+1,:,:)
    78           v_p(nzb,:,:) = -v_p(nzb+1,:,:)
    79        ELSE
     78!--    Bottom boundary
     79       IF ( ibc_uv_b == 1 )  THEN
    8080          u_p(nzb,:,:) = u_p(nzb+1,:,:)
    8181          v_p(nzb,:,:) = v_p(nzb+1,:,:)
    8282       ENDIF
    83        DO  i = nxl-1, nxr+1
    84           DO  j = nys-1, nyn+1
     83       DO  i = nxlg, nxrg
     84          DO  j = nysg, nyng
    8585             w_p(nzb_w_inner(j,i),j,i) = 0.0
    8686          ENDDO
    … …  
    9090!--    Top boundary
    9191       IF ( ibc_uv_t == 0 )  THEN
    92           u_p(nzt+1,:,:) = ug(nzt+1)
    93           v_p(nzt+1,:,:) = vg(nzt+1)
     92           u_p(nzt+1,:,:) = ug(nzt+1)
     93           v_p(nzt+1,:,:) = vg(nzt+1)
    9494       ELSE
    95           u_p(nzt+1,:,:) = u_p(nzt,:,:)
    96           v_p(nzt+1,:,:) = v_p(nzt,:,:)
     95           u_p(nzt+1,:,:) = u_p(nzt,:,:)
     96           v_p(nzt+1,:,:) = v_p(nzt,:,:)
    9797       ENDIF
    9898       w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
    … …  
    103103!--    the sea surface temperature of the coupled ocean model.
    104104       IF ( ibc_pt_b == 0 )  THEN
    105           DO  i = nxl-1, nxr+1
    106              DO  j = nys-1, nyn+1
     105          DO  i = nxlg, nxrg
     106             DO  j = nysg, nyng
    107107                pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
    108108             ENDDO
    109109          ENDDO
    110110       ELSEIF ( ibc_pt_b == 1 )  THEN
    111           DO  i = nxl-1, nxr+1
    112              DO  j = nys-1, nyn+1
     111          DO  i = nxlg, nxrg
     112             DO  j = nysg, nyng
    113113                pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
    114114             ENDDO
    … …  
    119119!--    Temperature at top boundary
    120120       IF ( ibc_pt_t == 0 )  THEN
    121           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
     121           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
    122122       ELSEIF ( ibc_pt_t == 1 )  THEN
    123           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
     123           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
    124124       ELSEIF ( ibc_pt_t == 2 )  THEN
    125           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
     125           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
    126126       ENDIF
    127127
    … …  
    130130!--    Generally Neumann conditions with de/dz=0 are assumed
    131131       IF ( .NOT. constant_diffusion )  THEN
    132           DO  i = nxl-1, nxr+1
    133              DO  j = nys-1, nyn+1
     132          DO  i = nxlg, nxrg
     133             DO  j = nysg, nyng
    134134                e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
    135135             ENDDO
    … …  
    144144!--       Bottom boundary: Neumann condition because salinity flux is always
    145145!--       given
    146           DO  i = nxl-1, nxr+1
    147              DO  j = nys-1, nyn+1
     146          DO  i = nxlg, nxrg
     147             DO  j = nysg, nyng
    148148                sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
    149149             ENDDO
    … …  
    153153!--       Top boundary: Dirichlet or Neumann
    154154          IF ( ibc_sa_t == 0 )  THEN
    155              sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
     155              sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
    156156          ELSEIF ( ibc_sa_t == 1 )  THEN
    157              sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
     157              sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
    158158          ENDIF
    159159
    … …  
    167167!--       Surface conditions for constant_humidity_flux
    168168          IF ( ibc_q_b == 0 ) THEN
    169              DO  i = nxl-1, nxr+1
    170                 DO  j = nys-1, nyn+1
     169             DO  i = nxlg, nxrg
     170                DO  j = nysg, nyng
    171171                   q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
    172172                ENDDO
    173173             ENDDO
    174174          ELSE
    175              DO  i = nxl-1, nxr+1
    176                 DO  j = nys-1, nyn+1
     175             DO  i = nxlg, nxrg
     176                DO  j = nysg, nyng
    177177                   q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
    178178                ENDDO
    … …  
    182182!--       Top boundary
    183183          q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
     184
     185
    184186       ENDIF
    185187
    … …  
    226228       c_max = dy / dt_3d
    227229
    228        DO i = nxl-1, nxr+1
     230       DO i = nxlg, nxrg
    229231          DO k = nzb+1, nzt+1
    230232
    … …  
    299301!--    Bottom boundary at the outflow
    300302       IF ( ibc_uv_b == 0 )  THEN
    301           u_p(nzb,-1,:) = -u_p(nzb+1,-1,:)
    302           v_p(nzb,0,:)  = -v_p(nzb+1,0,:) 
     303          u_p(nzb,-1,:) = 0.0
     304          v_p(nzb,0,:)  = 0.0 
    303305       ELSE                   
    304306          u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
    … …  
    324326       c_max = dy / dt_3d
    325327
    326        DO i = nxl-1, nxr+1
     328       DO i = nxlg, nxrg
    327329          DO k = nzb+1, nzt+1
    328330
    … …  
    397399!--    Bottom boundary at the outflow
    398400       IF ( ibc_uv_b == 0 )  THEN
    399           u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:)
    400           v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:) 
     401          u_p(nzb,ny+1,:) = 0.0
     402          v_p(nzb,ny+1,:) = 0.0  
    401403       ELSE                   
    402404          u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
    … …  
    422424       c_max = dx / dt_3d
    423425
    424        DO j = nys-1, nyn+1
     426       DO j = nysg, nyng
    425427          DO k = nzb+1, nzt+1
    426428
    … …  
    495497!--    Bottom boundary at the outflow
    496498       IF ( ibc_uv_b == 0 )  THEN
    497           u_p(nzb,:,-1) = -u_p(nzb+1,:,-1)
    498           v_p(nzb,:,-1) = -v_p(nzb+1,:,-1) 
     499          u_p(nzb,:,0)  = 0.0
     500          v_p(nzb,:,-1) = 0.0
    499501       ELSE                   
    500           u_p(nzb,:,-1) =  u_p(nzb+1,:,-1)
     502          u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
    501503          v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
    502504       ENDIF
    … …  
    520522       c_max = dx / dt_3d
    521523
    522        DO j = nys-1, nyn+1
     524       DO j = nysg, nyng
    523525          DO k = nzb+1, nzt+1
    524526
    … …  
    593595!--    Bottom boundary at the outflow
    594596       IF ( ibc_uv_b == 0 )  THEN
    595           u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1)
    596           v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1)
     597          u_p(nzb,:,nx+1) = 0.0
     598          v_p(nzb,:,nx+1) = 0.0
    597599       ELSE                   
    598600          u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)

  • palm/trunk/SOURCE/calc_liquid_water_content.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    77!
    88! Former revisions:
    … …  
    4747    REAL :: alpha, e_s, q_s, t_l
    4848
    49     DO  i = nxl-1, nxr+1
    50        DO  j = nys-1, nyn+1
     49    DO  i = nxlg, nxrg
     50       DO  j = nysg, nyng
    5151          DO  k = nzb_2d(j,i)+1, nzt
    5252

  • palm/trunk/SOURCE/calc_precipitation.f90

    r484 r667  
    88! Former revisions:
    99! -----------------
    10 ! $Id$
     10! $Id: calc_precipitation.f90 484 2010-02-05 07:36:54Z raasch
    1111!
    1212! 403 2009-10-22 13:57:16Z franke

  • palm/trunk/SOURCE/calc_spectra.f90

    r392 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
     7! of tend
    78!
    89! Former revisions:
    … …  
    152153    IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
    153154       DEALLOCATE( tend )
    154        ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     155       ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    155156    ENDIF
    156157

  • palm/trunk/SOURCE/check_for_restart.f90

    r623 r667  
    55! -----------------
    66!
     7! Exchange of terminate_coupled between ocean and atmosphere by PE0
    78!
    89! Former revisions:
    … …  
    9394       terminate_coupled = 3
    9495#if defined( __parallel )
    95        CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,          &
    96                           target_id, 0,                                      &
    97                           terminate_coupled_remote, 1, MPI_INTEGER,          &
    98                           target_id, 0,                                      &
    99                           comm_inter, status, ierr )
     96       IF ( myid == 0 ) THEN
     97          CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,          &
     98                             target_id, 0,                                      &
     99                             terminate_coupled_remote, 1, MPI_INTEGER,          &
     100                             target_id, 0,                                      &
     101                             comm_inter, status, ierr )
     102       ENDIF
     103       CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)
    100104#endif
    101105    ENDIF
    … …  
    140144             ENDIF
    141145#if defined( __parallel )
    142              CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,    &
    143                                 target_id,  0,                               &
    144                                 terminate_coupled_remote, 1, MPI_INTEGER,    &
    145                                 target_id,  0,                               &
    146                                 comm_inter, status, ierr )
     146             IF ( myid == 0 ) THEN
     147                CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER,    &
     148                                   target_id,  0,                               &
     149                                   terminate_coupled_remote, 1, MPI_INTEGER,    &
     150                                   target_id,  0,                               &
     151                                   comm_inter, status, ierr )   
     152             ENDIF
     153             CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr) 
     154           
    147155#endif
    148156          ENDIF

  • palm/trunk/SOURCE/check_open.f90

    r601 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! Output of total array size was adapted to nbgp.
    77!
    88! Former revisions:
    … …  
    278278!--          Output for combine_plot_fields
    279279             IF ( data_output_2d_on_each_pe  .AND.  myid_char /= '' )  THEN
    280                 WRITE (21)  -1, nx+1, -1, ny+1    ! total array size
     280                WRITE (21)  -nbgp, nx+nbgp, -nbgp, ny+nbgp    ! total array size
    281281                WRITE (21)   0, nx+1,  0, ny+1    ! output part
    282282             ENDIF
    … …  
    319319!--          Output for combine_plot_fields
    320320             IF ( data_output_2d_on_each_pe  .AND.  myid_char /= '' )  THEN
    321                 WRITE (22)  -1, nx+1, 0, nz+1    ! total array size
     321                WRITE (22)  -nbgp, nx+nbgp, 0, nz+1    ! total array size
    322322                WRITE (22)   0, nx+1, 0, nz+1    ! output part
    323323             ENDIF
    … …  
    357357!--          Output for combine_plot_fields
    358358             IF ( data_output_2d_on_each_pe  .AND.  myid_char /= '' )  THEN
    359                 WRITE (23)  -1, ny+1, 0, nz+1    ! total array size
     359                WRITE (23)  -nbgp, ny+nbgp, 0, nz+1    ! total array size
    360360                WRITE (23)   0, ny+1, 0, nz+1    ! output part
    361361             ENDIF
    … …  
    392392!--          Specifications for combine_plot_fields
    393393             IF ( .NOT. do3d_compress )  THEN
    394                 WRITE ( 30 )  -1,nx+1,-1,ny+1,0,nz_do3d
     394                WRITE ( 30 )  -nbgp,nx+nbgp,-nbgp,ny+nbgp, 0 ,nz_do3d
    395395                WRITE ( 30 )  0,nx+1,0,ny+1,0,nz_do3d
    396396             ENDIF
    … …  
    503503                openfile(33)%opened = .TRUE.
    504504                WRITE ( 33, 3300 )  TRIM( avs_coor_file ), &
    505                                     TRIM( avs_coor_file ), (nx+2)*4, &
    506                                     TRIM( avs_coor_file ), (nx+2)*4+(ny+2)*4
     505                                    TRIM( avs_coor_file ), (nx+2*nbgp)*4, &
     506                                    TRIM( avs_coor_file ), (nx+2*nbgp)*4+(ny+2*nbgp)*4
    507507           
    508508
    … …  
    548548!--       corresponding partial array of a PE only once at the top of the file
    549549          IF ( avs_output  .AND.  do3d_compress )  THEN
    550              WRITE ( 30 )  nxl-1, nxr+1, nys-1, nyn+1, nzb, nz_do3d
     550             WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb, nz_do3d
    551551          ENDIF
    552552
    … …  
    929929#endif
    930930             ENDIF
    931 
     931             
    932932             CALL handle_netcdf_error( 'check_open', 26 )
    933933!
    … …  
    13251325#endif
    13261326             ENDIF
    1327 
     1327             
    13281328             CALL handle_netcdf_error( 'check_open', 43 ) 
    13291329

  • palm/trunk/SOURCE/check_parameters.f90

    r601 r667  
    44! Current revisions:
    55! -----------------
     6!
    67!
     8! Exchange of parameters between ocean and atmosphere via PE0
     9! Check for illegal combination of ws-scheme and timestep scheme.
     10! Check for topography and ws-scheme.
     11! Check for not cyclic boundary conditions in combination with ws-scheme and
     12! loop_optimization = 'vector'.
     13! Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.
     14!
     15! Different processor/grid topology in atmosphere and ocean is now allowed!
     16!
     17! Bugfixes in checking for conserve_volume_flow_mode
     18!
     19! Adapt error messages.
    720!
    821! Former revisions:
    … …  
    180193!
    181194!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
    182     IF ( coupling_mode /= 'uncoupled' )  THEN
     195    IF ( coupling_mode /= 'uncoupled')  THEN
    183196
    184197       IF ( dt_coupling == 9999999.9 )  THEN
    … …  
    189202
    190203#if defined( __parallel )
    191        CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
    192                       ierr )
    193        CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &
    194                       status, ierr )
     204       IF ( myid == 0 ) THEN
     205          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
     206                         ierr )
     207          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &
     208                         status, ierr )
     209       ENDIF
     210       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
     211       
    195212       IF ( dt_coupling /= remote )  THEN
    196213          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
    … …  
    200217       ENDIF
    201218       IF ( dt_coupling <= 0.0 )  THEN
    202           CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
    203           CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &
    204                          status, ierr )
     219          IF ( myid == 0  ) THEN
     220             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
     221             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &
     222                            status, ierr )
     223          ENDIF   
     224          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
     225         
    205226          dt_coupling = MAX( dt_max, remote )
    206227          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
    … …  
    209230          CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
    210231       ENDIF
    211 
    212        CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
    213                       ierr )
    214        CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &
    215                       status, ierr )
     232       IF ( myid == 0 ) THEN
     233          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
     234                         ierr )
     235          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &
     236                         status, ierr )
     237       ENDIF
     238       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
     239     
    216240       IF ( restart_time /= remote )  THEN
    217241          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
    … …  
    220244          CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
    221245       ENDIF
    222 
    223        CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
    224                       ierr )
    225        CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &
    226                       status, ierr )
     246       IF ( myid == 0 ) THEN
     247          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
     248                         ierr )
     249          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &
     250                         status, ierr )
     251       ENDIF   
     252       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
     253     
    227254       IF ( dt_restart /= remote )  THEN
    228255          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
    … …  
    233260
    234261       simulation_time_since_reference = end_time - coupling_start_time
    235        CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
    236                       14, comm_inter, ierr )
    237        CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
    238                       status, ierr )
     262       IF  ( myid == 0 ) THEN
     263          CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
     264                         14, comm_inter, ierr )
     265          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
     266                         status, ierr )   
     267       ENDIF
     268       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
     269     
    239270       IF ( simulation_time_since_reference /= remote )  THEN
    240271          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
    … …  
    245276       ENDIF
    246277
    247        CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
    248        CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &
    249                       status, ierr )
    250        IF ( dx /= remote )  THEN
    251           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
    252                  '":  dx = ', dx, '& is not equal to dx_remote = ', remote
    253           CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
    254        ENDIF
    255 
    256        CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
    257        CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &
    258                       status, ierr )
    259        IF ( dy /= remote )  THEN
    260           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
    261                  '":  dy = ', dy, '& is not equal to dy_remote = ', remote
    262           CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
    263        ENDIF
    264 
    265        CALL MPI_SEND( nx, 1, MPI_INTEGER, target_id, 17, comm_inter, ierr )
    266        CALL MPI_RECV( iremote, 1, MPI_INTEGER, target_id, 17, comm_inter, &
    267                       status, ierr )
    268        IF ( nx /= iremote )  THEN
    269           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
    270                  '": nx = ', nx, '& is not equal to nx_remote = ', iremote
    271           CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
    272        ENDIF
    273 
    274        CALL MPI_SEND( ny, 1, MPI_INTEGER, target_id, 18, comm_inter, ierr )
    275        CALL MPI_RECV( iremote, 1, MPI_INTEGER, target_id, 18, comm_inter, &
    276                       status, ierr )
    277        IF ( ny /= iremote )  THEN
    278           WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
    279                  '": ny = ', ny, '& is not equal to ny_remote = ', iremote
    280           CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
     278 
     279       IF ( myid == 0 ) THEN
     280          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
     281          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &
     282                                                             status, ierr )
     283       ENDIF
     284       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
     285
     286
     287       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
     288
     289          IF ( dx < remote ) THEN
     290             WRITE( message_string, * ) 'coupling mode "', &
     291                   TRIM( coupling_mode ),                  &
     292           '": dx in Atmosphere is not equal to or not larger then dx in ocean'
     293             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
     294          ENDIF
     295
     296          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
     297             WRITE( message_string, * ) 'coupling mode "', &
     298                    TRIM( coupling_mode ), &
     299             '": Domain size in x-direction is not equal in ocean and atmosphere'
     300             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
     301          ENDIF
     302
     303       ENDIF
     304
     305       IF ( myid == 0) THEN
     306          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
     307          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &
     308                         status, ierr )
     309       ENDIF
     310       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
     311
     312       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
     313
     314          IF ( dy < remote )  THEN
     315             WRITE( message_string, * ) 'coupling mode "', &
     316                    TRIM( coupling_mode ), &
     317                 '": dy in Atmosphere is not equal to or not larger then dy in ocean'
     318             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
     319          ENDIF
     320
     321          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
     322             WRITE( message_string, * ) 'coupling mode "', &
     323                   TRIM( coupling_mode ), &
     324             '": Domain size in y-direction is not equal in ocean and atmosphere'
     325             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
     326          ENDIF
     327
     328          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
     329             WRITE( message_string, * ) 'coupling mode "', &
     330                   TRIM( coupling_mode ), &
     331             '": nx+1 in ocean is not divisible without remainder with nx+1 in', &
     332             ' atmosphere'
     333             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
     334          ENDIF
     335
     336          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
     337             WRITE( message_string, * ) 'coupling mode "', &
     338                   TRIM( coupling_mode ), &
     339             '": ny+1 in ocean is not divisible without remainder with ny+1 in', &
     340             ' atmosphere'
     341             CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
     342          ENDIF
     343
    281344       ENDIF
    282345#else
    … …  
    290353!
    291354!-- Exchange via intercommunicator
    292     IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
     355    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
    293356       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, &
    294357                      ierr )
    295     ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
     358    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
    296359       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, &
    297360                      comm_inter, status, ierr )
    298361    ENDIF
     362    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
     363   
    299364#endif
    300365
    … …  
    372437          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
    373438       ENDIF
     439       IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' ) &
     440       THEN
     441          message_string = 'topography is still not allowed with ' // &
     442                           'momentum_advec = "' // TRIM( momentum_advec ) //  &
     443                           '"or scalar_advec = "' // TRIM( scalar_advec ) //'"'
     444   ! message number still needs modification
     445           CALL message( 'check_parameters', 'PA0341', 1, 2, 0, 6, 0 )
     446       END IF
     447         
    374448!
    375449!--    In case of non-flat topography, check whether the convention how to
    … …  
    492566       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
    493567    ENDIF
    494 
     568   
     569    IF( momentum_advec == 'ws-scheme' .AND. &
     570        call_psolver_at_all_substeps == .FALSE. ) THEN
     571        message_string = 'psolver must be called at each RK3 substep when "'//&
     572                      TRIM(momentum_advec) // ' "is used for momentum_advec'
     573        CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
     574    END IF
    495575!
    496576!-- Advection schemes:
    497     IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ups-scheme' ) &
    498     THEN
     577    IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' .AND. &
     578         momentum_advec /= 'ups-scheme' ) THEN
    499579       message_string = 'unknown advection scheme: momentum_advec = "' // &
    500580                        TRIM( momentum_advec ) // '"'
    501581       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
    502582    ENDIF
    503     IF ( ( momentum_advec == 'ups-scheme'  .OR.  scalar_advec == 'ups-scheme' )&
    504                                       .AND.  timestep_scheme /= 'euler' )  THEN
    505        message_string = 'momentum_advec = "' // TRIM( momentum_advec ) // &
    506                         '" is not allowed with timestep_scheme = "' //    &
    507                         TRIM( timestep_scheme ) // '"'
     583    IF ((( momentum_advec == 'ups-scheme'  .OR.  scalar_advec == 'ups-scheme' )&
     584           .AND.  timestep_scheme /= 'euler' ) .OR. (( momentum_advec == 'ws-scheme'&
     585           .OR.  scalar_advec == 'ws-scheme') .AND. (timestep_scheme == 'euler' .OR. &
     586           timestep_scheme == 'leapfrog+euler' .OR. timestep_scheme == 'leapfrog'    &
     587           .OR. timestep_scheme == 'runge-kutta-2'))) THEN
     588       message_string = 'momentum_advec or scalar_advec = "' &
     589         // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &
     590         TRIM( timestep_scheme ) // '"'
    508591       CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
    509592    ENDIF
    510 
    511     IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'bc-scheme'  .AND.&
    512          scalar_advec /= 'ups-scheme' )  THEN
     593    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
     594        scalar_advec /= 'bc-scheme'  .AND.  scalar_advec /= 'ups-scheme' )  THEN
    513595       message_string = 'unknown advection scheme: scalar_advec = "' // &
    514596                        TRIM( scalar_advec ) // '"'
    … …  
    563645    ENDIF
    564646
    565     IF ( momentum_advec /= 'pw-scheme' .AND. timestep_scheme(1:5) == 'runge' ) &
    566     THEN
     647    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme') &
     648         .AND. timestep_scheme(1:5) == 'runge' ) THEN
    567649       message_string = 'momentum advection scheme "' // &
    568650                        TRIM( momentum_advec ) // '" & does not work with ' // &
    … …  
    712794          ug_vertical_gradient_level_ind(1) = nzt+1
    713795          ug(nzt+1) = ug_surface
    714           DO  k = nzt, 0, -1
     796          DO  k = nzt, nzb, -1
    715797             IF ( i < 11 ) THEN
    716798                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND. &
    … …  
    778860          vg_vertical_gradient_level_ind(1) = nzt+1
    779861          vg(nzt+1) = vg_surface
    780           DO  k = nzt, 0, -1
     862          DO  k = nzt, nzb, -1
    781863             IF ( i < 11 ) THEN
    782864                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND. &
    … …  
    10201102 
    10211103             
     1104
    10221105!
    10231106!-- Compute Coriolis parameter
    … …  
    11591242!
    11601243!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
    1161 !-- Willimas advection scheme. Several schemes and tools do not work with
    1162 !-- non-cyclic boundary conditions.
     1244!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
     1245!-- and tools do not work with non-cyclic boundary conditions.
    11631246    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
    11641247       IF ( psolver /= 'multigrid' )  THEN
    … …  
    11671250          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
    11681251       ENDIF
    1169        IF ( momentum_advec /= 'pw-scheme' )  THEN
     1252       IF ( momentum_advec /= 'pw-scheme' .AND. &
     1253            momentum_advec /= 'ws-scheme')  THEN
    11701254          message_string = 'non-cyclic lateral boundaries do not allow ' // &
    11711255                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
    11721256          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
    11731257       ENDIF
    1174        IF ( scalar_advec /= 'pw-scheme' )  THEN
     1258       IF ( scalar_advec /= 'pw-scheme' .AND. &
     1259            scalar_advec /= 'ws-scheme' )  THEN
    11751260          message_string = 'non-cyclic lateral boundaries do not allow ' // &
    11761261                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
    11771262          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
    11781263       ENDIF
     1264       IF ( (scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' ) &
     1265          .AND. loop_optimization == 'vector' ) THEN
     1266          message_string = 'non-cyclic lateral boundaries do not allow ' // &
     1267                           'loop_optimization = vector and ' //  &
     1268                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
     1269  ! The error message number still needs modification.
     1270          CALL message( 'check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
     1271       END IF
    11791272       IF ( galilei_transformation )  THEN
    11801273          message_string = 'non-cyclic lateral boundaries do not allow ' // &
    … …  
    14071500                        TRIM( bc_uv_b ) // '"'
    14081501       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
     1502    ENDIF
     1503!
     1504!-- In case of coupled simulations u and v at the ground in atmosphere will be
     1505!-- assigned with the u and v values of the ocean surface
     1506    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
     1507       ibc_uv_b = 2
    14091508    ENDIF
    14101509
    … …  
    21092208             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
    21102209
     2210
    21112211          CASE ( 'u"pt"' )
    21122212             dopr_index(i) = 58
    … …  
    22442344
    22452345       END SELECT
     2346
    22462347!
    22472348!--    Check to which of the predefined coordinate systems the profile belongs
    … …  
    25842685!-- Upper plot limit (grid point value) for 1D profiles
    25852686    IF ( z_max_do1d == -1.0 )  THEN
     2687
    25862688       nz_do1d = nzt+1
     2689
    25872690    ELSE
    25882691       DO  k = nzb+1, nzt+1
    … …  
    27372840
    27382841!
     2842
    27392843!-- Check netcdf precison
    27402844    ldum = .FALSE.
    … …  
    30703174    IF ( conserve_volume_flow )  THEN
    30713175       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
    3072           IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
    3073              conserve_volume_flow_mode = 'inflow_profile'
    3074           ELSE
    3075              conserve_volume_flow_mode = 'initial_profiles'
    3076           ENDIF
     3176
     3177          conserve_volume_flow_mode = 'initial_profiles'
     3178
    30773179       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
    30783180            TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' .AND.  &
    … …  
    30823184          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
    30833185       ENDIF
    3084        IF ( ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' ) .AND. &
    3085             TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' )  THEN
    3086           WRITE( message_string, * )  'noncyclic boundary conditions ', &
    3087                'require & conserve_volume_flow_mode = ''inflow_profile'''
     3186       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic') .AND. &
     3187          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
     3188          WRITE( message_string, * )  'non-cyclic boundary conditions ', &
     3189               'require  conserve_volume_flow_mode = ''initial_profiles'''
    30883190          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
    30893191       ENDIF
    … …  
    30913193            TRIM( conserve_volume_flow_mode ) == 'inflow_profile' )  THEN
    30923194          WRITE( message_string, * )  'cyclic boundary conditions ', &
    3093                'require & conserve_volume_flow_mode = ''initial_profiles''', &
     3195               'require conserve_volume_flow_mode = ''initial_profiles''', &
    30943196               ' or ''bulk_velocity'''
    30953197          CALL message( 'check_parameters', 'PA0156', 1, 2, 0, 6, 0 )
    … …  
    31003202         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
    31013203       WRITE( message_string, * )  'nonzero bulk velocity requires ', &
    3102             'conserve_volume_flow = .T. and & ', &
     3204            'conserve_volume_flow = .T. and ', &
    31033205            'conserve_volume_flow_mode = ''bulk_velocity'''
    31043206       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
    … …  
    31393241
    31403242
     3243
    31413244 END SUBROUTINE check_parameters

  • palm/trunk/SOURCE/data_output_2d.f90

    r623 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
     7! allocation of arrays local_2d and total_2d.
     8! Calls of exchange_horiz are modiefied.
    79!
    810! Former revisions:
    … …  
    112114
    113115       CASE ( 'xy' )
    114 
    115116          s = 1
    116           ALLOCATE( level_z(0:nzt+1), local_2d(nxl-1:nxr+1,nys-1:nyn+1) )
     117          ALLOCATE( level_z(nzb:nzt+1), local_2d(nxlg:nxrg,nysg:nyng) )
    117118
    118119!
    … …  
    130131                IF ( iso2d_output )  CALL check_open( 21 )
    131132#if defined( __parallel )
    132                 ALLOCATE( total_2d(-1:nx+1,-1:ny+1) )
     133                ALLOCATE( total_2d(-nbgp:nx+nbgp,-nbgp:ny+nbgp) )
    133134#endif
    134135             ENDIF
    … …  
    136137
    137138       CASE ( 'xz' )
    138 
    139139          s = 2
    140           ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) )
     140          ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
    141141
    142142!
    … …  
    154154                IF ( iso2d_output )  CALL check_open( 22 )
    155155#if defined( __parallel )
    156                 ALLOCATE( total_2d(-1:nx+1,nzb:nzt+1) )
     156                ALLOCATE( total_2d(-nbgp:nx+nbgp,nzb:nzt+1) )
    157157#endif
    158158             ENDIF
    … …  
    162162
    163163          s = 3
    164           ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) )
     164          ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
    165165
    166166!
    … …  
    178178                IF ( iso2d_output )  CALL check_open( 23 )
    179179#if defined( __parallel )
    180                 ALLOCATE( total_2d(-1:ny+1,nzb:nzt+1) )
     180                ALLOCATE( total_2d(-nbgp:ny+nbgp,nzb:nzt+1) )
    181181#endif
    182182             ENDIF
    … …  
    192192!
    193193!-- Allocate a temporary array for resorting (kji -> ijk).
    194     ALLOCATE( local_pf(nxl-1:nxr+1,nys-1:nyn+1,nzb:nzt+1) )
     194    ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nzt+1) )
    195195
    196196!
    … …  
    219219             CASE ( 'lwp*_xy' )        ! 2d-array
    220220                IF ( av == 0 )  THEN
    221                    DO  i = nxl-1, nxr+1
    222                       DO  j = nys-1, nyn+1
     221                   DO  i = nxlg, nxrg
     222                      DO  j = nysg, nyng
    223223                         local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) * &
    224224                                                    dzw(1:nzt+1) )
    … …  
    226226                   ENDDO
    227227                ELSE
    228                    DO  i = nxl-1, nxr+1
    229                       DO  j = nys-1, nyn+1
     228                   DO  i = nxlg, nxrg
     229                      DO  j = nysg, nyng
    230230                         local_pf(i,j,nzb+1) = lwp_av(j,i)
    231231                      ENDDO
    … …  
    248248                   IF ( simulated_time >= particle_advection_start )  THEN
    249249                      tend = prt_count
    250                       CALL exchange_horiz( tend )
     250                      CALL exchange_horiz( tend, nbgp )
    251251                   ELSE
    252252                      tend = 0.0
    253253                   ENDIF
    254                    DO  i = nxl-1, nxr+1
    255                       DO  j = nys-1, nyn+1
     254                   DO  i = nxlg, nxrg
     255                      DO  j = nysg, nyng
    256256                         DO  k = nzb, nzt+1
    257257                            local_pf(i,j,k) = tend(k,j,i)
    … …  
    261261                   resorted = .TRUE.
    262262                ELSE
    263                    CALL exchange_horiz( pc_av )
     263                   CALL exchange_horiz( pc_av, nbgp )
    264264                   to_be_resorted => pc_av
    265265                ENDIF
    … …  
    287287                         ENDDO
    288288                      ENDDO
    289                       CALL exchange_horiz( tend )
     289                      CALL exchange_horiz( tend, nbgp )
    290290                   ELSE
    291291                      tend = 0.0
    292                    ENDIF
    293                    DO  i = nxl-1, nxr+1
    294                       DO  j = nys-1, nyn+1
     292                   END IF
     293                   DO  i = nxlg, nxrg
     294                      DO  j = nysg, nyng
    295295                         DO  k = nzb, nzt+1
    296296                            local_pf(i,j,k) = tend(k,j,i)
    … …  
    300300                   resorted = .TRUE.
    301301                ELSE
    302                    CALL exchange_horiz( pr_av )
     302                   CALL exchange_horiz( pr_av, nbgp )
    303303                   to_be_resorted => pr_av
    304304                ENDIF
    … …  
    306306             CASE ( 'pra*_xy' )        ! 2d-array / integral quantity => no av
    307307                CALL exchange_horiz_2d( precipitation_amount )
    308                 DO  i = nxl-1, nxr+1
    309                    DO  j = nys-1, nyn+1
     308                   DO  i = nxlg, nxrg
     309                      DO  j = nysg, nyng
    310310                      local_pf(i,j,nzb+1) =  precipitation_amount(j,i)
    311311                   ENDDO
    … …  
    319319                IF ( av == 0 )  THEN
    320320                   CALL exchange_horiz_2d( precipitation_rate )
    321                    DO  i = nxl-1, nxr+1
    322                       DO  j = nys-1, nyn+1
     321                   DO  i = nxlg, nxrg
     322                      DO  j = nysg, nyng
    323323                         local_pf(i,j,nzb+1) =  precipitation_rate(j,i)
    324324                      ENDDO
    … …  
    326326                ELSE
    327327                   CALL exchange_horiz_2d( precipitation_rate_av )
    328                    DO  i = nxl-1, nxr+1
    329                       DO  j = nys-1, nyn+1
     328                   DO  i = nxlg, nxrg
     329                      DO  j = nysg, nyng
    330330                         local_pf(i,j,nzb+1) =  precipitation_rate_av(j,i)
    331331                      ENDDO
    … …  
    341341                      to_be_resorted => pt
    342342                   ELSE
    343                       DO  i = nxl-1, nxr+1
    344                          DO  j = nys-1, nyn+1
     343                   DO  i = nxlg, nxrg
     344                      DO  j = nysg, nyng
    345345                            DO  k = nzb, nzt+1
    346346                               local_pf(i,j,k) = pt(k,j,i) + l_d_cp *    &
    … …  
    399399             CASE ( 'qsws*_xy' )        ! 2d-array
    400400                IF ( av == 0 ) THEN
    401                    DO  i = nxl-1, nxr+1
    402                       DO  j = nys-1, nyn+1
     401                   DO  i = nxlg, nxrg
     402                      DO  j = nysg, nyng
    403403                         local_pf(i,j,nzb+1) =  qsws(j,i)
    404404                      ENDDO
    405405                   ENDDO
    406406                ELSE
    407                    DO  i = nxl-1, nxr+1
    408                       DO  j = nys-1, nyn+1
     407                   DO  i = nxlg, nxrg
     408                      DO  j = nysg, nyng
    409409                         local_pf(i,j,nzb+1) =  qsws_av(j,i)
    410410                      ENDDO
    … …  
    417417             CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
    418418                IF ( av == 0 )  THEN
    419                    DO  i = nxl-1, nxr+1
    420                       DO  j = nys-1, nyn+1
     419                   DO  i = nxlg, nxrg
     420                      DO  j = nysg, nyng
    421421                         DO  k = nzb, nzt+1
    422422                            local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
    … …  
    453453             CASE ( 'shf*_xy' )        ! 2d-array
    454454                IF ( av == 0 ) THEN
    455                    DO  i = nxl-1, nxr+1
    456                       DO  j = nys-1, nyn+1
     455                   DO  i = nxlg, nxrg
     456                      DO  j = nysg, nyng
    457457                         local_pf(i,j,nzb+1) =  shf(j,i)
    458458                      ENDDO
    459459                   ENDDO
    460460                ELSE
    461                    DO  i = nxl-1, nxr+1
    462                       DO  j = nys-1, nyn+1
     461                   DO  i = nxlg, nxrg
     462                      DO  j = nysg, nyng
    463463                         local_pf(i,j,nzb+1) =  shf_av(j,i)
    464464                      ENDDO
    … …  
    471471             CASE ( 't*_xy' )        ! 2d-array
    472472                IF ( av == 0 )  THEN
    473                    DO  i = nxl-1, nxr+1
    474                       DO  j = nys-1, nyn+1
     473                   DO  i = nxlg, nxrg
     474                      DO  j = nysg, nyng
    475475                         local_pf(i,j,nzb+1) = ts(j,i)
    476476                      ENDDO
    477477                   ENDDO
    478478                ELSE
    479                    DO  i = nxl-1, nxr+1
    480                       DO  j = nys-1, nyn+1
     479                   DO  i = nxlg, nxrg
     480                      DO  j = nysg, nyng
    481481                         local_pf(i,j,nzb+1) = ts_av(j,i)
    482482                      ENDDO
    … …  
    503503             CASE ( 'u*_xy' )        ! 2d-array
    504504                IF ( av == 0 )  THEN
    505                    DO  i = nxl-1, nxr+1
    506                       DO  j = nys-1, nyn+1
     505                   DO  i = nxlg, nxrg
     506                      DO  j = nysg, nyng
    507507                         local_pf(i,j,nzb+1) = us(j,i)
    508508                      ENDDO
    509509                   ENDDO
    510510                ELSE
    511                    DO  i = nxl-1, nxr+1
    512                       DO  j = nys-1, nyn+1
     511                   DO  i = nxlg, nxrg
     512                      DO  j = nysg, nyng
    513513                         local_pf(i,j,nzb+1) = us_av(j,i)
    514514                      ENDDO
    … …  
    551551             CASE ( 'z0*_xy' )        ! 2d-array
    552552                IF ( av == 0 ) THEN
    553                    DO  i = nxl-1, nxr+1
    554                       DO  j = nys-1, nyn+1
     553                   DO  i = nxlg, nxrg
     554                      DO  j = nysg, nyng
    555555                         local_pf(i,j,nzb+1) =  z0(j,i)
    556556                      ENDDO
    557557                   ENDDO
    558558                ELSE
    559                    DO  i = nxl-1, nxr+1
    560                       DO  j = nys-1, nyn+1
     559                   DO  i = nxlg, nxrg
     560                      DO  j = nysg, nyng
    561561                         local_pf(i,j,nzb+1) =  z0_av(j,i)
    562562                      ENDDO
    … …  
    593593!--       Resort the array to be output, if not done above
    594594          IF ( .NOT. resorted )  THEN
    595              DO  i = nxl-1, nxr+1
    596                 DO  j = nys-1, nyn+1
     595             DO  i = nxlg, nxrg
     596                DO  j = nysg, nyng
    597597                   DO  k = nzb, nzt+1
    598598                      local_pf(i,j,k) = to_be_resorted(k,j,i)
    … …  
    647647!--                   Carry out the averaging (all data are on the PE)
    648648                      DO  k = nzb, nzt+1
    649                          DO  j = nys-1, nyn+1
    650                             DO  i = nxl-1, nxr+1
     649                         DO  j = nysg, nyng
     650                            DO  i = nxlg, nxrg
    651651                               local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
    652652                            ENDDO
    … …  
    654654                      ENDDO
    655655
    656                       local_2d = local_2d / ( nzt -nzb + 2.0 )
     656                      local_2d = local_2d / ( nzt -nzb + 2.0)
    657657
    658658                   ELSE
    … …  
    723723                         ENDIF
    724724#endif
    725                          WRITE ( 21 )  nxl-1, nxr+1, nys-1, nyn+1
     725                         WRITE ( 21 )  nxlg, nxrg, nysg, nyng
    726726                         WRITE ( 21 )  local_2d
    727727
    … …  
    734734                         CALL MPI_BARRIER( comm2d, ierr )
    735735
    736                          ngp = ( nxr-nxl+3 ) * ( nyn-nys+3 )
     736                         ngp = ( nxrg-nxlg+1 ) * ( nyng-nysg+1 )
    737737                         IF ( myid == 0 )  THEN
    738738!
    739739!--                         Local array can be relocated directly.
    740                             total_2d(nxl-1:nxr+1,nys-1:nyn+1) = local_2d
     740                            total_2d(nxlg:nxrg,nysg:nyng) = local_2d
    741741!
    742742!--                         Receive data from all other PEs.
    … …  
    760760!--                         Output of the total cross-section.
    761761                            IF ( iso2d_output )  THEN
    762                                WRITE (21)  total_2d(0:nx+1,0:ny+1)
     762                               WRITE (21)  total_2d(-nbgp:nx+nbgp,-nbgp:ny+nbgp)
    763763                            ENDIF
    764764!
    765765!--                         Relocate the local array for the next loop increment
    766766                            DEALLOCATE( local_2d )
    767                             ALLOCATE( local_2d(nxl-1:nxr+1,nys-1:nyn+1) )
     767                            ALLOCATE( local_2d(nxlg:nxrg,nysg:nyng) )
    768768
    769769#if defined( __netcdf )
    … …  
    789789!
    790790!--                         First send the local index limits to PE0
    791                             ind(1) = nxl-1; ind(2) = nxr+1
    792                             ind(3) = nys-1; ind(4) = nyn+1
     791                            ind(1) = nxlg; ind(2) = nxrg
     792                            ind(3) = nysg; ind(4) = nyng
    793793                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
    794794                                           comm2d, ierr )
    795795!
    796796!--                         Send data to PE0
    797                             CALL MPI_SEND( local_2d(nxl-1,nys-1), ngp, &
     797                            CALL MPI_SEND( local_2d(nxlg,nysg), ngp, &
    798798                                           MPI_REAL, 0, 1, comm2d, ierr )
    799799                         ENDIF
    … …  
    882882
    883883                   ENDIF
     884
    884885!
    885886!--                If required, carry out averaging along y
    886887                   IF ( section(is,s) == -1 )  THEN
    887888
    888                       ALLOCATE( local_2d_l(nxl-1:nxr+1,nzb:nzt+1) )
     889                      ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
    889890                      local_2d_l = 0.0
    890                       ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
     891                      ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
    891892!
    892893!--                   First local averaging on the PE
    893894                      DO  k = nzb, nzt+1
    894895                         DO  j = nys, nyn
    895                             DO  i = nxl-1, nxr+1
     896                            DO  i = nxlg, nxrg
    896897                               local_2d_l(i,k) = local_2d_l(i,k) + &
    897898                                                 local_pf(i,j,k)
    … …  
    903904!--                   Now do the averaging over all PEs along y
    904905                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    905                       CALL MPI_ALLREDUCE( local_2d_l(nxl-1,nzb),              &
    906                                           local_2d(nxl-1,nzb), ngp, MPI_REAL, &
     906                      CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb),              &
     907                                          local_2d(nxlg,nzb), ngp, MPI_REAL, &
    907908                                          MPI_SUM, comm1dy, ierr )
    908909#else
    … …  
    936937!--                   BEGIN WORKAROUND---------------------------------------
    937938                      IF ( npey /= 1  .AND.  section(is,s) /= -1)  THEN
    938                          ALLOCATE( local_2d_l(nxl-1:nxr+1,nzb:nzt+1) )
     939                         ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
    939940                         local_2d_l = 0.0
    940941                         IF ( section(is,s) >= nys .AND. section(is,s) <= nyn )&
    … …  
    945946!
    946947!--                      Distribute data over all PEs along y
    947                          ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
     948                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
    948949                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
    949                          CALL MPI_ALLREDUCE( local_2d_l(nxl-1,nzb),            &
    950                                              local_2d(nxl-1,nzb), ngp,         &
     950                         CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb),            &
     951                                             local_2d(nxlg,nzb), ngp,         &
    951952                                             MPI_REAL, MPI_SUM, comm1dy, ierr )
    952953#else
    … …  
    10221023                              ( section(is,s) == -1  .AND.  nys-1 == -1 ) )  &
    10231024                         THEN
    1024                             WRITE (22)  nxl-1, nxr+1, nzb, nzt+1
     1025                            WRITE (22)  nxlg, nxrg, nzb, nzt+1
    10251026                            WRITE (22)  local_2d
    10261027                         ELSE
    … …  
    10361037                         CALL MPI_BARRIER( comm2d, ierr )
    10371038
    1038                          ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
     1039                         ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
    10391040                         IF ( myid == 0 )  THEN
    10401041!
    … …  
    10441045                                 ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
    10451046                            THEN
    1046                                total_2d(nxl-1:nxr+1,nzb:nzt+1) = local_2d
     1047                               total_2d(nxlg:nxrg,nzb:nzt+1) = local_2d
    10471048                            ENDIF
    10481049!
    … …  
    10731074!--                         Output of the total cross-section.
    10741075                            IF ( iso2d_output )  THEN
    1075                                WRITE (22)  total_2d(0:nx+1,nzb:nzt+1)
     1076                               WRITE (22)  total_2d(-nbgp:nx+nbgp,nzb:nzt+1)
    10761077                            ENDIF
    10771078!
    10781079!--                         Relocate the local array for the next loop increment
    10791080                            DEALLOCATE( local_2d )
    1080                             ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) )
     1081                            ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
    10811082
    10821083#if defined( __netcdf )
    … …  
    10991100                                 ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
    11001101                            THEN
    1101                                ind(1) = nxl-1; ind(2) = nxr+1
     1102                               ind(1) = nxlg; ind(2) = nxrg
    11021103                               ind(3) = nzb;   ind(4) = nzt+1
    11031104                            ELSE
    … …  
    11101111!--                         If applicable, send data to PE0.
    11111112                            IF ( ind(1) /= -9999 )  THEN
    1112                                CALL MPI_SEND( local_2d(nxl-1,nzb), ngp, &
     1113                               CALL MPI_SEND( local_2d(nxlg,nzb), ngp, &
    11131114                                              MPI_REAL, 0, 1, comm2d, ierr )
    11141115                            ENDIF
    … …  
    11871188                   IF ( section(is,s) == -1 )  THEN
    11881189
    1189                       ALLOCATE( local_2d_l(nys-1:nyn+1,nzb:nzt+1) )
     1190                      ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
    11901191                      local_2d_l = 0.0
    1191                       ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
     1192                      ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
    11921193!
    11931194!--                   First local averaging on the PE
    11941195                      DO  k = nzb, nzt+1
    1195                          DO  j = nys-1, nyn+1
     1196                         DO  j = nysg, nyng
    11961197                            DO  i = nxl, nxr
    11971198                               local_2d_l(j,k) = local_2d_l(j,k) + &
    … …  
    12041205!--                   Now do the averaging over all PEs along x
    12051206                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1206                       CALL MPI_ALLREDUCE( local_2d_l(nys-1,nzb),              &
    1207                                           local_2d(nys-1,nzb), ngp, MPI_REAL, &
     1207                      CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb),              &
     1208                                          local_2d(nysg,nzb), ngp, MPI_REAL, &
    12081209                                          MPI_SUM, comm1dx, ierr )
    12091210#else
    … …  
    12371238!--                   BEGIN WORKAROUND---------------------------------------
    12381239                      IF ( npex /= 1  .AND.  section(is,s) /= -1)  THEN
    1239                          ALLOCATE( local_2d_l(nys-1:nyn+1,nzb:nzt+1) )
     1240                         ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
    12401241                         local_2d_l = 0.0
    12411242                         IF ( section(is,s) >= nxl .AND. section(is,s) <= nxr )&
    … …  
    12461247!
    12471248!--                      Distribute data over all PEs along x
    1248                          ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
     1249                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 )
    12491250                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
    1250                          CALL MPI_ALLREDUCE( local_2d_l(nys-1,nzb),            &
    1251                                              local_2d(nys-1,nzb), ngp,         &
     1251                         CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb),            &
     1252                                             local_2d(nysg,nzb), ngp,         &
    12521253                                             MPI_REAL, MPI_SUM, comm1dx, ierr )
    12531254#else
    … …  
    13231324                              ( section(is,s) == -1  .AND.  nxl-1 == -1 ) )  &
    13241325                         THEN
    1325                             WRITE (23)  nys-1, nyn+1, nzb, nzt+1
     1326                            WRITE (23)  nysg, nyng, nzb, nzt+1
    13261327                            WRITE (23)  local_2d
    13271328                         ELSE
    … …  
    13371338                         CALL MPI_BARRIER( comm2d, ierr )
    13381339
    1339                          ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
     1340                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
    13401341                         IF ( myid == 0 )  THEN
    13411342!
    … …  
    13451346                                 ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
    13461347                            THEN
    1347                                total_2d(nys-1:nyn+1,nzb:nzt+1) = local_2d
     1348                               total_2d(nysg:nyng,nzb:nzt+1) = local_2d
    13481349                            ENDIF
    13491350!
    … …  
    13791380!--                         Relocate the local array for the next loop increment
    13801381                            DEALLOCATE( local_2d )
    1381                             ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) )
     1382                            ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
    13821383
    13831384#if defined( __netcdf )
    … …  
    14001401                                 ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
    14011402                            THEN
    1402                                ind(1) = nys-1; ind(2) = nyn+1
     1403                               ind(1) = nysg; ind(2) = nyng
    14031404                               ind(3) = nzb;   ind(4) = nzt+1
    14041405                            ELSE
    … …  
    14111412!--                         If applicable, send data to PE0.
    14121413                            IF ( ind(1) /= -9999 )  THEN
    1413                                CALL MPI_SEND( local_2d(nys-1,nzb), ngp, &
     1414                               CALL MPI_SEND( local_2d(nysg,nzb), ngp, &
    14141415                                              MPI_REAL, 0, 1, comm2d, ierr )
    14151416                            ENDIF

  • palm/trunk/SOURCE/data_output_3d.f90

    r647 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
     7! allocation of arrays.  Calls of exchange_horiz are modified.
     8! Skip-value skip_do_avs changed to a dynamic adaption of ghost points.
    79!
    810! Former revisions:
    … …  
    102104!
    103105!-- Allocate a temporary array with the desired output dimensions.
    104     ALLOCATE( local_pf(nxl-1:nxr+1,nys-1:nyn+1,nzb:nz_do3d) )
     106    ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) )
    105107
    106108!
    … …  
    157159             IF ( av == 0 )  THEN
    158160                tend = prt_count
    159                 CALL exchange_horiz( tend )
    160                 DO  i = nxl-1, nxr+1
    161                    DO  j = nys-1, nyn+1
     161                CALL exchange_horiz( tend, nbgp )
     162                DO  i = nxlg, nxrg
     163                   DO  j = nysg, nyng
    162164                      DO  k = nzb, nz_do3d
    163165                         local_pf(i,j,k) = tend(k,j,i)
    … …  
    167169                resorted = .TRUE.
    168170             ELSE
    169                 CALL exchange_horiz( pc_av )
     171                CALL exchange_horiz( pc_av, nbgp )
    170172                to_be_resorted => pc_av
    171173             ENDIF
    … …  
    192194                   ENDDO
    193195                ENDDO
    194                 CALL exchange_horiz( tend )
    195                 DO  i = nxl-1, nxr+1
    196                    DO  j = nys-1, nyn+1
     196                CALL exchange_horiz( tend, nbgp )
     197                DO  i = nxlg, nxrg
     198                   DO  j = nysg, nyng
    197199                      DO  k = nzb, nzt+1
    198200                         local_pf(i,j,k) = tend(k,j,i)
    … …  
    202204                resorted = .TRUE.
    203205             ELSE
    204                 CALL exchange_horiz( pr_av )
     206                CALL exchange_horiz( pr_av, nbgp )
    205207                to_be_resorted => pr_av
    206208             ENDIF
    … …  
    211213                   to_be_resorted => pt
    212214                ELSE
    213                    DO  i = nxl-1, nxr+1
    214                       DO  j = nys-1, nyn+1
     215                   DO  i = nxlg, nxrg
     216                      DO  j = nysg, nyng
    215217                         DO  k = nzb, nz_do3d
    216218                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp *    &
    … …  
    263265          CASE ( 'qv' )
    264266             IF ( av == 0 )  THEN
    265                 DO  i = nxl-1, nxr+1
    266                    DO  j = nys-1, nyn+1
     267                DO  i = nxlg, nxrg
     268                   DO  j = nysg, nyng
    267269                      DO  k = nzb, nz_do3d
    268270                         local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
    … …  
    342344!--    Resort the array to be output, if not done above
    343345       IF ( .NOT. resorted )  THEN
    344           DO  i = nxl-1, nxr+1
    345              DO  j = nys-1, nyn+1
     346          DO  i = nxlg, nxrg
     347             DO  j = nysg, nyng
    346348                DO  k = nzb, nz_do3d
    347349                   local_pf(i,j,k) = to_be_resorted(k,j,i)
    … …  
    376378!--       Determine the Skip-value for the next array. Record end and start
    377379!--       require 4 byte each.
    378           skip_do_avs = skip_do_avs + ( ((nx+2)*(ny+2)*(nz_do3d+1)) * 4 + 8 )
     380          skip_do_avs = skip_do_avs + ( ((nx+2*nbgp)*(ny+2*nbgp)*(nz_do3d+1)) * 4 + 8 )
    379381       ENDIF
    380382
    … …  
    386388!--       of compressed data.
    387389          CALL write_compressed( local_pf, 30, 33, myid, nxl, nxr, nyn, nys, &
    388                                  nzb, nz_do3d, prec )
     390                                 nzb, nz_do3d, prec, nbgp )
    389391       ELSE
    390392!
    … …  
    400402                   WRITE ( 30 )  simulated_time, do3d_time_count(av), av
    401403                ENDIF
    402                 WRITE ( 30 )  nxl-1, nxr+1, nys-1, nyn+1, nzb, nz_do3d
     404                WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb, nz_do3d
    403405                WRITE ( 30 )  local_pf
    404406
    … …  
    411413                IF ( nxr == nx  .AND.  nyn /= ny )  THEN
    412414                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
    413                                   local_pf(nxl:nxr+1,nys:nyn,nzb:nz_do3d),    &
     415                                  local_pf(nxl:nxrg,nys:nyn,nzb:nz_do3d),    &
    414416                      start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
    415                       count = (/ nxr-nxl+2, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
     417                      count = (/ nxr-nxl+1+nbgp, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
    416418                ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
    417419                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
    418                                   local_pf(nxl:nxr,nys:nyn+1,nzb:nz_do3d),    &
     420                                  local_pf(nxl:nxr,nys:nyng,nzb:nz_do3d),    &
    419421                      start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
    420                       count = (/ nxr-nxl+1, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
     422                      count = (/ nxr-nxl+1, nyn-nys+1+nbgp, nz_do3d-nzb+1, 1 /) )
    421423                ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
    422424                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
    423                                   local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),  &
     425                                  local_pf(nxl:nxrg,nys:nyng,nzb:nz_do3d),  &
    424426                      start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
    425                       count = (/ nxr-nxl+2, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
     427                      count = (/ nxr-nxl+1+nbgp, nyn-nys+1+nbgp, nz_do3d-nzb+1, 1 /) )
    426428                ELSE
    427429                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &

  • palm/trunk/SOURCE/data_output_mask.f90

    r565 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! Calls of exchange_horiz are modified.
    77!
    88! Former revisions:
    … …  
    123123             IF ( av == 0 )  THEN
    124124                tend = prt_count
    125                 CALL exchange_horiz( tend )
     125                CALL exchange_horiz( tend, nbgp )
    126126                DO  i = 1, mask_size_l(mid,1)
    127127                   DO  j = 1, mask_size_l(mid,2)
    … …  
    134134                resorted = .TRUE.
    135135             ELSE
    136                 CALL exchange_horiz( pc_av )
     136                CALL exchange_horiz( pc_av, nbgp )
    137137                to_be_resorted => pc_av
    138138             ENDIF
    … …  
    159159                   ENDDO
    160160                ENDDO
    161                 CALL exchange_horiz( tend )
     161                CALL exchange_horiz( tend, nbgp )
    162162                DO  i = 1, mask_size_l(mid,1)
    163163                   DO  j = 1, mask_size_l(mid,2)
    … …  
    170170                resorted = .TRUE.
    171171             ELSE
    172                 CALL exchange_horiz( pr_av )
     172                CALL exchange_horiz( pr_av, nbgp )
    173173                to_be_resorted => pr_av
    174174             ENDIF
    … …  
    439439
    440440       if = if + 1
     441
    441442    ENDDO
    442443

  • palm/trunk/SOURCE/diffusion_e.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    77!
    88! Former revisions:
    … …  
    6868       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
    6969       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
    70                            l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1)
    71        REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
     70                           l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
     71       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
    7272       REAL, DIMENSION(:,:), POINTER   ::  rif
    7373       REAL, DIMENSION(:,:,:), POINTER ::  e, km, var
    … …  
    289289       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
    290290       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
    291                            l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1)
    292        REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
     291                           l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
     292       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
    293293       REAL, DIMENSION(:,:), POINTER   ::  rif
    294294       REAL, DIMENSION(:,:,:), POINTER ::  e, km, var

  • palm/trunk/SOURCE/diffusion_s.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    77!
    88! Former revisions:
    … …  
    6464       REAL    ::  vertical_gridspace
    6565       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    66        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     66       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    6767       REAL    ::  wall_s_flux(0:4)
    6868       REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t
    … …  
    176176       REAL    ::  vertical_gridspace
    177177       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1)
    178        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     178       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    179179       REAL    ::  wall_s_flux(0:4)
    180180       REAL, DIMENSION(:,:),   POINTER ::  s_flux_b, s_flux_t

  • palm/trunk/SOURCE/diffusion_u.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    77!
    88! Former revisions:
    … …  
    7272       INTEGER ::  i, j, k
    7373       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
    74        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
    75        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     74       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)
     75       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    7676       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
    7777       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
    … …  
    177177                      &   - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k) &
    178178                      &            + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx     &
    179                       &            )                                         &
     179                      &            )                                          &
    180180                      &   ) * ddzw(k)
    181181             ENDDO
    … …  
    248248       INTEGER ::  i, j, k
    249249       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
    250        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
    251        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     250       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)
     251       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    252252       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
    253253       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst

  • palm/trunk/SOURCE/diffusion_v.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    77!
    88! Former revisions:
    … …  
    7070       INTEGER ::  i, j, k
    7171       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
    72        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
    73        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     72       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)
     73       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    7474       REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
    7575       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
    … …  
    246246       INTEGER ::  i, j, k
    247247       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
    248        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
    249        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     248       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)
     249       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    250250       REAL, DIMENSION(nzb:nzt+1)      ::  vsus
    251251       REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst

  • palm/trunk/SOURCE/diffusion_w.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    77!
    88! Former revisions:
    … …  
    6464       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
    6565                   kmyp_z
    66        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
    67                    km_damp_y(nys-1:nyn+1)
    68        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     66       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        &
     67                   km_damp_y(nysg:nyng)
     68       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    6969       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
    7070       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus, wsvs
    … …  
    211211       REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
    212212                   kmyp_z
    213        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
    214                    km_damp_y(nys-1:nyn+1)
    215        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     213       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg),        &
     214                   km_damp_y(nysg:nyng)
     215       REAL    ::  tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    216216       REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
    217217       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w

  • palm/trunk/SOURCE/diffusivities.f90

    r510 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
    77!
    88! Former revisions:
    … …  
    5454    REAL, SAVE ::  phi_m = 1.0
    5555
    56     REAL    ::  var(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     56    REAL    ::  var(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    5757
    5858    REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e
    … …  
    7373
    7474    !$OMP DO
    75     DO  i = nxl-1, nxr+1
    76        DO  j = nys-1, nyn+1
     75    DO  i = nxlg, nxrg
     76       DO  j = nysg, nyng
    7777
    7878!
    … …  
    157157
    158158    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
    159                                                 ! data output
     159                                                  ! data output
    160160
    161161    !$OMP END PARALLEL
    … …  
    167167!-- values of the diffusivities are not needed
    168168    !$OMP PARALLEL DO
    169     DO  i = nxl-1, nxr+1
    170        DO  j = nys-1, nyn+1
     169    DO  i = nxlg, nxrg
     170       DO  j = nysg, nyng
    171171          km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i)
    172172          km(nzt+1,j,i)            = km(nzt,j,i)

  • palm/trunk/SOURCE/disturb_field.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
     7! Calls of exchange_horiz are modified.
    78!
    89! Former revisions:
    … …  
    4445
    4546    INTEGER ::  i, j, k
    46     INTEGER ::  nzb_uv_inner(nys-1:nyn+1,nxl-1:nxr+1)
     47    INTEGER ::  nzb_uv_inner(nysg:nyng,nxlg:nxrg)
    4748
    4849    REAL    ::  randomnumber,                             &
    49                 dist1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    50                 field(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     50                dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
     51                field(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    5152    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  dist2
    5253
    … …  
    5758!-- Create an additional temporary array and initialize the arrays needed
    5859!-- to store the disturbance
    59     ALLOCATE( dist2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     60    ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    6061    dist1 = 0.0
    6162    dist2 = 0.0
    … …  
    102103!-- Exchange of ghost points for the random perturbation
    103104
    104     CALL exchange_horiz( dist1 )
    105 
     105    CALL exchange_horiz( dist1, nbgp )
    106106!
    107107!-- Applying the Shuman filter in order to smooth the perturbations.
    … …  
    128128!-- Exchange of ghost points for the filtered perturbation.
    129129!-- Afterwards, filter operation and exchange of ghost points are repeated.
    130     CALL exchange_horiz( dist2 )
     130    CALL exchange_horiz( dist2, nbgp )
    131131
    132132    DO  i = nxl, nxr
    … …  
    141141    ENDDO
    142142
    143     CALL exchange_horiz( dist1 )
     143    CALL exchange_horiz( dist1, nbgp )
    144144
    145145!
    … …  
    148148!-- (diffusion criterion))
    149149    IF ( TRIM( topography ) /= 'flat' )  THEN
    150        DO  i = nxl-1, nxr+1
    151           DO  j = nys-1, nyn+1
     150       DO  i = nxlg, nxrg
     151          DO  j = nysg, nyng
    152152             dist1(nzb:nzb_uv_inner(j,i)+1,j,i) = 0.0
    153153          ENDDO
    … …  
    157157!
    158158!-- Random perturbation is added to the array to be disturbed.
    159     DO  i = nxl-1, nxr+1
    160        DO  j = nys-1, nyn+1
     159    DO  i = nxlg, nxrg
     160       DO  j = nysg, nyng
    161161          DO  k = disturbance_level_ind_b-2, disturbance_level_ind_t+2
    162162             field(k,j,i) = field(k,j,i) + dist1(k,j,i)

  • palm/trunk/SOURCE/exchange_horiz.f90

    r484 r667  
    1  SUBROUTINE exchange_horiz( ar )
     1 SUBROUTINE exchange_horiz( ar, nbgp_local)
    22
    33!------------------------------------------------------------------------------!
    44! Current revisions:
    55! -----------------
    6 !
     6! Dynamic exchange of ghost points with nbgp_local to ensure that no useless
     7! ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0)
     8! used for normal grid, the remaining types used for the several grid levels.
     9! Exchange is done via MPI-Vectors with a dynamic value of ghost points which
     10! depend on the advection scheme. Exchange of left and right PEs is 10% faster
     11! with MPI-Vectors than without.
    712!
    813! Former revisions:
    … …  
    4146    INTEGER, DIMENSION(MPI_STATUS_SIZE,4) ::  wait_stat
    4247#endif
    43 
    44     REAL ::  ar(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    45 
     48    INTEGER :: i,nbgp_local
     49    REAL, DIMENSION(nzb:nzt+1,nys-nbgp_local:nyn+nbgp_local, &
     50                    nxl-nbgp_local:nxr+nbgp_local) ::  ar
    4651
    4752    CALL cpu_log( log_point_s(2), 'exchange_horiz', 'start' )
    4853
     54    IF ( exchange_mg == .TRUE. ) THEN
     55      i = grid_level
     56    ELSE
     57      i = 0
     58    END IF
    4959#if defined( __parallel )
    5060
    … …  
    5666!--    within the PE memory
    5767       IF ( bc_lr == 'cyclic' )  THEN
    58           ar(:,nys:nyn,nxl-1) = ar(:,nys:nyn,nxr)
    59           ar(:,nys:nyn,nxr+1) = ar(:,nys:nyn,nxl)
     68          ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)
     69          ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1)
    6070       ENDIF
    6171
    … …  
    6575!
    6676!--    Send left boundary, receive right one
    67        CALL MPI_ISEND(                                                     &
    68                ar(nzb,nys-1,nxl), ngp_yz(grid_level), MPI_REAL, pleft,  0, &
    69                           comm2d, req(1), ierr )
    70        CALL MPI_IRECV(                                                       &
    71                ar(nzb,nys-1,nxr+1), ngp_yz(grid_level), MPI_REAL, pright, 0, &
    72                           comm2d, req(2), ierr )
     77       CALL MPI_ISEND(ar(nzb,nys-nbgp_local,nxl),1,type_yz(i),pleft,0,comm2d,&
     78                      req(1),ierr)
     79       CALL MPI_IRECV(ar(nzb,nys-nbgp_local,nxr+1),1,type_yz(i),pright,0,&
     80                     comm2d,req(2),ierr)
    7381!
    7482!--    Send right boundary, receive left one
    75        CALL MPI_ISEND(                                                     &
    76                ar(nzb,nys-1,nxr), ngp_yz(grid_level), MPI_REAL, pright, 1, &
    77                           comm2d, req(3), ierr )
    78        CALL MPI_IRECV(                                                       &
    79                ar(nzb,nys-1,nxl-1), ngp_yz(grid_level), MPI_REAL, pleft,  1, &
    80                           comm2d, req(4), ierr )
     83
     84
     85       CALL MPI_ISEND(ar(nzb,nys-nbgp_local,nxr+1-nbgp_local),1,type_yz(i),pright, 1,  &
     86                      comm2d, req(3), ierr )
     87       CALL MPI_IRECV(ar(nzb,nys-nbgp_local,nxl-nbgp_local),1,type_yz(i),pleft,1,&
     88                      comm2d,req(4), ierr)
     89
    8190       CALL MPI_WAITALL( 4, req, wait_stat, ierr )
    8291
    … …  
    8998!--    within the PE memory
    9099       IF ( bc_ns == 'cyclic' )  THEN
    91           ar(:,nys-1,:) = ar(:,nyn,:)
    92           ar(:,nyn+1,:) = ar(:,nys,:)
     100          ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)
     101          ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:)
    93102       ENDIF
    94103
    … …  
    98107!
    99108!--    Send front boundary, receive rear one
    100        CALL MPI_ISEND( ar(nzb,nys,nxl-1),   1, type_xz(grid_level), psouth, 0, &
     109!--    MPI_ISEND initial send adress changed, type_xz() is sendet nbgp times
     110
     111       CALL MPI_ISEND( ar(nzb,nys,nxl-nbgp_local),1, type_xz(i), psouth, 0, &
    101112                       comm2d, req(1), ierr )
    102        CALL MPI_IRECV( ar(nzb,nyn+1,nxl-1), 1, type_xz(grid_level), pnorth, 0, &
     113       CALL MPI_IRECV( ar(nzb,nyn+1,nxl-nbgp_local),1, type_xz(i), pnorth, 0, &
    103114                       comm2d, req(2), ierr )
    104115!
    105116!--    Send rear boundary, receive front one
    106        CALL MPI_ISEND( ar(nzb,nyn,nxl-1),   1, type_xz(grid_level), pnorth, 1, &
     117       CALL MPI_ISEND( ar(nzb,nyn-nbgp_local+1,nxl-nbgp_local),1, type_xz(i), pnorth, 1, &
    107118                       comm2d, req(3), ierr )
    108        CALL MPI_IRECV( ar(nzb,nys-1,nxl-1), 1, type_xz(grid_level), psouth, 1, &
     119       CALL MPI_IRECV( ar(nzb,nys-nbgp_local,nxl-nbgp_local),1, type_xz(i), psouth, 1, &
    109120                       comm2d, req(4), ierr )
    110121       call MPI_WAITALL( 4, req, wait_stat, ierr )
    111122
    112123    ENDIF
    113 
    114124
    115125#else
    … …  
    118128!-- Lateral boundary conditions in the non-parallel case
    119129    IF ( bc_lr == 'cyclic' )  THEN
    120        ar(:,nys:nyn,nxl-1) = ar(:,nys:nyn,nxr)
    121        ar(:,nys:nyn,nxr+1) = ar(:,nys:nyn,nxl)
     130        ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)
     131        ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1)
    122132    ENDIF
    123133
    124134    IF ( bc_ns == 'cyclic' )  THEN
    125        ar(:,nys-1,:) = ar(:,nyn,:)
    126        ar(:,nyn+1,:) = ar(:,nys,:)
     135        ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)
     136        ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:)
    127137    ENDIF
    128138
    129139#endif
    130 
    131140    CALL cpu_log( log_point_s(2), 'exchange_horiz', 'stop' )
    132141
     142
    133143 END SUBROUTINE exchange_horiz

  • palm/trunk/SOURCE/exchange_horiz_2d.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! Dynamic exchange of ghost points with nbgp, which depends on the advection
     7! scheme. Exchange between left and right PEs is now done with MPI-vectors.
    78!
    89! Former revisions:
    … …  
    3738    IMPLICIT NONE
    3839
    39     REAL ::  ar(nys-1:nyn+1,nxl-1:nxr+1)
     40    REAL ::  ar(nysg:nyng,nxlg:nxrg)
     41    INTEGER :: i
    4042
    4143
    … …  
    5153!--    One-dimensional decomposition along y, boundary values can be exchanged
    5254!--    within the PE memory
    53        ar(nys:nyn,nxl-1) = ar(nys:nyn,nxr)
    54        ar(nys:nyn,nxr+1) = ar(nys:nyn,nxl)
     55       ar(:,nxl-nbgp:nxl-1) = ar(:,nxr-nbgp+1:nxr)
     56       ar(:,nxr+1:nxr+nbgp) = ar(:,nxl:nxl+nbgp-1)
    5557
    5658    ELSE
    5759!
    5860!--    Send left boundary, receive right one
    59        CALL MPI_SENDRECV( ar(nys,nxl),   ngp_y, MPI_REAL, pleft,  0, &
    60                           ar(nys,nxr+1), ngp_y, MPI_REAL, pright, 0, &
     61
     62       CALL MPI_SENDRECV( ar(nysg,nxl), 1, type_y, pleft,  0, &
     63                          ar(nysg,nxr+1), 1, type_y, pright, 0, &
    6164                          comm2d, status, ierr )
    6265!
    6366!--    Send right boundary, receive left one
    64        CALL MPI_SENDRECV( ar(nys,nxr),   ngp_y, MPI_REAL, pright,  1, &
    65                           ar(nys,nxl-1), ngp_y, MPI_REAL, pleft,   1, &
     67       CALL MPI_SENDRECV( ar(nysg,nxr+1-nbgp), 1, type_y, pright,  1, &
     68                          ar(nysg,nxlg), 1, type_y, pleft,   1, &
    6669                          comm2d, status, ierr )
    6770    ENDIF
    … …  
    7174!--    One-dimensional decomposition along x, boundary values can be exchanged
    7275!--    within the PE memory
    73        ar(nys-1,:) = ar(nyn,:)
    74        ar(nyn+1,:) = ar(nys,:)
     76       ar(nys-nbgp:nys-1,:) = ar(nyn-nbgp+1:nyn,:)
     77       ar(nyn+1:nyn+nbgp,:) = ar(nys:nys+nbgp-1,:)
    7578
    7679    ELSE
    7780!
    7881!--    Send front boundary, receive rear one
    79        CALL MPI_SENDRECV( ar(nys,nxl-1),   1, type_x, psouth, 0, &
    80                           ar(nyn+1,nxl-1), 1, type_x, pnorth, 0, &
     82
     83       CALL MPI_SENDRECV( ar(nys,nxlg), 1, type_x, psouth, 0, &        !replace number of sended elements from
     84                          ar(nyn+1,nxlg), 1, type_x, pnorth, 0, &      ! nbgp to 1
    8185                          comm2d, status, ierr )
    8286!
    8387!--    Send rear boundary, receive front one
    84        CALL MPI_SENDRECV( ar(nyn,nxl-1),   1, type_x, pnorth, 1, &
    85                           ar(nys-1,nxl-1), 1, type_x, psouth, 1, &
    86                           comm2d, status, ierr )
     88       CALL MPI_SENDRECV( ar(nyn+1-nbgp,nxlg), 1, type_x, pnorth, 1, &
     89                          ar(nysg,nxlg), 1, type_x, psouth, 1, &
     90                          comm2d, status, ierr )
     91
    8792    ENDIF
    8893
    … …  
    9297!-- Lateral boundary conditions in the non-parallel case
    9398    IF ( bc_lr == 'cyclic' )  THEN
    94        ar(nys:nyn,nxl-1) = ar(nys:nyn,nxr)
    95        ar(nys:nyn,nxr+1) = ar(nys:nyn,nxl)
     99       ar(:,nxl-nbgp:nxl-1) = ar(:,nxr-nbgp+1:nxr)
     100       ar(:,nxr+1:nxr+nbgp) = ar(:,nxl:nxl+nbgp-1)
    96101    ENDIF
    97102
    98103    IF ( bc_ns == 'cyclic' )  THEN
    99        ar(nys-1,:) = ar(nyn,:)
    100        ar(nyn+1,:) = ar(nys,:)
    101     ENDIF
     104       ar(nys-nbgp:nys-1,:) = ar(nyn-nbgp+1:nyn,:)
     105       ar(nyn+1:nyn+nbgp,:) = ar(nys:nys+nbgp-1,:)
     106    ENDIF
     107
    102108
    103109#endif
    … …  
    106112!-- Neumann-conditions at inflow/outflow in case of non-cyclic boundary
    107113!-- conditions
    108     IF ( inflow_l .OR. outflow_l )  ar(:,nxl-1) = ar(:,nxl)
    109     IF ( inflow_r .OR. outflow_r )  ar(:,nxr+1) = ar(:,nxr)
    110     IF ( inflow_s .OR. outflow_s )  ar(nys-1,:) = ar(nys,:)
    111     IF ( inflow_n .OR. outflow_n )  ar(nyn+1,:) = ar(nyn,:)
    112 
     114    IF ( inflow_l .OR. outflow_l )  THEN
     115       DO i=nbgp, 1, -1
     116         ar(:,nxl-i) = ar(:,nxl)
     117       END DO
     118    END IF
     119    IF ( inflow_r .OR. outflow_r )  THEN
     120       DO i=1, nbgp
     121          ar(:,nxr+i) = ar(:,nxr)
     122       END DO
     123    END IF
     124    IF ( inflow_s .OR. outflow_s )  THEN
     125       DO i=nbgp, 1, -1
     126         ar(nys-i,:) = ar(nys,:)
     127       END DO
     128    END IF
     129    IF ( inflow_n .OR. outflow_n )  THEN
     130       DO i=1, nbgp
     131         ar(nyn+i,:) = ar(nyn,:)
     132       END DO
     133    END IF
    113134    CALL cpu_log( log_point_s(13), 'exchange_horiz_2d', 'stop' )
    114135
    … …  
    134155    IMPLICIT NONE
    135156
    136     INTEGER ::  ar(nys-1:nyn+1,nxl-1:nxr+1)
     157    REAL ::  ar(nysg:nyng,nxlg:nxrg)
     158    INTEGER :: i
     159
    137160
    138161
    … …  
    154177!
    155178!--    Send left boundary, receive right one
    156        CALL MPI_SENDRECV( ar(nys,nxl),   ngp_y, MPI_INTEGER, pleft,  0, &
    157                           ar(nys,nxr+1), ngp_y, MPI_INTEGER, pright, 0, &
     179       CALL MPI_SENDRECV( ar(nysg,nxl), 1, type_y_int, pleft,  0, &
     180                          ar(nysg,nxr+1), 1, type_y_int, pright, 0, &
    158181                          comm2d, status, ierr )
    159182!
    160183!--    Send right boundary, receive left one
    161        CALL MPI_SENDRECV( ar(nys,nxr),   ngp_y, MPI_INTEGER, pright,  1, &
    162                           ar(nys,nxl-1), ngp_y, MPI_INTEGER, pleft,   1, &
    163                           comm2d, status, ierr )
     184       CALL MPI_SENDRECV( ar(nysg,nxr+1-nbgp), 1, type_y_int, pright,  1, &
     185                          ar(nysg,nxlg), 1, type_y_int, pleft,   1, &
     186                          comm2d, status, ierr )
     187
    164188    ENDIF
    165189
    … …  
    168192!--    One-dimensional decomposition along x, boundary values can be exchanged
    169193!--    within the PE memory
    170        ar(nys-1,:) = ar(nyn,:)
    171        ar(nyn+1,:) = ar(nys,:)
     194       ar(nysg:nys-1,:) = ar(nyn+1-nbgp:nyn,:)
     195       ar(nyn+1:nyng,:) = ar(nys:nys-1+nbgp,:)
     196
    172197
    173198    ELSE
    174199!
    175200!--    Send front boundary, receive rear one
    176        CALL MPI_SENDRECV( ar(nys,nxl-1),   1, type_x_int, psouth, 0, &
    177                           ar(nyn+1,nxl-1), 1, type_x_int, pnorth, 0, &
    178                           comm2d, status, ierr )
     201       CALL MPI_SENDRECV( ar(nys,nxlg),1, type_x_int, psouth, 0, &
     202                          ar(nyn+1,nxlg),1, type_x_int, pnorth, 0, &
     203                          comm2d, status, ierr )
     204
    179205!
    180206!--    Send rear boundary, receive front one
    181        CALL MPI_SENDRECV( ar(nyn,nxl-1),   1, type_x_int, pnorth, 1, &
    182                           ar(nys-1,nxl-1), 1, type_x_int, psouth, 1, &
    183                           comm2d, status, ierr )
     207       CALL MPI_SENDRECV( ar(nyn+1-nbgp,nxlg), nbgp, type_x_int, pnorth, 1, &
     208                          ar(nysg,nxlg), nbgp, type_x_int, psouth, 1, &
     209                          comm2d, status, ierr )
     210
    184211    ENDIF
    185212
    … …  
    194221
    195222    IF ( bc_ns == 'cyclic' )  THEN
    196        ar(nys-1,:) = ar(nyn,:)
    197        ar(nyn+1,:) = ar(nys,:)
     223       ar(nysg:nys-1,:) = ar(nyn+1-nbgp:nyn,:)
     224       ar(nyn+1:nyng,:) = ar(nys:nys-1+nbgp,:)
    198225    ENDIF
    199226

  • palm/trunk/SOURCE/flow_statistics.f90

    r625 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! When advection is computed with ws-scheme, turbulent fluxes are already
     7! computed in the respective advection routines and buffered in arrays
     8! sums_xx_ws_l(). This is due to a consistent treatment of statistics with the
     9! numerics and to avoid unphysical kinks near the surface.
     10! So some if requests has to be done to dicern between fluxes from ws-scheme
     11! other advection schemes.
     12! Furthermore the computation of z_i is only done if the heat flux exceeds a
     13! minimum value. This affects only simulations of a neutral boundary layer and
     14! is due to reasons of computations in the advection scheme.
     15!
    716!
    817! Former revisions:
    … …  
    102111    REAL    ::  sums_ll(nzb:nzt+1,2)
    103112
    104 
    105113    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
    106114
    … …  
    133141       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
    134142
     143!
     144!--    Copy the turbulent quantities, evaluated in the advection routines to
     145!--    the local array sums_l() for further computations
     146       IF ( ws_scheme_mom )  THEN
     147!       
     148!--       Boundary condition for u'u' and v'v', because below the surface no
     149!--       computation for these quantities is done.
     150          DO  i = nxl, nxr
     151             DO  j =  nys, nyn
     152                sums_us2_ws_l(nzb_u_inner(j,i),sr) =                          &
     153                    sums_us2_ws_l(nzb_u_inner(j,i)+1,sr)
     154                sums_vs2_ws_l(nzb_v_inner(j,i),sr) =                          & 
     155                    sums_vs2_ws_l(nzb_v_inner(j,i)+1,sr)
     156             ENDDO
     157          ENDDO
     158!         
     159!--       Swap the turbulent quantities evaluated in advec_ws.
     160          sums_l(:,13,0) = sums_wsus_ws_l(:,sr)       ! w*u*
     161          sums_l(:,15,0) = sums_wsvs_ws_l(:,sr)       ! w*v*
     162          sums_l(:,30,0) = sums_us2_ws_l(:,sr)        ! u*2
     163          sums_l(:,31,0) = sums_vs2_ws_l(:,sr)        ! v*2
     164          sums_l(:,32,0) = sums_ws2_ws_l(:,sr)        ! w*2
     165          sums_l(:,34,0) = sums_l(:,34,0) + 0.5 *                             &
     166                (sums_us2_ws_l(:,sr) + sums_vs2_ws_l(:,sr)                    &
     167                + sums_ws2_ws_l(:,sr))                      ! e*
     168          DO  k = nzb, nzt
     169             sums_l(nzb+5,pr_palm,0) = sums_l(nzb+5,pr_palm,0) + 0.5 * (      &
     170                sums_us2_ws_l(k,sr) + sums_vs2_ws_l(k,sr) +                   &
     171                sums_ws2_ws_l(k,sr))
     172          ENDDO
     173       ENDIF
     174       IF ( ws_scheme_sca )  THEN
     175          sums_l(:,17,0) = sums_wspts_ws_l(:,sr)      ! w*pts* from advec_s_ws
     176          IF ( ocean ) sums_l(:,66,0) = sums_wssas_ws_l(:,sr) ! w*sa*
     177          IF ( humidity  .OR.  passive_scalar ) sums_l(:,49,0) =              &
     178                                                   sums_wsqs_ws_l(:,sr) !w*q*
     179       ENDIF
    135180!
    136181!--    Horizontally averaged profiles of horizontal velocities and temperature.
    … …  
    138183!--    for other horizontal averages.
    139184       tn = 0
     185
    140186       !$OMP PARALLEL PRIVATE( i, j, k, tn )
    141187#if defined( __intel_openmp_bug )
    … …  
    215261       ENDIF
    216262       !$OMP END PARALLEL
    217 
    218263!
    219264!--    Summation of thread sums
    … …  
    304349       hom(:,1,2,sr) = sums(:,2)             ! v
    305350       hom(:,1,4,sr) = sums(:,4)             ! pt
     351
    306352
    307353!
    … …  
    354400          DO  j =  nys, nyn
    355401             sums_l_etot = 0.0
    356              sums_l_eper = 0.0
    357402             DO  k = nzb_s_inner(j,i), nzt+1
    358                 u2   = u(k,j,i)**2
    359                 v2   = v(k,j,i)**2
    360                 w2   = w(k,j,i)**2
    361                 ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
    362                 vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
    363403!
    364404!--             Prognostic and diagnostic variables
    … …  
    369409                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
    370410
    371 !
    372 !--             Variances
    373                 sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
    374                 sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
    375                 sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
    376411                sums_l(k,33,tn) = sums_l(k,33,tn) + &
    377412                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
    … …  
    381416                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
    382417                ENDIF
    383 !
    384 !--             Higher moments
    385 !--             (Computation of the skewness of w further below)
    386                 sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
    387                                                     rmask(j,i,sr)
    388 !
    389 !--             Perturbation energy
    390                 sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 * ( ust2 + vst2 + w2 ) &
    391                                                     * rmask(j,i,sr)
     418
    392419                sums_l_etot  = sums_l_etot + &
    393                                         0.5 * ( u2 + v2 + w2 ) * rmask(j,i,sr)
    394                 sums_l_eper  = sums_l_eper + &
    395                                         0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
     420                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
     421                                        w(k,j,i)**2 ) * rmask(j,i,sr)
    396422             ENDDO
    397423!
    … …  
    401427!--          allow vectorization of that loop.
    402428             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
    403              sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) + sums_l_eper
    404429!
    405430!--          2D-arrays (being collected in the last column of sums_l)
    … …  
    420445
    421446!
     447!--    Computation of statistics when ws-scheme is not used. Else these
     448!--    quantities are evaluated in the advection routines.
     449       IF ( .NOT. ws_scheme_mom )  THEN
     450          !$OMP DO
     451          DO  i = nxl, nxr
     452             DO  j =  nys, nyn
     453                sums_l_eper = 0.0
     454                DO  k = nzb_s_inner(j,i), nzt+1
     455                   u2   = u(k,j,i)**2
     456                   v2   = v(k,j,i)**2
     457                   w2   = w(k,j,i)**2
     458                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
     459                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
     460
     461                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
     462                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
     463                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
     464!
     465!--   Higher moments
     466!--  (Computation of the skewness of w further below)
     467                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
     468                                                    rmask(j,i,sr)
     469!
     470!--             Perturbation energy
     471
     472                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
     473                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
     474                   sums_l_eper  = sums_l_eper + &
     475                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
     476
     477                ENDDO
     478                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
     479                     + sums_l_eper
     480             ENDDO
     481          ENDDO
     482       ELSE
     483          !$OMP DO
     484          DO  i = nxl, nxr
     485             DO  j =  nys, nyn
     486                DO  k = nzb_s_inner(j,i), nzt + 1
     487                   w2   = w(k,j,i)**2
     488!
     489!--                Higher moments
     490!--                (Computation of the skewness of w further below)
     491                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
     492                                                    rmask(j,i,sr)
     493                ENDDO
     494             ENDDO
     495          ENDDO
     496       ENDIF
     497
     498!
    422499!--    Horizontally averaged profiles of the vertical fluxes
     500
    423501       !$OMP DO
    424502       DO  i = nxl, nxr
    … …  
    587665                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
    588666                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
    589 !
    590 !--             Momentum flux w*u*
    591                 sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                     &
    592                                                     ( w(k,j,i-1) + w(k,j,i) ) &
    593                                                     * ust * rmask(j,i,sr)
    594 !
    595 !--             Momentum flux w*v*
    596                 sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                     &
    597                                                     ( w(k,j-1,i) + w(k,j,i) ) &
    598                                                     * vst * rmask(j,i,sr)
    599 !
    600 !--             Heat flux w*pt*
    601 !--             The following formula (comment line, not executed) does not
    602 !--             work if applied to subregions
    603 !                sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 *                     &
    604 !                                                    ( pt(k,j,i)+pt(k+1,j,i) ) &
    605 !                                                    * w(k,j,i) * rmask(j,i,sr)
    606                 sums_l(k,17,tn) = sums_l(k,17,tn) + pts * w(k,j,i) * &
    607                                                     rmask(j,i,sr)
    608 !
     667
    609668!--             Higher moments
    610669                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
    … …  
    617676!--             but so far there is no other suitable place to calculate)
    618677                IF ( ocean )  THEN
    619                    pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
     678                   IF( .NOT. ws_scheme_sca )  THEN
     679                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
    620680                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
    621                    sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
     681                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
    622682                                                       rmask(j,i,sr)
     683                   ENDIF
    623684                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
    624685                                                       rmask(j,i,sr)
    … …  
    635696                   sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
    636697                                                       rmask(j,i,sr)
    637                    pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
    638                                  q(k+1,j,i) - hom(k+1,1,41,sr) )
    639                    sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
    640                                                        rmask(j,i,sr)
     698
    641699                   IF ( cloud_physics  .OR.  cloud_droplets )  THEN
    642700                      pts = 0.5 *                                            &
    … …  
    652710!
    653711!--             Passive scalar flux
    654                 IF ( passive_scalar )  THEN
     712                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca ))  THEN
    655713                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
    656714                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
    … …  
    661719!
    662720!--             Energy flux w*e*
    663                 sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *           &
    664                                               ( ust**2 + vst**2 + w(k,j,i)**2 )&
    665                                               * rmask(j,i,sr)
    666          
     721!--             has to be adjusted
     722                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
     723                                             ( ust**2 + vst**2 + w(k,j,i)**2 )&
     724                                             * rmask(j,i,sr)
    667725             ENDDO
    668726          ENDDO
    669727       ENDDO
     728!-     for reasons of speed optimization the loop is splitted, to avoid if-else
     729!-     statements inside the loop
     730!-     Fluxes which have been computed in part directly inside the advection routines
     731!-     treated seperatly.
     732!-     First treat the momentum fluxes
     733       IF ( .NOT. ws_scheme_mom )  THEN
     734         !$OMP DO
     735         DO  i = nxl, nxr
     736            DO  j = nys, nyn
     737               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
     738                  ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
     739                              u(k+1,j,i) - hom(k+1,1,1,sr) )
     740                  vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
     741                              v(k+1,j,i) - hom(k+1,1,2,sr) )
     742!                             
     743!--               Momentum flux w*u*
     744                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
     745                                                    ( w(k,j,i-1) + w(k,j,i) ) &
     746                                                    * ust * rmask(j,i,sr)
     747!
     748!--               Momentum flux w*v*
     749                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
     750                                                    ( w(k,j-1,i) + w(k,j,i) ) &
     751                                                    * vst * rmask(j,i,sr)
     752               ENDDO
     753            ENDDO
     754         ENDDO
     755
     756       ENDIF
     757       IF ( .NOT. ws_scheme_sca )  THEN
     758         !$OMP DO
     759         DO  i = nxl, nxr
     760            DO  j = nys, nyn
     761               DO  k = nzb_diff_s_inner(j,i) - 1, nzt_diff
     762!-                vertical heat flux
     763                  sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
     764                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
     765                           pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
     766                           * w(k,j,i) * rmask(j,i,sr)
     767                  IF ( humidity )  THEN
     768                     pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
     769                                q(k+1,j,i) - hom(k+1,1,41,sr) )
     770                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
     771                                                      rmask(j,i,sr)
     772                  ENDIF
     773               ENDDO
     774            ENDDO
     775         ENDDO
     776
     777       ENDIF
     778
    670779
    671780!
    … …  
    814923
    815924#if defined( __parallel )
     925
    816926!
    817927!--    Compute total sum from local sums
    … …  
    843953          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
    844954       ENDDO
     955
    845956!--    Upstream-parts
    846957       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
    … …  
    868979          ENDDO
    869980       ENDIF
    870 
    871981!
    872982!--    Collect horizontal average in hom.
    … …  
    9341044                                       ! upstream-parts u_x, u_y, u_z, v_x,
    9351045                                       ! v_y, usw. (in last but one profile)
    936        hom(:,1,pr_palm,sr) =   sums(:,pr_palm) 
     1046       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
    9371047                                       ! u*, w'u', w'v', t* (in last profile)
    9381048
    … …  
    9511061       z_i(1) = 0.0
    9521062       first = .TRUE.
     1063
    9531064       IF ( ocean )  THEN
    9541065          DO  k = nzt, nzb+1, -1
    955              IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
     1066             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
     1067                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
    9561068                first = .FALSE.
    9571069                height = zw(k)
    9581070             ENDIF
    9591071             IF ( hom(k,1,18,sr) < 0.0  .AND. &
     1072                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
    9601073                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
    9611074                IF ( zw(k) < 1.5 * height )  THEN
    … …  
    9691082       ELSE
    9701083          DO  k = nzb, nzt-1
    971              IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
     1084             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
     1085                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
    9721086                first = .FALSE.
    9731087                height = zw(k)
    9741088             ENDIF
    975              IF ( hom(k,1,18,sr) < 0.0  .AND. &
     1089             IF ( hom(k,1,18,sr) < 0.0  .AND. &
     1090                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
    9761091                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
    9771092                IF ( zw(k) < 1.5 * height )  THEN
    … …  
    10261141!--    the characteristic convective boundary layer temperature.
    10271142!--    The horizontal average at nzb+1 is input for the average temperature.
    1028        IF ( hom(nzb,1,18,sr) > 0.0  .AND.  z_i(1) /= 0.0 )  THEN
     1143       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
     1144           .AND.  z_i(1) /= 0.0 )  THEN
    10291145          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
    10301146                                       hom(nzb,1,18,sr) *      &

  • palm/trunk/SOURCE/global_min_max.f90

    r623 r667  
    55! Current revisions:
    66! -----------------
    7 !
     7! Adapting of the index arrays, because MINLOC assumes lowerbound at 1 and not
     8! at nbgp.
    89!
    910! Former revisions:
    … …  
    5960!--    Determine the local minimum
    6061       fmin_ijk_l = MINLOC( ar )
    61        fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - 1    ! MINLOC assumes lowerbound = 1
    62        fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - 1
     62       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - nbgp    ! MINLOC assumes lowerbound = 1
     63       fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
    6364       fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - 1
    6465       fmin_l(1)  = ar(fmin_ijk_l(1),fmin_ijk_l(2),fmin_ijk_l(3))
    … …  
    100101!--    Determine the local maximum
    101102       fmax_ijk_l = MAXLOC( ar )
    102        fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - 1    ! MAXLOC assumes lowerbound = 1
    103        fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - 1
     103       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - nbgp    ! MAXLOC assumes lowerbound = 1
     104       fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
    104105       fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - 1
    105106       fmax_l(1) = ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3))
    … …  
    221222          IF ( fmax_ijk(1) < 0 )  THEN
    222223             value        = -value
    223              value_ijk(1) = -value_ijk(1) - 10
     224             value_ijk(1) = -value_ijk(1) - 10         !???
    224225          ENDIF
    225226
    … …  
    228229!
    229230!-- Limit index values to the range 0..nx, 0..ny
    230     IF ( value_ijk(3) ==   -1 )  value_ijk(3) = nx
    231     IF ( value_ijk(3) == nx+1 )  value_ijk(3) =  0
    232     IF ( value_ijk(2) ==   -1 )  value_ijk(2) = ny
    233     IF ( value_ijk(2) == ny+1 )  value_ijk(2) =  0
     231    IF ( value_ijk(3) < 0  ) value_ijk(3) = nx +1 + value_ijk(3)
     232    IF ( value_ijk(3) > nx ) value_ijk(3) = value_ijk(3) - (nx+1)
     233    IF ( value_ijk(2) < 0  ) value_ijk(2) = ny +1 + value_ijk(2)
     234    IF ( value_ijk(2) > ny ) value_ijk(2) = value_ijk(2) - (ny+1)
    234235
    235236

  • palm/trunk/SOURCE/header.f90

    r581 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! Output of advection scheme.
     7! Modified output of Prandtl-layer height.
    78!
    89! Former revisions:
    … …  
    250251    IF ( momentum_advec == 'pw-scheme' )  THEN
    251252       WRITE ( io, 113 )
    252     ELSE
     253    ELSEIF (momentum_advec == 'ws-scheme' ) THEN
     254       WRITE ( io, 503 )
     255    ELSEIF (momentum_advec == 'ups-scheme' ) THEN
    253256       WRITE ( io, 114 )
    254257       IF ( cut_spline_overshoot )  WRITE ( io, 124 )
    … …  
    267270    IF ( scalar_advec == 'pw-scheme' )  THEN
    268271       WRITE ( io, 116 )
     272    ELSEIF ( scalar_advec == 'ws-scheme' )  THEN
     273       WRITE ( io, 504 )
    269274    ELSEIF ( scalar_advec == 'ups-scheme' )  THEN
    270275       WRITE ( io, 117 )
    … …  
    575580    ELSEIF( ibc_pt_t == 2 )  THEN
    576581       roben  = TRIM( roben  ) // ' pt(nzt+1) = pt(nzt) + dpt/dz_ini'
     582
    577583    ENDIF
    578584
    … …  
    662668
    663669    IF ( prandtl_layer )  THEN
    664        WRITE ( io, 305 )  0.5 * (zu(1)-zu(0)), roughness_length, kappa, &
     670       WRITE ( io, 305 )  (zu(1)-zu(0)), roughness_length, kappa, &
    665671                          rif_min, rif_max
    666672       IF ( .NOT. constant_heatflux )  WRITE ( io, 308 )
    … …  
    19811987            '    Dissipation calculation:           ',A/)
    19821988502 FORMAT ('    Damping layer starts from ',F7.1,' m (GP ',I4,')'/)
     1989503 FORMAT (' --> Momentum advection via Wicker-Skamarock-Scheme 5th order')
     1990504 FORMAT (' --> Scalar advection via Wicker-Skamarock-Scheme 5th order')
    19831991
    19841992

  • palm/trunk/SOURCE/inflow_turbulence.f90

    r623 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! Using nbgp recycling planes for a better resolution of the turbulent flow
     7! near the inflow.
    78!
    89! Former revisions:
    … …  
    3536    IMPLICIT NONE
    3637
    37     INTEGER ::  i, imax, j, k, ngp_ifd, ngp_pr
     38    INTEGER ::  i, imax, j, k, l, ngp_ifd, ngp_pr
    3839
    3940    REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
    40     REAL, DIMENSION(nzb:nzt+1,5) ::  avpr, avpr_l
    41     REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,5) ::  inflow_dist
     41    REAL, DIMENSION(nzb:nzt+1,5,nbgp) ::  avpr, avpr_l
     42    REAL, DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) ::  inflow_dist
    4243
    4344    CALL cpu_log( log_point(40), 'inflow_turbulence', 'start' )
    4445
    4546!
    46 !-- Carry out horizontal averaging in the recycling plane
     47!-- Carry out spanwise averaging in the recycling plane
    4748    avpr_l = 0.0
    48     ngp_pr = ( nzt - nzb + 2 ) * 5
    49     ngp_ifd = ngp_pr * ( nyn - nys + 3 )
     49    ngp_pr = ( nzt - nzb + 2 ) * 5 * nbgp
     50    ngp_ifd = ngp_pr * ( nyn - nys + 1 + 2 * nbgp )
    5051
    5152!
    5253!-- First, local averaging within the recycling domain
    53     IF ( recycling_plane >= nxl )  THEN
    54 
    55        imax = MIN( nxr, recycling_plane )
    56 
    57        DO  i = nxl, imax
     54
     55    i = recycling_plane
     56
     57#if defined( __parallel )
     58    IF ( myidx == id_recycling )  THEN
     59       
     60       DO  l = 1, nbgp
    5861          DO  j = nys, nyn
    59              DO  k = nzb, nzt+1
    60 
    61                 avpr_l(k,1) = avpr_l(k,1) + u(k,j,i)
    62                 avpr_l(k,2) = avpr_l(k,2) + v(k,j,i)
    63                 avpr_l(k,3) = avpr_l(k,3) + w(k,j,i)
    64                 avpr_l(k,4) = avpr_l(k,4) + pt(k,j,i)
    65                 avpr_l(k,5) = avpr_l(k,5) + e(k,j,i)
     62             DO  k = nzb, nzt + 1
     63
     64                avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
     65                avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
     66                avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
     67                avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
     68                avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
    6669
    6770             ENDDO
    6871          ENDDO
    69        ENDDO
    70 
    71     ENDIF
    72 
    73 !    WRITE (9,*) '*** averaged profiles avpr_l'
    74 !    DO  k = nzb, nzt+1
    75 !       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)
    76 !    ENDDO
    77 !    WRITE (9,*) ' '
    78 
    79 #if defined( __parallel )
     72          i = i + 1
     73       ENDDO
     74
     75    ENDIF
    8076!
    8177!-- Now, averaging over all PEs
    8278    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    83     CALL MPI_ALLREDUCE( avpr_l(nzb,1), avpr(nzb,1), ngp_pr, MPI_REAL, MPI_SUM, &
    84                         comm2d, ierr )
     79    CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr,  &
     80                    MPI_REAL, MPI_SUM, comm2d, ierr )
     81
    8582#else
     83    DO  l = 1, nbgp
     84       DO  j = nys, nyn
     85          DO  k = nzb, nzt + 1
     86
     87             avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
     88             avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
     89             avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
     90             avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
     91             avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
     92
     93          ENDDO
     94       ENDDO
     95       i = i + 1
     96    ENDDO
     97   
    8698    avpr = avpr_l
    8799#endif
    88100
    89     avpr = avpr / ( ( ny + 1 ) * ( recycling_plane + 1 ) )
    90 
    91 !    WRITE (9,*) '*** averaged profiles'
    92 !    DO  k = nzb, nzt+1
    93 !       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)
    94 !    ENDDO
    95 !    WRITE (9,*) ' '
    96 
     101    avpr = avpr / ( ny + 1 )
    97102!
    98103!-- Calculate the disturbances at the recycling plane
    … …  
    101106#if defined( __parallel )
    102107    IF ( myidx == id_recycling )  THEN
    103 
    104        DO  j = nys-1, nyn+1
     108       DO  l = 1, nbgp
     109          DO  j = nysg, nyng
     110             DO  k = nzb, nzt + 1
     111
     112                inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
     113                inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
     114                inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
     115                inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
     116                inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
     117             
     118            ENDDO
     119          ENDDO
     120          i = i + 1
     121       ENDDO
     122
     123    ENDIF
     124#else
     125    DO  l = 1, nbgp
     126       DO  j = nysg, nyng
    105127          DO  k = nzb, nzt+1
    106128
    107               inflow_dist(k,j,1) = u(k,j,i+1) - avpr(k,1)
    108               inflow_dist(k,j,2) = v(k,j,i)   - avpr(k,2)
    109               inflow_dist(k,j,3) = w(k,j,i)   - avpr(k,3)
    110               inflow_dist(k,j,4) = pt(k,j,i)  - avpr(k,4)
    111               inflow_dist(k,j,5) = e(k,j,i)   - avpr(k,5)
    112 
    113           ENDDO
    114        ENDDO
    115 
    116     ENDIF
    117 #else
    118     DO  j = nys-1, nyn+1
    119        DO  k = nzb, nzt+1
    120 
    121           inflow_dist(k,j,1) = u(k,j,i+1) - avpr(k,1)
    122           inflow_dist(k,j,2) = v(k,j,i)   - avpr(k,2)
    123           inflow_dist(k,j,3) = w(k,j,i)   - avpr(k,3)
    124           inflow_dist(k,j,4) = pt(k,j,i)  - avpr(k,4)
    125           inflow_dist(k,j,5) = e(k,j,i)   - avpr(k,5)
    126 
    127        ENDDO
     129             inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
     130             inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
     131             inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
     132             inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
     133             inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
     134             
     135          ENDDO
     136       ENDDO
     137       i = i + 1
    128138    ENDDO
    129139#endif
    … …  
    134144    IF ( myidx == id_recycling  .AND.  myidx /= id_inflow )  THEN
    135145
    136        CALL MPI_SEND( inflow_dist(nzb,nys-1,1), ngp_ifd, MPI_REAL, &
     146       CALL MPI_SEND( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
    137147                      id_inflow, 1, comm1dx, ierr )
    138148
    … …  
    140150
    141151       inflow_dist = 0.0
    142        CALL MPI_RECV( inflow_dist(nzb,nys-1,1), ngp_ifd, MPI_REAL, &
     152       CALL MPI_RECV( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
    143153                      id_recycling, 1, comm1dx, status, ierr )
    144154
    … …  
    150160    IF ( nxl == 0 )  THEN
    151161
    152        DO  j = nys-1, nyn+1
    153           DO  k = nzb, nzt+1
    154 
    155 !              WRITE (9,*) 'j=',j,' k=',k
    156 !              WRITE (9,*) 'mean_u = ', mean_inflow_profiles(k,1), ' dist_u = ',&
    157 !                          inflow_dist(k,j,1)
    158 !              WRITE (9,*) 'mean_v = ', mean_inflow_profiles(k,2), ' dist_v = ',&
    159 !                          inflow_dist(k,j,2)
    160 !              WRITE (9,*) 'mean_w = 0.0', ' dist_w = ',&
    161 !                          inflow_dist(k,j,3)
    162 !              WRITE (9,*) 'mean_pt = ', mean_inflow_profiles(k,4), ' dist_pt = ',&
    163 !                          inflow_dist(k,j,4)
    164 !              WRITE (9,*) 'mean_e = ', mean_inflow_profiles(k,5), ' dist_e = ',&
    165 !                          inflow_dist(k,j,5)
    166               u(k,j,0)   = mean_inflow_profiles(k,1) + &
    167                            inflow_dist(k,j,1) * inflow_damping_factor(k)
    168               v(k,j,-1)  = mean_inflow_profiles(k,2) + &
    169                            inflow_dist(k,j,2) * inflow_damping_factor(k)
    170               w(k,j,-1)  = inflow_dist(k,j,3) * inflow_damping_factor(k)
    171               pt(k,j,-1) = mean_inflow_profiles(k,4) + &
    172                            inflow_dist(k,j,4) * inflow_damping_factor(k)
    173               e(k,j,-1)  = mean_inflow_profiles(k,5) + &
    174                            inflow_dist(k,j,5) * inflow_damping_factor(k)
    175               e(k,j,-1)  = MAX( e(k,j,-1), 0.0 )
     162       DO  j = nysg, nyng
     163          DO  k = nzb, nzt + 1
     164
     165              u(k,j,-nbgp+1:0)   = mean_inflow_profiles(k,1) + &
     166                           inflow_dist(k,j,1,1:nbgp) * inflow_damping_factor(k)
     167              v(k,j,-nbgp:-1)  = mean_inflow_profiles(k,2) + &
     168                           inflow_dist(k,j,2,1:nbgp) * inflow_damping_factor(k)
     169              w(k,j,-nbgp:-1)  = inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)
     170              pt(k,j,-nbgp:-1) = mean_inflow_profiles(k,4) + &
     171                           inflow_dist(k,j,4,1:nbgp) * inflow_damping_factor(k)
     172              e(k,j,-nbgp:-1)  = mean_inflow_profiles(k,5) + &
     173                           inflow_dist(k,j,5,1:nbgp) * inflow_damping_factor(k)
     174              e(k,j,-nbgp:-1)  = MAX( e(k,j,-nbgp:-1), 0.0 )
    176175
    177176          ENDDO

  • palm/trunk/SOURCE/init_1d_model.f90

    r392 r667  
    6363              l1d_m(nzb:nzt+1),  rif1d(nzb:nzt+1),   te_e(nzb:nzt+1),  &
    6464              te_em(nzb:nzt+1),  te_u(nzb:nzt+1),    te_um(nzb:nzt+1), &
    65               te_v(nzb:nzt+1),   te_vm(nzb:nzt+1),   u1d(nzb:nzt+1),   &
     65              te_v(nzb:nzt+1),   te_vm(nzb:nzt+1),    u1d(nzb:nzt+1),   &
    6666              u1d_m(nzb:nzt+1),  u1d_p(nzb:nzt+1),   v1d(nzb:nzt+1),   &
    6767              v1d_m(nzb:nzt+1),  v1d_p(nzb:nzt+1) )
    … …  
    385385!--       boundary condition applies to u and v.
    386386!--       The boundary condition for e is set further below ( (u*/cm)**2 ).
    387           u1d_p(nzb) = -u1d_p(nzb+1)
    388           v1d_p(nzb) = -v1d_p(nzb+1)
     387         ! u1d_p(nzb) = -u1d_p(nzb+1)
     388         ! v1d_p(nzb) = -v1d_p(nzb+1)
     389
     390          u1d_p(nzb) = 0.0
     391          v1d_p(nzb) = 0.0
    389392
    390393!

  • palm/trunk/SOURCE/init_3d_model.f90

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

  • palm/trunk/SOURCE/init_advec.f90

    r484 r667  
    88! Former revisions:
    99! -----------------
    10 ! $Id$
     10! $Id: init_advec.f90 484 2010-02-05 07:36:54Z raasch
    1111! RCS Log replace by Id keyword, revision history cleaned up
    1212!

  • palm/trunk/SOURCE/init_coupling.f90

    r484 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6!
     7! determination of target_id's moved to init_pegrid
     8!
    79!
    810! Former revisions:
    … …  
    2426    USE pegrid
    2527    USE control_parameters
     28    USE indices
    2629
    2730    IMPLICIT NONE
    … …  
    4851!-- ATTENTION: numprocs will be reset according to the new communicators
    4952#if defined ( __parallel )
     53
     54!myid_absolut = myid
    5055    IF ( myid == 0 )  THEN
    5156       READ (*,*,ERR=10,END=10)  coupling_mode, bc_data(1), bc_data(2)
    … …  
    9499
    95100       IF ( myid < bc_data(1) ) THEN
    96           inter_color = 0
    97           numprocs = bc_data(1)
     101          inter_color     = 0
     102          numprocs        = bc_data(1)
     103          coupling_mode   = 'atmosphere_to_ocean'
    98104       ELSE
    99           inter_color = 1
    100           numprocs = bc_data(2)
     105          inter_color     = 1
     106          numprocs        = bc_data(2)
     107          coupling_mode   = 'ocean_to_atmosphere'
    101108       ENDIF
    102 !
    103 !--    Calculate the target PE for coupling and set up the new communicators.
    104 !--    Currently only 1:1 topologies are supported.
    105        target_id = myid - ISIGN( numprocs, inter_color - 1 )
    106        IF ( inter_color == 0 ) THEN
    107           coupling_mode = 'atmosphere_to_ocean'
    108        ELSE
    109           coupling_mode = 'ocean_to_atmosphere'
    110        ENDIF
     109
    111110       CALL MPI_COMM_SPLIT( MPI_COMM_WORLD, inter_color, 0, comm_palm, ierr )
    112111       comm2d = comm_palm

  • palm/trunk/SOURCE/init_grid.f90

    r559 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
     7! Furthermore the allocation of arrays and steering of loops is done with these
     8! parameters. Call of exchange_horiz are modified.
     9! In case of dirichlet bounday condition at the bottom zu(0)=0.0
     10! dzu_mg has to be set explicitly for a equally spaced grid near bottom.
     11! ddzu_pres added to use a equally spaced grid near bottom.
    712!
    813! Former revisions:
    … …  
    7681
    7782    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  distance
    78 
     83   
     84!
     85!   Computation of the array bounds.
     86    nxlg = nxl - nbgp
     87    nxrg = nxr + nbgp
     88    nysg = nys - nbgp
     89    nyng = nyn + nbgp
    7990!
    8091!-- Allocate grid arrays
    8192    ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), &
    82               dzw(1:nzt+1), l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1) )
     93              dzw(1:nzt+1), l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) )
    8394
    8495!
    … …  
    97108!
    98109!--    Grid for atmosphere with surface at z=0 (k=0, w-grid).
    99 !--    Since the w-level lies on the surface, the first u-level (staggered!)
    100 !--    lies below the surface (used for "mirror" boundary condition).
    101110!--    The first u-level above the surface corresponds to the top of the
    102111!--    Prandtl-layer.
    103        zu(0) = - dz * 0.5
     112
     113       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
     114          zu(0) = 0.0
     115      !    zu(0) = - dz * 0.5
     116       ELSE
     117          zu(0) = - dz * 0.5
     118       ENDIF
    104119       zu(1) =   dz * 0.5
    105120
    … …  
    173188       dd2zu(k) = 1.0 / ( dzu(k) + dzu(k+1) )
    174189    ENDDO
     190   
     191!   
     192!-- In case of FFT method or SOR swap inverse grid lenght ddzu to ddzu_fft and
     193!-- modify the lowest entry. This is necessary for atmosphere runs, because
     194!-- zu(0) and so ddzu(1) changed. Accompanied with this modified arrays
     195!-- pressure solver uses wrong values and this causes kinks in the profiles
     196!-- of turbulent quantities. 
     197    IF ( psolver /= 'multigrid' ) THEN
     198       ALLOCATE( ddzu_pres(1:nzt+1) )
     199       ddzu_pres = ddzu
     200       IF( .NOT. ocean ) ddzu_pres(1) = ddzu_pres(2)
     201    ENDIF   
    175202
    176203!
    … …  
    187214
    188215       dzu_mg(:,maximum_grid_level) = dzu
     216!       
     217!--    To ensure a equally spaced grid. For ocean runs this is not necessary,
     218!--    because zu(0) does not changed so far. Also this would cause errors
     219!--    if a vertical stretching for ocean runs is used.
     220       IF ( .NOT. ocean ) dzu_mg(1,maximum_grid_level) = dzu(2)
    189221       dzw_mg(:,maximum_grid_level) = dzw
    190222       nzt_l = nzt
    … …  
    239271!-- the flag arrays needed for the multigrid method
    240272    gls = 2**( maximum_grid_level )
     273
    241274    ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr),       &
    242275              corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr),       &
    243               nzb_local(-gls:ny+gls,-gls:nx+gls), nzb_tmp(-1:ny+1,-1:nx+1), &
     276              nzb_local(-gls:ny+gls,-gls:nx+gls),                                   &
     277              nzb_tmp(-nbgp:ny+nbgp,-nbgp:nx+nbgp),                         &
    244278              wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr),             &
    245279              wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
    246     ALLOCATE( fwxm(nys-1:nyn+1,nxl-1:nxr+1), fwxp(nys-1:nyn+1,nxl-1:nxr+1), &
    247               fwym(nys-1:nyn+1,nxl-1:nxr+1), fwyp(nys-1:nyn+1,nxl-1:nxr+1), &
    248               fxm(nys-1:nyn+1,nxl-1:nxr+1), fxp(nys-1:nyn+1,nxl-1:nxr+1),   &
    249               fym(nys-1:nyn+1,nxl-1:nxr+1), fyp(nys-1:nyn+1,nxl-1:nxr+1),   &
    250               nzb_s_inner(nys-1:nyn+1,nxl-1:nxr+1),                         &
    251               nzb_s_outer(nys-1:nyn+1,nxl-1:nxr+1),                         &
    252               nzb_u_inner(nys-1:nyn+1,nxl-1:nxr+1),                         &
    253               nzb_u_outer(nys-1:nyn+1,nxl-1:nxr+1),                         &
    254               nzb_v_inner(nys-1:nyn+1,nxl-1:nxr+1),                         &
    255               nzb_v_outer(nys-1:nyn+1,nxl-1:nxr+1),                         &
    256               nzb_w_inner(nys-1:nyn+1,nxl-1:nxr+1),                         &
    257               nzb_w_outer(nys-1:nyn+1,nxl-1:nxr+1),                         &
    258               nzb_diff_s_inner(nys-1:nyn+1,nxl-1:nxr+1),                    &
    259               nzb_diff_s_outer(nys-1:nyn+1,nxl-1:nxr+1),                    &
    260               nzb_diff_u(nys-1:nyn+1,nxl-1:nxr+1),                          &
    261               nzb_diff_v(nys-1:nyn+1,nxl-1:nxr+1),                          &
    262               nzb_2d(nys-1:nyn+1,nxl-1:nxr+1),                              &
    263               wall_e_x(nys-1:nyn+1,nxl-1:nxr+1),                            &
    264               wall_e_y(nys-1:nyn+1,nxl-1:nxr+1),                            &
    265               wall_u(nys-1:nyn+1,nxl-1:nxr+1),                              &
    266               wall_v(nys-1:nyn+1,nxl-1:nxr+1),                              &
    267               wall_w_x(nys-1:nyn+1,nxl-1:nxr+1),                            &
    268               wall_w_y(nys-1:nyn+1,nxl-1:nxr+1) )
    269 
    270     ALLOCATE( l_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     280    ALLOCATE( fwxm(nysg:nyng,nxlg:nxrg), fwxp(nysg:nyng,nxlg:nxrg),         &
     281              fwym(nysg:nyng,nxlg:nxrg), fwyp(nysg:nyng,nxlg:nxrg),         &
     282              fxm(nysg:nyng,nxlg:nxrg), fxp(nysg:nyng,nxlg:nxrg),           &
     283              fym(nysg:nyng,nxlg:nxrg), fyp(nysg:nyng,nxlg:nxrg),           &
     284              nzb_s_inner(nysg:nyng,nxlg:nxrg),                             &
     285              nzb_s_outer(nysg:nyng,nxlg:nxrg),                             &
     286              nzb_u_inner(nysg:nyng,nxlg:nxrg),                             &
     287              nzb_u_outer(nysg:nyng,nxlg:nxrg),                             &
     288              nzb_v_inner(nysg:nyng,nxlg:nxrg),                             &
     289              nzb_v_outer(nysg:nyng,nxlg:nxrg),                             &
     290              nzb_w_inner(nysg:nyng,nxlg:nxrg),                             &
     291              nzb_w_outer(nysg:nyng,nxlg:nxrg),                             &
     292              nzb_diff_s_inner(nysg:nyng,nxlg:nxrg),                        &
     293              nzb_diff_s_outer(nysg:nyng,nxlg:nxrg),                        &
     294              nzb_diff_u(nysg:nyng,nxlg:nxrg),                              &
     295              nzb_diff_v(nysg:nyng,nxlg:nxrg),                              &
     296              nzb_2d(nysg:nyng,nxlg:nxrg),                                  &
     297              wall_e_x(nysg:nyng,nxlg:nxrg),                                &
     298              wall_e_y(nysg:nyng,nxlg:nxrg),                                &
     299              wall_u(nysg:nyng,nxlg:nxrg),                                  &
     300              wall_v(nysg:nyng,nxlg:nxrg),                                  &
     301              wall_w_x(nysg:nyng,nxlg:nxrg),                                &
     302              wall_w_y(nysg:nyng,nxlg:nxrg) )
     303
     304
     305
     306    ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    271307
    272308    nzb_s_inner = nzb;  nzb_s_outer = nzb
    … …  
    327363    vertical_influence(0) = vertical_influence(1)
    328364
    329     DO  i = nxl-1, nxr+1
    330        DO  j = nys-1, nyn+1
     365    DO  i = nxlg, nxrg
     366       DO  j = nysg, nyng
    331367          DO  k = nzb_s_inner(j,i) + 1, &
    332368                  nzb_s_inner(j,i) + vertical_influence(nzb_s_inner(j,i))
    … …  
    477513          nzb_local(:,-gls:-1)        = nzb_local(:,nx-gls+1:nx)
    478514          nzb_local(:,nx+1:nx+gls)    = nzb_local(:,0:gls-1)
     515
     516
    479517     
    480518          GOTO 12
    … …  
    588626!--    nzb_s_outer:
    589627!--    extend nzb_local east-/westwards first, then north-/southwards
    590        nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
     628       nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
    591629       DO  j = -1, ny + 1
    592630          DO  i = 0, nx
    … …  
    620658!--    nzb_u_inner:
    621659!--    extend nzb_local rightwards only
    622        nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
     660       nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
    623661       DO  j = -1, ny + 1
    624662          DO  i = 0, nx + 1
    … …  
    626664          ENDDO
    627665       ENDDO
    628        nzb_u_inner = nzb_tmp(nys-1:nyn+1,nxl-1:nxr+1)
     666       nzb_u_inner = nzb_tmp(nysg:nyng,nxlg:nxrg)
    629667
    630668!
    … …  
    652690!--    nzb_v_inner:
    653691!--    extend nzb_local northwards only
    654        nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
     692       nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
    655693       DO  i = -1, nx + 1
    656694          DO  j = 0, ny + 1
    … …  
    658696          ENDDO
    659697       ENDDO
    660        nzb_v_inner = nzb_tmp(nys-1:nyn+1,nxl-1:nxr+1)
     698       nzb_v_inner = nzb_tmp(nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp)
    661699
    662700!
    … …  
    10961134!
    10971135!-- Need to set lateral boundary conditions for l_wall
    1098     CALL exchange_horiz( l_wall )
     1136
     1137    CALL exchange_horiz( l_wall, nbgp )
    10991138
    11001139    DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, &

  • palm/trunk/SOURCE/init_particles.f90

    r623 r667  
    44! Current revisions:
    55! -----------------
    6 !
     6! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
     7! of arrays.
    78!
    89! Former revisions:
    … …  
    185186       ENDIF
    186187
    187        ALLOCATE( prt_count(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),       &
    188                  prt_start_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
     188       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
     189                 prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    189190                 particle_mask(maximum_number_of_particles),         &
    190191                 particles(maximum_number_of_particles) )
    … …  
    209210!--    particles, which can be also periodically released at later times.
    210211!--    Also allocate array for particle tail coordinates, if needed.
    211        ALLOCATE( prt_count(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),       &
    212                  prt_start_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
     212       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
     213                 prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    213214                 particle_mask(maximum_number_of_particles),         &
    214215                 particles(maximum_number_of_particles) )

  • palm/trunk/SOURCE/init_pegrid.f90

    r647 r667  
    44! Current revisions:
    55! -----------------
     6!
     7! Moved determination of target_id's from init_coupling
     8!
     9! Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
     10! with different grid/processor-topology in ocean and atmosphere
     11!
     12!
     13! Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
     14! The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
     15! maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
     16! This distinction is due to reasons of data exchange and performance for the
     17! normal grid and grids in poismg.
     18! The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
     19! New MPI-Vectors for data exchange between left and right boundaries added.
     20! This is due to reasons of performance (10% faster).
    621!
    722! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
    … …  
    7994
    8095
     96
    8197    IMPLICIT NONE
    8298
    … …  
    88104
    89105    INTEGER, DIMENSION(:), ALLOCATABLE ::  ind_all, nxlf, nxrf, nynf, nysf
     106
     107    INTEGER, DIMENSION(2) :: pdims_remote
    90108
    91109    LOGICAL ::  found
    … …  
    103121
    104122#if defined( __parallel )
     123
    105124!
    106125!-- Determine the processor topology or check it, if prescribed by the user
    … …  
    624643#endif
    625644
     645!
     646!-- Determine the number of ghost points
     647    IF (scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme') THEN
     648       nbgp = 3
     649    ELSE
     650       nbgp = 1
     651    END IF
     652
    626653!
    627654!-- In case of coupled runs, create a new MPI derived datatype for the
    628655!-- exchange of surface (xy) data .
    629656!-- Gridpoint number for the exchange of ghost points (xy-plane)
    630     ngp_xy  = ( nxr - nxl + 3 ) * ( nyn - nys + 3 )
     657
     658    ngp_xy  = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
    631659
    632660!
    … …  
    635663    CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
    636664    CALL MPI_TYPE_COMMIT( type_xy, ierr )
     665
     666
     667    IF ( TRIM( coupling_mode ) .NE. 'uncoupled' ) THEN
     668   
     669!
     670!--    Pass the number of grid points of the atmosphere model to
     671!--    the ocean model and vice versa
     672       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
     673
     674          nx_a = nx
     675          ny_a = ny
     676
     677          IF ( myid == 0 ) THEN
     678             CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, &
     679                            comm_inter, ierr )
     680             CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, &
     681                            comm_inter, ierr )
     682             CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, &
     683                            comm_inter, ierr )
     684             CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, &
     685                            comm_inter, status, ierr )
     686             CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, &
     687                            comm_inter, status, ierr )
     688             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &
     689                            comm_inter, status, ierr )
     690          ENDIF
     691
     692          CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr)
     693          CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr)
     694          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
     695       
     696       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
     697
     698          nx_o = nx
     699          ny_o = ny
     700
     701          IF ( myid == 0 ) THEN
     702             CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, &
     703                            comm_inter, status, ierr )
     704             CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, &
     705                            comm_inter, status, ierr )
     706             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, &
     707                            comm_inter, status, ierr )
     708             CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, &
     709                            comm_inter, ierr )
     710             CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, &
     711                            comm_inter, ierr )
     712             CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, &
     713                            comm_inter, ierr )
     714          ENDIF
     715
     716          CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
     717          CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr)
     718          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
     719
     720       ENDIF
     721 
     722       ngp_a = (nx_a+1+2*nbgp)*(ny_a+1+2*nbgp)
     723       ngp_o = (nx_o+1+2*nbgp)*(ny_o+1+2*nbgp)
     724
     725!
     726!--    determine if the horizontal grid and the number of PEs
     727!--    in ocean and atmosphere is same or not
     728!--    (different number of PEs still not implemented)
     729       IF ( nx_o == nx_a .AND. ny_o == ny_a .AND.  &
     730            pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
     731       THEN
     732          coupling_topology = 0
     733       ELSE
     734          coupling_topology = 1
     735       ENDIF
     736
     737!
     738!--    Determine the target PEs for the exchange between ocean and
     739!--    atmosphere (comm2d)
     740       IF ( coupling_topology == 0) THEN
     741          IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
     742             target_id = myid + numprocs
     743          ELSE
     744             target_id = myid
     745          ENDIF
     746
     747       ELSE
     748
     749!
     750!--       In case of nonequivalent topology in ocean and atmosphere only for
     751!--       PE0 in ocean and PE0 in atmosphere a target_id is needed, since
     752!--       data echxchange between ocean and atmosphere will be done only by
     753!--       those PEs.   
     754          IF ( myid == 0 ) THEN
     755             IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
     756                target_id = numprocs
     757             ELSE
     758                target_id = 0
     759             ENDIF
     760 print*, coupling_mode, myid, " -> ", target_id, "numprocs: ", numprocs
     761          ENDIF
     762       ENDIF
     763
     764    ENDIF
     765
     766
    637767#endif
    638768
    … …  
    854984    ELSE
    855985
    856        maximum_grid_level = 1
     986       maximum_grid_level = 0
    857987
    858988    ENDIF
    … …  
    863993!
    864994!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
    865     ngp_y  = nyn - nys + 1
     995    ngp_y  = nyn - nys + 1 + 2 * nbgp
    866996
    867997!
    868998!-- Define a new MPI derived datatype for the exchange of ghost points in
    869999!-- y-direction for 2D-arrays (line)
    870     CALL MPI_TYPE_VECTOR( nxr-nxl+3, 1, ngp_y+2, MPI_REAL, type_x, ierr )