source: palm/trunk/SOURCE/data_output_spectra.f90 @ 1205

 SUBROUTINE data_output_spectra

!--------------------------------------------------------------------------------!
! 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-2012  Leibniz University Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
!
!
! Former revisions:
! -----------------
! $Id: data_output_spectra.f90 1037 2012-10-22 14:10:22Z kanani $
!
! 1036 2012-10-22 13:43:42Z raasch
! code put under GPL (PALM 3.9)
!
! 964 2012-07-26 09:14:24Z raasch
! code for profil-output removed
!
! 291 2009-04-16 12:07:26Z raasch
! simulated_time in NetCDF output replaced by time_since_reference_point.
! Output of NetCDF messages with aid of message handling routine.
! Output of messages replaced by message handling routine.
!
! 189 2008-08-13 17:09:26Z letzel
! allow 100 spectra levels instead of 10 for consistency with
! define_netcdf_header, +user-defined spectra
!
! February 2007
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.7  2006/04/11 14:56:38  raasch
! pl_spectra renamed data_output_sp
!
! Revision 1.1  2001/01/05 15:14:20  raasch
! Initial revision
!
!
! Description:
! ------------
! Writing spectra data on file, using a special format which allows
! plotting of these data with PROFIL-graphic-software
!------------------------------------------------------------------------------!
#if defined( __spectra )

    USE arrays_3d
    USE control_parameters
    USE cpulog
    USE interfaces
    USE netcdf_control
    USE pegrid
    USE spectrum
    USE statistics


    IMPLICIT NONE

    INTEGER :: m, pr, cranz_x, cranz_y
    LOGICAL :: frame_x, frame_y

    CALL cpu_log( log_point(31), 'data_output_spectra', 'start' )

!
!-- Output is only performed on PE0
    IF ( myid == 0 )  THEN

!
!--    Open file for spectra output in NetCDF format
       IF ( netcdf_output )  CALL check_open( 107 )

!
!--    Increment the counter for number of output times
       dosp_time_count = dosp_time_count + 1

#if defined( __netcdf )
!
!--    Update the spectra time axis
       nc_stat = NF90_PUT_VAR( id_set_sp, id_var_time_sp,        &
                               (/ time_since_reference_point /), &
                               start = (/ dosp_time_count /), count = (/ 1 /) )
       CALL handle_netcdf_error( 'data_output_spectra', 47 )
#endif

!
!--    If necessary, calculate time average and reset average counter
       IF ( average_count_sp == 0 )  THEN
           message_string = 'no spectra data available'
           CALL message( 'data_output_spectra', 'PA0186', 0, 0, 0, 6, 0 )
       ENDIF
       IF ( average_count_sp /= 1 )  THEN
          spectrum_x = spectrum_x / REAL( average_count_sp )
          spectrum_y = spectrum_y / REAL( average_count_sp )
          average_count_sp = 0
       ENDIF

!
!--    Loop over all spectra defined by the user
       m = 1
       DO WHILE ( data_output_sp(m) /= ' '  .AND.  m <= 10 )

          SELECT CASE ( TRIM( data_output_sp(m) ) )

             CASE ( 'u' )
                pr = 1

             CASE ( 'v' )
                pr = 2

             CASE ( 'w' )
                pr = 3

             CASE ( 'pt' )
                pr = 4

             CASE ( 'q' )
                pr = 5

             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( 'data_output', m, pr )

          END SELECT

!
!--       Output of spectra in NetCDF format
          IF ( netcdf_output )  THEN
!
!--          Output of x-spectra
             IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN
                CALL output_spectra_netcdf( m, 'x' )
             ENDIF
!
!--          Output of y-spectra
             IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN
                CALL output_spectra_netcdf( m, 'y' )
             ENDIF
          ENDIF

!
!--       Increase counter for next spectrum
          m = m + 1

       ENDDO

!
!--    Reset spectra values
       spectrum_x = 0.0; spectrum_y = 0.0

    ENDIF

    CALL cpu_log( log_point(31), 'data_output_spectra', 'stop' )

#if defined( __parallel )
!    CALL MPI_BARRIER( comm2d, ierr )  ! really necessary
#endif

#endif
 END SUBROUTINE data_output_spectra


 SUBROUTINE output_spectra_netcdf( nsp, direction )
#if defined( __netcdf )

    USE constants
    USE control_parameters
    USE grid_variables
    USE indices
    USE netcdf_control
    USE spectrum
    USE statistics

    IMPLICIT NONE

    CHARACTER (LEN=1), INTENT(IN) ::  direction

    INTEGER, INTENT(IN) ::  nsp

    INTEGER ::  i, k

    REAL ::  frequency

    REAL, DIMENSION(nx/2) ::  netcdf_data_x
    REAL, DIMENSION(ny/2) ::  netcdf_data_y


    IF ( direction == 'x' )  THEN

       DO  k = 1, n_sp_x

          DO  i = 1, nx/2
             frequency = 2.0 * pi * i / ( dx * ( nx + 1 ) )
             netcdf_data_x(i) = frequency * spectrum_x(i,k,nsp)
          ENDDO

          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_dospx(nsp), netcdf_data_x, &
                                  start = (/ 1, k, dosp_time_count /), &
                                  count = (/ nx/2, 1, 1 /) )
          CALL handle_netcdf_error( 'data_output_spectra', 348 )

       ENDDO

    ENDIF

    IF ( direction == 'y' )  THEN

       DO  k = 1, n_sp_y

          DO  i = 1, ny/2
             frequency = 2.0 * pi * i / ( dy * ( ny + 1 ) )
             netcdf_data_y(i) = frequency * spectrum_y(i,k,nsp)
          ENDDO

          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_dospy(nsp), netcdf_data_y, &
                                  start = (/ 1, k, dosp_time_count /), &
                                  count = (/ ny/2, 1, 1 /) )
          CALL handle_netcdf_error( 'data_output_spectra', 349 )

       ENDDO

    ENDIF

#endif
 END SUBROUTINE output_spectra_netcdf


#if defined( __spectra )
 SUBROUTINE data_output_spectra_x( m, cranz, pr, frame_written )

    USE arrays_3d
    USE constants
    USE control_parameters
    USE grid_variables
    USE indices
    USE pegrid
    USE singleton
    USE spectrum
    USE statistics
    USE transpose_indices

    IMPLICIT NONE

    CHARACTER (LEN=30) ::  atext
    INTEGER            ::  i, j, k, m, pr
    LOGICAL            ::  frame_written
    REAL               ::  frequency = 0.0

!
!-- Variables needed for PROFIL-namelist
    INTEGER                  :: cranz, labforx = 3, labfory = 3, legpos = 3, &
                                timodex = 1
    INTEGER, DIMENSION(1:100):: cucol = 1, klist = 999999, lstyle = 0
    LOGICAL                  :: datleg = .TRUE., grid = .TRUE., &
                                lclose = .TRUE., rand = .TRUE., &
                                swap = .TRUE., twoxa = .TRUE.,  &
                                xlog = .TRUE., ylog = .TRUE.
    CHARACTER (LEN=80)       :: rtext, utext, xtext = 'k in m>->1', ytext
    REAL                     :: gwid = 0.1, rlegfak = 0.7, uxmin, uxmax, &
                                uymin, uymax
    REAL, DIMENSION(1:100)   :: lwid = 0.6
    REAL, DIMENSION(100)     :: uyma, uymi

    NAMELIST /RAHMEN/  cranz, datleg, rtext, swap
    NAMELIST /CROSS/   rand, cucol, grid, gwid, klist, labforx, labfory,      &
                       legpos, lclose, lstyle, lwid, rlegfak, timodex, utext, &
                       uxmin, uxmax, uymin, uymax, twoxa, xlog, xtext, ylog,  &
                       ytext


    rtext = '\0.5 ' // run_description_header

!
!-- Open parameter- and data-file
    CALL check_open( 81 )
    CALL check_open( 82 )

!
!-- Write file header,
!-- write RAHMEN-parameters (pr=3: w-array is on zw, other arrys on zu,
!-- pr serves as an index for output of strings (axis-labels) of the
!-- different quantities u, v, w, pt and q)
    DO  k = 1, n_sp_x
       IF ( k < 100 )  THEN
          IF ( pr == 3 )  THEN
             WRITE ( 82, 100 )  '#', k, header_char( pr ),        &
                                INT( zw(comp_spectra_level(k)) ), &
                                simulated_time_chr
          ELSE
             WRITE ( 82, 100 )  '#', k, header_char( pr ),        &
                                INT( zu(comp_spectra_level(k)) ), &
                                simulated_time_chr
          ENDIF
       ELSE
          IF ( pr == 3 )  THEN
             WRITE ( 82, 101 )  '#', k, header_char( pr ),        &
                                INT( zw(comp_spectra_level(k)) ), &
                                simulated_time_chr
          ELSE
             WRITE ( 82, 101 )  '#', k, header_char( pr ),        &
                                INT( zu(comp_spectra_level(k)) ), &
                                simulated_time_chr
          ENDIF
       ENDIF
    ENDDO

    IF ( .NOT. frame_written )  THEN
       WRITE ( 81, RAHMEN )
       frame_written = .TRUE.
    ENDIF

!
!-- Write all data and calculate uymi and uyma. They serve to calculate
!-- the CROSS-parameters uymin and uymax
    uymi = 999.999; uyma = -999.999
    DO  i = 1, nx/2
       frequency = 2.0 * pi * i / ( dx * ( nx + 1 ) )
       WRITE ( 82, 102 )  frequency, ( frequency * spectrum_x(i,k,m), k = 1, &
                          n_sp_x )
       DO  k = 1, n_sp_x
          uymi(k) = MIN( uymi(k), frequency * spectrum_x(i,k,m) )
          uyma(k) = MAX( uyma(k), frequency * spectrum_x(i,k,m) )
       ENDDO
    ENDDO

!
!-- Determine CROSS-parameters
    cucol(1:n_sp_x)  = (/ ( k, k = 1, n_sp_x ) /)
    lstyle(1:n_sp_x) = (/ ( lstyles(k), k = 1, n_sp_x ) /)

!
!-- Calculate klist-values from the available comp_spectra_level values
    i = 1; k = 1
    DO WHILE ( i <= 100  .AND.  plot_spectra_level(i) /= 999999 )
       DO WHILE ( k <= n_sp_x  .AND. &
                  plot_spectra_level(i) >= comp_spectra_level(k) )
          IF ( plot_spectra_level(i) == comp_spectra_level(k) )  THEN
             klist(i) = k + klist_x
          ELSE
             uymi(k) =  999.999
             uyma(k) = -999.999
          ENDIF
          k = k + 1
       ENDDO
       i = i + 1
    ENDDO
    uymi(k:n_sp_x) =  999.999
    uyma(k:n_sp_x) = -999.999
    utext = 'x'//utext_char( pr )
    IF ( averaging_interval_sp /= 0.0 ) THEN
       WRITE ( atext, 104 )  averaging_interval_sp
       utext = TRIM(utext) // ',  ' // TRIM( atext )
    ENDIF
    uxmin = 0.8 * 2.0 * pi        / ( dx * ( nx + 1 ) )
    uxmax = 1.2 * 2.0 * pi * nx/2 / ( dx * ( nx + 1 ) )
    uymin = 0.8 * MIN (  999.999, MINVAL ( uymi ) )
    uymax = 1.2 * MAX ( -999.999, MAXVAL ( uyma ) )
    ytext = ytext_char( pr )

!
!-- Output of CROSS-parameters
    WRITE ( 81, CROSS )

!
!-- Increase counter by the number of profiles written in the actual block
    klist_x = klist_x + n_sp_x

!
!-- Write end-mark
    WRITE ( 82, 103 )

!
!-- Close parameter- and data-file
    CALL close_file( 81 )
    CALL close_file( 82 )

!
!-- Formats
100 FORMAT (A,I1,1X,A,1X,I4,'m ',A)
101 FORMAT (A,I2,1X,A,1X,I4,'m ',A)
102 FORMAT (E15.7,100(1X,E15.7))
103 FORMAT ('NEXT')
104 FORMAT ('time averaged over',F7.1,' s')

 END SUBROUTINE data_output_spectra_x


 SUBROUTINE data_output_spectra_y( m, cranz, pr, frame_written )

    USE arrays_3d
    USE constants
    USE control_parameters
    USE grid_variables
    USE indices
    USE pegrid
    USE singleton
    USE spectrum
    USE statistics
    USE transpose_indices

    IMPLICIT NONE

    CHARACTER (LEN=30) ::  atext
    INTEGER            :: i, j, k, m, pr
    LOGICAL            :: frame_written
    REAL               :: frequency = 0.0

!
!-- Variables needed for PROFIL-namelist
    INTEGER                  :: cranz, labforx = 3, labfory = 3, legpos = 3, &
                                timodex = 1
    INTEGER, DIMENSION(1:100):: cucol = 1, klist = 999999, lstyle = 0
    LOGICAL                  :: datleg = .TRUE., grid = .TRUE., &
                                lclose = .TRUE., rand = .TRUE., &
                                swap = .TRUE., twoxa = .TRUE.,  &
                                xlog = .TRUE., ylog = .TRUE.
    CHARACTER (LEN=80)       :: rtext, utext, xtext = 'k in m>->1', ytext
    REAL                     :: gwid = 0.1, rlegfak = 0.7, uxmin, uxmax, &
                                uymin, uymax
    REAL, DIMENSION(1:100)   :: lwid = 0.6
    REAL, DIMENSION(100)     :: uyma, uymi

    NAMELIST /RAHMEN/  cranz, datleg, rtext, swap
    NAMELIST /CROSS/   rand, cucol, grid, gwid, klist, labforx, labfory,      &
                       legpos, lclose, lstyle, lwid, rlegfak, timodex, utext, &
                       uxmin, uxmax, uymin, uymax, twoxa, xlog, xtext, ylog,  &
                       ytext


    rtext = '\0.5 ' // run_description_header

!
!-- Open parameter- and data-file
    CALL check_open( 83 )
    CALL check_open( 84 )

!
!-- Write file header,
!-- write RAHMEN-parameters (pr=3: w-array is on zw, other arrys on zu,
!-- pr serves as an index for output of strings (axis-labels) of the
!-- different quantities u, v, w, pt and q)
    DO  k = 1, n_sp_y
       IF ( k < 100 )  THEN
          IF ( pr == 3 ) THEN
             WRITE ( 84, 100 )  '#', k, header_char( pr ),        &
                                INT( zw(comp_spectra_level(k)) ), &
                                simulated_time_chr
          ELSE
             WRITE ( 84, 100 )  '#', k, header_char( pr ),        &
                                INT( zu(comp_spectra_level(k)) ), &
                                simulated_time_chr
          ENDIF
       ELSE
          IF ( pr == 3 )  THEN
             WRITE ( 84, 101 )  '#', k, header_char( pr ),        &
                                INT( zw(comp_spectra_level(k)) ), &
                                simulated_time_chr
          ELSE
             WRITE ( 84, 101 )  '#', k, header_char( pr ),        &
                                INT( zu(comp_spectra_level(k)) ), &
                                simulated_time_chr
          ENDIF
       ENDIF
    ENDDO

    IF ( .NOT. frame_written )  THEN
       WRITE ( 83, RAHMEN )
       frame_written = .TRUE.
    ENDIF

!
!-- Write all data and calculate uymi and uyma. They serve to calculate
!-- the CROSS-parameters uymin and uymax
    uymi = 999.999; uyma = -999.999
    DO  j = 1, ny/2
       frequency = 2.0 * pi * j / ( dy * ( ny + 1 ) )
       WRITE ( 84, 102 ) frequency, ( frequency * spectrum_y(j,k,m), &
                                      k = 1, n_sp_y ) 
       DO k = 1, n_sp_y
          uymi(k) = MIN( uymi(k), frequency * spectrum_y(j,k,m) )
          uyma(k) = MAX( uyma(k), frequency * spectrum_y(j,k,m) )
       ENDDO
    ENDDO

!
!-- Determine CROSS-parameters
    cucol(1:n_sp_y)  = (/ ( k, k = 1, n_sp_y ) /)
    lstyle(1:n_sp_y) = (/ ( lstyles(k), k = 1, n_sp_y ) /)

!
!-- Calculate klist-values from the available comp_spectra_level values
    j = 1; k = 1
    DO WHILE ( j <= 100  .AND.  plot_spectra_level(j) /= 999999 )
       DO WHILE ( k <= n_sp_y  .AND. &
                  plot_spectra_level(j) >= comp_spectra_level(k) )
          IF ( plot_spectra_level(j) == comp_spectra_level(k) )  THEN
             klist(j) = k + klist_y
          ELSE
             uymi(k) =  999.999
             uyma(k) = -999.999
          ENDIF
          k = k + 1
       ENDDO
       j = j + 1
    ENDDO
    uymi(k:n_sp_y) =  999.999
    uyma(k:n_sp_y) = -999.999
    utext = 'y'//utext_char( pr )
    IF ( averaging_interval_sp /= 0.0 )  THEN
       WRITE ( atext, 104 )  averaging_interval_sp
       utext = TRIM(utext) // ',  ' // TRIM( atext )
    ENDIF
    uxmin = 0.8 * 2.0 * pi        / ( dy * ( ny + 1 ) )
    uxmax = 1.2 * 2.0 * pi * ny/2 / ( dy * ( ny + 1 ) )
    uymin = 0.8 * MIN (  999.999, MINVAL ( uymi ) )
    uymax = 1.2 * MAX ( -999.999, MAXVAL ( uyma ) )
    ytext = ytext_char( pr )

!
!-- Output CROSS-parameters
    WRITE ( 83, CROSS )

!
!-- Increase counter by the number of profiles written in the actual block
    klist_y = klist_y + n_sp_y

!
!-- Write end-mark
    WRITE ( 84, 103 )  

!
!-- Close parameter- and data-file
    CALL close_file( 83 )
    CALL close_file( 84 )

!
!-- Formats
100 FORMAT (A,I1,1X,A,1X,I4,'m ',A)
101 FORMAT (A,I2,1X,A,1X,I4,'m ',A)
102 FORMAT (E15.7,100(1X,E15.7))
103 FORMAT ('NEXT')
104 FORMAT ('time averaged over',F7.1,' s')

 END SUBROUTINE data_output_spectra_y
#endif