Home


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

Changeset 1786


Ignore:
Timestamp:
Mar 8, 2016 5:49:27 AM (9 years ago)
Author:
raasch
Message:

pmc-change in server-client get-put, spectra-directives removed, spectra-package modularized

Location:
palm/trunk/SOURCE
Files:
12 edited
1 moved

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/Makefile

    r1784 r1786  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# rename calc_spectra to spectrum + modularization of spectrum
    2323#
    2424# Former revisions:
    … …  
    231231        advec_w_pw.f90 advec_w_up.f90 average_3d_data.f90 boundary_conds.f90 \
    232232        buoyancy.f90 calc_liquid_water_content.f90 calc_mean_profile.f90 \
    233         calc_precipitation.f90 calc_radiation.f90 calc_spectra.f90 \
     233        calc_precipitation.f90 calc_radiation.f90 \
    234234        check_for_restart.f90 check_open.f90 check_parameters.f90 \
    235235        close_file.f90 compute_vpt.f90 coriolis.f90 cpulog.f90 \
    … …  
    268268        random_function.f90 random_gauss.f90 random_generator_parallel.f90 \
    269269        read_3d_binary.f90 read_var_list.f90 run_control.f90 \
    270         set_slicer_attributes_dvrp.f90 singleton.f90 sor.f90 \
     270        set_slicer_attributes_dvrp.f90 singleton.f90 sor.f90 spectrum.f90 \
    271271        subsidence.f90 sum_up_3d_data.f90 \
    272272        surface_coupler.f90 surface_layer_fluxes.f90 swap_timelevel.f90 temperton_fft.f90 \
    … …  
    336336calc_precipitation.o: modules.o mod_kinds.o
    337337calc_radiation.o: modules.o mod_kinds.o
    338 calc_spectra.o: modules.o cpulog.o fft_xy.o mod_kinds.o
    339338check_for_restart.o: modules.o mod_kinds.o
    340339check_open.o: modules.o mod_kinds.o mod_particle_attributes.o netcdf_interface.o
    … …  
    355354data_output_ptseries.o: modules.o cpulog.o mod_kinds.o \
    356355   netcdf_interface.o mod_particle_attributes.o
    357 data_output_spectra.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o
     356data_output_spectra.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o \
     357   spectrum.o
    358358data_output_tseries.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o
    359359data_output_2d.o: modules.o cpulog.o mod_kinds.o mod_particle_attributes.o\
    … …  
    378378global_min_max.o: modules.o mod_kinds.o
    379379header.o: modules.o cpulog.o mod_kinds.o netcdf_interface.o land_surface_model.o\
    380           plant_canopy_model.o pmc_interface.o radiation_model.o subsidence.o
     380   plant_canopy_model.o pmc_interface.o radiation_model.o spectrum.o \
     381   subsidence.o
    381382impact_of_latent_heat.o: modules.o mod_kinds.o
    382383inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
    … …  
    440441mod_kinds.o: mod_kinds.f90
    441442mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
    442 netcdf_interface.o: netcdf_interface.f90 modules.o mod_kinds.o land_surface_model.o
     443netcdf_interface.o: netcdf_interface.f90 modules.o mod_kinds.o \
     444   land_surface_model.o spectrum.o
    443445nudging.o: modules.o cpulog.o mod_kinds.o
    444 package_parin.o: modules.o mod_kinds.o land_surface_model.o\
    445                  plant_canopy_model.o radiation_model.o
     446package_parin.o: modules.o mod_kinds.o land_surface_model.o \
     447   plant_canopy_model.o radiation_model.o spectrum.o
    446448palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o\
    447449        pmc_interface.o surface_layer_fluxes.o
    … …  
    483485singleton.o: mod_kinds.o singleton.f90
    484486sor.o: modules.o mod_kinds.o
     487spectrum.o: spectrum.f90 modules.o mod_kinds.o cpulog.o fft_xy.o
    485488subsidence.o: modules.o mod_kinds.o
    486489sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o land_surface_model.o\
    … …  
    495498        ls_forcing.o mod_kinds.o nudging.o pmc_interface.o production_e.o \
    496499        prognostic_equations.o progress_bar.o radiation_model.o \
    497         user_actions.o surface_layer_fluxes.o
     500        spectrum.o user_actions.o surface_layer_fluxes.o
    498501time_to_string.o: mod_kinds.o
    499502timestep.o: modules.o cpulog.o mod_kinds.o
    … …  
    528531user_parin.o: modules.o mod_kinds.o user_module.o
    529532user_read_restart_data.o: modules.o mod_kinds.o user_module.o
    530 user_spectra.o: modules.o mod_kinds.o user_module.o
     533user_spectra.o: modules.o mod_kinds.o spectrum.o user_module.o
    531534user_statistics.o: modules.o mod_kinds.o netcdf_interface.o user_module.o
    532535wall_fluxes.o: modules.o mod_kinds.o

  • palm/trunk/SOURCE/check_parameters.f90

    r1784 r1786  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! cpp-direktives for spectra removed, check of spectra level removed
    2222!
    2323! Former revisions:
    … …  
    326326    USE profil_parameter
    327327    USE radiation_model_mod
    328     USE spectrum
    329328    USE statistics
    330329    USE subsidence_mod
    … …  
    37523751    END SELECT
    37533752
    3754 #if defined( __spectra )
    3755 !
    3756 !-- Check the number of spectra level to be output
    3757     i = 1
    3758     DO WHILE ( comp_spectra_level(i) /= 999999  .AND.  i <= 100 )
    3759        i = i + 1
    3760     ENDDO
    3761     i = i - 1
    3762     IF ( i == 0 )  THEN
    3763        WRITE( message_string, * )  'no spectra levels given'
    3764        CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
    3765     ENDIF
    3766 #endif
    3767 
    37683753!
    37693754!-- Check mask conditions

  • palm/trunk/SOURCE/data_output_spectra.f90

    r1784 r1786  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! cpp-directives for spectra removed, immediate return if no spectra levels are
     22! given
    2223!
    2324! Former revisions:
    … …  
    7374 
    7475#if defined( __netcdf )
    75 #if defined( __spectra )
    76 
    7776    USE control_parameters,                                                    &
    7877        ONLY:  average_count_sp, averaging_interval_sp, dosp_time_count,       &
    … …  
    9392
    9493    USE spectrum,                                                              &
    95         ONLY:  data_output_sp, spectra_direction
     94        ONLY:  comp_spectra_level, data_output_sp, spectra_direction
    9695
    9796    USE statistics,                                                            &
    … …  
    110109
    111110    CALL cpu_log( log_point(31), 'data_output_spectra', 'start' )
     111
     112!
     113!-- Check if user gave any levels for spectra to be calculated
     114    IF ( comp_spectra_level(1) == 999999 )  RETURN
    112115
    113116!
    … …  
    206209
    207210#endif
    208 #endif
    209211 END SUBROUTINE data_output_spectra
    210212
    … …  
    299301
    300302
    301 #if defined( __spectra )
    302303!------------------------------------------------------------------------------!
    303304! Description:
    … …  
    714715
    715716 END SUBROUTINE data_output_spectra_y
    716 #endif

  • palm/trunk/SOURCE/header.f90

    r1784 r1786  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! cpp-direktives for spectra removed
    2222!
    2323! Former revisions:
    … …  
    17331733#endif
    17341734
    1735 #if defined( __spectra )
    17361735!
    17371736!-- Spectra output
    … …  
    17531752                          averaging_interval_sp, dt_averaging_input_pr
    17541753    ENDIF
    1755 #endif
    17561754
    17571755    WRITE ( io, 99 )
    … …  
    22342232367 FORMAT ('       Polygon reduction for topography: cluster_size = ', I1)
    22352233#endif
    2236 #if defined( __spectra )
    22372234370 FORMAT ('    Spectra:')
    22382235371 FORMAT ('       Output every ',F7.1,' s'/)
    … …  
    22532250            '       Profiles for the time averaging are taken every ', &
    22542251                    F6.1,' s')
    2255 #endif
    22562252400 FORMAT (//' Physical quantities:'/ &
    22572253              ' -------------------'/)

  • palm/trunk/SOURCE/modules.f90

    r1784 r1786  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! module spectrum moved to new separate module
    2222!
    2323! Former revisions:
    … …  
    12771277 END MODULE profil_parameter
    12781278
    1279 
    1280 !------------------------------------------------------------------------------!
    1281 ! Description:
    1282 ! ------------
    1283 !> Definition of quantities used for computing spectra.
    1284 !------------------------------------------------------------------------------!
    1285  MODULE spectrum
    1286 
    1287     USE kinds
    1288 
    1289     CHARACTER (LEN=6),  DIMENSION(1:5) ::  header_char = (/ 'PS(u) ', 'PS(v) ',&
    1290                                            'PS(w) ', 'PS(pt)', 'PS(q) ' /)
    1291     CHARACTER (LEN=2),  DIMENSION(10)  ::  spectra_direction = 'x'
    1292     CHARACTER (LEN=10), DIMENSION(10)  ::  data_output_sp  = ' '
    1293     CHARACTER (LEN=25), DIMENSION(1:5) ::  utext_char =                    &
    1294                                            (/ '-power spectrum of u     ', &
    1295                                               '-power spectrum of v     ', &
    1296                                               '-power spectrum of w     ', &
    1297                                               '-power spectrum of ^1185 ', &
    1298                                               '-power spectrum of q     ' /)
    1299     CHARACTER (LEN=39), DIMENSION(1:5) ::  ytext_char =                        &
    1300                                  (/ 'k ^2236 ^2566^2569<u(k) in m>2s>->2    ', &
    1301                                     'k ^2236 ^2566^2569<v(k) in m>2s>->2    ', &
    1302                                     'k ^2236 ^2566^2569<w(k) in m>2s>->2    ', &
    1303                                     'k ^2236 ^2566^2569<^1185(k) in m>2s>->2', &
    1304                                     'k ^2236 ^2566^2569<q(k) in m>2s>->2    ' /)
    1305 
    1306     INTEGER(iwp) ::  klist_x = 0, klist_y = 0, n_sp_x = 0, n_sp_y = 0
    1307 
    1308     INTEGER(iwp) ::  comp_spectra_level(100) = 999999,                   &
    1309                      lstyles(100) = (/ 0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
    1310                                        0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
    1311                                        0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
    1312                                        0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
    1313                                        0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
    1314                                        0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
    1315                                        0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
    1316                                        0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
    1317                                        0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
    1318                                        0, 7, 3, 10, 1, 4, 9, 2, 6, 8 /), &
    1319                      plot_spectra_level(100) = 999999
    1320 
    1321     REAL(wp)    ::  time_to_start_sp = 0.0_wp
    1322 
    1323     SAVE
    1324 
    1325  END MODULE spectrum
    1326 
    1327 
    13281279!------------------------------------------------------------------------------!
    13291280! Description:
    … …  
    13651316
    13661317
     1318
    13671319!------------------------------------------------------------------------------!
    13681320! Description:

  • palm/trunk/SOURCE/netcdf_interface.f90

    r1784 r1786  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Bugfix: id_var_time_sp made public
    2222!
    2323! Former revisions:
    … …  
    296296            id_var_dopr, id_var_dopts, id_var_dospx, id_var_dospy,             &
    297297            id_var_dots, id_var_do2d, id_var_do3d, id_var_norm_dopr,           &
    298             id_var_time_mask, id_var_time_pr, id_var_time_pts, id_var_time_ts, &
    299             id_var_time_xy, id_var_time_xz, id_var_time_yz,id_var_time_3d,     &
    300             nc_stat, netcdf_data_format, netcdf_data_format_string,            &
    301             netcdf_deflate, netcdf_precision, output_for_t0
     298            id_var_time_mask, id_var_time_pr, id_var_time_pts, id_var_time_sp, &
     299            id_var_time_ts, id_var_time_xy, id_var_time_xz, id_var_time_yz,    &
     300            id_var_time_3d, nc_stat, netcdf_data_format,                       &
     301            netcdf_data_format_string, netcdf_deflate, netcdf_precision,       &
     302            output_for_t0
    302303
    303304    SAVE

  • palm/trunk/SOURCE/package_parin.f90

    r1758 r1786  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! cpp-direktives for spectra removed
    2222!
    2323! Former revisions:
    … …  
    318318 30 CONTINUE
    319319
    320 
    321 #if defined( __spectra )
    322320    REWIND ( 11 )
    323321    line = ' '
    … …  
    337335
    338336 40 CONTINUE
    339 #endif
    340337
    341338!

  • palm/trunk/SOURCE/pmc_client.f90

    r1784 r1786  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! change in client-server data transfer: server now gets data from client
     23! instead that client put's it to the server
    2324!
    2425! Former revisions:
    … …  
    361362       do i=1,me%inter_npes
    362363          aPE => me%PEs(i)
    363           ar  => aPE%array_list(next_array_in_list)    !actual array is last array in list
     364          ar  => aPE%array_list(next_array_in_list)
    364365          ar%NrDims    = NrDims
    365366          ar%A_dim     = dims
    … …  
    407408      IMPLICIT none
    408409
    409       INTEGER                                 :: i, ierr, j
     410!--   naming convention:  appending       _sc  -> server to client transfer
     411!--                                       _cs  -> client to server transfer
     412!--                                       Recv -> server to client transfer
     413!--                                       Send -> client to server transfer
     414
     415      INTEGER                                 :: i, istat, ierr, j
     416      INTEGER,PARAMETER                       :: NoINdex=-1
     417      INTEGER                                 :: rcount
    410418      INTEGER                                 :: arlen, myIndex, tag
    411419      INTEGER(idp)                            :: bufsize                   ! Size of MPI data Window
    412420      TYPE(PeDef),POINTER                     :: aPE
    413421      TYPE(ArrayDef),POINTER                  :: ar
     422      INTEGER,DIMENSION(1024)                 :: req
    414423      character(len=DA_Namelen)               :: myName
    415424      Type(c_ptr)                             :: base_ptr
    416       REAL(kind=wp),DIMENSION(:),POINTER      :: base_array
     425      REAL(kind=wp),DIMENSION(:),POINTER,save :: base_array_sc             !Base array
     426      REAL(kind=wp),DIMENSION(:),POINTER,save :: base_array_cs             !Base array
    417427      INTEGER(KIND=MPI_ADDRESS_KIND)          :: WinSize
    418428
    … …  
    420430      bufsize = 8
    421431
    422       ! First stride, Compute size and set index
     432!--   Server to client direction
     433
     434!--   First stride, Compute size and set index
    423435
    424436      do i=1,me%inter_npes
    … …  
    433445            CALL MPI_Recv (myIndex, 1, MPI_INTEGER, i-1, tag, me%inter_comm, MPI_STATUS_IGNORE, ierr)
    434446
    435             if(ar%NrDims == 3) then                    ! PALM has k in first dimension
     447            if(ar%NrDims == 3) then
    436448               bufsize = max(bufsize,ar%A_dim(1)*ar%A_dim(2)*ar%A_dim(3))    ! determine max, because client buffer is allocated only once
    437449            else
    438450               bufsize = max(bufsize,ar%A_dim(1)*ar%A_dim(2))
    439451            end if
    440             ar%BufIndex = myIndex
    441           end do
     452            ar%RecvIndex = myIndex
     453
     454           end do
    442455      end do
    443456
    444       ! Create RMA (One Sided Communication) window for data buffer
    445 
    446       CALL PMC_Alloc_mem (base_array, bufsize, base_ptr)
    447       me%TotalBufferSize = bufsize*wp                          !Total buffer size in Byte
    448 
    449       WinSize = me%TotalBufferSize
    450 !      write(9,'(a,8i7)') 'PMC_S_SetInd_and_Mem ',m_model_rank,me%inter_npes,WinSize,ar%A_dim
    451       CALL MPI_Win_create (base_array, WinSize, wp, MPI_INFO_NULL, me%intra_comm, me%BufWin, ierr);
    452       CALL MPI_Win_fence (0, me%BufWin, ierr);                    !  Open Window to set data
    453       CALL MPI_Barrier(me%intra_comm, ierr)
     457
     458!--   Create RMA (One Sided Communication) data buffer
     459!--   The buffer for MPI_Get can be PE local, i.e. it can but must not be part of the MPI RMA window
     460
     461      CALL PMC_Alloc_mem (base_array_sc, bufsize, base_ptr)
     462      me%TotalBufferSize = bufsize*wp                          ! Total buffer size in Byte
     463
     464!--   Second stride, Set Buffer pointer
    454465
    455466      do i=1,me%inter_npes
    … …  
    458469         do j=1,aPE%Nr_arrays
    459470            ar  => aPE%array_list(j)
    460             ar%SendBuf = base_ptr
     471            ar%RecvBuf = base_ptr
     472         end do
     473      end do
     474
     475!--   Client to server direction
     476
     477      myIndex = 1
     478      rCount  = 0
     479      bufsize = 8
     480
     481      do i=1,me%inter_npes
     482         aPE => me%PEs(i)
     483         tag = 300
     484         do j=1,aPE%Nr_arrays
     485            ar  => aPE%array_list(j)
     486            if(ar%NrDims == 2) then
     487               arlen     = aPE%NrEle                        ! 2D
     488            else if(ar%NrDims == 3) then
     489               arlen     = aPE%NrEle*ar%A_dim(1)            ! 3D
     490            end if
     491
     492            tag    = tag+1
     493            rCount = rCount+1
     494            if(aPE%NrEle > 0)  then
     495               CALL MPI_Isend (myIndex, 1, MPI_INTEGER, i-1, tag, me%inter_comm, req(rCount),ierr)
     496               ar%SendIndex = myIndex
     497            else
     498               CALL MPI_Isend (NoIndex, 1, MPI_INTEGER, i-1, tag, me%inter_comm, req(rCount),ierr)
     499               ar%SendIndex = NoIndex
     500            end if
     501
     502            if(rCount == 1024) then                                  ! Maximum of 1024 outstanding requests
     503               CALL MPI_Waitall (rCount, req, MPI_STATUSES_IGNORE, ierr)
     504               rCount = 0;
     505            end if
     506
     507            if(aPE%NrEle > 0)  then
     508               ar%SendSize  = arlen
     509               myIndex     = myIndex+arlen
     510               bufsize     = bufsize+arlen
     511            end if
     512          end do
     513         if(rCount > 0) then                                          ! Wait for all send completed
     514            CALL MPI_Waitall (rCount, req, MPI_STATUSES_IGNORE, ierr)
     515         end if
     516      end do
     517
     518!--   Create RMA (One Sided Communication) window for data buffer client to server transfer
     519!--   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
     520!--   Only one RMA window is required to prepare the data for server -> client transfer on the server side and
     521!--                                                       for client -> server transfer on the client side
     522
     523      CALL PMC_Alloc_mem (base_array_cs, bufsize)
     524      me%TotalBufferSize = bufsize*wp                          !Total buffer size in Byte
     525
     526      WinSize = me%TotalBufferSize
     527      CALL MPI_Win_create (base_array_cs, WinSize, wp, MPI_INFO_NULL, me%intra_comm, me%win_server_client, ierr);
     528      CALL MPI_Win_fence (0, me%win_server_client, ierr);                    !  Open Window to set data
     529      CALL MPI_Barrier(me%intra_comm, ierr)
     530
     531!--   Second stride, Set Buffer pointer
     532
     533      do i=1,me%inter_npes
     534         aPE => me%PEs(i)
     535
     536         do j=1,aPE%Nr_arrays
     537            ar  => aPE%array_list(j)
     538            if(aPE%NrEle > 0)  then
     539              ar%SendBuf = c_loc(base_array_cs(ar%SendIndex))
     540              if(ar%SendIndex+ar%SendSize > bufsize) then
     541                 write(0,'(a,i4,4i7,1x,a)') 'Client Buffer too small ',i,ar%SendIndex,ar%SendSize,ar%SendIndex+ar%SendSize,bufsize,trim(ar%name)
     542                 CALL MPI_Abort (MPI_COMM_WORLD, istat, ierr)
     543              end if
     544            end if
    461545         end do
    462546      end do
    … …  
    504588
    505589            buf_shape(1) = nr
    506             CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
     590            CALL c_f_pointer(ar%RecvBuf, buf, buf_shape)
    507591!
    508592!--         MPI passive target RMA
    509593            if(nr > 0)   then
    510                target_disp = (ar%BufIndex-1)
    511                CALL MPI_Win_lock (MPI_LOCK_SHARED , ip-1, 0, me%BufWin, ierr)
    512                CALL MPI_Get (buf, nr, MPI_REAL, ip-1, target_disp, nr, MPI_REAL, me%BufWin, ierr)
    513                CALL MPI_Win_unlock (ip-1, me%BufWin, ierr)
     594               target_disp = (ar%RecvIndex-1)
     595               CALL MPI_Win_lock (MPI_LOCK_SHARED , ip-1, 0, me%win_server_client, ierr)
     596               CALL MPI_Get (buf, nr, MPI_REAL, ip-1, target_disp, nr, MPI_REAL, me%win_server_client, ierr)
     597               CALL MPI_Win_unlock (ip-1, me%win_server_client, ierr)
    514598            end if
    515599
    … …  
    555639      INTEGER(kind=MPI_ADDRESS_KIND)    ::  target_disp
    556640
     641      t1 = PMC_Time()
     642      CALL MPI_Barrier(me%intra_comm, ierr)              ! Wait for empty buffer
     643      t2 = PMC_Time()
     644      if(present(WaitTime)) WaitTime = t2-t1
    557645
    558646      do ip=1,me%inter_npes
    … …  
    561649         do j=1,aPE%Nr_arrays
    562650            ar  => aPE%array_list(j)
     651            myIndex=1
    563652            if(ar%NrDims == 2) then
    564                nr = aPE%NrEle
    565             else if(ar%NrDims == 3) then
    566                nr = aPE%NrEle*ar%A_dim(1)
    567             end if
    568 
    569             buf_shape(1) = nr
    570             CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
    571 
    572             myIndex = 1
    573             if(ar%NrDims == 2) then
     653               buf_shape(1) = aPE%NrEle
     654               CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
    574655               CALL c_f_pointer(ar%data, data_2d, ar%A_dim(1:2))
    575656               do ij=1,aPE%NrEle
    … …  
    578659               end do
    579660            else if(ar%NrDims == 3) then
     661               buf_shape(1) = aPE%NrEle*ar%A_dim(1)
     662               CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
    580663               CALL c_f_pointer(ar%data, data_3d, ar%A_dim(1:3))
    581664               do ij=1,aPE%NrEle
    … …  
    584667               end do
    585668            end if
    586 !
    587 !--         MPI passiv target RMA
    588             if(nr > 0)   then
    589                target_disp = (ar%BufIndex-1)
    590                CALL MPI_Win_lock (MPI_LOCK_SHARED , ip-1, 0, me%BufWin, ierr)
    591                CALL MPI_Put (buf, nr, MPI_REAL, ip-1, target_disp, nr, MPI_REAL, me%BufWin, ierr)
    592                CALL MPI_Win_unlock (ip-1, me%BufWin, ierr)
    593             end if
    594 
    595669          end do
    596670      end do
    597671
    598672
    599       t1 = PMC_Time()
    600       CALL MPI_Barrier(me%intra_comm, ierr)                         ! Wait for server to fill buffer
    601       t2 = PMC_Time()
    602       if(present(WaitTime)) WaitTime = t2-t1
    603 
    604       CALL MPI_Barrier(me%intra_comm, ierr)                         ! Wait for buffer is filled
    605 
     673!      CALL MPI_Win_fence (0, me%win_server_client, ierr)      ! Fence might do it, test later
     674      CALL MPI_Barrier(me%intra_comm, ierr)                   ! buffer is filled
    606675
    607676      return

  • palm/trunk/SOURCE/pmc_general.f90

    r1780 r1786  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! change in client-server data transfer: server now gets data from client
     23! instead that client put's it to the server
    2324!
    2425! Former revisions:
    … …  
    8283      TYPE(c_ptr)                   :: data                        ! Pointer of data in server space
    8384      TYPE(c_ptr), DIMENSION(2)     :: po_data                     ! Base Pointers, PMC_S_Set_Active_data_array sets active pointer
    84       INTEGER(idp)                  :: BufIndex                    ! index in Send Buffer
    85       INTEGER                       :: BufSize                     ! size in Send Buffer
     85      INTEGER(idp)                  :: SendIndex                   ! index in Send Buffer
     86      INTEGER(idp)                  :: RecvIndex                   ! index in Receive Buffer
     87      INTEGER                       :: SendSize                    ! size in Send Buffer
     88      INTEGER                       :: RecvSize                    ! size in Receiver Buffer
    8689      TYPE (c_ptr)                  :: SendBuf                     ! Pointer of Data in Send buffer
     90      TYPE (c_ptr)                  :: RecvBuf                     ! Pointer of Data in ReceiveS buffer
    8791      CHARACTER(len=8)              :: Name                        ! Name of Array
    8892
    … …  
    108112      INTEGER                       :: inter_npes                  ! Number of PEs client model
    109113      INTEGER                       :: intra_rank                  ! rank within intra_comm
    110       INTEGER                       :: BufWin                      ! MPI RMA windows
     114      INTEGER                       :: win_server_client           ! MPI RMA for preparing data on server AND client side
    111115      TYPE (PeDef), DIMENSION(:), POINTER      :: PEs              ! List of all Client PEs
    112116   END TYPE ClientDef

  • palm/trunk/SOURCE/pmc_handle_communicator.f90

    r1780 r1786  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Bugfix: nesting_mode is broadcast now
    2323!
    2424! Former revisions:
    … …  
    211211         CALL MPI_BCAST( m_couplers(i)%lower_left_x, 1, MPI_REAL,    0, MPI_COMM_WORLD, istat )
    212212         CALL MPI_BCAST( m_couplers(i)%lower_left_y, 1, MPI_REAL,    0, MPI_COMM_WORLD, istat )
     213         CALL MPI_BCAST( nesting_mode, LEN( nesting_mode ), MPI_CHARACTER, 0, MPI_COMM_WORLD, istat )
    213214      ENDDO
    214215

  • palm/trunk/SOURCE/pmc_server.f90

    r1780 r1786  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! change in client-server data transfer: server now gets data from client
     23! instead that client put's it to the server
    2324!
    2425! Former revisions:
    … …  
    336337      IMPLICIT none
    337338
     339!
     340!--   Naming convention for appendices:   _sc  -> server to client transfer
     341!--                                       _cs  -> client to server transfer
     342!--                                       Send -> Server to client transfer
     343!--                                       Recv -> client to server transfer
    338344      INTEGER,INTENT(IN)                      :: ClientId
    339345
    … …  
    347353      INTEGER,DIMENSION(1024)                 :: req
    348354      Type(c_ptr)                             :: base_ptr
    349       REAL(kind=wp),DIMENSION(:),POINTER      :: base_array
     355      REAL(wp),DIMENSION(:),POINTER,SAVE :: base_array_sc  !Base array for server to client transfer
     356      REAL(wp),DIMENSION(:),POINTER,SAVE :: base_array_cs  !Base array for client to server transfer
    350357      INTEGER(KIND=MPI_ADDRESS_KIND)          :: WinSize
    351358
     359!
     360!--   Server to client direction
    352361      myIndex = 1
    353362      rCount  = 0
    354363      bufsize = 8
    355364
    356       ! First stride, Compute size and set index
    357 
     365!
     366!--   First stride: Compute size and set index
    358367      do i=1,Clients(ClientId)%inter_npes
    359368         aPE => Clients(ClientId)%PEs(i)
    … …  
    362371            ar  => aPE%array_list(j)
    363372            if(ar%NrDims == 2) then
    364                arlen     = aPE%NrEle;                             ! 2D
     373               arlen     = aPE%NrEle                              ! 2D
    365374            else if(ar%NrDims == 3) then
    366                arlen     = aPE%NrEle * ar%A_dim(4);               ! PALM 3D
     375               arlen     = aPE%NrEle * ar%A_dim(4);               ! 3D
    367376            else
    368377               arlen     = -1
    369378            end if
    370             ar%BufIndex = myIndex
     379            ar%SendIndex = myIndex
    371380
    372381            tag    = tag+1
    … …  
    381390            myIndex = myIndex+arlen
    382391            bufsize = bufsize+arlen
    383             ar%BufSize = arlen
     392            ar%SendSize = arlen
    384393
    385394         end do
    386          if(rCount > 0) then                                          ! Wait for all send completed
     395         if(rCount > 0) then                       ! Wait for all send completed
    387396            CALL MPI_Waitall (rCount, req, MPI_STATUSES_IGNORE, ierr)
    388397         end if
    389398      end do
    390399
    391       ! Create RMA (One Sided Communication) window for data buffer
    392 
    393       CALL PMC_Alloc_mem (base_array, bufsize, base_ptr)
    394       Clients(ClientId)%TotalBufferSize = bufsize*wp       !Total buffer size in Byte
     400!
     401!--   Create RMA (One Sided Communication) window for data buffer server to
     402!--   client transfer.
     403!--   The buffer of MPI_Get (counter part of transfer) can be PE-local, i.e.
     404!--   it can but must not be part of the MPI RMA window.
     405!--   Only one RMA window is required to prepare the data
     406!--   for server -> client transfer on the server side and
     407!--   for client -> server transfer on the client side
     408      CALL PMC_Alloc_mem (base_array_sc, bufsize)
     409      Clients(ClientId)%TotalBufferSize = bufsize*wp   !Total buffer size in Byte
    395410
    396411      WinSize = bufsize*wp
    397 !      write(9,*) 'PMC_S_SetInd_and_Mem ',m_model_rank,Clients(ClientId)%inter_npes,WinSize,bufsize
    398       CALL MPI_Win_create (base_array, WinSize, wp, MPI_INFO_NULL, Clients(ClientId)%intra_comm, Clients(ClientId)%BufWin, ierr);
    399       CALL MPI_Win_fence (0, Clients(ClientId)%BufWin, ierr);                    !  Open Window to set data
    400       CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)
    401 
    402       ! Second stride, Set Buffer pointer
    403 
     412      CALL MPI_Win_create (base_array_sc, WinSize, wp, MPI_INFO_NULL,          &
     413                  Clients(ClientId)%intra_comm,  Clients(ClientId)%win_server_client, ierr)
     414      CALL MPI_Win_fence (0, Clients(ClientId)%win_server_client, ierr);        !  Open Window to set data
     415!
     416!--   Second stride: Set Buffer pointer
    404417      do i=1,Clients(ClientId)%inter_npes
    405418         aPE => Clients(ClientId)%PEs(i)
    406419         do j=1,aPE%Nr_arrays
    407420            ar  => aPE%array_list(j)
    408 !--         TO_DO:  Adressrechnung ueberlegen?
    409             ar%SendBuf = c_loc(base_array(ar%BufIndex))                         !kk Adressrechnung ueberlegen
    410             if(ar%BufIndex+ar%BufSize > bufsize) then
    411 !--            TO_DO: can this error really happen, and what can be the reason?
    412                write(0,'(a,i4,4i7,1x,a)') 'Buffer too small ',i,ar%BufIndex,ar%BufSize,ar%BufIndex+ar%BufSize,bufsize,trim(ar%name)
     421            ar%SendBuf = c_loc(base_array_sc(ar%SendIndex))
     422            if(ar%SendIndex+ar%SendSize > bufsize) then
     423               write(0,'(a,i4,4i7,1x,a)') 'Server Buffer too small ',i,ar%SendIndex,ar%SendSize,ar%SendIndex+ar%SendSize,bufsize,trim(ar%name)
    413424               CALL MPI_Abort (MPI_COMM_WORLD, istat, ierr)
    414425            end if
     426         end do
     427      end do
     428
     429!--   Client to server direction
     430
     431      bufsize  = 8
     432
     433!--   First stride, Compute size and set index
     434
     435      do i=1,Clients(ClientId)%inter_npes
     436         aPE => Clients(ClientId)%PEs(i)
     437         tag = 300
     438
     439         do j=1,aPE%Nr_arrays
     440            ar  => aPE%array_list(j)
     441
     442            ! Receive Index from client
     443            tag = tag+1
     444            CALL MPI_Recv (myIndex, 1, MPI_INTEGER, i-1, tag, Clients(ClientId)%inter_comm, MPI_STATUS_IGNORE, ierr)
     445
     446            if(ar%NrDims == 3) then
     447               bufsize = max(bufsize,aPE%NrEle * ar%A_dim(4))               ! 3D
     448            else
     449               bufsize = max(bufsize,aPE%NrEle)                             ! 2D
     450            end if
     451            ar%RecvIndex = myIndex
     452          end do
     453
     454      end do
     455
     456!--   Create RMA (One Sided Communication) data buffer
     457!--   The buffer for MPI_Get can be PE local, i.e. it can but must not be part of the MPI RMA window
     458
     459      CALL PMC_Alloc_mem (base_array_cs, bufsize, base_ptr)
     460      Clients(ClientId)%TotalBufferSize = bufsize*wp       !Total buffer size in Byte
     461
     462      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)
     463
     464!--   Second stride, Set Buffer pointer
     465
     466      do i=1,Clients(ClientId)%inter_npes
     467         aPE => Clients(ClientId)%PEs(i)
     468
     469         do j=1,aPE%Nr_arrays
     470            ar  => aPE%array_list(j)
     471            ar%RecvBuf = base_ptr
    415472         end do
    416473      end do
    … …  
    465522      end do
    466523
    467       CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)    ! buffer is full
     524      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)    ! buffer is filled
    468525
    469526      return
    … …  
    480537      INTEGER                             ::  ip,ij,istat,ierr,j
    481538      INTEGER                             ::  myIndex
     539      INTEGER                             ::  nr
    482540      REAL(wp)                            ::  t1,t2
    483541      TYPE(PeDef),POINTER                 ::  aPE
    … …  
    488546      REAL(wp), POINTER, DIMENSION(:,:)   ::  data_2d
    489547      REAL(wp), POINTER, DIMENSION(:,:,:) ::  data_3d
     548      INTEGER                             ::  target_pe
     549      INTEGER(kind=MPI_ADDRESS_KIND)      ::  target_disp
    490550
    491551      t1 = PMC_Time()
    492       CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)              ! Wait for MPI_Put from client
    493       t2 = PMC_Time()
    494       if(present(WaitTime)) WaitTime = t2-t1
     552      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)                         ! Wait for client to fill buffer
     553      t2 = PMC_Time()-t1
     554      if(present(WaitTime)) WaitTime = t2
     555
     556!      CALL MPI_Win_fence (0, Clients(ClientId)%win_server_client, ierr)            ! Fence might do it, test later
     557      CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)                         ! Wait for buffer is filled
    495558
    496559      do ip=1,Clients(ClientId)%inter_npes
    … …  
    498561         do j=1,aPE%Nr_arrays
    499562            ar  => aPE%array_list(j)
     563
     564            if(ar%RecvIndex < 0)  CYCLE
     565
     566            if(ar%NrDims == 2) then
     567               nr = aPE%NrEle
     568            else if(ar%NrDims == 3) then
     569               nr = aPE%NrEle*ar%A_dim(4)
     570            end if
     571
     572            buf_shape(1) = nr
     573            CALL c_f_pointer(ar%RecvBuf, buf, buf_shape)
     574!
     575!--         MPI passive target RMA
     576
     577            if(nr > 0)   then
     578               target_disp = ar%RecvIndex-1
     579               target_pe = ip-1+m_model_npes                         ! client PEs are located behind server PEs
     580               CALL MPI_Win_lock (MPI_LOCK_SHARED , target_pe, 0, Clients(ClientId)%win_server_client, ierr)
     581               CALL MPI_Get (buf, nr, MPI_REAL, target_pe, target_disp, nr, MPI_REAL, Clients(ClientId)%win_server_client, ierr)
     582               CALL MPI_Win_unlock (target_pe, Clients(ClientId)%win_server_client, ierr)
     583            end if
     584
    500585            myIndex=1
    501586            if(ar%NrDims == 2) then
    502                buf_shape(1) = aPE%NrEle
    503                CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
    504587               CALL c_f_pointer(ar%data, data_2d, ar%A_dim(1:2))
    505588               do ij=1,aPE%NrEle
    … …  
    508591               end do
    509592            else if(ar%NrDims == 3) then
    510                buf_shape(1) = aPE%NrEle*ar%A_dim(4)
    511                CALL c_f_pointer(ar%SendBuf, buf, buf_shape)
    512593               CALL c_f_pointer(ar%data, data_3d, ar%A_dim(1:3))
    513594               do ij=1,aPE%NrEle
    … …  
    519600      end do
    520601
    521       CALL MPI_Barrier(Clients(ClientId)%intra_comm, ierr)              ! data copy finished, buffer is free for use agein
    522 
    523         return
    524602   END SUBROUTINE PMC_S_GetData_from_Buffer
    525603

  • palm/trunk/SOURCE/spectrum.f90

    r1785 r1786  
    1 !> @file calc_spectra.f90
     1!> @file spectrum.f90
    22!--------------------------------------------------------------------------------!
    33! This file is part of PALM.
    … …  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! routine is modularized, filename renamed from calc_spectra to spectrum,
     22! privious data module spectrum moved from modules.f90 to here,
     23! cpp-direktives for spectra removed, immediate return if no spectra levels are
     24! given
    2225!
    2326! Former revisions:
    … …  
    7982!>            transpose_zyd needs modification).
    8083!------------------------------------------------------------------------------!
    81  SUBROUTINE calc_spectra
    82  
    83 
    84 #if defined( __spectra )
    85     USE arrays_3d,                                                             &
    86         ONLY:  d, tend
    87 
    88     USE control_parameters,                                                    &
    89         ONLY:  average_count_sp, bc_lr_cyc, bc_ns_cyc, message_string, psolver
    90 
    91     USE cpulog,                                                                &
    92         ONLY:  cpu_log, log_point
    93 
    94     USE fft_xy,                                                                &
    95         ONLY:  fft_init
    96 
    97     USE indices,                                                               &
    98         ONLY:  nxl, nxr, nyn, nys, nzb, nzt
     84 MODULE spectrum
    9985
    10086    USE kinds
    10187
    102     USE pegrid
    103 
    104     USE spectrum,                                                              &
    105         ONLY:  data_output_sp, spectra_direction
    106 
    107 
    108     IMPLICIT NONE
    109 
    110     INTEGER(iwp) ::  m  !<
    111     INTEGER(iwp) ::  pr !<
    112 
    113 
    114     CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
    115 
    116 !
    117 !-- Initialize ffts
    118     CALL fft_init
    119 
    120 !
    121 !-- Reallocate array d in required size
    122     IF ( psolver(1:9) == 'multigrid' )  THEN
    123        DEALLOCATE( d )
    124        ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
    125     ENDIF
    126 
    127     m = 1
    128     DO WHILE ( data_output_sp(m) /= ' '  .AND.  m <= 10 )
    129 !
    130 !--    Transposition from z --> x  ( y --> x in case of a 1d-decomposition
    131 !--    along x)
    132        IF ( INDEX( spectra_direction(m), 'x' ) /= 0 )  THEN
    133 
    134 !
    135 !--       Calculation of spectra works for cyclic boundary conditions only
    136           IF ( .NOT. bc_lr_cyc )  THEN
    137 
    138              message_string = 'non-cyclic lateral boundaries along x do not'// &
    139                               '& allow calculation of spectra along x'
    140              CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
     88    PRIVATE
     89
     90    CHARACTER (LEN=6),  DIMENSION(1:5) ::  header_char = (/ 'PS(u) ', 'PS(v) ',&
     91                                           'PS(w) ', 'PS(pt)', 'PS(q) ' /)
     92    CHARACTER (LEN=2),  DIMENSION(10)  ::  spectra_direction = 'x'
     93    CHARACTER (LEN=10), DIMENSION(10)  ::  data_output_sp  = ' '
     94    CHARACTER (LEN=25), DIMENSION(1:5) ::  utext_char =                    &
     95                                           (/ '-power spectrum of u     ', &
     96                                              '-power spectrum of v     ', &
     97                                              '-power spectrum of w     ', &
     98                                              '-power spectrum of ^1185 ', &
     99                                              '-power spectrum of q     ' /)
     100    CHARACTER (LEN=39), DIMENSION(1:5) ::  ytext_char =                        &
     101                                 (/ 'k ^2236 ^2566^2569<u(k) in m>2s>->2    ', &
     102                                    'k ^2236 ^2566^2569<v(k) in m>2s>->2    ', &
     103                                    'k ^2236 ^2566^2569<w(k) in m>2s>->2    ', &
     104                                    'k ^2236 ^2566^2569<^1185(k) in m>2s>->2', &
     105                                    'k ^2236 ^2566^2569<q(k) in m>2s>->2    ' /)
     106
     107    INTEGER(iwp) ::  klist_x = 0, klist_y = 0, n_sp_x = 0, n_sp_y = 0
     108
     109    INTEGER(iwp) ::  comp_spectra_level(100) = 999999,                   &
     110                     lstyles(100) = (/ 0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
     111                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
     112                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
     113                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
     114                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
     115                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
     116                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
     117                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
     118                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8,    &
     119                                       0, 7, 3, 10, 1, 4, 9, 2, 6, 8 /), &
     120                     plot_spectra_level(100) = 999999
     121
     122    REAL(wp)    ::  time_to_start_sp = 0.0_wp
     123
     124    PUBLIC  comp_spectra_level, data_output_sp, header_char, klist_x, klist_y, &
     125            lstyles, n_sp_x, n_sp_y, plot_spectra_level, spectra_direction,    &
     126            utext_char, ytext_char
     127
     128    SAVE
     129
     130    INTERFACE calc_spectra
     131       MODULE PROCEDURE calc_spectra
     132    END INTERFACE calc_spectra
     133
     134    INTERFACE preprocess_spectra
     135       MODULE PROCEDURE preprocess_spectra
     136    END INTERFACE preprocess_spectra
     137
     138    INTERFACE calc_spectra_x
     139       MODULE PROCEDURE calc_spectra_x
     140    END INTERFACE calc_spectra_x
     141
     142    INTERFACE calc_spectra_y
     143       MODULE PROCEDURE calc_spectra_y
     144    END INTERFACE calc_spectra_y
     145
     146    PUBLIC calc_spectra
     147
     148
     149 CONTAINS
     150
     151    SUBROUTINE calc_spectra
     152
     153       USE arrays_3d,                                                          &
     154           ONLY:  d, tend
     155
     156       USE control_parameters,                                                 &
     157           ONLY:  average_count_sp, bc_lr_cyc, bc_ns_cyc, message_string, psolver
     158
     159       USE cpulog,                                                             &
     160           ONLY:  cpu_log, log_point
     161
     162       USE fft_xy,                                                             &
     163           ONLY:  fft_init
     164
     165       USE indices,                                                            &
     166           ONLY:  nxl, nxr, nyn, nys, nzb, nzt
     167
     168       USE kinds
     169
     170       USE pegrid,                                                             &
     171           ONLY:  myid, pdims
     172
     173       IMPLICIT NONE
     174
     175       INTEGER(iwp) ::  m  !<
     176       INTEGER(iwp) ::  pr !<
     177
     178
     179!
     180!--    Check if user gave any levels for spectra to be calculated
     181       IF ( comp_spectra_level(1) == 999999 )  RETURN
     182
     183       CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
     184
     185!
     186!--    Initialize ffts
     187       CALL fft_init
     188
     189!
     190!--    Reallocate array d in required size
     191       IF ( psolver(1:9) == 'multigrid' )  THEN
     192          DEALLOCATE( d )
     193          ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
     194       ENDIF
     195
     196       m = 1
     197       DO WHILE ( data_output_sp(m) /= ' '  .AND.  m <= 10 )
     198!
     199!--       Transposition from z --> x  ( y --> x in case of a 1d-decomposition
     200!--       along x)
     201          IF ( INDEX( spectra_direction(m), 'x' ) /= 0 )  THEN
     202
     203!
     204!--          Calculation of spectra works for cyclic boundary conditions only
     205             IF ( .NOT. bc_lr_cyc )  THEN
     206
     207                message_string = 'non-cyclic lateral boundaries along x do'//  &
     208                                 ' not & allow calculation of spectra along x'
     209                CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
     210             ENDIF
     211
     212             CALL preprocess_spectra( m, pr )
     213
     214#if defined( __parallel )
     215             IF ( pdims(2) /= 1 )  THEN
     216                CALL resort_for_zx( d, tend )
     217                CALL transpose_zx( tend, d )
     218             ELSE
     219                CALL transpose_yxd( d, d )
     220             ENDIF
     221             CALL calc_spectra_x( d, pr, m )
     222#else
     223             message_string = 'sorry, calculation of spectra in non paral' //  &
     224                              'lel mode& is still not realized'
     225             CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
     226#endif
     227
    141228          ENDIF
    142229
    143           CALL preprocess_spectra( m, pr )
     230!
     231!--       Transposition from z --> y (d is rearranged only in case of a
     232!--       1d-decomposition along x)
     233          IF ( INDEX( spectra_direction(m), 'y' ) /= 0 )  THEN
     234
     235!
     236!--          Calculation of spectra works for cyclic boundary conditions only
     237             IF ( .NOT. bc_ns_cyc )  THEN
     238                IF ( myid == 0 )  THEN
     239                   message_string = 'non-cyclic lateral boundaries along y' // &
     240                                    ' do not & allow calculation of spectr' // &
     241                                    'a along y'
     242                   CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
     243                ENDIF
     244                CALL local_stop
     245             ENDIF
     246
     247             CALL preprocess_spectra( m, pr )
    144248
    145249#if defined( __parallel )
    146           IF ( pdims(2) /= 1 )  THEN
    147              CALL resort_for_zx( d, tend )
    148              CALL transpose_zx( tend, d )
    149           ELSE
    150              CALL transpose_yxd( d, d )
     250             CALL transpose_zyd( d, d )
     251             CALL calc_spectra_y( d, pr, m )
     252#else
     253             message_string = 'sorry, calculation of spectra in non paral' //  &
     254                              'lel mode& is still not realized'
     255             CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
     256#endif
     257
    151258          ENDIF
    152           CALL calc_spectra_x( d, pr, m )
    153 #else
    154           message_string = 'sorry, calculation of spectra in non parallel ' // &
    155                            'mode& is still not realized'
    156           CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )     
    157 #endif
    158 
    159        ENDIF
    160 
    161 !
    162 !--    Transposition from z --> y (d is rearranged only in case of a
    163 !--    1d-decomposition along x)
    164        IF ( INDEX( spectra_direction(m), 'y' ) /= 0 )  THEN
    165 
    166 !
    167 !--       Calculation of spectra works for cyclic boundary conditions only
    168           IF ( .NOT. bc_ns_cyc )  THEN
    169              IF ( myid == 0 )  THEN
    170                 message_string = 'non-cyclic lateral boundaries along y do' // &
    171                                  ' not & allow calculation of spectra along y'
    172                 CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
    173              ENDIF
    174              CALL local_stop
    175           ENDIF
    176 
    177           CALL preprocess_spectra( m, pr )
    178 
    179 #if defined( __parallel )
    180           CALL transpose_zyd( d, d )
    181           CALL calc_spectra_y( d, pr, m )
    182 #else
    183           message_string = 'sorry, calculation of spectra in non parallel' //  &
    184                            'mode& is still not realized'
    185           CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
    186 #endif
    187 
    188        ENDIF
    189 
    190 !
    191 !--    Increase counter for next spectrum
    192        m = m + 1
     259
     260!
     261!--       Increase counter for next spectrum
     262          m = m + 1
    193263         
    194     ENDDO
    195 
    196 !
    197 !-- Increase counter for averaging process in routine plot_spectra
    198     average_count_sp = average_count_sp + 1
    199 
    200     CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
    201 
    202 #endif
    203  END SUBROUTINE calc_spectra
    204 
    205 
    206 #if defined( __spectra )
     264       ENDDO
     265
     266!
     267!--    Increase counter for averaging process in routine plot_spectra
     268       average_count_sp = average_count_sp + 1
     269
     270       CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
     271
     272    END SUBROUTINE calc_spectra
     273
     274
    207275!------------------------------------------------------------------------------!
    208276! Description:
    … …  
    210278!> @todo Missing subroutine description.
    211279!------------------------------------------------------------------------------!
    212  SUBROUTINE preprocess_spectra( m, pr )
    213 
    214     USE arrays_3d,                                                             &
    215         ONLY:  d, pt, q, u, v, w
    216 
    217     USE indices,                                                               &
    218         ONLY:  ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt
    219 
    220     USE kinds
    221 
    222     USE pegrid
    223 
    224     USE spectrum,                                                              &
    225         ONLY:  data_output_sp
    226 
    227     USE statistics,                                                            &
    228         ONLY:  hom, var_d
    229 
    230 
    231     IMPLICIT NONE
    232 
    233     INTEGER(iwp) :: i  !<
    234     INTEGER(iwp) :: j  !<
    235     INTEGER(iwp) :: k  !<
    236     INTEGER(iwp) :: m  !<
    237     INTEGER(iwp) :: pr !<
    238 
    239     REAL(wp), DIMENSION(nzb:nzt+1) :: var_d_l
    240 
    241     SELECT CASE ( TRIM( data_output_sp(m) ) )
     280    SUBROUTINE preprocess_spectra( m, pr )
     281
     282       USE arrays_3d,                                                          &
     283           ONLY:  d, pt, q, u, v, w
     284
     285       USE indices,                                                            &
     286           ONLY:  ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt
     287
     288       USE kinds
     289
     290#if defined( __lc )
     291       USE MPI
     292#else
     293       INCLUDE "mpif.h"
     294#endif
     295       USE pegrid,                                                             &
     296           ONLY:  collective_wait, comm2d, ierr
     297
     298       USE statistics,                                                         &
     299           ONLY:  hom, var_d
     300
     301
     302       IMPLICIT NONE
     303
     304       INTEGER(iwp) :: i  !<
     305       INTEGER(iwp) :: j  !<
     306       INTEGER(iwp) :: k  !<
     307       INTEGER(iwp) :: m  !<
     308       INTEGER(iwp) :: pr !<
     309
     310       REAL(wp), DIMENSION(nzb:nzt+1) :: var_d_l
     311
     312       SELECT CASE ( TRIM( data_output_sp(m) ) )
    242313         
    243     CASE ( 'u' )
    244        pr = 1
    245        d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
     314       CASE ( 'u' )
     315          pr = 1
     316          d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
    246317       
    247     CASE ( 'v' )
    248        pr = 2
    249        d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
     318       CASE ( 'v' )
     319          pr = 2
     320          d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
    250321       
    251     CASE ( 'w' )
    252        pr = 3
    253        d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
     322       CASE ( 'w' )
     323          pr = 3
     324          d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
    254325       
    255     CASE ( 'pt' )
    256        pr = 4
    257        d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
     326       CASE ( 'pt' )
     327          pr = 4
     328          d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
    258329       
    259     CASE ( 'q' )
    260        pr = 41
    261        d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
     330       CASE ( 'q' )
     331          pr = 41
     332          d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
    262333       
    263     CASE DEFAULT
    264 !
    265 !--    The DEFAULT case is reached either if the parameter data_output_sp(m)
    266 !--    contains a wrong character string or if the user has coded a special
    267 !--    case in the user interface. There, the subroutine user_spectra
    268 !--    checks which of these two conditions applies.
    269        CALL user_spectra( 'preprocess', m, pr )
     334       CASE DEFAULT
     335!
     336!--       The DEFAULT case is reached either if the parameter data_output_sp(m)
     337!--       contains a wrong character string or if the user has coded a special
     338!--       case in the user interface. There, the subroutine user_spectra
     339!--       checks which of these two conditions applies.
     340          CALL user_spectra( 'preprocess', m, pr )
    270341         
    271     END SELECT
    272 
    273 !
    274 !-- Subtract horizontal mean from the array, for which spectra have to be
    275 !-- calculated
    276     var_d_l(:) = 0.0_wp
    277     DO  i = nxl, nxr
    278        DO  j = nys, nyn
    279           DO  k = nzb+1, nzt
    280              d(k,j,i)   = d(k,j,i) - hom(k,1,pr,0)
    281              var_d_l(k) = var_d_l(k) + d(k,j,i) * d(k,j,i)
     342       END SELECT
     343
     344!
     345!--    Subtract horizontal mean from the array, for which spectra have to be
     346!--    calculated
     347       var_d_l(:) = 0.0_wp
     348       DO  i = nxl, nxr
     349          DO  j = nys, nyn
     350             DO  k = nzb+1, nzt
     351                d(k,j,i)   = d(k,j,i) - hom(k,1,pr,0)
     352                var_d_l(k) = var_d_l(k) + d(k,j,i) * d(k,j,i)
     353             ENDDO
    282354          ENDDO
    283355       ENDDO
    284     ENDDO
    285 !
    286 !-- Compute total variance from local variances
    287     var_d(:) = 0.0_wp
     356!
     357!--    Compute total variance from local variances
     358       var_d(:) = 0.0_wp
    288359#if defined( __parallel )
    289     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    290     CALL MPI_ALLREDUCE( var_d_l(0), var_d(0), nzt+1-nzb, MPI_REAL,  &
    291                         MPI_SUM, comm2d, ierr )
    292 #else
    293     var_d(:) = var_d_l(:)
    294 #endif
    295     var_d(:) = var_d(:) / ngp_2dh(0)
    296 
    297  END SUBROUTINE preprocess_spectra
     360       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     361       CALL MPI_ALLREDUCE( var_d_l(0), var_d(0), nzt+1-nzb, MPI_REAL, MPI_SUM, &
     362                          comm2d, ierr )
     363#else
     364       var_d(:) = var_d_l(:)
     365#endif
     366       var_d(:) = var_d(:) / ngp_2dh(0)
     367
     368    END SUBROUTINE preprocess_spectra
    298369
    299370
    … …  
    303374!> @todo Missing subroutine description.
    304375!------------------------------------------------------------------------------!
    305  SUBROUTINE calc_spectra_x( ddd, pr, m )
    306 
    307     USE arrays_3d,                                                             &
    308         ONLY: 
    309 
    310     USE control_parameters,                                                    &
    311         ONLY:  fft_method
    312 
    313     USE fft_xy,                                                                &
    314         ONLY:  fft_x_1d
    315 
    316     USE grid_variables,                                                        &
    317         ONLY:  dx
    318 
    319     USE indices,                                                               &
    320         ONLY:  nx, ny
    321 
    322     USE kinds
    323 
    324     USE pegrid
    325 
    326     USE spectrum,                                                              &
    327         ONLY:  comp_spectra_level, n_sp_x
    328 
    329     USE statistics,                                                            &
    330         ONLY:  spectrum_x, var_d
    331 
    332     USE transpose_indices,                                                     &
    333         ONLY:  nyn_x, nys_x, nzb_x, nzt_x
    334 
    335 
    336     IMPLICIT NONE
    337 
    338     INTEGER(iwp) ::  i         !<
    339     INTEGER(iwp) ::  ishape(1) !<
    340     INTEGER(iwp) ::  j         !<
    341     INTEGER(iwp) ::  k         !<
    342     INTEGER(iwp) ::  m         !<
    343     INTEGER(iwp) ::  n         !<
    344     INTEGER(iwp) ::  pr        !<
    345 
    346     REAL(wp) ::  exponent     !<
    347     REAL(wp) ::  sum_spec_dum !< wavenumber-integrated spectrum
    348    
    349     REAL(wp), DIMENSION(0:nx) ::  work !<
    350    
    351     REAL(wp), DIMENSION(0:nx/2) ::  sums_spectra_l !<
    352    
    353     REAL(wp), DIMENSION(0:nx/2,100) ::  sums_spectra !<
    354    
    355     REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ddd !<
    356 
    357 !
    358 !-- Exponent for geometric average
    359     exponent = 1.0_wp / ( ny + 1.0_wp )
    360 
    361 !
    362 !-- Loop over all levels defined by the user
    363     n = 1
    364     DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
    365 
    366        k = comp_spectra_level(n)
    367 
    368 !
    369 !--    Calculate FFT only if the corresponding level is situated on this PE
    370        IF ( k >= nzb_x  .AND.  k <= nzt_x )  THEN
     376    SUBROUTINE calc_spectra_x( ddd, pr, m )
     377
     378       USE control_parameters,                                                 &
     379           ONLY:  fft_method
     380
     381       USE fft_xy,                                                             &
     382           ONLY:  fft_x_1d
     383
     384       USE grid_variables,                                                     &
     385           ONLY:  dx
     386
     387       USE indices,                                                            &
     388           ONLY:  nx, ny
     389
     390       USE kinds
     391
     392#if defined( __lc )
     393       USE MPI
     394#else
     395       INCLUDE "mpif.h"
     396#endif
     397       USE pegrid,                                                             &
     398           ONLY:  comm2d, ierr, myid
     399
     400       USE statistics,                                                         &
     401           ONLY:  spectrum_x, var_d
     402
     403       USE transpose_indices,                                                  &
     404           ONLY:  nyn_x, nys_x, nzb_x, nzt_x
     405
     406
     407       IMPLICIT NONE
     408
     409       INTEGER(iwp) ::  i         !<
     410       INTEGER(iwp) ::  ishape(1) !<
     411       INTEGER(iwp) ::  j         !<
     412       INTEGER(iwp) ::  k         !<
     413       INTEGER(iwp) ::  m         !<
     414       INTEGER(iwp) ::  n         !<
     415       INTEGER(iwp) ::  pr        !<
     416
     417       REAL(wp) ::  exponent     !<
     418       REAL(wp) ::  sum_spec_dum !< wavenumber-integrated spectrum
     419   
     420       REAL(wp), DIMENSION(0:nx) ::  work !<
     421   
     422       REAL(wp), DIMENSION(0:nx/2) ::  sums_spectra_l !<
     423   
     424       REAL(wp), DIMENSION(0:nx/2,100) ::  sums_spectra !<
     425   
     426       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ddd !<
     427
     428!
     429!--    Exponent for geometric average
     430       exponent = 1.0_wp / ( ny + 1.0_wp )
     431
     432!
     433!--    Loop over all levels defined by the user
     434       n = 1
     435       DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
     436
     437          k = comp_spectra_level(n)
     438
     439!
     440!--       Calculate FFT only if the corresponding level is situated on this PE
     441          IF ( k >= nzb_x  .AND.  k <= nzt_x )  THEN
    371442         
    372           DO  j = nys_x, nyn_x
    373 
    374              work = ddd(0:nx,j,k)
    375              CALL fft_x_1d( work, 'forward' )
    376 
    377              ddd(0,j,k) = dx * work(0)**2
    378              DO  i = 1, nx/2
    379                 ddd(i,j,k) = dx * ( work(i)**2 + work(nx+1-i)**2 )
     443             DO  j = nys_x, nyn_x
     444
     445                work = ddd(0:nx,j,k)
     446                CALL fft_x_1d( work, 'forward' )
     447
     448                ddd(0,j,k) = dx * work(0)**2
     449                DO  i = 1, nx/2
     450                   ddd(i,j,k) = dx * ( work(i)**2 + work(nx+1-i)**2 )
     451                ENDDO
     452
    380453             ENDDO
    381454
     455!
     456!--          Local sum and geometric average of these spectra
     457!--          (WARNING: no global sum should be performed, because floating
     458!--          point overflow may occur)
     459             DO  i = 0, nx/2
     460
     461                sums_spectra_l(i) = 1.0_wp
     462                DO  j = nys_x, nyn_x
     463                   sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
     464                ENDDO
     465
     466             ENDDO
     467         
     468          ELSE
     469
     470             sums_spectra_l = 1.0_wp
     471
     472          ENDIF
     473
     474!
     475!--       Global sum of spectra on PE0 (from where they are written on file)
     476          sums_spectra(:,n) = 0.0_wp
     477#if defined( __parallel )   
     478          CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
     479          CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1,       &
     480                           MPI_REAL, MPI_PROD, 0, comm2d, ierr )
     481#else
     482          sums_spectra(:,n) = sums_spectra_l
     483#endif
     484!
     485!--       Normalize spectra by variance
     486          sum_spec_dum = SUM( sums_spectra(:,n) )
     487          IF ( sum_spec_dum /= 0.0_wp )  THEN
     488             sums_spectra(:,n) = sums_spectra(:,n) * var_d(k) / sum_spec_dum
     489          ENDIF
     490          n = n + 1
     491
     492       ENDDO
     493       n = n - 1
     494
     495       IF ( myid == 0 )  THEN
     496!
     497!--       Sum of spectra for later averaging (see routine data_output_spectra)
     498          DO  i = 1, nx/2
     499             DO k = 1, n
     500                spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k)
     501             ENDDO
    382502          ENDDO
    383503
    384 !
    385 !--       Local sum and geometric average of these spectra
    386 !--       (WARNING: no global sum should be performed, because floating
    387 !--       point overflow may occur)
    388           DO  i = 0, nx/2
    389 
    390              sums_spectra_l(i) = 1.0_wp
    391              DO  j = nys_x, nyn_x
    392                 sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
    393              ENDDO
    394 
    395           ENDDO
    396          
    397        ELSE
    398 
    399           sums_spectra_l = 1.0_wp
    400 
    401504       ENDIF
    402 
    403 !
    404 !--    Global sum of spectra on PE0 (from where they are written on file)
    405        sums_spectra(:,n) = 0.0_wp
    406 #if defined( __parallel )   
    407        CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
    408        CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1,          &
    409                         MPI_REAL, MPI_PROD, 0, comm2d, ierr )
    410 #else
    411        sums_spectra(:,n) = sums_spectra_l
    412 #endif
    413 !
    414 !--    Normalize spectra by variance
    415        sum_spec_dum = SUM( sums_spectra(:,n) )
    416        IF ( sum_spec_dum /= 0.0_wp )  THEN
    417           sums_spectra(:,n) = sums_spectra(:,n) * var_d(k) / sum_spec_dum
    418        ENDIF
    419        n = n + 1
    420 
    421     ENDDO
    422     n = n - 1
    423 
    424     IF ( myid == 0 )  THEN
    425 !
    426 !--    Sum of spectra for later averaging (see routine data_output_spectra)
    427        DO  i = 1, nx/2
    428           DO k = 1, n
    429              spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k)
    430           ENDDO
    431        ENDDO
    432 
    433     ENDIF
    434 !
    435 !-- n_sp_x is needed by data_output_spectra_x
    436     n_sp_x = n
    437 
    438  END SUBROUTINE calc_spectra_x
     505!
     506!--    n_sp_x is needed by data_output_spectra_x
     507       n_sp_x = n
     508
     509    END SUBROUTINE calc_spectra_x
    439510
    440511
    … …  
    444515!> @todo Missing subroutine description.
    445516!------------------------------------------------------------------------------!
    446  SUBROUTINE calc_spectra_y( ddd, pr, m )
    447 
    448     USE arrays_3d,                                                             &
    449         ONLY: 
    450 
    451     USE control_parameters,                                                    &
    452         ONLY:  fft_method
    453 
    454     USE fft_xy,                                                                &
    455         ONLY:  fft_y_1d
    456 
    457     USE grid_variables,                                                        &
    458         ONLY:  dy
    459 
    460     USE indices,                                                               &
    461         ONLY:  nx, ny
    462 
    463     USE kinds
    464 
    465     USE pegrid
    466 
    467     USE spectrum,                                                              &
    468         ONLY:  comp_spectra_level, n_sp_y
    469 
    470     USE statistics,                                                            &
    471         ONLY:  spectrum_y, var_d
    472 
    473     USE transpose_indices,                                                     &
    474         ONLY:  nxl_yd, nxr_yd, nzb_yd, nzt_yd
    475 
    476 
    477     IMPLICIT NONE
    478 
    479     INTEGER(iwp) ::  i         !<
    480     INTEGER(iwp) ::  j         !<
    481     INTEGER(iwp) ::  jshape(1) !<
    482     INTEGER(iwp) ::  k         !<
    483     INTEGER(iwp) ::  m         !<
    484     INTEGER(iwp) ::  n         !<
    485     INTEGER(iwp) ::  pr        !<
    486 
    487     REAL(wp) ::  exponent !<
    488     REAL(wp) ::  sum_spec_dum !< wavenumber-integrated spectrum
    489    
    490     REAL(wp), DIMENSION(0:ny) ::  work !<
    491    
    492     REAL(wp), DIMENSION(0:ny/2) ::  sums_spectra_l !<
    493    
    494     REAL(wp), DIMENSION(0:ny/2,100) ::  sums_spectra !<
    495    
    496     REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !<
    497 
    498 
    499 !
    500 !-- Exponent for geometric average
    501     exponent = 1.0_wp / ( nx + 1.0_wp )
    502 
    503 !
    504 !-- Loop over all levels defined by the user
    505     n = 1
    506     DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
    507 
    508        k = comp_spectra_level(n)
    509 
    510 !
    511 !--    Calculate FFT only if the corresponding level is situated on this PE
    512        IF ( k >= nzb_yd  .AND.  k <= nzt_yd )  THEN
     517    SUBROUTINE calc_spectra_y( ddd, pr, m )
     518
     519       USE control_parameters,                                                 &
     520           ONLY:  fft_method
     521
     522       USE fft_xy,                                                             &
     523           ONLY:  fft_y_1d
     524
     525       USE grid_variables,                                                     &
     526           ONLY:  dy
     527
     528       USE indices,                                                            &
     529           ONLY:  nx, ny
     530
     531       USE kinds
     532
     533#if defined( __lc )
     534       USE MPI
     535#else
     536       INCLUDE "mpif.h"
     537#endif
     538       USE pegrid,                                                             &
     539           ONLY:  comm2d, ierr, myid
     540
     541       USE statistics,                                                         &
     542           ONLY:  spectrum_y, var_d
     543
     544       USE transpose_indices,                                                  &
     545           ONLY:  nxl_yd, nxr_yd, nzb_yd, nzt_yd
     546
     547
     548       IMPLICIT NONE
     549
     550       INTEGER(iwp) ::  i         !<
     551       INTEGER(iwp) ::  j         !<
     552       INTEGER(iwp) ::  jshape(1) !<
     553       INTEGER(iwp) ::  k         !<
     554       INTEGER(iwp) ::  m         !<
     555       INTEGER(iwp) ::  n         !<
     556       INTEGER(iwp) ::  pr        !<
     557
     558       REAL(wp) ::  exponent !<
     559       REAL(wp) ::  sum_spec_dum !< wavenumber-integrated spectrum
     560   
     561       REAL(wp), DIMENSION(0:ny) ::  work !<
     562   
     563       REAL(wp), DIMENSION(0:ny/2) ::  sums_spectra_l !<
     564   
     565       REAL(wp), DIMENSION(0:ny/2,100) ::  sums_spectra !<
     566   
     567       REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !<
     568
     569
     570!
     571!--    Exponent for geometric average
     572       exponent = 1.0_wp / ( nx + 1.0_wp )
     573
     574!
     575!--    Loop over all levels defined by the user
     576       n = 1
     577       DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
     578
     579          k = comp_spectra_level(n)
     580
     581!
     582!--       Calculate FFT only if the corresponding level is situated on this PE
     583          IF ( k >= nzb_yd  .AND.  k <= nzt_yd )  THEN
    513584         
    514           DO  i = nxl_yd, nxr_yd
    515 
    516              work = ddd(0:ny,i,k)
    517              CALL fft_y_1d( work, 'forward' )
    518 
    519              ddd(0,i,k) = dy * work(0)**2
    520              DO  j = 1, ny/2
    521                 ddd(j,i,k) = dy * ( work(j)**2 + work(ny+1-j)**2 )
     585             DO  i = nxl_yd, nxr_yd
     586
     587                work = ddd(0:ny,i,k)
     588                CALL fft_y_1d( work, 'forward' )
     589
     590                ddd(0,i,k) = dy * work(0)**2
     591                DO  j = 1, ny/2
     592                   ddd(j,i,k) = dy * ( work(j)**2 + work(ny+1-j)**2 )
     593                ENDDO
     594
    522595             ENDDO
    523596
     597!
     598!--          Local sum and geometric average of these spectra
     599!--          (WARNING: no global sum should be performed, because floating
     600!--          point overflow may occur)
     601             DO  j = 0, ny/2
     602
     603                sums_spectra_l(j) = 1.0_wp
     604                DO  i = nxl_yd, nxr_yd
     605                   sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
     606                ENDDO
     607
     608             ENDDO
     609         
     610          ELSE
     611
     612             sums_spectra_l = 1.0_wp
     613
     614          ENDIF
     615
     616!
     617!--       Global sum of spectra on PE0 (from where they are written on file)
     618          sums_spectra(:,n) = 0.0_wp
     619#if defined( __parallel )   
     620          CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
     621          CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1,       &
     622                           MPI_REAL, MPI_PROD, 0, comm2d, ierr )
     623#else
     624          sums_spectra(:,n) = sums_spectra_l
     625#endif
     626!
     627!--       Normalize spectra by variance
     628          sum_spec_dum = SUM( sums_spectra(:,n) )
     629          IF ( SUM(sums_spectra(:,n)) /= 0.0_wp )  THEN
     630             sums_spectra(:,n) = sums_spectra(:,n) *                           &
     631                                 var_d(k) / SUM(sums_spectra(:,n))
     632          ENDIF
     633          n = n + 1
     634
     635       ENDDO
     636       n = n - 1
     637
     638
     639       IF ( myid == 0 )  THEN
     640!
     641!--       Sum of spectra for later averaging (see routine data_output_spectra)
     642          DO  j = 1, ny/2
     643             DO k = 1, n
     644                spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k)
     645             ENDDO
    524646          ENDDO
    525647
    526 !
    527 !--       Local sum and geometric average of these spectra
    528 !--       (WARNING: no global sum should be performed, because floating
    529 !--       point overflow may occur)
    530           DO  j = 0, ny/2
    531 
    532              sums_spectra_l(j) = 1.0_wp
    533              DO  i = nxl_yd, nxr_yd
    534                 sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
    535              ENDDO
    536 
    537           ENDDO
    538          
    539        ELSE
    540 
    541           sums_spectra_l = 1.0_wp
    542 
    543648       ENDIF
    544649
    545650!
    546 !--    Global sum of spectra on PE0 (from where they are written on file)
    547        sums_spectra(:,n) = 0.0_wp
    548 #if defined( __parallel )   
    549        CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
    550        CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1,          &
    551                         MPI_REAL, MPI_PROD, 0, comm2d, ierr )
    552 #else
    553        sums_spectra(:,n) = sums_spectra_l
    554 #endif
    555 !
    556 !--    Normalize spectra by variance
    557        sum_spec_dum = SUM( sums_spectra(:,n) )
    558        IF ( SUM(sums_spectra(:,n)) /= 0.0_wp )  THEN
    559           sums_spectra(:,n) = sums_spectra(:,n) * var_d(k) / SUM(sums_spectra(:,n))
    560        ENDIF
    561        n = n + 1
    562 
    563     ENDDO
    564     n = n - 1
    565 
    566 
    567     IF ( myid == 0 )  THEN
    568 !
    569 !--    Sum of spectra for later averaging (see routine data_output_spectra)
    570        DO  j = 1, ny/2
    571           DO k = 1, n
    572              spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k)
    573           ENDDO
    574        ENDDO
    575 
    576     ENDIF
    577 
    578 !
    579 !-- n_sp_y is needed by data_output_spectra_y
    580     n_sp_y = n
    581 
    582  END SUBROUTINE calc_spectra_y
    583 #endif
     651!--    n_sp_y is needed by data_output_spectra_y
     652       n_sp_y = n
     653
     654    END SUBROUTINE calc_spectra_y
     655
     656 END MODULE spectrum

  • palm/trunk/SOURCE/time_integration.f90

    r1784 r1786  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! +module spectrum
    2222!
    2323! Former revisions:
    … …  
    271271              skip_time_do_radiation, time_radiation
    272272
     273    USE spectrum,                                                              &
     274        ONLY: calc_spectra
     275
    273276    USE statistics,                                                            &
    274277        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |