SUBROUTINE calc_spectra !------------------------------------------------------------------------------! ! Actual revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: calc_spectra.f90 4 2007-02-13 11:33:16Z knoop $ ! 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 ! !-- 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 IF ( myid == 0 ) THEN PRINT*, '+++ calc_spectra:' PRINT*, ' non-cyclic lateral boundaries along x do not ', & 'allow calculation of spectra along x' ENDIF CALL local_stop ENDIF CALL preprocess_spectra( m, pr ) #if defined( __parallel ) IF ( pdims(2) /= 1 ) THEN CALL transpose_zx( d, tend, d, tend, d ) ELSE CALL transpose_yxd( d, tend, d, tend, d ) ENDIF CALL calc_spectra_x( d, pr, m ) #else PRINT*, '+++ calc_spectra: sorry, calculation of spectra ', & 'in non parallel mode' PRINT*, ' is still not realized' CALL local_stop #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 PRINT*, '+++ calc_spectra:' PRINT*, ' non-cyclic lateral boundaries along y do not ', & 'allow calculation of spectra along y' ENDIF CALL local_stop ENDIF CALL preprocess_spectra( m, pr ) #if defined( __parallel ) CALL transpose_zyd( d, tend, d, tend, d ) CALL calc_spectra_y( d, pr, m ) #else PRINT*, '+++ calc_spectra: sorry, calculation of spectra', & 'in non parallel mode' PRINT*, ' still not realized' CALL local_stop #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 PRINT*, '+++ preprocess_spectra: Spectra of ', & TRIM( data_output_sp(m) ), ' can not be calculated' 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,10) :: sums_spectra REAL, DIMENSION(0:nxa,nys_x:nyn_xa,nzb_x:nzt_xa) :: 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 <= 10 ) 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 plot_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 plot_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,10) :: 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 <= 10 ) 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 plot_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 plot_spectra_y n_sp_y = n END SUBROUTINE calc_spectra_y #endif