Changeset 1786 for palm/trunk/SOURCE/spectrum.f90
- Timestamp:
- Mar 8, 2016 5:49:27 AM (8 years ago)
- File:
-
- 1 moved
-
palm/trunk/SOURCE/spectrum.f90 (moved) (moved from palm/trunk/SOURCE/calc_spectra.f90) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/spectrum.f90
r1785 r1786 1 !> @file calc_spectra.f901 !> @file spectrum.f90 2 2 !--------------------------------------------------------------------------------! 3 3 ! This file is part of PALM. … … 19 19 ! Current revisions: 20 20 ! ----------------- 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 22 25 ! 23 26 ! Former revisions: … … 79 82 !> transpose_zyd needs modification). 80 83 !------------------------------------------------------------------------------! 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 99 85 100 86 USE kinds 101 87 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 141 228 ENDIF 142 229 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 ) 144 248 145 249 #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 151 258 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 193 263 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 207 275 !------------------------------------------------------------------------------! 208 276 ! Description: … … 210 278 !> @todo Missing subroutine description. 211 279 !------------------------------------------------------------------------------! 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) ) ) 242 313 243 CASE ( 'u' )244 pr = 1245 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) 246 317 247 CASE ( 'v' )248 pr = 2249 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) 250 321 251 CASE ( 'w' )252 pr = 3253 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) 254 325 255 CASE ( 'pt' )256 pr = 4257 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) 258 329 259 CASE ( 'q' )260 pr = 41261 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) 262 333 263 CASE DEFAULT264 ! 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 special267 !-- case in the user interface. There, the subroutine user_spectra268 !-- 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 ) 270 341 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 282 354 ENDDO 283 355 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 288 359 #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_spectra360 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 298 369 299 370 … … 303 374 !> @todo Missing subroutine description. 304 375 !------------------------------------------------------------------------------! 305 SUBROUTINE calc_spectra_x( ddd, pr, m )306 307 USE arrays_3d,&308 ONLY:309 310 USE control_parameters,&311 ONLY: fft_method312 313 USE fft_xy,&314 ONLY: fft_x_1d315 316 USE grid_variables,&317 ONLY: dx318 319 USE indices, &320 ONLY: nx, ny 321 322 USE kinds323 324 USE pegrid325 326 USE spectrum,&327 ONLY: comp_spectra_level, n_sp_x328 329 USE statistics,&330 ONLY: spectrum_x, var_d331 332 USE transpose_indices,&333 ONLY: nyn_x, nys_x, nzb_x, nzt_x334 335 336 IMPLICIT NONE337 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 spectrum348 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 average359 exponent = 1.0_wp / ( ny + 1.0_wp )360 361 ! 362 !-- Loop over all levels defined by the user363 n = 1364 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 PE370 IF ( k >= nzb_x .AND. k <= nzt_x ) THEN376 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 371 442 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 380 453 ENDDO 381 454 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 382 502 ENDDO 383 503 384 !385 !-- Local sum and geometric average of these spectra386 !-- (WARNING: no global sum should be performed, because floating387 !-- point overflow may occur)388 DO i = 0, nx/2389 390 sums_spectra_l(i) = 1.0_wp391 DO j = nys_x, nyn_x392 sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent393 ENDDO394 395 ENDDO396 397 ELSE398 399 sums_spectra_l = 1.0_wp400 401 504 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 439 510 440 511 … … 444 515 !> @todo Missing subroutine description. 445 516 !------------------------------------------------------------------------------! 446 SUBROUTINE calc_spectra_y( ddd, pr, m )447 448 USE arrays_3d,&449 ONLY:450 451 USE control_parameters,&452 ONLY: fft_method453 454 USE fft_xy,&455 ONLY: fft_y_1d456 457 USE grid_variables,&458 ONLY: dy459 460 USE indices, &461 ONLY: nx, ny 462 463 USE kinds464 465 USE pegrid466 467 USE spectrum,&468 ONLY: comp_spectra_level, n_sp_y469 470 USE statistics,&471 ONLY: spectrum_y, var_d472 473 USE transpose_indices,&474 ONLY: nxl_yd, nxr_yd, nzb_yd, nzt_yd475 476 477 IMPLICIT NONE478 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 spectrum489 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 average501 exponent = 1.0_wp / ( nx + 1.0_wp )502 503 ! 504 !-- Loop over all levels defined by the user505 n = 1506 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 PE512 IF ( k >= nzb_yd .AND. k <= nzt_yd ) THEN517 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 513 584 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 522 595 ENDDO 523 596 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 524 646 ENDDO 525 647 526 !527 !-- Local sum and geometric average of these spectra528 !-- (WARNING: no global sum should be performed, because floating529 !-- point overflow may occur)530 DO j = 0, ny/2531 532 sums_spectra_l(j) = 1.0_wp533 DO i = nxl_yd, nxr_yd534 sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent535 ENDDO536 537 ENDDO538 539 ELSE540 541 sums_spectra_l = 1.0_wp542 543 648 ENDIF 544 649 545 650 ! 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
Note: See TracChangeset
for help on using the changeset viewer.
