Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 1785)
+++ /palm/trunk/SOURCE/Makefile	(revision 1786)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# rename calc_spectra to spectrum + modularization of spectrum
 # 
 # Former revisions:
@@ -231,5 +231,5 @@
 	advec_w_pw.f90 advec_w_up.f90 average_3d_data.f90 boundary_conds.f90 \
 	buoyancy.f90 calc_liquid_water_content.f90 calc_mean_profile.f90 \
-	calc_precipitation.f90 calc_radiation.f90 calc_spectra.f90 \
+	calc_precipitation.f90 calc_radiation.f90 \
 	check_for_restart.f90 check_open.f90 check_parameters.f90 \
 	close_file.f90 compute_vpt.f90 coriolis.f90 cpulog.f90 \
@@ -268,5 +268,5 @@
 	random_function.f90 random_gauss.f90 random_generator_parallel.f90 \
 	read_3d_binary.f90 read_var_list.f90 run_control.f90 \
-	set_slicer_attributes_dvrp.f90 singleton.f90 sor.f90 \
+	set_slicer_attributes_dvrp.f90 singleton.f90 sor.f90 spectrum.f90 \
 	subsidence.f90 sum_up_3d_data.f90 \
 	surface_coupler.f90 surface_layer_fluxes.f90 swap_timelevel.f90 temperton_fft.f90 \
@@ -336,5 +336,4 @@
 calc_precipitation.o: modules.o mod_kinds.o
 calc_radiation.o: modules.o mod_kinds.o
-calc_spectra.o: modules.o cpulog.o fft_xy.o mod_kinds.o
 check_for_restart.o: modules.o mod_kinds.o
 check_open.o: modules.o mod_kinds.o mod_particle_attributes.o netcdf_interface.o
@@ -355,5 +354,6 @@
 data_output_ptseries.o: modules.o cpulog.o mod_kinds.o \
    netcdf_interface.o mod_particle_attributes.o
-data_output_spectra.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o
+data_output_spectra.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o \
+   spectrum.o
 data_output_tseries.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o
 data_output_2d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\
@@ -378,5 +378,6 @@
 global_min_max.o: modules.o mod_kinds.o
 header.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o land_surface_model.o\
-          plant_canopy_model.o pmc_interface.o radiation_model.o subsidence.o 
+   plant_canopy_model.o pmc_interface.o radiation_model.o spectrum.o \
+   subsidence.o 
 impact_of_latent_heat.o: modules.o mod_kinds.o
 inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
@@ -440,8 +441,9 @@
 mod_kinds.o: mod_kinds.f90
 mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
-netcdf_interface.o: netcdf_interface.f90 modules.o mod_kinds.o land_surface_model.o
+netcdf_interface.o: netcdf_interface.f90 modules.o mod_kinds.o \
+   land_surface_model.o spectrum.o
 nudging.o: modules.o cpulog.o mod_kinds.o
-package_parin.o: modules.o mod_kinds.o land_surface_model.o\
-                 plant_canopy_model.o radiation_model.o
+package_parin.o: modules.o mod_kinds.o land_surface_model.o \
+   plant_canopy_model.o radiation_model.o spectrum.o
 palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o\
 	pmc_interface.o surface_layer_fluxes.o
@@ -483,4 +485,5 @@
 singleton.o: mod_kinds.o singleton.f90
 sor.o: modules.o mod_kinds.o
+spectrum.o: spectrum.f90 modules.o mod_kinds.o cpulog.o fft_xy.o
 subsidence.o: modules.o mod_kinds.o
 sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
@@ -495,5 +498,5 @@
 	ls_forcing.o mod_kinds.o nudging.o pmc_interface.o production_e.o \
 	prognostic_equations.o progress_bar.o radiation_model.o \
-	 user_actions.o surface_layer_fluxes.o
+	spectrum.o user_actions.o surface_layer_fluxes.o
 time_to_string.o: mod_kinds.o
 timestep.o: modules.o cpulog.o mod_kinds.o
@@ -528,5 +531,5 @@
 user_parin.o: modules.o mod_kinds.o user_module.o
 user_read_restart_data.o: modules.o mod_kinds.o user_module.o
-user_spectra.o: modules.o mod_kinds.o user_module.o
+user_spectra.o: modules.o mod_kinds.o spectrum.o user_module.o
 user_statistics.o: modules.o mod_kinds.o netcdf_interface.o user_module.o
 wall_fluxes.o: modules.o mod_kinds.o
Index: lm/trunk/SOURCE/calc_spectra.f90
===================================================================
--- /palm/trunk/SOURCE/calc_spectra.f90	(revision 1785)
+++ 	(revision )
@@ -1,583 +1,0 @@
-!> @file calc_spectra.f90
-!--------------------------------------------------------------------------------!
-! This file is part of PALM.
-!
-! PALM is free software: you can redistribute it and/or modify it under the terms
-! of the GNU General Public License as published by the Free Software Foundation,
-! either version 3 of the License, or (at your option) any later version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
-! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with
-! PALM. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 1997-2014 Leibniz Universitaet Hannover
-!--------------------------------------------------------------------------------!
-!
-! Current revisions:
-! -----------------
-! 
-! 
-! Former revisions:
-! -----------------
-! $Id$
-!
-! 1682 2015-10-07 23:56:08Z knoop
-! Code annotations made doxygen readable 
-!
-! 1575 2015-03-27 09:56:27Z raasch
-! adjustments for psolver-queries
-!
-! 1511 2014-12-16 15:54:16Z suehring
-! Bugfix concerning spectra normalization
-!
-! 1431 2014-07-15 14:47:17Z suehring
-! Wavenumber-integrated spectra coincide with respective variance. 
-! 
-! 1342 2014-03-26 17:04:47Z kanani
-! REAL constants defined as wp-kinds
-!
-! 1324 2014-03-21 09:13:16Z suehring
-! Bugfix: nzb_x, nzb_yd, nyn_x, nyn_x, nzt_x, nzt_yd belong to transpose_indices 
-!
-! 1320 2014-03-20 08:40:49Z raasch
-! ONLY-attribute added to USE-statements,
-! kind-parameters added to all INTEGER and REAL declaration statements,
-! kinds are defined in new module kinds,
-! revision history before 2012 removed,
-! comment fields (!:) to be used for variable explanations added to
-! all variable declaration statements
-!
-! 1318 2014-03-17 13:35:16Z raasch
-! module interfaces removed
-!
-! 1216 2013-08-26 09:31:42Z raasch
-! resorting of array moved to separate routine resort_for_zx,
-! one argument removed from the transpose_..d routines
-!
-! 1120 2013-04-05 15:11:35Z raasch
-! bugfix: calls of fft_x|y replaced by fft_x|y_1d
-!
-! 1036 2012-10-22 13:43:42Z raasch
-! code put under GPL (PALM 3.9)
-!
-! 1003 2012-09-14 14:35:53Z raasch
-! adjustment of array tend for cases with unequal subdomain sizes removed
-!
-! Revision 1.1  2001/01/05 15:08:07  raasch
-! Initial revision
-!
-!
-! Description:
-! ------------
-!> Calculate horizontal spectra along x and y. 
-!> ATTENTION: 1d-decomposition along y still needs improvement, because in that
-!>            case the gridpoint number along z still depends on the PE number
-!>            because transpose_xz has to be used (and possibly also
-!>            transpose_zyd needs modification).
-!------------------------------------------------------------------------------!
- SUBROUTINE calc_spectra
- 
-
-#if defined( __spectra )
-    USE arrays_3d,                                                             &
-        ONLY:  d, tend
-
-    USE control_parameters,                                                    &
-        ONLY:  average_count_sp, bc_lr_cyc, bc_ns_cyc, message_string, psolver
-
-    USE cpulog,                                                                &
-        ONLY:  cpu_log, log_point
-
-    USE fft_xy,                                                                &
-        ONLY:  fft_init
-
-    USE indices,                                                               &
-        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
-
-    USE kinds
-
-    USE pegrid
-
-    USE spectrum,                                                              &
-        ONLY:  data_output_sp, spectra_direction
-
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  m  !<
-    INTEGER(iwp) ::  pr !<
-
-
-    CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
-
-!
-!-- Initialize ffts
-    CALL fft_init
-
-!
-!-- Reallocate array d in required size
-    IF ( psolver(1:9) == 'multigrid' )  THEN
-       DEALLOCATE( d )
-       ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
-    ENDIF
-
-    m = 1
-    DO WHILE ( data_output_sp(m) /= ' '  .AND.  m <= 10 )
-!
-!--    Transposition from z --> x  ( y --> x in case of a 1d-decomposition
-!--    along x)
-       IF ( INDEX( spectra_direction(m), 'x' ) /= 0 )  THEN
-
-!
-!--       Calculation of spectra works for cyclic boundary conditions only
-          IF ( .NOT. bc_lr_cyc )  THEN
-
-             message_string = 'non-cyclic lateral boundaries along x do not'// &
-                              '& allow calculation of spectra along x'
-             CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
-          ENDIF
-
-          CALL preprocess_spectra( m, pr )
-
-#if defined( __parallel )
-          IF ( pdims(2) /= 1 )  THEN
-             CALL resort_for_zx( d, tend )
-             CALL transpose_zx( tend, d )
-          ELSE
-             CALL transpose_yxd( d, d )
-          ENDIF
-          CALL calc_spectra_x( d, pr, m )
-#else
-          message_string = 'sorry, calculation of spectra in non parallel ' // &
-                           'mode& is still not realized'
-          CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )      
-#endif
-
-       ENDIF
-
-!
-!--    Transposition from z --> y (d is rearranged only in case of a
-!--    1d-decomposition along x)
-       IF ( INDEX( spectra_direction(m), 'y' ) /= 0 )  THEN
-
-!
-!--       Calculation of spectra works for cyclic boundary conditions only
-          IF ( .NOT. bc_ns_cyc )  THEN
-             IF ( myid == 0 )  THEN
-                message_string = 'non-cyclic lateral boundaries along y do' // &
-                                 ' not & allow calculation of spectra along y' 
-                CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
-             ENDIF
-             CALL local_stop
-          ENDIF
-
-          CALL preprocess_spectra( m, pr )
-
-#if defined( __parallel )
-          CALL transpose_zyd( d, d )
-          CALL calc_spectra_y( d, pr, m )
-#else
-          message_string = 'sorry, calculation of spectra in non parallel' //  &
-                           'mode& is still not realized'
-          CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
-#endif
-
-       ENDIF
-
-!
-!--    Increase counter for next spectrum
-       m = m + 1
-          
-    ENDDO
-
-!
-!-- Increase counter for averaging process in routine plot_spectra
-    average_count_sp = average_count_sp + 1
-
-    CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
-
-#endif
- END SUBROUTINE calc_spectra
-
-
-#if defined( __spectra )
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> @todo Missing subroutine description.
-!------------------------------------------------------------------------------!
- SUBROUTINE preprocess_spectra( m, pr )
-
-    USE arrays_3d,                                                             &
-        ONLY:  d, pt, q, u, v, w
-
-    USE indices,                                                               &
-        ONLY:  ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt
-
-    USE kinds
-
-    USE pegrid
-
-    USE spectrum,                                                              &
-        ONLY:  data_output_sp
-
-    USE statistics,                                                            &
-        ONLY:  hom, var_d
-
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) :: i  !<
-    INTEGER(iwp) :: j  !<
-    INTEGER(iwp) :: k  !<
-    INTEGER(iwp) :: m  !<
-    INTEGER(iwp) :: pr !<
-
-    REAL(wp), DIMENSION(nzb:nzt+1) :: var_d_l
-
-    SELECT CASE ( TRIM( data_output_sp(m) ) )
-          
-    CASE ( 'u' )
-       pr = 1
-       d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
-       
-    CASE ( 'v' )
-       pr = 2
-       d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
-       
-    CASE ( 'w' )
-       pr = 3
-       d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
-       
-    CASE ( 'pt' )
-       pr = 4
-       d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
-       
-    CASE ( 'q' )
-       pr = 41
-       d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
-       
-    CASE DEFAULT
-!
-!--    The DEFAULT case is reached either if the parameter data_output_sp(m)
-!--    contains a wrong character string or if the user has coded a special
-!--    case in the user interface. There, the subroutine user_spectra 
-!--    checks which of these two conditions applies.
-       CALL user_spectra( 'preprocess', m, pr )
-          
-    END SELECT
-
-!
-!-- Subtract horizontal mean from the array, for which spectra have to be
-!-- calculated
-    var_d_l(:) = 0.0_wp
-    DO  i = nxl, nxr
-       DO  j = nys, nyn
-          DO  k = nzb+1, nzt
-             d(k,j,i)   = d(k,j,i) - hom(k,1,pr,0)
-             var_d_l(k) = var_d_l(k) + d(k,j,i) * d(k,j,i)
-          ENDDO
-       ENDDO
-    ENDDO
-!
-!-- Compute total variance from local variances
-    var_d(:) = 0.0_wp
-#if defined( __parallel ) 
-    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( var_d_l(0), var_d(0), nzt+1-nzb, MPI_REAL,   &
-                        MPI_SUM, comm2d, ierr )
-#else
-    var_d(:) = var_d_l(:)
-#endif
-    var_d(:) = var_d(:) / ngp_2dh(0)
-
- END SUBROUTINE preprocess_spectra
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> @todo Missing subroutine description.
-!------------------------------------------------------------------------------!
- SUBROUTINE calc_spectra_x( ddd, pr, m )
-
-    USE arrays_3d,                                                             &
-        ONLY:  
-
-    USE control_parameters,                                                    &
-        ONLY:  fft_method
-
-    USE fft_xy,                                                                &
-        ONLY:  fft_x_1d
-
-    USE grid_variables,                                                        &
-        ONLY:  dx
-
-    USE indices,                                                               &
-        ONLY:  nx, ny
-
-    USE kinds
-
-    USE pegrid
-
-    USE spectrum,                                                              &
-        ONLY:  comp_spectra_level, n_sp_x
-
-    USE statistics,                                                            &
-        ONLY:  spectrum_x, var_d
-
-    USE transpose_indices,                                                     &
-        ONLY:  nyn_x, nys_x, nzb_x, nzt_x
-
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  i         !<
-    INTEGER(iwp) ::  ishape(1) !<
-    INTEGER(iwp) ::  j         !<
-    INTEGER(iwp) ::  k         !<
-    INTEGER(iwp) ::  m         !<
-    INTEGER(iwp) ::  n         !<
-    INTEGER(iwp) ::  pr        !<
-
-    REAL(wp) ::  exponent     !<
-    REAL(wp) ::  sum_spec_dum !< wavenumber-integrated spectrum
-    
-    REAL(wp), DIMENSION(0:nx) ::  work !<
-    
-    REAL(wp), DIMENSION(0:nx/2) ::  sums_spectra_l !<
-    
-    REAL(wp), DIMENSION(0:nx/2,100) ::  sums_spectra !<
-    
-    REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ddd !<
-
-!
-!-- Exponent for geometric average
-    exponent = 1.0_wp / ( ny + 1.0_wp )
-
-!
-!-- Loop over all levels defined by the user
-    n = 1
-    DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
-
-       k = comp_spectra_level(n)
-
-!
-!--    Calculate FFT only if the corresponding level is situated on this PE
-       IF ( k >= nzb_x  .AND.  k <= nzt_x )  THEN
-          
-          DO  j = nys_x, nyn_x
-
-             work = ddd(0:nx,j,k)
-             CALL fft_x_1d( work, 'forward' )
-
-             ddd(0,j,k) = dx * work(0)**2
-             DO  i = 1, nx/2
-                ddd(i,j,k) = dx * ( work(i)**2 + work(nx+1-i)**2 )
-             ENDDO
-
-          ENDDO
-
-!
-!--       Local sum and geometric average of these spectra
-!--       (WARNING: no global sum should be performed, because floating
-!--       point overflow may occur)
-          DO  i = 0, nx/2
-
-             sums_spectra_l(i) = 1.0_wp
-             DO  j = nys_x, nyn_x
-                sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
-             ENDDO
-
-          ENDDO
-          
-       ELSE
-
-          sums_spectra_l = 1.0_wp
-
-       ENDIF
-
-!
-!--    Global sum of spectra on PE0 (from where they are written on file) 
-       sums_spectra(:,n) = 0.0_wp
-#if defined( __parallel )   
-       CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
-       CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1,          &
-                        MPI_REAL, MPI_PROD, 0, comm2d, ierr )
-#else
-       sums_spectra(:,n) = sums_spectra_l
-#endif
-!
-!--    Normalize spectra by variance
-       sum_spec_dum = SUM( sums_spectra(:,n) )
-       IF ( sum_spec_dum /= 0.0_wp )  THEN
-          sums_spectra(:,n) = sums_spectra(:,n) * var_d(k) / sum_spec_dum
-       ENDIF
-       n = n + 1
-
-    ENDDO
-    n = n - 1
-
-    IF ( myid == 0 )  THEN
-!
-!--    Sum of spectra for later averaging (see routine data_output_spectra)
-       DO  i = 1, nx/2
-          DO k = 1, n
-             spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k)
-          ENDDO
-       ENDDO
-
-    ENDIF
-!
-!-- n_sp_x is needed by data_output_spectra_x
-    n_sp_x = n
-
- END SUBROUTINE calc_spectra_x
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> @todo Missing subroutine description.
-!------------------------------------------------------------------------------!
- SUBROUTINE calc_spectra_y( ddd, pr, m )
-
-    USE arrays_3d,                                                             &
-        ONLY:  
-
-    USE control_parameters,                                                    &
-        ONLY:  fft_method
-
-    USE fft_xy,                                                                &
-        ONLY:  fft_y_1d
-
-    USE grid_variables,                                                        &
-        ONLY:  dy
-
-    USE indices,                                                               &
-        ONLY:  nx, ny
-
-    USE kinds
-
-    USE pegrid
-
-    USE spectrum,                                                              &
-        ONLY:  comp_spectra_level, n_sp_y
-
-    USE statistics,                                                            &
-        ONLY:  spectrum_y, var_d
-
-    USE transpose_indices,                                                     &
-        ONLY:  nxl_yd, nxr_yd, nzb_yd, nzt_yd
-
-
-    IMPLICIT NONE
-
-    INTEGER(iwp) ::  i         !<
-    INTEGER(iwp) ::  j         !<
-    INTEGER(iwp) ::  jshape(1) !<
-    INTEGER(iwp) ::  k         !<
-    INTEGER(iwp) ::  m         !<
-    INTEGER(iwp) ::  n         !<
-    INTEGER(iwp) ::  pr        !<
-
-    REAL(wp) ::  exponent !<
-    REAL(wp) ::  sum_spec_dum !< wavenumber-integrated spectrum
-    
-    REAL(wp), DIMENSION(0:ny) ::  work !<
-    
-    REAL(wp), DIMENSION(0:ny/2) ::  sums_spectra_l !<
-    
-    REAL(wp), DIMENSION(0:ny/2,100) ::  sums_spectra !<
-    
-    REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !<
-
-
-!
-!-- Exponent for geometric average
-    exponent = 1.0_wp / ( nx + 1.0_wp )
-
-!
-!-- Loop over all levels defined by the user
-    n = 1
-    DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
-
-       k = comp_spectra_level(n)
-
-!
-!--    Calculate FFT only if the corresponding level is situated on this PE
-       IF ( k >= nzb_yd  .AND.  k <= nzt_yd )  THEN
-          
-          DO  i = nxl_yd, nxr_yd
-
-             work = ddd(0:ny,i,k)
-             CALL fft_y_1d( work, 'forward' )
-
-             ddd(0,i,k) = dy * work(0)**2
-             DO  j = 1, ny/2
-                ddd(j,i,k) = dy * ( work(j)**2 + work(ny+1-j)**2 )
-             ENDDO
-
-          ENDDO
-
-!
-!--       Local sum and geometric average of these spectra
-!--       (WARNING: no global sum should be performed, because floating
-!--       point overflow may occur)
-          DO  j = 0, ny/2
-
-             sums_spectra_l(j) = 1.0_wp
-             DO  i = nxl_yd, nxr_yd
-                sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
-             ENDDO
-
-          ENDDO
-          
-       ELSE
-
-          sums_spectra_l = 1.0_wp
-
-       ENDIF
-
-!
-!--    Global sum of spectra on PE0 (from where they are written on file) 
-       sums_spectra(:,n) = 0.0_wp
-#if defined( __parallel )   
-       CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
-       CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1,          &
-                        MPI_REAL, MPI_PROD, 0, comm2d, ierr )
-#else
-       sums_spectra(:,n) = sums_spectra_l
-#endif
-!
-!--    Normalize spectra by variance
-       sum_spec_dum = SUM( sums_spectra(:,n) )
-       IF ( SUM(sums_spectra(:,n)) /= 0.0_wp )  THEN
-          sums_spectra(:,n) = sums_spectra(:,n) * var_d(k) / SUM(sums_spectra(:,n))
-       ENDIF
-       n = n + 1
-
-    ENDDO
-    n = n - 1
-
-
-    IF ( myid == 0 )  THEN
-!
-!--    Sum of spectra for later averaging (see routine data_output_spectra)
-       DO  j = 1, ny/2
-          DO k = 1, n
-             spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k)
-          ENDDO
-       ENDDO
-
-    ENDIF
-
-!
-!-- n_sp_y is needed by data_output_spectra_y
-    n_sp_y = n
-
- END SUBROUTINE calc_spectra_y
-#endif
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 1785)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 1786)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-!
+! cpp-direktives for spectra removed, check of spectra level removed
 !
 ! Former revisions:
@@ -326,5 +326,4 @@
     USE profil_parameter
     USE radiation_model_mod
-    USE spectrum
     USE statistics
     USE subsidence_mod
@@ -3752,18 +3751,4 @@
     END SELECT
 
-#if defined( __spectra )
-!
-!-- Check the number of spectra level to be output
-    i = 1
-    DO WHILE ( comp_spectra_level(i) /= 999999  .AND.  i <= 100 )
-       i = i + 1
-    ENDDO
-    i = i - 1
-    IF ( i == 0 )  THEN
-       WRITE( message_string, * )  'no spectra levels given'
-       CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
-    ENDIF
-#endif
-
 !
 !-- Check mask conditions
Index: /palm/trunk/SOURCE/data_output_spectra.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_spectra.f90	(revision 1785)
+++ /palm/trunk/SOURCE/data_output_spectra.f90	(revision 1786)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! cpp-directives for spectra removed, immediate return if no spectra levels are
+! given
 !
 ! Former revisions:
@@ -73,6 +74,4 @@
  
 #if defined( __netcdf )
-#if defined( __spectra )
-
     USE control_parameters,                                                    &
         ONLY:  average_count_sp, averaging_interval_sp, dosp_time_count,       &
@@ -93,5 +92,5 @@
 
     USE spectrum,                                                              &
-        ONLY:  data_output_sp, spectra_direction
+        ONLY:  comp_spectra_level, data_output_sp, spectra_direction
 
     USE statistics,                                                            &
@@ -110,4 +109,8 @@
 
     CALL cpu_log( log_point(31), 'data_output_spectra', 'start' )
+
+!
+!-- Check if user gave any levels for spectra to be calculated
+    IF ( comp_spectra_level(1) == 999999 )  RETURN
 
 !
@@ -206,5 +209,4 @@
 
 #endif
-#endif
  END SUBROUTINE data_output_spectra
 
@@ -299,5 +301,4 @@
 
 
-#if defined( __spectra )
 !------------------------------------------------------------------------------!
 ! Description:
@@ -714,3 +715,2 @@
 
  END SUBROUTINE data_output_spectra_y
-#endif
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 1785)
+++ /palm/trunk/SOURCE/header.f90	(revision 1786)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-!
+! cpp-direktives for spectra removed
 !
 ! Former revisions:
@@ -1733,5 +1733,4 @@
 #endif
 
-#if defined( __spectra )
 !
 !-- Spectra output
@@ -1753,5 +1752,4 @@
                           averaging_interval_sp, dt_averaging_input_pr
     ENDIF
-#endif
 
     WRITE ( io, 99 )
@@ -2234,5 +2232,4 @@
 367 FORMAT ('       Polygon reduction for topography: cluster_size = ', I1)
 #endif
-#if defined( __spectra )
 370 FORMAT ('    Spectra:')
 371 FORMAT ('       Output every ',F7.1,' s'/)
@@ -2253,5 +2250,4 @@
             '       Profiles for the time averaging are taken every ', &
                     F6.1,' s')
-#endif
 400 FORMAT (//' Physical quantities:'/ &
               ' -------------------'/)
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 1785)
+++ /palm/trunk/SOURCE/modules.f90	(revision 1786)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-!
+! module spectrum moved to new separate module
 !
 ! Former revisions:
@@ -1277,53 +1277,4 @@
  END MODULE profil_parameter
 
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Definition of quantities used for computing spectra.
-!------------------------------------------------------------------------------!
- MODULE spectrum
-
-    USE kinds
-
-    CHARACTER (LEN=6),  DIMENSION(1:5) ::  header_char = (/ 'PS(u) ', 'PS(v) ',&
-                                           'PS(w) ', 'PS(pt)', 'PS(q) ' /)
-    CHARACTER (LEN=2),  DIMENSION(10)  ::  spectra_direction = 'x'
-    CHARACTER (LEN=10), DIMENSION(10)  ::  data_output_sp  = ' '
-    CHARACTER (LEN=25), DIMENSION(1:5) ::  utext_char =                    &
-                                           (/ '-power spectrum of u     ', &
-                                              '-power spectrum of v     ', &
-                                              '-power spectrum of w     ', &
-                                              '-power spectrum of ^1185 ', &
-                                              '-power spectrum of q     ' /)
-    CHARACTER (LEN=39), DIMENSION(1:5) ::  ytext_char =                        &
-                                 (/ 'k ^2236 ^2566^2569<u(k) in m>2s>->2    ', &
-                                    'k ^2236 ^2566^2569<v(k) in m>2s>->2    ', &
-                                    'k ^2236 ^2566^2569<w(k) in m>2s>->2    ', &
-                                    'k ^2236 ^2566^2569<^1185(k) in m>2s>->2', &
-                                    'k ^2236 ^2566^2569<q(k) in m>2s>->2    ' /)
-
-    INTEGER(iwp) ::  klist_x = 0, klist_y = 0, n_sp_x = 0, n_sp_y = 0
-
-    INTEGER(iwp) ::  comp_spectra_level(100) = 999999,                   &
-                     lstyles(100) = (/ 0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
-                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
-                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
-                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
-                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
-                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
-                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
-                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
-                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
-                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8 /), &
-                     plot_spectra_level(100) = 999999
-
-    REAL(wp)    ::  time_to_start_sp = 0.0_wp
-
-    SAVE
-
- END MODULE spectrum
-
-
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1365,4 +1316,5 @@
 
 
+
 !------------------------------------------------------------------------------!
 ! Description:
Index: /palm/trunk/SOURCE/netcdf_interface.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface.f90	(revision 1785)
+++ /palm/trunk/SOURCE/netcdf_interface.f90	(revision 1786)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Bugfix: id_var_time_sp made public
 !
 ! Former revisions:
@@ -296,8 +296,9 @@
             id_var_dopr, id_var_dopts, id_var_dospx, id_var_dospy,             &
             id_var_dots, id_var_do2d, id_var_do3d, id_var_norm_dopr,           &
-            id_var_time_mask, id_var_time_pr, id_var_time_pts, id_var_time_ts, &
-            id_var_time_xy, id_var_time_xz, id_var_time_yz,id_var_time_3d,     &
-            nc_stat, netcdf_data_format, netcdf_data_format_string,            &
-            netcdf_deflate, netcdf_precision, output_for_t0
+            id_var_time_mask, id_var_time_pr, id_var_time_pts, id_var_time_sp, &
+            id_var_time_ts, id_var_time_xy, id_var_time_xz, id_var_time_yz,    &
+            id_var_time_3d, nc_stat, netcdf_data_format,                       &
+            netcdf_data_format_string, netcdf_deflate, netcdf_precision,       &
+            output_for_t0
 
     SAVE
Index: /palm/trunk/SOURCE/package_parin.f90
===================================================================
--- /palm/trunk/SOURCE/package_parin.f90	(revision 1785)
+++ /palm/trunk/SOURCE/package_parin.f90	(revision 1786)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! cpp-direktives for spectra removed
 ! 
 ! Former revisions:
@@ -318,6 +318,4 @@
  30 CONTINUE
 
-
-#if defined( __spectra )
     REWIND ( 11 )
     line = ' '
@@ -337,5 +335,4 @@
 
  40 CONTINUE
-#endif
 
 !
Index: /palm/trunk/SOURCE/pmc_client.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_client.f90	(revision 1785)
+++ /palm/trunk/SOURCE/pmc_client.f90	(revision 1786)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! change in client-server data transfer: server now gets data from client
+! instead that client put's it to the server
 !
 ! Former revisions:
@@ -361,5 +362,5 @@
        do i=1,me%inter_npes
           aPE => me%PEs(i)
-          ar  => aPE%array_list(next_array_in_list)    !actual array is last array in list
+          ar  => aPE%array_list(next_array_in_list)
           ar%NrDims    = NrDims
           ar%A_dim     = dims
@@ -407,12 +408,21 @@
       IMPLICIT none
 
-      INTEGER                                 :: i, ierr, j
+!--   naming convention:  appending       _sc  -> server to client transfer
+!--                                       _cs  -> client to server transfer
+!--                                       Recv -> server to client transfer
+!--                                       Send -> client to server transfer
+
+      INTEGER                                 :: i, istat, ierr, j
+      INTEGER,PARAMETER                       :: NoINdex=-1
+      INTEGER                                 :: rcount
       INTEGER                                 :: arlen, myIndex, tag
       INTEGER(idp)                            :: bufsize                   ! Size of MPI data Window
       TYPE(PeDef),POINTER                     :: aPE
       TYPE(ArrayDef),POINTER                  :: ar
+      INTEGER,DIMENSION(1024)                 :: req
       character(len=DA_Namelen)               :: myName
       Type(c_ptr)                             :: base_ptr
-      REAL(kind=wp),DIMENSION(:),POINTER      :: base_array
+      REAL(kind=wp),DIMENSION(:),POINTER,save :: base_array_sc             !Base array
+      REAL(kind=wp),DIMENSION(:),POINTER,save :: base_array_cs             !Base array
       INTEGER(KIND=MPI_ADDRESS_KIND)          :: WinSize
 
@@ -420,5 +430,7 @@
       bufsize = 8
 
-      ! First stride, Compute size and set index
+!--   Server to client direction
+
+!--   First stride, Compute size and set index
 
       do i=1,me%inter_npes
@@ -433,23 +445,22 @@
             CALL MPI_Recv (myIndex, 1, MPI_INTEGER, i-1, tag, me%inter_comm, MPI_STATUS_IGNORE, ierr)
 
-            if(ar%NrDims == 3) then                    ! PALM has k in first dimension
+            if(ar%NrDims == 3) then
                bufsize = max(bufsize,ar%A_dim(1)*ar%A_dim(2)*ar%A_dim(3))    ! determine max, because client buffer is allocated only once
             else
                bufsize = max(bufsize,ar%A_dim(1)*ar%A_dim(2))
             end if
-            ar%BufIndex = myIndex
-          end do
+            ar%RecvIndex = myIndex
+
+           end do
       end do
 
-      ! Create RMA (One Sided Communication) window for data buffer
-
-      CALL PMC_Alloc_mem (base_array, bufsize, base_ptr)
-      me%TotalBufferSize = bufsize*wp                          !Total buffer size in Byte
-
-      WinSize = me%TotalBufferSize
-!      write(9,'(a,8i7)') 'PMC_S_SetInd_and_Mem ',m_model_rank,me%inter_npes,WinSize,ar%A_dim
-      CALL MPI_Win_create (base_array, WinSize, wp, MPI_INFO_NULL, me%intra_comm, me%BufWin, ierr);
-      CALL MPI_Win_fence (0, me%BufWin, ierr);                    !  Open Window to set data
-      CALL MPI_Barrier(me%intra_comm, ierr)
+
+!--   Create RMA (One Sided Communication) data buffer
+!--   The buffer for MPI_Get can be PE local, i.e. it can but must not be part of the MPI RMA window
+
+      CALL PMC_Alloc_mem (base_array_sc, bufsize, base_ptr)
+      me%TotalBufferSize = bufsize*wp                          ! Total buffer size in Byte
+
+!--   Second stride, Set Buffer pointer
 
       do i=1,me%inter_npes
@@ -458,5 +469,78 @@
          do j=1,aPE%Nr_arrays
             ar  => aPE%array_list(j)
-            ar%SendBuf = base_ptr
+            ar%RecvBuf = base_ptr
+         end do
+      end do
+
+!--   Client to server direction
+
+      myIndex = 1
+      rCount  = 0
+      bufsize = 8
+
+      do i=1,me%inter_npes
+         aPE => me%PEs(i)
+         tag = 300
+         do j=1,aPE%Nr_arrays
+            ar  => aPE%array_list(j)
+            if(ar%NrDims == 2) then
+               arlen     = aPE%NrEle                        ! 2D
+            else if(ar%NrDims == 3) then
+               arlen     = aPE%NrEle*ar%A_dim(1)            ! 3D
+            end if
+
+            tag    = tag+1
+            rCount = rCount+1
+            if(aPE%NrEle > 0)  then
+               CALL MPI_Isend (myIndex, 1, MPI_INTEGER, i-1, tag, me%inter_comm, req(rCount),ierr)
+               ar%SendIndex = myIndex
+            else
+               CALL MPI_Isend (NoIndex, 1, MPI_INTEGER, i-1, tag, me%inter_comm, req(rCount),ierr)
+               ar%SendIndex = NoIndex
+            end if
+
+            if(rCount == 1024) then                                  ! Maximum of 1024 outstanding requests
+               CALL MPI_Waitall (rCount, req, MPI_STATUSES_IGNORE, ierr)
+               rCount = 0;
+            end if
+
+            if(aPE%NrEle > 0)  then
+               ar%SendSize  = arlen
+               myIndex     = myIndex+arlen
+               bufsize     = bufsize+arlen
+            end if
+          end do
+         if(rCount > 0) then                                          ! Wait for all send completed
+            CALL MPI_Waitall (rCount, req, MPI_STATUSES_IGNORE, ierr)
+         end if
+      end do
+
+!--   Create RMA (One Sided Communication) window for data buffer client to server transfer
+!--   The buffer of MPI_Get (counter part of transfer) can be PE-local, i.e. it can but must not be part of the MPI RMA window
+!--   Only one RMA window is required to prepare the data for server -> client transfer on the server side and
+!--                                                       for client -> server transfer on the client side
+
+      CALL PMC_Alloc_mem (base_array_cs, bufsize)
+      me%TotalBufferSize = bufsize*wp                          !Total buffer size in Byte
+
+      WinSize = me%TotalBufferSize
+      CALL MPI_Win_create (base_array_cs, WinSize, wp, MPI_INFO_NULL, me%intra_comm, me%win_server_client, ierr);
+      CALL MPI_Win_fence (0, me%win_server_client, ierr);                    !  Open Window to set data
+      CALL MPI_Barrier(me%intra_comm, ierr)
+
+!--   Second stride, Set Buffer pointer
+
+      do i=1,me%inter_npes
+         aPE => me%PEs(i)
+
+         do j=1,aPE%Nr_arrays
+            ar  => aPE%array_list(j)
+            if(aPE%NrEle > 0)  then
+              ar%SendBuf = c_loc(base_array_cs(ar%SendIndex))
+              if(ar%SendIndex+ar%SendSize > bufsize) then
+                 write(0,'(a,i4,4i7,1x,a)') 'Client Buffer too small ',i,ar%SendIndex,ar%SendSize,ar%SendIndex+ar%SendSize,bufsize,trim(ar%name)
+                 CALL MPI_Abort (MPI_COMM_WORLD, istat, ierr)
+              end if
+            end if
          end do
       end do
@@ -504,12 +588,12 @@
 
             buf_shape(1) = nr
-            CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
+            CALL c_f_pointer(ar%RecvBuf, buf, buf_shape)
 !
 !--         MPI passive target RMA
             if(nr > 0)   then
-               target_disp = (ar%BufIndex-1)
-               CALL MPI_Win_lock (MPI_LOCK_SHARED , ip-1, 0, me%BufWin, ierr)
-               CALL MPI_Get (buf, nr, MPI_REAL, ip-1, target_disp, nr, MPI_REAL, me%BufWin, ierr)
-               CALL MPI_Win_unlock (ip-1, me%BufWin, ierr)
+               target_disp = (ar%RecvIndex-1)
+               CALL MPI_Win_lock (MPI_LOCK_SHARED , ip-1, 0, me%win_server_client, ierr)
+               CALL MPI_Get (buf, nr, MPI_REAL, ip-1, target_disp, nr, MPI_REAL, me%win_server_client, ierr)
+               CALL MPI_Win_unlock (ip-1, me%win_server_client, ierr)
             end if
 
@@ -555,4 +639,8 @@
       INTEGER(kind=MPI_ADDRESS_KIND)    ::  target_disp
 
+      t1 = PMC_Time()
+      CALL MPI_Barrier(me%intra_comm, ierr)              ! Wait for empty buffer
+      t2 = PMC_Time()
+      if(present(WaitTime)) WaitTime = t2-t1
 
       do ip=1,me%inter_npes
@@ -561,15 +649,8 @@
          do j=1,aPE%Nr_arrays
             ar  => aPE%array_list(j)
+            myIndex=1
             if(ar%NrDims == 2) then
-               nr = aPE%NrEle
-            else if(ar%NrDims == 3) then
-               nr = aPE%NrEle*ar%A_dim(1)
-            end if
-
-            buf_shape(1) = nr
-            CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
-
-            myIndex = 1
-            if(ar%NrDims == 2) then
+               buf_shape(1) = aPE%NrEle
+               CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
                CALL c_f_pointer(ar%data, data_2d, ar%A_dim(1:2))
                do ij=1,aPE%NrEle
@@ -578,4 +659,6 @@
                end do
             else if(ar%NrDims == 3) then
+               buf_shape(1) = aPE%NrEle*ar%A_dim(1)
+               CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
                CALL c_f_pointer(ar%data, data_3d, ar%A_dim(1:3))
                do ij=1,aPE%NrEle
@@ -584,24 +667,10 @@
                end do
             end if
-!
-!--         MPI passiv target RMA
-            if(nr > 0)   then
-               target_disp = (ar%BufIndex-1)
-               CALL MPI_Win_lock (MPI_LOCK_SHARED , ip-1, 0, me%BufWin, ierr)
-               CALL MPI_Put (buf, nr, MPI_REAL, ip-1, target_disp, nr, MPI_REAL, me%BufWin, ierr)
-               CALL MPI_Win_unlock (ip-1, me%BufWin, ierr)
-            end if
-
           end do
       end do
 
 
-      t1 = PMC_Time()
-      CALL MPI_Barrier(me%intra_comm, ierr)                         ! Wait for server to fill buffer
-      t2 = PMC_Time()
-      if(present(WaitTime)) WaitTime = t2-t1
-
-      CALL MPI_Barrier(me%intra_comm, ierr)                         ! Wait for buffer is filled
-
+!      CALL MPI_Win_fence (0, me%win_server_client, ierr)      ! Fence might do it, test later
+      CALL MPI_Barrier(me%intra_comm, ierr)                   ! buffer is filled
 
       return
Index: /palm/trunk/SOURCE/pmc_general.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_general.f90	(revision 1785)
+++ /palm/trunk/SOURCE/pmc_general.f90	(revision 1786)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! change in client-server data transfer: server now gets data from client
+! instead that client put's it to the server
 !
 ! Former revisions:
@@ -82,7 +83,10 @@
       TYPE(c_ptr)                   :: data                        ! Pointer of data in server space
       TYPE(c_ptr), DIMENSION(2)     :: po_data                     ! Base Pointers, PMC_S_Set_Active_data_array sets active pointer
-      INTEGER(idp)                  :: BufIndex                    ! index in Send Buffer
-      INTEGER                       :: BufSize                     ! size in Send Buffer
+      INTEGER(idp)                  :: SendIndex                   ! index in Send Buffer
+      INTEGER(idp)                  :: RecvIndex                   ! index in Receive Buffer
+      INTEGER                       :: SendSize                    ! size in Send Buffer
+      INTEGER                       :: RecvSize                    ! size in Receiver Buffer
       TYPE (c_ptr)                  :: SendBuf                     ! Pointer of Data in Send buffer
+      TYPE (c_ptr)                  :: RecvBuf                     ! Pointer of Data in ReceiveS buffer
       CHARACTER(len=8)              :: Name                        ! Name of Array
 
@@ -108,5 +112,5 @@
       INTEGER                       :: inter_npes                  ! Number of PEs client model
       INTEGER                       :: intra_rank                  ! rank within intra_comm
-      INTEGER                       :: BufWin                      ! MPI RMA windows
+      INTEGER                       :: win_server_client           ! MPI RMA for preparing data on server AND client side
       TYPE (PeDef), DIMENSION(:), POINTER      :: PEs              ! List of all Client PEs
    END TYPE ClientDef
Index: /palm/trunk/SOURCE/pmc_handle_communicator.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_handle_communicator.f90	(revision 1785)
+++ /palm/trunk/SOURCE/pmc_handle_communicator.f90	(revision 1786)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-!
+! Bugfix: nesting_mode is broadcast now
 !
 ! Former revisions:
@@ -211,4 +211,5 @@
          CALL MPI_BCAST( m_couplers(i)%lower_left_x, 1, MPI_REAL,    0, MPI_COMM_WORLD, istat )
          CALL MPI_BCAST( m_couplers(i)%lower_left_y, 1, MPI_REAL,    0, MPI_COMM_WORLD, istat )
+         CALL MPI_BCAST( nesting_mode, LEN( nesting_mode ), MPI_CHARACTER, 0, MPI_COMM_WORLD, istat )
       ENDDO
 
Index: /palm/trunk/SOURCE/pmc_server.f90
===================================================================
--- /palm/trunk/SOURCE/pmc_server.f90	(revision 1785)
+++ /palm/trunk/SOURCE/pmc_server.f90	(revision 1786)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! change in client-server data transfer: server now gets data from client
+! instead that client put's it to the server
 !
 ! Former revisions:
@@ -336,4 +337,9 @@
       IMPLICIT none
 
+!
+!--   Naming convention for appendices:   _sc  -> server to client transfer
+!--                                       _cs  -> client to server transfer
+!--                                       Send -> Server to client transfer
+!--                                       Recv -> client to server transfer
       INTEGER,INTENT(IN)                      :: ClientId
 
@@ -347,13 +353,16 @@
       INTEGER,DIMENSION(1024)                 :: req
       Type(c_ptr)                             :: base_ptr
-      REAL(kind=wp),DIMENSION(:),POINTER      :: base_array
+      REAL(wp),DIMENSION(:),POINTER,SAVE :: base_array_sc  !Base array for server to client transfer
+      REAL(wp),DIMENSION(:),POINTER,SAVE :: base_array_cs  !Base array for client to server transfer
       INTEGER(KIND=MPI_ADDRESS_KIND)          :: WinSize
 
+!
+!--   Server to client direction
       myIndex = 1
       rCount  = 0
       bufsize = 8
 
-      ! First stride, Compute size and set index
-
+!
+!--   First stride: Compute size and set index
       do i=1,Clients(ClientId)%inter_npes
          aPE => Clients(ClientId)%PEs(i)
@@ -362,11 +371,11 @@
             ar  => aPE%array_list(j)
             if(ar%NrDims == 2) then
-               arlen     = aPE%NrEle;                             ! 2D
+               arlen     = aPE%NrEle                              ! 2D
             else if(ar%NrDims == 3) then
-               arlen     = aPE%NrEle * ar%A_dim(4);               ! PALM 3D
+               arlen     = aPE%NrEle * ar%A_dim(4);               ! 3D
             else
                arlen     = -1
             end if
-            ar%BufIndex = myIndex
+            ar%SendIndex = myIndex
 
             tag    = tag+1
@@ -381,36 +390,84 @@
             myIndex = myIndex+arlen
             bufsize = bufsize+arlen
-            ar%BufSize = arlen
+            ar%SendSize = arlen
 
          end do
-         if(rCount > 0) then                                          ! Wait for all send completed
+         if(rCount > 0) then                       ! Wait for all send completed
             CALL MPI_Waitall (rCount, req, MPI_STATUSES_IGNORE, ierr)
          end if
       end do
 
-      ! Create RMA (One Sided Communication) window for data buffer
-
-      CALL PMC_Alloc_mem (base_array, bufsize, base_ptr)
-      Clients(ClientId)%TotalBufferSize = bufsize*wp       !Total buffer size in Byte
+!
+!--   Create RMA (One Sided Communication) window for data buffer server to
+!--   client transfer.
+!--   The buffer of MPI_Get (counter part of transfer) can be PE-local, i.e.
+!--   it can but must not be part of the MPI RMA window.
+!--   Only one RMA window is required to prepare the data
+!--   for server -> client transfer on the server side and
+!--   for client -> server transfer on the client side
+      CALL PMC_Alloc_mem (base_array_sc, bufsize)
+      Clients(ClientId)%TotalBufferSize = bufsize*wp   !Total buffer size in Byte
 
       WinSize = bufsize*wp
-!      write(9,*) 'PMC_S_SetInd_and_Mem ',m_model_rank,Clients(ClientId)%inter_npes,WinSize,bufsize
-      CALL MPI_Win_create (base_array, WinSize, wp, MPI_INFO_NULL, Clients(ClientId)%intra_comm, Clients(ClientId)%BufWin, ierr);
-      CALL MPI_Win_fence (0, Clients(ClientId)%BufWin, ierr);                    !  Open Window to set data
-      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)
-
-      ! Second stride, Set Buffer pointer
-
+      CALL MPI_Win_create (base_array_sc, WinSize, wp, MPI_INFO_NULL,          &
+                  Clients(ClientId)%intra_comm,  Clients(ClientId)%win_server_client, ierr)
+      CALL MPI_Win_fence (0, Clients(ClientId)%win_server_client, ierr);        !  Open Window to set data
+!
+!--   Second stride: Set Buffer pointer
       do i=1,Clients(ClientId)%inter_npes
          aPE => Clients(ClientId)%PEs(i)
          do j=1,aPE%Nr_arrays
             ar  => aPE%array_list(j)
-!--         TO_DO:  Adressrechnung ueberlegen?
-            ar%SendBuf = c_loc(base_array(ar%BufIndex))                         !kk Adressrechnung ueberlegen
-            if(ar%BufIndex+ar%BufSize > bufsize) then
-!--            TO_DO: can this error really happen, and what can be the reason?
-               write(0,'(a,i4,4i7,1x,a)') 'Buffer too small ',i,ar%BufIndex,ar%BufSize,ar%BufIndex+ar%BufSize,bufsize,trim(ar%name)
+            ar%SendBuf = c_loc(base_array_sc(ar%SendIndex))
+            if(ar%SendIndex+ar%SendSize > bufsize) then
+               write(0,'(a,i4,4i7,1x,a)') 'Server Buffer too small ',i,ar%SendIndex,ar%SendSize,ar%SendIndex+ar%SendSize,bufsize,trim(ar%name)
                CALL MPI_Abort (MPI_COMM_WORLD, istat, ierr)
             end if
+         end do
+      end do
+
+!--   Client to server direction
+
+      bufsize  = 8
+
+!--   First stride, Compute size and set index
+
+      do i=1,Clients(ClientId)%inter_npes
+         aPE => Clients(ClientId)%PEs(i)
+         tag = 300
+
+         do j=1,aPE%Nr_arrays
+            ar  => aPE%array_list(j)
+
+            ! Receive Index from client
+            tag = tag+1
+            CALL MPI_Recv (myIndex, 1, MPI_INTEGER, i-1, tag, Clients(ClientId)%inter_comm, MPI_STATUS_IGNORE, ierr)
+
+            if(ar%NrDims == 3) then
+               bufsize = max(bufsize,aPE%NrEle * ar%A_dim(4))               ! 3D
+            else
+               bufsize = max(bufsize,aPE%NrEle)                             ! 2D
+            end if
+            ar%RecvIndex = myIndex
+          end do
+
+      end do
+
+!--   Create RMA (One Sided Communication) data buffer
+!--   The buffer for MPI_Get can be PE local, i.e. it can but must not be part of the MPI RMA window
+
+      CALL PMC_Alloc_mem (base_array_cs, bufsize, base_ptr)
+      Clients(ClientId)%TotalBufferSize = bufsize*wp       !Total buffer size in Byte
+
+      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)
+
+!--   Second stride, Set Buffer pointer
+
+      do i=1,Clients(ClientId)%inter_npes
+         aPE => Clients(ClientId)%PEs(i)
+
+         do j=1,aPE%Nr_arrays
+            ar  => aPE%array_list(j)
+            ar%RecvBuf = base_ptr
          end do
       end do
@@ -465,5 +522,5 @@
       end do
 
-      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)    ! buffer is full
+      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)    ! buffer is filled
 
       return
@@ -480,4 +537,5 @@
       INTEGER                             ::  ip,ij,istat,ierr,j
       INTEGER                             ::  myIndex
+      INTEGER                             ::  nr
       REAL(wp)                            ::  t1,t2
       TYPE(PeDef),POINTER                 ::  aPE
@@ -488,9 +546,14 @@
       REAL(wp), POINTER, DIMENSION(:,:)   ::  data_2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  data_3d
+      INTEGER                             ::  target_pe
+      INTEGER(kind=MPI_ADDRESS_KIND)      ::  target_disp
 
       t1 = PMC_Time()
-      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)              ! Wait for MPI_Put from client
-      t2 = PMC_Time()
-      if(present(WaitTime)) WaitTime = t2-t1
+      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)                         ! Wait for client to fill buffer
+      t2 = PMC_Time()-t1
+      if(present(WaitTime)) WaitTime = t2
+
+!      CALL MPI_Win_fence (0, Clients(ClientId)%win_server_client, ierr)            ! Fence might do it, test later
+      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)                         ! Wait for buffer is filled
 
       do ip=1,Clients(ClientId)%inter_npes
@@ -498,8 +561,28 @@
          do j=1,aPE%Nr_arrays
             ar  => aPE%array_list(j)
+
+            if(ar%RecvIndex < 0)  CYCLE
+
+            if(ar%NrDims == 2) then
+               nr = aPE%NrEle
+            else if(ar%NrDims == 3) then
+               nr = aPE%NrEle*ar%A_dim(4)
+            end if
+
+            buf_shape(1) = nr
+            CALL c_f_pointer(ar%RecvBuf, buf, buf_shape)
+!
+!--         MPI passive target RMA
+
+            if(nr > 0)   then
+               target_disp = ar%RecvIndex-1
+               target_pe = ip-1+m_model_npes                         ! client PEs are located behind server PEs
+               CALL MPI_Win_lock (MPI_LOCK_SHARED , target_pe, 0, Clients(ClientId)%win_server_client, ierr)
+               CALL MPI_Get (buf, nr, MPI_REAL, target_pe, target_disp, nr, MPI_REAL, Clients(ClientId)%win_server_client, ierr)
+               CALL MPI_Win_unlock (target_pe, Clients(ClientId)%win_server_client, ierr)
+            end if
+
             myIndex=1
             if(ar%NrDims == 2) then
-               buf_shape(1) = aPE%NrEle
-               CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
                CALL c_f_pointer(ar%data, data_2d, ar%A_dim(1:2))
                do ij=1,aPE%NrEle
@@ -508,6 +591,4 @@
                end do
             else if(ar%NrDims == 3) then
-               buf_shape(1) = aPE%NrEle*ar%A_dim(4)
-               CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
                CALL c_f_pointer(ar%data, data_3d, ar%A_dim(1:3))
                do ij=1,aPE%NrEle
@@ -519,7 +600,4 @@
       end do
 
-      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)              ! data copy finished, buffer is free for use agein
-
-        return
    END SUBROUTINE PMC_S_GetData_from_Buffer
 
Index: /palm/trunk/SOURCE/spectrum.f90
===================================================================
--- /palm/trunk/SOURCE/spectrum.f90	(revision 1786)
+++ /palm/trunk/SOURCE/spectrum.f90	(revision 1786)
@@ -0,0 +1,656 @@
+!> @file spectrum.f90
+!--------------------------------------------------------------------------------!
+! This file is part of PALM.
+!
+! PALM is free software: you can redistribute it and/or modify it under the terms
+! of the GNU General Public License as published by the Free Software Foundation,
+! either version 3 of the License, or (at your option) any later version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2014 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! routine is modularized, filename renamed from calc_spectra to spectrum,
+! privious data module spectrum moved from modules.f90 to here,
+! cpp-direktives for spectra removed, immediate return if no spectra levels are
+! given
+! 
+! Former revisions:
+! -----------------
+! $Id$
+!
+! 1682 2015-10-07 23:56:08Z knoop
+! Code annotations made doxygen readable 
+!
+! 1575 2015-03-27 09:56:27Z raasch
+! adjustments for psolver-queries
+!
+! 1511 2014-12-16 15:54:16Z suehring
+! Bugfix concerning spectra normalization
+!
+! 1431 2014-07-15 14:47:17Z suehring
+! Wavenumber-integrated spectra coincide with respective variance. 
+! 
+! 1342 2014-03-26 17:04:47Z kanani
+! REAL constants defined as wp-kinds
+!
+! 1324 2014-03-21 09:13:16Z suehring
+! Bugfix: nzb_x, nzb_yd, nyn_x, nyn_x, nzt_x, nzt_yd belong to transpose_indices 
+!
+! 1320 2014-03-20 08:40:49Z raasch
+! ONLY-attribute added to USE-statements,
+! kind-parameters added to all INTEGER and REAL declaration statements,
+! kinds are defined in new module kinds,
+! revision history before 2012 removed,
+! comment fields (!:) to be used for variable explanations added to
+! all variable declaration statements
+!
+! 1318 2014-03-17 13:35:16Z raasch
+! module interfaces removed
+!
+! 1216 2013-08-26 09:31:42Z raasch
+! resorting of array moved to separate routine resort_for_zx,
+! one argument removed from the transpose_..d routines
+!
+! 1120 2013-04-05 15:11:35Z raasch
+! bugfix: calls of fft_x|y replaced by fft_x|y_1d
+!
+! 1036 2012-10-22 13:43:42Z raasch
+! code put under GPL (PALM 3.9)
+!
+! 1003 2012-09-14 14:35:53Z raasch
+! adjustment of array tend for cases with unequal subdomain sizes removed
+!
+! Revision 1.1  2001/01/05 15:08:07  raasch
+! Initial revision
+!
+!
+! Description:
+! ------------
+!> Calculate horizontal spectra along x and y. 
+!> ATTENTION: 1d-decomposition along y still needs improvement, because in that
+!>            case the gridpoint number along z still depends on the PE number
+!>            because transpose_xz has to be used (and possibly also
+!>            transpose_zyd needs modification).
+!------------------------------------------------------------------------------!
+ MODULE spectrum
+
+    USE kinds
+
+    PRIVATE
+
+    CHARACTER (LEN=6),  DIMENSION(1:5) ::  header_char = (/ 'PS(u) ', 'PS(v) ',&
+                                           'PS(w) ', 'PS(pt)', 'PS(q) ' /)
+    CHARACTER (LEN=2),  DIMENSION(10)  ::  spectra_direction = 'x'
+    CHARACTER (LEN=10), DIMENSION(10)  ::  data_output_sp  = ' '
+    CHARACTER (LEN=25), DIMENSION(1:5) ::  utext_char =                    &
+                                           (/ '-power spectrum of u     ', &
+                                              '-power spectrum of v     ', &
+                                              '-power spectrum of w     ', &
+                                              '-power spectrum of ^1185 ', &
+                                              '-power spectrum of q     ' /)
+    CHARACTER (LEN=39), DIMENSION(1:5) ::  ytext_char =                        &
+                                 (/ 'k ^2236 ^2566^2569<u(k) in m>2s>->2    ', &
+                                    'k ^2236 ^2566^2569<v(k) in m>2s>->2    ', &
+                                    'k ^2236 ^2566^2569<w(k) in m>2s>->2    ', &
+                                    'k ^2236 ^2566^2569<^1185(k) in m>2s>->2', &
+                                    'k ^2236 ^2566^2569<q(k) in m>2s>->2    ' /)
+
+    INTEGER(iwp) ::  klist_x = 0, klist_y = 0, n_sp_x = 0, n_sp_y = 0
+
+    INTEGER(iwp) ::  comp_spectra_level(100) = 999999,                   &
+                     lstyles(100) = (/ 0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
+                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
+                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
+                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
+                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
+                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
+                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
+                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
+                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
+                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8 /), &
+                     plot_spectra_level(100) = 999999
+
+    REAL(wp)    ::  time_to_start_sp = 0.0_wp
+
+    PUBLIC  comp_spectra_level, data_output_sp, header_char, klist_x, klist_y, &
+            lstyles, n_sp_x, n_sp_y, plot_spectra_level, spectra_direction,    &
+            utext_char, ytext_char
+
+    SAVE
+
+    INTERFACE calc_spectra
+       MODULE PROCEDURE calc_spectra
+    END INTERFACE calc_spectra
+
+    INTERFACE preprocess_spectra
+       MODULE PROCEDURE preprocess_spectra
+    END INTERFACE preprocess_spectra
+
+    INTERFACE calc_spectra_x
+       MODULE PROCEDURE calc_spectra_x
+    END INTERFACE calc_spectra_x
+
+    INTERFACE calc_spectra_y
+       MODULE PROCEDURE calc_spectra_y
+    END INTERFACE calc_spectra_y
+
+    PUBLIC calc_spectra
+
+
+ CONTAINS
+
+    SUBROUTINE calc_spectra
+
+       USE arrays_3d,                                                          &
+           ONLY:  d, tend
+
+       USE control_parameters,                                                 &
+           ONLY:  average_count_sp, bc_lr_cyc, bc_ns_cyc, message_string, psolver
+
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point
+
+       USE fft_xy,                                                             &
+           ONLY:  fft_init
+
+       USE indices,                                                            &
+           ONLY:  nxl, nxr, nyn, nys, nzb, nzt
+
+       USE kinds
+
+       USE pegrid,                                                             &
+           ONLY:  myid, pdims
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  m  !<
+       INTEGER(iwp) ::  pr !<
+
+
+!
+!--    Check if user gave any levels for spectra to be calculated
+       IF ( comp_spectra_level(1) == 999999 )  RETURN
+
+       CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
+
+!
+!--    Initialize ffts
+       CALL fft_init
+
+!
+!--    Reallocate array d in required size
+       IF ( psolver(1:9) == 'multigrid' )  THEN
+          DEALLOCATE( d )
+          ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
+       ENDIF
+
+       m = 1
+       DO WHILE ( data_output_sp(m) /= ' '  .AND.  m <= 10 )
+!
+!--       Transposition from z --> x  ( y --> x in case of a 1d-decomposition
+!--       along x)
+          IF ( INDEX( spectra_direction(m), 'x' ) /= 0 )  THEN
+
+!
+!--          Calculation of spectra works for cyclic boundary conditions only
+             IF ( .NOT. bc_lr_cyc )  THEN
+
+                message_string = 'non-cyclic lateral boundaries along x do'//  &
+                                 ' not & allow calculation of spectra along x'
+                CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
+             ENDIF
+
+             CALL preprocess_spectra( m, pr )
+
+#if defined( __parallel )
+             IF ( pdims(2) /= 1 )  THEN
+                CALL resort_for_zx( d, tend )
+                CALL transpose_zx( tend, d )
+             ELSE
+                CALL transpose_yxd( d, d )
+             ENDIF
+             CALL calc_spectra_x( d, pr, m )
+#else
+             message_string = 'sorry, calculation of spectra in non paral' //  &
+                              'lel mode& is still not realized'
+             CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
+#endif
+
+          ENDIF
+
+!
+!--       Transposition from z --> y (d is rearranged only in case of a
+!--       1d-decomposition along x)
+          IF ( INDEX( spectra_direction(m), 'y' ) /= 0 )  THEN
+
+!
+!--          Calculation of spectra works for cyclic boundary conditions only
+             IF ( .NOT. bc_ns_cyc )  THEN
+                IF ( myid == 0 )  THEN
+                   message_string = 'non-cyclic lateral boundaries along y' // &
+                                    ' do not & allow calculation of spectr' // &
+                                    'a along y'
+                   CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
+                ENDIF
+                CALL local_stop
+             ENDIF
+
+             CALL preprocess_spectra( m, pr )
+
+#if defined( __parallel )
+             CALL transpose_zyd( d, d )
+             CALL calc_spectra_y( d, pr, m )
+#else
+             message_string = 'sorry, calculation of spectra in non paral' //  &
+                              'lel mode& is still not realized'
+             CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
+#endif
+
+          ENDIF
+
+!
+!--       Increase counter for next spectrum
+          m = m + 1
+          
+       ENDDO
+
+!
+!--    Increase counter for averaging process in routine plot_spectra
+       average_count_sp = average_count_sp + 1
+
+       CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
+
+    END SUBROUTINE calc_spectra
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> @todo Missing subroutine description.
+!------------------------------------------------------------------------------!
+    SUBROUTINE preprocess_spectra( m, pr )
+
+       USE arrays_3d,                                                          &
+           ONLY:  d, pt, q, u, v, w
+
+       USE indices,                                                            &
+           ONLY:  ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt
+
+       USE kinds
+
+#if defined( __lc )
+       USE MPI
+#else
+       INCLUDE "mpif.h"
+#endif
+       USE pegrid,                                                             &
+           ONLY:  collective_wait, comm2d, ierr
+
+       USE statistics,                                                         &
+           ONLY:  hom, var_d
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) :: i  !<
+       INTEGER(iwp) :: j  !<
+       INTEGER(iwp) :: k  !<
+       INTEGER(iwp) :: m  !<
+       INTEGER(iwp) :: pr !<
+
+       REAL(wp), DIMENSION(nzb:nzt+1) :: var_d_l
+
+       SELECT CASE ( TRIM( data_output_sp(m) ) )
+          
+       CASE ( 'u' )
+          pr = 1
+          d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
+       
+       CASE ( 'v' )
+          pr = 2
+          d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
+       
+       CASE ( 'w' )
+          pr = 3
+          d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
+       
+       CASE ( 'pt' )
+          pr = 4
+          d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
+       
+       CASE ( 'q' )
+          pr = 41
+          d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
+       
+       CASE DEFAULT
+!
+!--       The DEFAULT case is reached either if the parameter data_output_sp(m)
+!--       contains a wrong character string or if the user has coded a special
+!--       case in the user interface. There, the subroutine user_spectra
+!--       checks which of these two conditions applies.
+          CALL user_spectra( 'preprocess', m, pr )
+          
+       END SELECT
+
+!
+!--    Subtract horizontal mean from the array, for which spectra have to be
+!--    calculated
+       var_d_l(:) = 0.0_wp
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb+1, nzt
+                d(k,j,i)   = d(k,j,i) - hom(k,1,pr,0)
+                var_d_l(k) = var_d_l(k) + d(k,j,i) * d(k,j,i)
+             ENDDO
+          ENDDO
+       ENDDO
+!
+!--    Compute total variance from local variances
+       var_d(:) = 0.0_wp
+#if defined( __parallel ) 
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( var_d_l(0), var_d(0), nzt+1-nzb, MPI_REAL, MPI_SUM, &
+                           comm2d, ierr )
+#else
+       var_d(:) = var_d_l(:)
+#endif
+       var_d(:) = var_d(:) / ngp_2dh(0)
+
+    END SUBROUTINE preprocess_spectra
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> @todo Missing subroutine description.
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_spectra_x( ddd, pr, m )
+
+       USE control_parameters,                                                 &
+           ONLY:  fft_method
+
+       USE fft_xy,                                                             &
+           ONLY:  fft_x_1d
+
+       USE grid_variables,                                                     &
+           ONLY:  dx
+
+       USE indices,                                                            &
+           ONLY:  nx, ny
+
+       USE kinds
+
+#if defined( __lc )
+       USE MPI
+#else
+       INCLUDE "mpif.h"
+#endif
+       USE pegrid,                                                             &
+           ONLY:  comm2d, ierr, myid
+
+       USE statistics,                                                         &
+           ONLY:  spectrum_x, var_d
+
+       USE transpose_indices,                                                  &
+           ONLY:  nyn_x, nys_x, nzb_x, nzt_x
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i         !<
+       INTEGER(iwp) ::  ishape(1) !<
+       INTEGER(iwp) ::  j         !<
+       INTEGER(iwp) ::  k         !<
+       INTEGER(iwp) ::  m         !<
+       INTEGER(iwp) ::  n         !<
+       INTEGER(iwp) ::  pr        !<
+
+       REAL(wp) ::  exponent     !<
+       REAL(wp) ::  sum_spec_dum !< wavenumber-integrated spectrum
+    
+       REAL(wp), DIMENSION(0:nx) ::  work !<
+    
+       REAL(wp), DIMENSION(0:nx/2) ::  sums_spectra_l !<
+    
+       REAL(wp), DIMENSION(0:nx/2,100) ::  sums_spectra !<
+    
+       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ddd !<
+
+!
+!--    Exponent for geometric average
+       exponent = 1.0_wp / ( ny + 1.0_wp )
+
+!
+!--    Loop over all levels defined by the user
+       n = 1
+       DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
+
+          k = comp_spectra_level(n)
+
+!
+!--       Calculate FFT only if the corresponding level is situated on this PE
+          IF ( k >= nzb_x  .AND.  k <= nzt_x )  THEN
+          
+             DO  j = nys_x, nyn_x
+
+                work = ddd(0:nx,j,k)
+                CALL fft_x_1d( work, 'forward' )
+
+                ddd(0,j,k) = dx * work(0)**2
+                DO  i = 1, nx/2
+                   ddd(i,j,k) = dx * ( work(i)**2 + work(nx+1-i)**2 )
+                ENDDO
+
+             ENDDO
+
+!
+!--          Local sum and geometric average of these spectra
+!--          (WARNING: no global sum should be performed, because floating
+!--          point overflow may occur)
+             DO  i = 0, nx/2
+
+                sums_spectra_l(i) = 1.0_wp
+                DO  j = nys_x, nyn_x
+                   sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
+                ENDDO
+
+             ENDDO
+          
+          ELSE
+
+             sums_spectra_l = 1.0_wp
+
+          ENDIF
+
+!
+!--       Global sum of spectra on PE0 (from where they are written on file)
+          sums_spectra(:,n) = 0.0_wp
+#if defined( __parallel )   
+          CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
+          CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1,       &
+                           MPI_REAL, MPI_PROD, 0, comm2d, ierr )
+#else
+          sums_spectra(:,n) = sums_spectra_l
+#endif
+!
+!--       Normalize spectra by variance
+          sum_spec_dum = SUM( sums_spectra(:,n) )
+          IF ( sum_spec_dum /= 0.0_wp )  THEN
+             sums_spectra(:,n) = sums_spectra(:,n) * var_d(k) / sum_spec_dum
+          ENDIF
+          n = n + 1
+
+       ENDDO
+       n = n - 1
+
+       IF ( myid == 0 )  THEN
+!
+!--       Sum of spectra for later averaging (see routine data_output_spectra)
+          DO  i = 1, nx/2
+             DO k = 1, n
+                spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k)
+             ENDDO
+          ENDDO
+
+       ENDIF
+!
+!--    n_sp_x is needed by data_output_spectra_x
+       n_sp_x = n
+
+    END SUBROUTINE calc_spectra_x
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> @todo Missing subroutine description.
+!------------------------------------------------------------------------------!
+    SUBROUTINE calc_spectra_y( ddd, pr, m )
+
+       USE control_parameters,                                                 &
+           ONLY:  fft_method
+
+       USE fft_xy,                                                             &
+           ONLY:  fft_y_1d
+
+       USE grid_variables,                                                     &
+           ONLY:  dy
+
+       USE indices,                                                            &
+           ONLY:  nx, ny
+
+       USE kinds
+
+#if defined( __lc )
+       USE MPI
+#else
+       INCLUDE "mpif.h"
+#endif
+       USE pegrid,                                                             &
+           ONLY:  comm2d, ierr, myid
+
+       USE statistics,                                                         &
+           ONLY:  spectrum_y, var_d
+
+       USE transpose_indices,                                                  &
+           ONLY:  nxl_yd, nxr_yd, nzb_yd, nzt_yd
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i         !<
+       INTEGER(iwp) ::  j         !<
+       INTEGER(iwp) ::  jshape(1) !<
+       INTEGER(iwp) ::  k         !<
+       INTEGER(iwp) ::  m         !<
+       INTEGER(iwp) ::  n         !<
+       INTEGER(iwp) ::  pr        !<
+
+       REAL(wp) ::  exponent !<
+       REAL(wp) ::  sum_spec_dum !< wavenumber-integrated spectrum
+    
+       REAL(wp), DIMENSION(0:ny) ::  work !<
+    
+       REAL(wp), DIMENSION(0:ny/2) ::  sums_spectra_l !<
+    
+       REAL(wp), DIMENSION(0:ny/2,100) ::  sums_spectra !<
+    
+       REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !<
+
+
+!
+!--    Exponent for geometric average
+       exponent = 1.0_wp / ( nx + 1.0_wp )
+
+!
+!--    Loop over all levels defined by the user
+       n = 1
+       DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
+
+          k = comp_spectra_level(n)
+
+!
+!--       Calculate FFT only if the corresponding level is situated on this PE
+          IF ( k >= nzb_yd  .AND.  k <= nzt_yd )  THEN
+          
+             DO  i = nxl_yd, nxr_yd
+
+                work = ddd(0:ny,i,k)
+                CALL fft_y_1d( work, 'forward' )
+
+                ddd(0,i,k) = dy * work(0)**2
+                DO  j = 1, ny/2
+                   ddd(j,i,k) = dy * ( work(j)**2 + work(ny+1-j)**2 )
+                ENDDO
+
+             ENDDO
+
+!
+!--          Local sum and geometric average of these spectra
+!--          (WARNING: no global sum should be performed, because floating
+!--          point overflow may occur)
+             DO  j = 0, ny/2
+
+                sums_spectra_l(j) = 1.0_wp
+                DO  i = nxl_yd, nxr_yd
+                   sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
+                ENDDO
+
+             ENDDO
+          
+          ELSE
+
+             sums_spectra_l = 1.0_wp
+
+          ENDIF
+
+!
+!--       Global sum of spectra on PE0 (from where they are written on file)
+          sums_spectra(:,n) = 0.0_wp
+#if defined( __parallel )   
+          CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
+          CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1,       &
+                           MPI_REAL, MPI_PROD, 0, comm2d, ierr )
+#else
+          sums_spectra(:,n) = sums_spectra_l
+#endif
+!
+!--       Normalize spectra by variance
+          sum_spec_dum = SUM( sums_spectra(:,n) )
+          IF ( SUM(sums_spectra(:,n)) /= 0.0_wp )  THEN
+             sums_spectra(:,n) = sums_spectra(:,n) *                           &
+                                 var_d(k) / SUM(sums_spectra(:,n))
+          ENDIF
+          n = n + 1
+
+       ENDDO
+       n = n - 1
+
+
+       IF ( myid == 0 )  THEN
+!
+!--       Sum of spectra for later averaging (see routine data_output_spectra)
+          DO  j = 1, ny/2
+             DO k = 1, n
+                spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k)
+             ENDDO
+          ENDDO
+
+       ENDIF
+
+!
+!--    n_sp_y is needed by data_output_spectra_y
+       n_sp_y = n
+
+    END SUBROUTINE calc_spectra_y
+
+ END MODULE spectrum
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 1785)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 1786)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! ------------------
-!
+! +module spectrum
 !
 ! Former revisions:
@@ -271,4 +271,7 @@
               skip_time_do_radiation, time_radiation
 
+    USE spectrum,                                                              &
+        ONLY: calc_spectra
+
     USE statistics,                                                            &
         ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
