SUBROUTINE calc_spectra !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: calc_spectra.f90 392 2009-09-24 10:39:14Z maronga $ ! ! 274 2009-03-26 15:11:21Z heinze ! Output of messages replaced by message handling routine ! ! 225 2009-01-26 14:44:20Z raasch ! Bugfix: array d is reallocated in case that multigrid is used ! ! 192 2008-08-27 16:51:49Z letzel ! bugfix in calc_spectra_x: exponent = 1.0 / ( ny + 1.0 ) ! allow 100 spectra levels instead of 10 for consistency with ! define_netcdf_header ! user-defined spectra, arguments removed from transpose routines ! ! February 2007 ! RCS Log replace by Id keyword, revision history cleaned up ! ! Revision 1.9 2006/04/11 14:56:00 raasch ! pl_spectra renamed data_output_sp ! ! Revision 1.1 2001/01/05 15:08:07 raasch ! Initial revision ! ! ! Description: ! ------------ ! Calculate horizontal spectra along x and y. ! ATTENTION: 1d-decomposition along y still needs improvement, because in that ! case the gridpoint number along z still depends on the PE number ! because transpose_xz has to be used (and possibly also ! transpose_zyd needs modification). !------------------------------------------------------------------------------! #if defined( __spectra ) USE arrays_3d USE control_parameters USE cpulog USE fft_xy USE indices USE interfaces USE pegrid USE spectrum IMPLICIT NONE INTEGER :: m, pr CALL cpu_log( log_point(30), 'calc_spectra', 'start' ) ! !-- Initialize ffts CALL fft_init ! !-- Reallocate array d in required size IF ( psolver == 'multigrid' ) THEN DEALLOCATE( d ) ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra) ) ENDIF ! !-- Enlarge the size of tend, used as a working array for the transpositions IF ( nxra > nxr .OR. nyna > nyn .OR. nza > nz ) THEN DEALLOCATE( tend ) ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) ) ENDIF m = 1 DO WHILE ( data_output_sp(m) /= ' ' .AND. m <= 10 ) ! !-- Transposition from z --> x ( y --> x in case of a 1d-decomposition !-- along x) IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN ! !-- Calculation of spectra works for cyclic boundary conditions only IF ( bc_lr /= 'cyclic' ) THEN message_string = 'non-cyclic lateral boundaries along x do not'// & '& allow calculation of spectra along x' CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 ) ENDIF CALL preprocess_spectra( m, pr ) #if defined( __parallel ) IF ( pdims(2) /= 1 ) THEN CALL transpose_zx( d, tend, d ) ELSE CALL transpose_yxd( d, tend, d ) ENDIF CALL calc_spectra_x( d, pr, m ) #else message_string = 'sorry, calculation of spectra in non parallel ' // & 'mode& is still not realized' CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 ) #endif ENDIF ! !-- Transposition from z --> y (d is rearranged only in case of a !-- 1d-decomposition along x) IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN ! !-- Calculation of spectra works for cyclic boundary conditions only IF ( bc_ns /= 'cyclic' ) THEN IF ( myid == 0 ) THEN message_string = 'non-cyclic lateral boundaries along y do' // & ' not & allow calculation of spectra along y' CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 ) ENDIF CALL local_stop ENDIF CALL preprocess_spectra( m, pr ) #if defined( __parallel ) CALL transpose_zyd( d, tend, d ) CALL calc_spectra_y( d, pr, m ) #else message_string = 'sorry, calculation of spectra in non parallel' // & 'mode& is still not realized' CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 ) #endif ENDIF ! !-- Increase counter for next spectrum m = m + 1 ENDDO ! !-- Increase counter for averaging process in routine plot_spectra average_count_sp = average_count_sp + 1 ! !-- Resize tend to its normal size IF ( nxra > nxr .OR. nyna > nyn .OR. nza > nz ) THEN DEALLOCATE( tend ) ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) ENDIF CALL cpu_log( log_point(30), 'calc_spectra', 'stop' ) #endif END SUBROUTINE calc_spectra #if defined( __spectra ) SUBROUTINE preprocess_spectra( m, pr ) USE arrays_3d USE indices USE pegrid USE spectrum USE statistics IMPLICIT NONE INTEGER :: i, j, k, m, pr SELECT CASE ( TRIM( data_output_sp(m) ) ) CASE ( 'u' ) pr = 1 d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr) CASE ( 'v' ) pr = 2 d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr) CASE ( 'w' ) pr = 3 d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr) CASE ( 'pt' ) pr = 4 d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr) CASE ( 'q' ) pr = 41 d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr) CASE DEFAULT ! !-- The DEFAULT case is reached either if the parameter data_output_sp(m) !-- contains a wrong character string or if the user has coded a special !-- case in the user interface. There, the subroutine user_spectra !-- checks which of these two conditions applies. CALL user_spectra( 'preprocess', m, pr ) END SELECT ! !-- Subtract horizontal mean from the array, for which spectra have to be !-- calculated DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt d(k,j,i) = d(k,j,i) - sums(k,pr) ENDDO ENDDO ENDDO END SUBROUTINE preprocess_spectra SUBROUTINE calc_spectra_x( ddd, pr, m ) USE arrays_3d USE constants USE control_parameters USE fft_xy USE grid_variables USE indices USE pegrid USE spectrum USE statistics USE transpose_indices IMPLICIT NONE INTEGER :: i, ishape(1), j, k, m, n, pr REAL :: fac, exponent REAL, DIMENSION(0:nx) :: work REAL, DIMENSION(0:nx/2) :: sums_spectra_l REAL, DIMENSION(0:nx/2,100):: sums_spectra REAL, DIMENSION(0:nxa,nys_x:nyn_xa,nzb_x:nzt_xa) :: ddd ! !-- Exponent for geometric average exponent = 1.0 / ( ny + 1.0 ) ! !-- Loop over all levels defined by the user n = 1 DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 ) k = comp_spectra_level(n) ! !-- Calculate FFT only if the corresponding level is situated on this PE IF ( k >= nzb_x .AND. k <= nzt_x ) THEN DO j = nys_x, nyn_x work = ddd(0:nx,j,k) CALL fft_x( work, 'forward' ) ddd(0,j,k) = dx * work(0)**2 DO i = 1, nx/2 ddd(i,j,k) = dx * ( work(i)**2 + work(nx+1-i)**2 ) ENDDO ENDDO ! !-- Local sum and geometric average of these spectra !-- (WARNING: no global sum should be performed, because floating !-- point overflow may occur) DO i = 0, nx/2 sums_spectra_l(i) = 1.0 DO j = nys_x, nyn_x sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent ENDDO ENDDO ELSE sums_spectra_l = 1.0 ENDIF ! !-- Global sum of spectra on PE0 (from where they are written on file) sums_spectra(:,n) = 0.0 #if defined( __parallel ) CALL MPI_BARRIER( comm2d, ierr ) ! Necessary? CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1, & MPI_REAL, MPI_PROD, 0, comm2d, ierr ) #else sums_spectra(:,n) = sums_spectra_l #endif n = n + 1 ENDDO n = n - 1 IF ( myid == 0 ) THEN ! !-- Sum of spectra for later averaging (see routine data_output_spectra) !-- Temperton fft results need to be normalized IF ( fft_method == 'temperton-algorithm' ) THEN fac = nx + 1.0 ELSE fac = 1.0 ENDIF DO i = 1, nx/2 DO k = 1, n spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k) * fac ENDDO ENDDO ENDIF ! !-- n_sp_x is needed by data_output_spectra_x n_sp_x = n END SUBROUTINE calc_spectra_x SUBROUTINE calc_spectra_y( ddd, pr, m ) USE arrays_3d USE constants USE control_parameters USE fft_xy USE grid_variables USE indices USE pegrid USE spectrum USE statistics USE transpose_indices IMPLICIT NONE INTEGER :: i, j, jshape(1), k, m, n, pr REAL :: fac, exponent REAL, DIMENSION(0:ny) :: work REAL, DIMENSION(0:ny/2) :: sums_spectra_l REAL, DIMENSION(0:ny/2,100):: sums_spectra REAL, DIMENSION(0:nya,nxl_yd:nxr_yda,nzb_yd:nzt_yda) :: ddd ! !-- Exponent for geometric average exponent = 1.0 / ( nx + 1.0 ) ! !-- Loop over all levels defined by the user n = 1 DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 ) k = comp_spectra_level(n) ! !-- Calculate FFT only if the corresponding level is situated on this PE IF ( k >= nzb_yd .AND. k <= nzt_yd ) THEN DO i = nxl_yd, nxr_yd work = ddd(0:ny,i,k) CALL fft_y( work, 'forward' ) ddd(0,i,k) = dy * work(0)**2 DO j = 1, ny/2 ddd(j,i,k) = dy * ( work(j)**2 + work(ny+1-j)**2 ) ENDDO ENDDO ! !-- Local sum and geometric average of these spectra !-- (WARNING: no global sum should be performed, because floating !-- point overflow may occur) DO j = 0, ny/2 sums_spectra_l(j) = 1.0 DO i = nxl_yd, nxr_yd sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent ENDDO ENDDO ELSE sums_spectra_l = 1.0 ENDIF ! !-- Global sum of spectra on PE0 (from where they are written on file) sums_spectra(:,n) = 0.0 #if defined( __parallel ) CALL MPI_BARRIER( comm2d, ierr ) ! Necessary? CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1, & MPI_REAL, MPI_PROD, 0, comm2d, ierr ) #else sums_spectra(:,n) = sums_spectra_l #endif n = n + 1 ENDDO n = n - 1 IF ( myid == 0 ) THEN ! !-- Sum of spectra for later averaging (see routine data_output_spectra) !-- Temperton fft results need to be normalized IF ( fft_method == 'temperton-algorithm' ) THEN fac = ny + 1.0 ELSE fac = 1.0 ENDIF DO j = 1, ny/2 DO k = 1, n spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k) * fac ENDDO ENDDO ENDIF ! !-- n_sp_y is needed by data_output_spectra_y n_sp_y = n END SUBROUTINE calc_spectra_y #endif