Changeset 4591 for palm/trunk/SOURCE/spectra_mod.f90
- Timestamp:
- Jul 6, 2020 3:56:08 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/spectra_mod.f90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/spectra_mod.f90
r4429 r4591 1 1 !> @file spectra_mod.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 5 ! PALM is free software: you can redistribute it and/or modify it under the 6 ! terms of the GNU General Public License as published by the Free Software 7 ! Foundation, either version 3 of the License, or (at your option) any later 8 ! version. 9 ! 10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY 11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 12 ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. 13 ! 14 ! You should have received a copy of the GNU General Public License along with 15 ! PALM. If not, see <http://www.gnu.org/licenses/>. 5 ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General 6 ! Public License as published by the Free Software Foundation, either version 3 of the License, or 7 ! (at your option) any later version. 8 ! 9 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the 10 ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 ! Public License for more details. 12 ! 13 ! You should have received a copy of the GNU General Public License along with PALM. If not, see 14 ! <http://www.gnu.org/licenses/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------! 17 !--------------------------------------------------------------------------------------------------! 18 ! 19 19 ! 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 23 ! 22 ! 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- 26 26 ! $Id$ 27 ! bugfix: preprocessor directives rearranged for serial mode 28 ! 27 ! File re-formatted to follow the PALM coding standard 28 ! 29 ! 4429 2020-02-27 15:24:30Z raasch 30 ! Bugfix: preprocessor directives rearranged for serial mode 31 ! 29 32 ! 4360 2020-01-07 11:25:50Z suehring 30 33 ! Corrected "Former revisions" section 31 ! 34 ! 32 35 ! 3805 2019-03-20 15:26:35Z raasch 33 ! unused variable removed34 ! 36 ! Unused variable removed 37 ! 35 38 ! 3655 2019-01-07 16:51:22Z knoop 36 39 ! Renamed output variables … … 39 42 ! Initial revision 40 43 ! 41 ! 44 !--------------------------------------------------------------------------------------------------! 42 45 ! Description: 43 46 ! ------------ 44 !> Calculate horizontal spectra along x and y. 45 !> ATTENTION: 1d-decomposition along y still needs improvement, because in that 46 !> case the gridpoint number along z still depends on the PE number 47 !> because transpose_xz has to be used (and possibly also 48 !> transpose_zyd needs modification). 49 !------------------------------------------------------------------------------! 47 !> Calculate horizontal spectra along x and y. 48 !> ATTENTION: 1d-decomposition along y still needs improvement, because in that case the gridpoint 49 !> number along z still depends on the PE number because transpose_xz has to be used 50 !> (and possibly also transpose_zyd needs modification). 51 !--------------------------------------------------------------------------------------------------! 50 52 MODULE spectra_mod 51 53 … … 54 56 PRIVATE 55 57 56 CHARACTER (LEN=2), DIMENSION(10) :: spectra_direction = 'x' 57 CHARACTER (LEN=10), DIMENSION(10) :: data_output_sp = ' ' 58 59 INTEGER(iwp) :: average_count_sp = 0 60 INTEGER(iwp) :: dosp_time_count = 0 61 INTEGER(iwp) :: n_sp_x = 0, n_sp_y = 0 62 63 INTEGER(iwp) :: comp_spectra_level(100) = 999999 58 CHARACTER (LEN=10), DIMENSION(10) :: data_output_sp = ' ' !< 59 CHARACTER (LEN=2), DIMENSION(10) :: spectra_direction = 'x' !< 60 61 INTEGER(iwp) :: average_count_sp = 0 !< 62 INTEGER(iwp) :: comp_spectra_level(100) = 999999 !< 63 INTEGER(iwp) :: dosp_time_count = 0 !< 64 INTEGER(iwp) :: n_sp_x = 0, n_sp_y = 0 !< 64 65 65 66 LOGICAL :: calculate_spectra = .FALSE. !< internal switch that spectra are calculated … … 70 71 REAL(wp) :: skip_time_dosp = 9999999.9_wp !< no output of spectra data before this interval has passed 71 72 72 REAL(wp), DIMENSION(:), ALLOCATABLE :: var_d 73 74 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: spectrum_x, spectrum_y 73 REAL(wp), DIMENSION(:), ALLOCATABLE :: var_d !< 74 75 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: spectrum_x !< 76 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: spectrum_y !< 75 77 76 78 SAVE … … 110 112 END INTERFACE spectra_parin 111 113 112 PUBLIC average_count_sp, averaging_interval_sp, calc_spectra, & 113 calculate_spectra, comp_spectra_level, data_output_sp, & 114 dosp_time_count, dt_dosp, n_sp_x, n_sp_y, & 115 skip_time_dosp, spectra_check_parameters, spectra_direction, & 116 spectra_header, spectra_init, spectra_parin, spectrum_x, & 117 spectrum_y, var_d 114 PUBLIC average_count_sp, & 115 averaging_interval_sp, & 116 calc_spectra, & 117 calculate_spectra, & 118 comp_spectra_level, & 119 data_output_sp, & 120 dosp_time_count, & 121 dt_dosp, & 122 n_sp_x, & 123 n_sp_y, & 124 skip_time_dosp, & 125 spectra_check_parameters, & 126 spectra_direction, & 127 spectra_header, & 128 spectra_init, & 129 spectra_parin, & 130 spectrum_x, & 131 spectrum_y, & 132 var_d 118 133 119 134 120 135 CONTAINS 121 136 122 !------------------------------------------------------------------------------ !137 !--------------------------------------------------------------------------------------------------! 123 138 ! Description: 124 139 ! ------------ 125 140 !> Parin for &spectra_par for calculating spectra 126 !------------------------------------------------------------------------------! 127 SUBROUTINE spectra_parin 128 129 USE control_parameters, & 130 ONLY: dt_data_output, message_string 131 132 IMPLICIT NONE 133 134 CHARACTER (LEN=80) :: line !< dummy string that contains the current & 135 !< line of the parameter file 136 137 NAMELIST /spectra_par/ averaging_interval_sp, comp_spectra_level, & 138 data_output_sp, dt_dosp, skip_time_dosp, & 139 spectra_direction 140 141 NAMELIST /spectra_parameters/ & 142 averaging_interval_sp, comp_spectra_level, & 143 data_output_sp, dt_dosp, skip_time_dosp, & 144 spectra_direction 145 ! 146 !-- Position the namelist-file at the beginning (it was already opened in 147 !-- parin), search for the namelist-group of the package and position the 148 !-- file at this line. 149 line = ' ' 150 151 ! 152 !-- Try to find the spectra package 153 REWIND ( 11 ) 154 line = ' ' 155 DO WHILE ( INDEX( line, '&spectra_parameters' ) == 0 ) 156 READ ( 11, '(A)', END=12 ) line 157 ENDDO 158 BACKSPACE ( 11 ) 159 160 ! 161 !-- Read namelist 162 READ ( 11, spectra_parameters, ERR = 10 ) 163 164 ! 165 !-- Default setting of dt_dosp here (instead of check_parameters), because 166 !-- its current value is needed in init_pegrid 167 IF ( dt_dosp == 9999999.9_wp ) dt_dosp = dt_data_output 168 169 ! 170 !-- Set general switch that spectra shall be calculated 171 calculate_spectra = .TRUE. 172 173 GOTO 14 174 175 10 BACKSPACE( 11 ) 176 READ( 11 , '(A)') line 177 CALL parin_fail_message( 'spectra_parameters', line ) 178 ! 179 !-- Try to find the old namelist 180 12 REWIND ( 11 ) 181 line = ' ' 182 DO WHILE ( INDEX( line, '&spectra_par' ) == 0 ) 183 READ ( 11, '(A)', END=14 ) line 184 ENDDO 185 BACKSPACE ( 11 ) 186 187 ! 188 !-- Read namelist 189 READ ( 11, spectra_par, ERR = 13, END = 14 ) 190 191 192 message_string = 'namelist spectra_par is deprecated and will be ' // & 193 'removed in near future. Please use namelist ' // & 194 'spectra_parameters instead' 195 CALL message( 'spectra_parin', 'PA0487', 0, 1, 0, 6, 0 ) 196 ! 197 !-- Default setting of dt_dosp here (instead of check_parameters), because 198 !-- its current value is needed in init_pegrid 199 IF ( dt_dosp == 9999999.9_wp ) dt_dosp = dt_data_output 200 201 ! 202 !-- Set general switch that spectra shall be calculated 203 calculate_spectra = .TRUE. 204 205 GOTO 14 206 207 13 BACKSPACE( 11 ) 208 READ( 11 , '(A)') line 209 CALL parin_fail_message( 'spectra_par', line ) 210 211 212 14 CONTINUE 213 214 END SUBROUTINE spectra_parin 215 216 217 218 !------------------------------------------------------------------------------! 141 !--------------------------------------------------------------------------------------------------! 142 SUBROUTINE spectra_parin 143 144 USE control_parameters, & 145 ONLY: dt_data_output, & 146 message_string 147 148 IMPLICIT NONE 149 150 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 151 152 NAMELIST /spectra_par/ averaging_interval_sp, & 153 comp_spectra_level, & 154 data_output_sp, & 155 dt_dosp, & 156 skip_time_dosp, & 157 spectra_direction 158 159 NAMELIST /spectra_parameters/ & 160 averaging_interval_sp, & 161 comp_spectra_level, & 162 data_output_sp, & 163 dt_dosp, & 164 skip_time_dosp, & 165 spectra_direction 166 ! 167 !-- Position the namelist-file at the beginning (it was already opened in parin), search for the 168 !-- namelist-group of the package and position the file at this line. 169 line = ' ' 170 171 ! 172 !-- Try to find the spectra package 173 REWIND ( 11 ) 174 line = ' ' 175 DO WHILE ( INDEX( line, '&spectra_parameters' ) == 0 ) 176 READ ( 11, '(A)', END=12 ) line 177 ENDDO 178 BACKSPACE ( 11 ) 179 180 ! 181 !-- Read namelist 182 READ ( 11, spectra_parameters, ERR = 10 ) 183 184 ! 185 !-- Default setting of dt_dosp here (instead of check_parameters), because its current value is 186 !-- needed in init_pegrid 187 IF ( dt_dosp == 9999999.9_wp ) dt_dosp = dt_data_output 188 189 ! 190 !-- Set general switch that spectra shall be calculated 191 calculate_spectra = .TRUE. 192 193 GOTO 14 194 195 10 BACKSPACE( 11 ) 196 READ( 11 , '(A)') line 197 CALL parin_fail_message( 'spectra_parameters', line ) 198 ! 199 !-- Try to find the old namelist 200 12 REWIND ( 11 ) 201 line = ' ' 202 DO WHILE ( INDEX( line, '&spectra_par' ) == 0 ) 203 READ ( 11, '(A)', END=14 ) line 204 ENDDO 205 BACKSPACE ( 11 ) 206 207 ! 208 !-- Read namelist 209 READ ( 11, spectra_par, ERR = 13, END = 14 ) 210 211 212 message_string = 'namelist spectra_par is deprecated and will be removed in near future.' // & 213 ' Please use namelist spectra_parameters instead' 214 CALL message( 'spectra_parin', 'PA0487', 0, 1, 0, 6, 0 ) 215 ! 216 !-- Default setting of dt_dosp here (instead of check_parameters), because its current value is 217 !-- needed in init_pegrid 218 IF ( dt_dosp == 9999999.9_wp ) dt_dosp = dt_data_output 219 220 ! 221 !-- Set general switch that spectra shall be calculated 222 calculate_spectra = .TRUE. 223 224 GOTO 14 225 226 13 BACKSPACE( 11 ) 227 READ( 11 , '(A)') line 228 CALL parin_fail_message( 'spectra_par', line ) 229 230 231 14 CONTINUE 232 233 END SUBROUTINE spectra_parin 234 235 236 237 !--------------------------------------------------------------------------------------------------! 219 238 ! Description: 220 239 ! ------------ 221 240 !> Initialization of spectra related variables 222 !------------------------------------------------------------------------------! 223 SUBROUTINE spectra_init 224 225 USE indices, & 226 ONLY: nx, ny, nzb, nzt 227 228 IMPLICIT NONE 229 230 IF ( spectra_initialized ) RETURN 231 232 IF ( dt_dosp /= 9999999.9_wp ) THEN 233 234 IF ( .NOT. ALLOCATED( spectrum_x ) ) THEN 235 ALLOCATE( spectrum_x( 1:nx/2, 1:100, 1:10 ) ) 236 spectrum_x = 0.0_wp 237 ENDIF 238 239 IF ( .NOT. ALLOCATED( spectrum_y ) ) THEN 240 ALLOCATE( spectrum_y( 1:ny/2, 1:100, 1:10 ) ) 241 spectrum_y = 0.0_wp 242 ENDIF 243 244 ALLOCATE( var_d(nzb:nzt+1) ) 245 var_d = 0.0_wp 241 !--------------------------------------------------------------------------------------------------! 242 SUBROUTINE spectra_init 243 244 USE indices, & 245 ONLY: nx, & 246 ny, & 247 nzb, & 248 nzt 249 250 IMPLICIT NONE 251 252 IF ( spectra_initialized ) RETURN 253 254 IF ( dt_dosp /= 9999999.9_wp ) THEN 255 256 IF ( .NOT. ALLOCATED( spectrum_x ) ) THEN 257 ALLOCATE( spectrum_x( 1:nx/2, 1:100, 1:10 ) ) 258 spectrum_x = 0.0_wp 246 259 ENDIF 247 260 248 spectra_initialized = .TRUE. 249 250 END SUBROUTINE spectra_init 251 252 253 254 !------------------------------------------------------------------------------! 261 IF ( .NOT. ALLOCATED( spectrum_y ) ) THEN 262 ALLOCATE( spectrum_y( 1:ny/2, 1:100, 1:10 ) ) 263 spectrum_y = 0.0_wp 264 ENDIF 265 266 ALLOCATE( var_d(nzb:nzt+1) ) 267 var_d = 0.0_wp 268 ENDIF 269 270 spectra_initialized = .TRUE. 271 272 END SUBROUTINE spectra_init 273 274 275 276 !--------------------------------------------------------------------------------------------------! 255 277 ! Description: 256 278 ! ------------ 257 279 !> Check spectra related quantities 258 !------------------------------------------------------------------------------! 259 SUBROUTINE spectra_check_parameters 260 261 USE control_parameters, & 262 ONLY: averaging_interval, message_string, skip_time_data_output 263 264 IMPLICIT NONE 265 266 ! 267 !-- Check the average interval 268 IF ( averaging_interval_sp == 9999999.9_wp ) THEN 269 averaging_interval_sp = averaging_interval 270 ENDIF 271 272 IF ( averaging_interval_sp > dt_dosp ) THEN 273 WRITE( message_string, * ) 'averaging_interval_sp = ', & 274 averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp 275 CALL message( 'spectra_check_parameters', 'PA0087', 1, 2, 0, 6, 0 ) 276 ENDIF 277 278 ! 279 !-- Set the default skip time interval for data output, if necessary 280 IF ( skip_time_dosp == 9999999.9_wp ) & 281 skip_time_dosp = skip_time_data_output 282 283 END SUBROUTINE spectra_check_parameters 284 285 286 287 !------------------------------------------------------------------------------! 280 !--------------------------------------------------------------------------------------------------! 281 SUBROUTINE spectra_check_parameters 282 283 USE control_parameters, & 284 ONLY: averaging_interval, & 285 message_string, & 286 skip_time_data_output 287 288 IMPLICIT NONE 289 290 ! 291 !-- Check the average interval 292 IF ( averaging_interval_sp == 9999999.9_wp ) THEN 293 averaging_interval_sp = averaging_interval 294 ENDIF 295 296 IF ( averaging_interval_sp > dt_dosp ) THEN 297 WRITE( message_string, * ) 'averaging_interval_sp = ', averaging_interval_sp, & 298 ' must be <= dt_dosp = ', dt_dosp 299 CALL message( 'spectra_check_parameters', 'PA0087', 1, 2, 0, 6, 0 ) 300 ENDIF 301 302 ! 303 !-- Set the default skip time interval for data output, if necessary 304 IF ( skip_time_dosp == 9999999.9_wp ) skip_time_dosp = skip_time_data_output 305 306 END SUBROUTINE spectra_check_parameters 307 308 309 310 !--------------------------------------------------------------------------------------------------! 288 311 ! Description: 289 312 ! ------------ … … 291 314 !> 292 315 !> @todo Output of netcdf data format and compression level 293 !------------------------------------------------------------------------------! 294 SUBROUTINE spectra_header ( io ) 295 296 USE control_parameters, & 297 ONLY: dt_averaging_input_pr 298 299 ! USE netcdf_interface, & 300 ! ONLY: netcdf_data_format_string, netcdf_deflate 301 302 IMPLICIT NONE 303 304 ! CHARACTER (LEN=40) :: output_format !< internal string 305 306 INTEGER(iwp) :: i !< internal counter 307 INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file 308 309 ! 310 !-- Spectra output 311 IF ( dt_dosp /= 9999999.9_wp ) THEN 312 WRITE ( io, 1 ) 313 314 ! output_format = netcdf_data_format_string 315 ! IF ( netcdf_deflate == 0 ) THEN 316 ! WRITE ( io, 2 ) output_format 317 ! ELSE 318 ! WRITE ( io, 3 ) TRIM( output_format ), netcdf_deflate 319 ! ENDIF 320 WRITE ( io, 2 ) 'see profiles or other quantities' 321 WRITE ( io, 4 ) dt_dosp 322 IF ( skip_time_dosp /= 0.0_wp ) WRITE ( io, 5 ) skip_time_dosp 323 WRITE ( io, 6 ) ( data_output_sp(i), i = 1,10 ), & 324 ( spectra_direction(i), i = 1,10 ), & 325 ( comp_spectra_level(i), i = 1,100 ), & 326 averaging_interval_sp, dt_averaging_input_pr 327 ENDIF 328 329 1 FORMAT (' Spectra:') 330 2 FORMAT (' Output format: ',A/) 331 ! 3 FORMAT (' Output format: ',A, ' compressed with level: ',I1/) 332 4 FORMAT (' Output every ',F7.1,' s'/) 333 5 FORMAT (' No output during initial ',F8.2,' s') 334 6 FORMAT (' Arrays: ', 10(A5,',')/ & 335 ' Directions: ', 10(A5,',')/ & 336 ' height levels k = ', 20(I3,',')/ & 337 ' ', 20(I3,',')/ & 338 ' ', 20(I3,',')/ & 339 ' ', 20(I3,',')/ & 340 ' ', 19(I3,','),I3,'.'/ & 341 ' Time averaged over ', F7.1, ' s,' / & 342 ' Profiles for the time averaging are taken every ', & 343 F6.1,' s') 344 345 END SUBROUTINE spectra_header 346 347 348 349 SUBROUTINE calc_spectra 350 351 USE arrays_3d, & 352 ONLY: d 353 #if defined( __parallel ) 354 USE arrays_3d, & 355 ONLY: tend 356 #endif 357 358 USE control_parameters, & 359 ONLY: bc_lr_cyc, bc_ns_cyc, message_string, psolver 360 361 USE cpulog, & 362 ONLY: cpu_log, log_point 363 364 USE fft_xy, & 365 ONLY: fft_init 366 367 USE indices, & 368 ONLY: nxl, nxr, nyn, nys, nzb, nzt 369 370 USE kinds 371 372 USE pegrid, & 373 ONLY: myid 374 #if defined( __parallel ) 375 USE pegrid, & 376 ONLY: pdims 377 #endif 378 379 IMPLICIT NONE 380 381 INTEGER(iwp) :: m !< 382 INTEGER(iwp) :: pr !< 383 384 385 ! 386 !-- Check if user gave any levels for spectra to be calculated 387 IF ( comp_spectra_level(1) == 999999 ) RETURN 388 389 CALL cpu_log( log_point(30), 'calc_spectra', 'start' ) 390 391 ! 392 !-- Initialize spectra related quantities 393 CALL spectra_init 394 395 ! 396 !-- Initialize ffts 397 CALL fft_init 398 399 ! 400 !-- Reallocate array d in required size 401 IF ( psolver(1:9) == 'multigrid' ) THEN 402 DEALLOCATE( d ) 403 ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) ) 404 ENDIF 405 406 m = 1 407 DO WHILE ( data_output_sp(m) /= ' ' .AND. m <= 10 ) 408 ! 409 !-- Transposition from z --> x ( y --> x in case of a 1d-decomposition 410 !-- along x) 411 IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN 412 413 ! 414 !-- Calculation of spectra works for cyclic boundary conditions only 415 IF ( .NOT. bc_lr_cyc ) THEN 416 417 message_string = 'non-cyclic lateral boundaries along x do'// & 418 ' not & allow calculation of spectra along x' 419 CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 ) 420 ENDIF 421 422 CALL preprocess_spectra( m, pr ) 423 424 #if defined( __parallel ) 425 IF ( pdims(2) /= 1 ) THEN 426 CALL resort_for_zx( d, tend ) 427 CALL transpose_zx( tend, d ) 428 ELSE 429 CALL transpose_yxd( d, d ) 430 ENDIF 431 CALL calc_spectra_x( d, m ) 432 #else 433 message_string = 'sorry, calculation of spectra in non paral' // & 434 'lel mode& is still not realized' 435 CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 ) 436 #endif 437 438 ENDIF 439 440 ! 441 !-- Transposition from z --> y (d is rearranged only in case of a 442 !-- 1d-decomposition along x) 443 IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN 444 445 ! 446 !-- Calculation of spectra works for cyclic boundary conditions only 447 IF ( .NOT. bc_ns_cyc ) THEN 448 IF ( myid == 0 ) THEN 449 message_string = 'non-cyclic lateral boundaries along y' // & 450 ' do not & allow calculation of spectra' //& 451 ' along y' 452 CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 ) 453 ENDIF 454 CALL local_stop 455 ENDIF 456 457 CALL preprocess_spectra( m, pr ) 458 459 #if defined( __parallel ) 460 CALL transpose_zyd( d, d ) 461 CALL calc_spectra_y( d, m ) 462 #else 463 message_string = 'sorry, calculation of spectra in non paral' // & 464 'lel mode& is still not realized' 465 CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 ) 466 #endif 467 468 ENDIF 469 470 ! 471 !-- Increase counter for next spectrum 472 m = m + 1 473 474 ENDDO 475 476 ! 477 !-- Increase counter for averaging process in routine plot_spectra 478 average_count_sp = average_count_sp + 1 479 480 CALL cpu_log( log_point(30), 'calc_spectra', 'stop' ) 481 482 END SUBROUTINE calc_spectra 483 484 485 !------------------------------------------------------------------------------! 316 !--------------------------------------------------------------------------------------------------! 317 SUBROUTINE spectra_header ( io ) 318 319 USE control_parameters, & 320 ONLY: dt_averaging_input_pr 321 322 ! USE netcdf_interface, & 323 ! ONLY: netcdf_data_format_string, & 324 ! netcdf_deflate 325 326 IMPLICIT NONE 327 328 ! CHARACTER (LEN=40) :: output_format !< internal string 329 330 INTEGER(iwp) :: i !< internal counter 331 INTEGER(iwp), INTENT(IN) :: io !< unit of the output file 332 333 ! 334 !-- Spectra output 335 IF ( dt_dosp /= 9999999.9_wp ) THEN 336 WRITE ( io, 1 ) 337 338 ! output_format = netcdf_data_format_string 339 ! IF ( netcdf_deflate == 0 ) THEN 340 ! WRITE ( io, 2 ) output_format 341 ! ELSE 342 ! WRITE ( io, 3 ) TRIM( output_format ), netcdf_deflate 343 ! ENDIF 344 WRITE ( io, 2 ) 'see profiles or other quantities' 345 WRITE ( io, 4 ) dt_dosp 346 IF ( skip_time_dosp /= 0.0_wp ) WRITE ( io, 5 ) skip_time_dosp 347 WRITE ( io, 6 ) ( data_output_sp(i), i = 1,10 ), & 348 ( spectra_direction(i), i = 1,10 ), & 349 ( comp_spectra_level(i), i = 1,100 ), & 350 averaging_interval_sp, dt_averaging_input_pr 351 ENDIF 352 353 1 FORMAT (' Spectra:') 354 2 FORMAT (' Output format: ',A/) 355 ! 3 FORMAT (' Output format: ',A, ' compressed with level: ',I1/) 356 4 FORMAT (' Output every ',F7.1,' s'/) 357 5 FORMAT (' No output during initial ',F8.2,' s') 358 6 FORMAT (' Arrays: ', 10(A5,',')/ & 359 ' Directions: ', 10(A5,',')/ & 360 ' height levels k = ', 20(I3,',')/ & 361 ' ', 20(I3,',')/ & 362 ' ', 20(I3,',')/ & 363 ' ', 20(I3,',')/ & 364 ' ', 19(I3,','),I3,'.'/ & 365 ' Time averaged over ', F7.1, ' s,' / & 366 ' Profiles for the time averaging are taken every ', & 367 F6.1,' s') 368 369 END SUBROUTINE spectra_header 370 371 372 !--------------------------------------------------------------------------------------------------! 486 373 ! Description: 487 374 ! ------------ 488 375 !> @todo Missing subroutine description. 489 !------------------------------------------------------------------------------! 490 SUBROUTINE preprocess_spectra( m, pr ) 491 492 USE arrays_3d, & 493 ONLY: d, pt, q, s, u, v, w 494 495 USE indices, & 496 ONLY: ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt 497 498 USE kinds 499 500 #if defined( __parallel ) 501 #if !defined( __mpifh ) 502 USE MPI 503 #endif 504 505 USE pegrid, & 506 ONLY: collective_wait, comm2d, ierr 507 #endif 508 509 USE statistics, & 510 ONLY: hom 511 512 513 IMPLICIT NONE 514 515 #if defined( __parallel ) 516 #if defined( __mpifh ) 517 INCLUDE "mpif.h" 518 #endif 519 #endif 520 521 INTEGER(iwp) :: i !< 522 INTEGER(iwp) :: j !< 523 INTEGER(iwp) :: k !< 524 INTEGER(iwp) :: m !< 525 INTEGER(iwp) :: pr !< 526 527 REAL(wp), DIMENSION(nzb:nzt+1) :: var_d_l 528 529 SELECT CASE ( TRIM( data_output_sp(m) ) ) 530 531 CASE ( 'u' ) 532 pr = 1 533 d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr) 534 535 CASE ( 'v' ) 536 pr = 2 537 d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr) 538 539 CASE ( 'w' ) 540 pr = 3 541 d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr) 542 543 CASE ( 'theta' ) 544 pr = 4 545 d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr) 546 547 CASE ( 'q' ) 548 pr = 41 549 d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr) 550 551 CASE ( 's' ) 552 pr = 117 553 d(nzb+1:nzt,nys:nyn,nxl:nxr) = s(nzb+1:nzt,nys:nyn,nxl:nxr) 554 555 CASE DEFAULT 556 ! 557 !-- The DEFAULT case is reached either if the parameter data_output_sp(m) 558 !-- contains a wrong character string or if the user has coded a special 559 !-- case in the user interface. There, the subroutine user_spectra 560 !-- checks which of these two conditions applies. 561 CALL user_spectra( 'preprocess', m, pr ) 562 563 END SELECT 564 565 ! 566 !-- Subtract horizontal mean from the array, for which spectra have to be 567 !-- calculated. Moreover, calculate variance of the respective quantitiy, 568 !-- later used for normalizing spectra output. 569 var_d_l(:) = 0.0_wp 570 DO i = nxl, nxr 571 DO j = nys, nyn 572 DO k = nzb+1, nzt 573 d(k,j,i) = d(k,j,i) - hom(k,1,pr,0) 574 var_d_l(k) = var_d_l(k) + d(k,j,i) * d(k,j,i) 575 ENDDO 576 ENDDO 577 ENDDO 578 ! 579 !-- Compute total variance from local variances 580 var_d(:) = 0.0_wp 581 #if defined( __parallel ) 582 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 583 CALL MPI_ALLREDUCE( var_d_l(0), var_d(0), nzt+1-nzb, MPI_REAL, MPI_SUM, & 584 comm2d, ierr ) 376 !--------------------------------------------------------------------------------------------------! 377 378 SUBROUTINE calc_spectra 379 380 USE arrays_3d, & 381 ONLY: d 382 #if defined( __parallel ) 383 USE arrays_3d, & 384 ONLY: tend 385 #endif 386 387 USE control_parameters, & 388 ONLY: bc_lr_cyc, & 389 bc_ns_cyc, & 390 message_string, & 391 psolver 392 393 USE cpulog, & 394 ONLY: cpu_log, & 395 log_point 396 397 USE fft_xy, & 398 ONLY: fft_init 399 400 USE indices, & 401 ONLY: nxl, & 402 nxr, & 403 nyn, & 404 nys, & 405 nzb, & 406 nzt 407 408 USE kinds 409 410 USE pegrid, & 411 ONLY: myid 412 #if defined( __parallel ) 413 USE pegrid, & 414 ONLY: pdims 415 #endif 416 417 IMPLICIT NONE 418 419 INTEGER(iwp) :: m !< 420 INTEGER(iwp) :: pr !< 421 422 423 ! 424 !-- Check if user gave any levels for spectra to be calculated 425 IF ( comp_spectra_level(1) == 999999 ) RETURN 426 427 CALL cpu_log( log_point(30), 'calc_spectra', 'start' ) 428 429 ! 430 !-- Initialize spectra related quantities 431 CALL spectra_init 432 433 ! 434 !-- Initialize ffts 435 CALL fft_init 436 437 ! 438 !-- Reallocate array d in required size 439 IF ( psolver(1:9) == 'multigrid' ) THEN 440 DEALLOCATE( d ) 441 ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) ) 442 ENDIF 443 444 m = 1 445 DO WHILE ( data_output_sp(m) /= ' ' .AND. m <= 10 ) 446 ! 447 !-- Transposition from z --> x ( y --> x in case of a 1d-decomposition along x) 448 IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN 449 450 ! 451 !-- Calculation of spectra works for cyclic boundary conditions only 452 IF ( .NOT. bc_lr_cyc ) THEN 453 454 message_string = 'non-cyclic lateral boundaries along x do' // & 455 ' not & allow calculation of spectra along x' 456 CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 ) 457 ENDIF 458 459 CALL preprocess_spectra( m, pr ) 460 461 #if defined( __parallel ) 462 IF ( pdims(2) /= 1 ) THEN 463 CALL resort_for_zx( d, tend ) 464 CALL transpose_zx( tend, d ) 465 ELSE 466 CALL transpose_yxd( d, d ) 467 ENDIF 468 CALL calc_spectra_x( d, m ) 585 469 #else 586 var_d(:) = var_d_l(:) 587 #endif 588 var_d(:) = var_d(:) / ngp_2dh(0) 589 590 END SUBROUTINE preprocess_spectra 591 592 593 !------------------------------------------------------------------------------! 470 message_string = 'sorry, calculation of spectra in non parallel mode& is still not ' // & 471 'realized' 472 CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 ) 473 #endif 474 475 ENDIF 476 477 ! 478 !-- Transposition from z --> y (d is rearranged only in case of a 1d-decomposition along x) 479 IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN 480 481 ! 482 !-- Calculation of spectra works for cyclic boundary conditions only 483 IF ( .NOT. bc_ns_cyc ) THEN 484 IF ( myid == 0 ) THEN 485 message_string = 'non-cyclic lateral boundaries along y' // & 486 ' do not & allow calculation of spectra along y' 487 CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 ) 488 ENDIF 489 CALL local_stop 490 ENDIF 491 492 CALL preprocess_spectra( m, pr ) 493 494 #if defined( __parallel ) 495 CALL transpose_zyd( d, d ) 496 CALL calc_spectra_y( d, m ) 497 #else 498 message_string = 'sorry, calculation of spectra in non parallel mode& is still not ' // & 499 'realized' 500 CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 ) 501 #endif 502 503 ENDIF 504 505 ! 506 !-- Increase counter for next spectrum 507 m = m + 1 508 509 ENDDO 510 511 ! 512 !-- Increase counter for averaging process in routine plot_spectra 513 average_count_sp = average_count_sp + 1 514 515 CALL cpu_log( log_point(30), 'calc_spectra', 'stop' ) 516 517 END SUBROUTINE calc_spectra 518 519 520 !--------------------------------------------------------------------------------------------------! 594 521 ! Description: 595 522 ! ------------ 596 523 !> @todo Missing subroutine description. 597 !------------------------------------------------------------------------------! 598 #if defined( __parallel ) 599 SUBROUTINE calc_spectra_x( ddd, m ) 600 601 USE fft_xy, & 602 ONLY: fft_x_1d 603 604 USE grid_variables, & 605 ONLY: dx 606 607 USE indices, & 608 ONLY: nx, ny 609 610 USE kinds 611 524 !--------------------------------------------------------------------------------------------------! 525 SUBROUTINE preprocess_spectra( m, pr ) 526 527 USE arrays_3d, & 528 ONLY: d, & 529 pt, & 530 q, & 531 s, & 532 u, & 533 v, & 534 w 535 536 USE indices, & 537 ONLY: ngp_2dh, & 538 nxl, & 539 nxr, & 540 nyn, & 541 nys, & 542 nzb, & 543 nzt 544 545 USE kinds 546 547 #if defined( __parallel ) 612 548 #if !defined( __mpifh ) 613 USE MPI 614 #endif 615 616 USE pegrid, & 617 ONLY: comm2d, ierr, myid 618 619 USE transpose_indices, & 620 ONLY: nyn_x, nys_x, nzb_x, nzt_x 621 622 623 IMPLICIT NONE 624 549 USE MPI 550 #endif 551 552 USE pegrid, & 553 ONLY: collective_wait, & 554 comm2d, & 555 ierr 556 #endif 557 558 USE statistics, & 559 ONLY: hom 560 561 562 IMPLICIT NONE 563 564 #if defined( __parallel ) 625 565 #if defined( __mpifh ) 626 INCLUDE "mpif.h" 627 #endif 628 629 INTEGER(iwp) :: i !< 630 INTEGER(iwp) :: j !< 631 INTEGER(iwp) :: k !< 632 INTEGER(iwp) :: m !< 633 INTEGER(iwp) :: n !< 634 635 REAL(wp) :: exponent !< 636 REAL(wp) :: sum_spec_dum !< wavenumber-integrated spectrum 637 638 REAL(wp), DIMENSION(0:nx) :: work !< 639 640 REAL(wp), DIMENSION(0:nx/2) :: sums_spectra_l !< 641 642 REAL(wp), DIMENSION(0:nx/2,100) :: sums_spectra !< 643 644 REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ddd !< 645 646 ! 647 !-- Exponent for geometric average 648 exponent = 1.0_wp / ( ny + 1.0_wp ) 649 650 ! 651 !-- Loop over all levels defined by the user 652 n = 1 653 DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 ) 654 655 k = comp_spectra_level(n) 656 657 ! 658 !-- Calculate FFT only if the corresponding level is situated on this PE 659 IF ( k >= nzb_x .AND. k <= nzt_x ) THEN 660 661 DO j = nys_x, nyn_x 662 663 work = ddd(0:nx,j,k) 664 CALL fft_x_1d( work, 'forward' ) 665 666 ddd(0,j,k) = dx * work(0)**2 667 DO i = 1, nx/2 668 ddd(i,j,k) = dx * ( work(i)**2 + work(nx+1-i)**2 ) 669 ENDDO 670 671 ENDDO 672 673 ! 674 !-- Local sum and geometric average of these spectra 675 !-- (WARNING: no global sum should be performed, because floating 676 !-- point overflow may occur) 677 DO i = 0, nx/2 678 679 sums_spectra_l(i) = 1.0_wp 680 DO j = nys_x, nyn_x 681 sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent 682 ENDDO 683 684 ENDDO 685 686 ELSE 687 688 sums_spectra_l = 1.0_wp 689 690 ENDIF 691 692 ! 693 !-- Global sum of spectra on PE0 (from where they are written on file) 694 sums_spectra(:,n) = 0.0_wp 695 #if defined( __parallel ) 696 CALL MPI_BARRIER( comm2d, ierr ) ! Necessary? 697 CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1, & 698 MPI_REAL, MPI_PROD, 0, comm2d, ierr ) 566 INCLUDE "mpif.h" 567 #endif 568 #endif 569 570 INTEGER(iwp) :: i !< 571 INTEGER(iwp) :: j !< 572 INTEGER(iwp) :: k !< 573 INTEGER(iwp) :: m !< 574 INTEGER(iwp) :: pr !< 575 576 REAL(wp), DIMENSION(nzb:nzt+1) :: var_d_l !< 577 578 SELECT CASE ( TRIM( data_output_sp(m) ) ) 579 580 CASE ( 'u' ) 581 pr = 1 582 d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr) 583 584 CASE ( 'v' ) 585 pr = 2 586 d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr) 587 588 CASE ( 'w' ) 589 pr = 3 590 d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr) 591 592 CASE ( 'theta' ) 593 pr = 4 594 d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr) 595 596 CASE ( 'q' ) 597 pr = 41 598 d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr) 599 600 CASE ( 's' ) 601 pr = 117 602 d(nzb+1:nzt,nys:nyn,nxl:nxr) = s(nzb+1:nzt,nys:nyn,nxl:nxr) 603 604 CASE DEFAULT 605 ! 606 !-- The DEFAULT case is reached either if the parameter data_output_sp(m) contains a wrong 607 !-- character string or if the user has coded a special case in the user interface. There, the 608 !-- subroutine user_spectra checks which of these two conditions applies. 609 CALL user_spectra( 'preprocess', m, pr ) 610 611 END SELECT 612 613 ! 614 !-- Subtract horizontal mean from the array, for which spectra have to be calculated. Moreover, 615 !-- calculate variance of the respective quantitiy, later used for normalizing spectra output. 616 var_d_l(:) = 0.0_wp 617 DO i = nxl, nxr 618 DO j = nys, nyn 619 DO k = nzb+1, nzt 620 d(k,j,i) = d(k,j,i) - hom(k,1,pr,0) 621 var_d_l(k) = var_d_l(k) + d(k,j,i) * d(k,j,i) 622 ENDDO 623 ENDDO 624 ENDDO 625 ! 626 !-- Compute total variance from local variances 627 var_d(:) = 0.0_wp 628 #if defined( __parallel ) 629 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 630 CALL MPI_ALLREDUCE( var_d_l(0), var_d(0), nzt+1-nzb, MPI_REAL, MPI_SUM, comm2d, ierr ) 699 631 #else 700 sums_spectra(:,n) = sums_spectra_l 701 #endif 702 ! 703 !-- Normalize spectra by variance 704 sum_spec_dum = SUM( sums_spectra(1:nx/2,n) ) 705 706 IF ( sum_spec_dum /= 0.0_wp ) THEN 707 sums_spectra(1:nx/2,n) = sums_spectra(1:nx/2,n) * & 708 var_d(k) / sum_spec_dum 709 ENDIF 710 n = n + 1 711 712 ENDDO 713 n = n - 1 714 715 IF ( myid == 0 ) THEN 716 ! 717 !-- Sum of spectra for later averaging (see routine data_output_spectra) 718 DO i = 1, nx/2 719 DO k = 1, n 720 spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k) 721 ENDDO 722 ENDDO 723 724 ENDIF 725 ! 726 !-- n_sp_x is needed by data_output_spectra_x 727 n_sp_x = n 728 729 END SUBROUTINE calc_spectra_x 730 #endif 731 732 733 !------------------------------------------------------------------------------! 632 var_d(:) = var_d_l(:) 633 #endif 634 var_d(:) = var_d(:) / ngp_2dh(0) 635 636 END SUBROUTINE preprocess_spectra 637 638 639 !--------------------------------------------------------------------------------------------------! 734 640 ! Description: 735 641 ! ------------ 736 642 !> @todo Missing subroutine description. 737 !------------------------------------------------------------------------------! 738 #if defined( __parallel ) 739 SUBROUTINE calc_spectra_y( ddd, m ) 740 741 USE fft_xy, & 742 ONLY: fft_y_1d 743 744 USE grid_variables, & 745 ONLY: dy 746 747 USE indices, & 748 ONLY: nx, ny 749 750 USE kinds 643 !--------------------------------------------------------------------------------------------------! 644 #if defined( __parallel ) 645 SUBROUTINE calc_spectra_x( ddd, m ) 646 647 USE fft_xy, & 648 ONLY: fft_x_1d 649 650 USE grid_variables, & 651 ONLY: dx 652 653 USE indices, & 654 ONLY: nx, & 655 ny 656 657 USE kinds 751 658 752 659 #if !defined( __mpifh ) 753 USE MPI 754 #endif 755 756 USE pegrid, & 757 ONLY: comm2d, ierr, myid 758 759 USE transpose_indices, & 760 ONLY: nxl_yd, nxr_yd, nzb_yd, nzt_yd 761 762 763 IMPLICIT NONE 660 USE MPI 661 #endif 662 663 USE pegrid, & 664 ONLY: comm2d, & 665 ierr, & 666 myid 667 668 USE transpose_indices, & 669 ONLY: nyn_x, & 670 nys_x, & 671 nzb_x, & 672 nzt_x 673 674 675 IMPLICIT NONE 764 676 765 677 #if defined( __mpifh ) 766 INCLUDE "mpif.h" 767 #endif 768 769 INTEGER(iwp) :: i !< 770 INTEGER(iwp) :: j !< 771 INTEGER(iwp) :: k !< 772 INTEGER(iwp) :: m !< 773 INTEGER(iwp) :: n !< 774 775 REAL(wp) :: exponent !< 776 REAL(wp) :: sum_spec_dum !< wavenumber-integrated spectrum 777 778 REAL(wp), DIMENSION(0:ny) :: work !< 779 780 REAL(wp), DIMENSION(0:ny/2) :: sums_spectra_l !< 781 782 REAL(wp), DIMENSION(0:ny/2,100) :: sums_spectra !< 783 784 REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !< 785 786 787 ! 788 !-- Exponent for geometric average 789 exponent = 1.0_wp / ( nx + 1.0_wp ) 790 791 ! 792 !-- Loop over all levels defined by the user 793 n = 1 794 DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 ) 795 796 k = comp_spectra_level(n) 797 798 ! 799 !-- Calculate FFT only if the corresponding level is situated on this PE 800 IF ( k >= nzb_yd .AND. k <= nzt_yd ) THEN 801 678 INCLUDE "mpif.h" 679 #endif 680 681 INTEGER(iwp) :: i !< 682 INTEGER(iwp) :: j !< 683 INTEGER(iwp) :: k !< 684 INTEGER(iwp) :: m !< 685 INTEGER(iwp) :: n !< 686 687 REAL(wp) :: exponent !< 688 REAL(wp) :: sum_spec_dum !< wavenumber-integrated spectrum 689 690 REAL(wp), DIMENSION(0:nx) :: work !< 691 692 REAL(wp), DIMENSION(0:nx/2) :: sums_spectra_l !< 693 694 REAL(wp), DIMENSION(0:nx/2,100) :: sums_spectra !< 695 696 REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ddd !< 697 698 ! 699 !-- Exponent for geometric average 700 exponent = 1.0_wp / ( ny + 1.0_wp ) 701 702 ! 703 !-- Loop over all levels defined by the user 704 n = 1 705 DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 ) 706 707 k = comp_spectra_level(n) 708 709 ! 710 !-- Calculate FFT only if the corresponding level is situated on this PE 711 IF ( k >= nzb_x .AND. k <= nzt_x ) THEN 712 713 DO j = nys_x, nyn_x 714 715 work = ddd(0:nx,j,k) 716 CALL fft_x_1d( work, 'forward' ) 717 718 ddd(0,j,k) = dx * work(0)**2 719 DO i = 1, nx/2 720 ddd(i,j,k) = dx * ( work(i)**2 + work(nx+1-i)**2 ) 721 ENDDO 722 723 ENDDO 724 725 ! 726 !-- Local sum and geometric average of these spectra (WARNING: no global sum should be 727 !-- performed, because floating point overflow may occur) 728 DO i = 0, nx/2 729 730 sums_spectra_l(i) = 1.0_wp 731 DO j = nys_x, nyn_x 732 sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent 733 ENDDO 734 735 ENDDO 736 737 ELSE 738 739 sums_spectra_l = 1.0_wp 740 741 ENDIF 742 743 ! 744 !-- Global sum of spectra on PE0 (from where they are written on file) 745 sums_spectra(:,n) = 0.0_wp 746 #if defined( __parallel ) 747 CALL MPI_BARRIER( comm2d, ierr ) ! Necessary? 748 CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1, MPI_REAL, MPI_PROD, 0, & 749 comm2d, ierr ) 750 #else 751 sums_spectra(:,n) = sums_spectra_l 752 #endif 753 ! 754 !-- Normalize spectra by variance 755 sum_spec_dum = SUM( sums_spectra(1:nx/2,n) ) 756 757 IF ( sum_spec_dum /= 0.0_wp ) THEN 758 sums_spectra(1:nx/2,n) = sums_spectra(1:nx/2,n) * var_d(k) / sum_spec_dum 759 ENDIF 760 n = n + 1 761 762 ENDDO 763 n = n - 1 764 765 IF ( myid == 0 ) THEN 766 ! 767 !-- Sum of spectra for later averaging (see routine data_output_spectra) 768 DO i = 1, nx/2 769 DO k = 1, n 770 spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k) 771 ENDDO 772 ENDDO 773 774 ENDIF 775 ! 776 !-- n_sp_x is needed by data_output_spectra_x 777 n_sp_x = n 778 779 END SUBROUTINE calc_spectra_x 780 #endif 781 782 783 !--------------------------------------------------------------------------------------------------! 784 ! Description: 785 ! ------------ 786 !> @todo Missing subroutine description. 787 !--------------------------------------------------------------------------------------------------! 788 #if defined( __parallel ) 789 SUBROUTINE calc_spectra_y( ddd, m ) 790 791 USE fft_xy, & 792 ONLY: fft_y_1d 793 794 USE grid_variables, & 795 ONLY: dy 796 797 USE indices, & 798 ONLY: nx, & 799 ny 800 801 USE kinds 802 803 #if !defined( __mpifh ) 804 USE MPI 805 #endif 806 807 USE pegrid, & 808 ONLY: comm2d, & 809 ierr, & 810 myid 811 812 USE transpose_indices, & 813 ONLY: nxl_yd, & 814 nxr_yd, & 815 nzb_yd, & 816 nzt_yd 817 818 819 IMPLICIT NONE 820 821 #if defined( __mpifh ) 822 INCLUDE "mpif.h" 823 #endif 824 825 INTEGER(iwp) :: i !< 826 INTEGER(iwp) :: j !< 827 INTEGER(iwp) :: k !< 828 INTEGER(iwp) :: m !< 829 INTEGER(iwp) :: n !< 830 831 REAL(wp) :: exponent !< 832 REAL(wp) :: sum_spec_dum !< wavenumber-integrated spectrum 833 834 REAL(wp), DIMENSION(0:ny) :: work !< 835 836 REAL(wp), DIMENSION(0:ny/2) :: sums_spectra_l !< 837 838 REAL(wp), DIMENSION(0:ny/2,100) :: sums_spectra !< 839 840 REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !< 841 842 843 ! 844 !-- Exponent for geometric average 845 exponent = 1.0_wp / ( nx + 1.0_wp ) 846 847 ! 848 !-- Loop over all levels defined by the user 849 n = 1 850 DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 ) 851 852 k = comp_spectra_level(n) 853 854 ! 855 !-- Calculate FFT only if the corresponding level is situated on this PE 856 IF ( k >= nzb_yd .AND. k <= nzt_yd ) THEN 857 858 DO i = nxl_yd, nxr_yd 859 860 work = ddd(0:ny,i,k) 861 CALL fft_y_1d( work, 'forward' ) 862 863 ddd(0,i,k) = dy * work(0)**2 864 DO j = 1, ny/2 865 ddd(j,i,k) = dy * ( work(j)**2 + work(ny+1-j)**2 ) 866 ENDDO 867 868 ENDDO 869 870 ! 871 !-- Local sum and geometric average of these spectra (WARNING: no global sum should be 872 !-- performed, because floating point overflow may occur) 873 DO j = 0, ny/2 874 875 sums_spectra_l(j) = 1.0_wp 802 876 DO i = nxl_yd, nxr_yd 803 804 work = ddd(0:ny,i,k) 805 CALL fft_y_1d( work, 'forward' ) 806 807 ddd(0,i,k) = dy * work(0)**2 808 DO j = 1, ny/2 809 ddd(j,i,k) = dy * ( work(j)**2 + work(ny+1-j)**2 ) 810 ENDDO 811 877 sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent 812 878 ENDDO 813 879 814 ! 815 !-- Local sum and geometric average of these spectra 816 !-- (WARNING: no global sum should be performed, because floating 817 !-- point overflow may occur) 818 DO j = 0, ny/2 819 820 sums_spectra_l(j) = 1.0_wp 821 DO i = nxl_yd, nxr_yd 822 sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent 823 ENDDO 824 825 ENDDO 826 827 ELSE 828 829 sums_spectra_l = 1.0_wp 830 831 ENDIF 832 833 ! 834 !-- Global sum of spectra on PE0 (from where they are written on file) 835 sums_spectra(:,n) = 0.0_wp 836 #if defined( __parallel ) 837 CALL MPI_BARRIER( comm2d, ierr ) ! Necessary? 838 CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1, & 839 MPI_REAL, MPI_PROD, 0, comm2d, ierr ) 880 ENDDO 881 882 ELSE 883 884 sums_spectra_l = 1.0_wp 885 886 ENDIF 887 888 ! 889 !-- Global sum of spectra on PE0 (from where they are written on file) 890 sums_spectra(:,n) = 0.0_wp 891 #if defined( __parallel ) 892 CALL MPI_BARRIER( comm2d, ierr ) ! Necessary? 893 CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1, MPI_REAL, MPI_PROD, 0, & 894 comm2d, ierr ) 840 895 #else 841 sums_spectra(:,n) = sums_spectra_l 842 #endif 843 ! 844 !-- Normalize spectra by variance 845 sum_spec_dum = SUM( sums_spectra(1:ny/2,n) ) 846 IF ( sum_spec_dum /= 0.0_wp ) THEN 847 sums_spectra(1:ny/2,n) = sums_spectra(1:ny/2,n) * & 848 var_d(k) / sum_spec_dum 849 ENDIF 850 n = n + 1 851 896 sums_spectra(:,n) = sums_spectra_l 897 #endif 898 ! 899 !-- Normalize spectra by variance 900 sum_spec_dum = SUM( sums_spectra(1:ny/2,n) ) 901 IF ( sum_spec_dum /= 0.0_wp ) THEN 902 sums_spectra(1:ny/2,n) = sums_spectra(1:ny/2,n) * var_d(k) / sum_spec_dum 903 ENDIF 904 n = n + 1 905 906 ENDDO 907 n = n - 1 908 909 910 IF ( myid == 0 ) THEN 911 ! 912 !-- Sum of spectra for later averaging (see routine data_output_spectra) 913 DO j = 1, ny/2 914 DO k = 1, n 915 spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k) 916 ENDDO 852 917 ENDDO 853 n = n - 1 854 855 856 IF ( myid == 0 ) THEN 857 ! 858 !-- Sum of spectra for later averaging (see routine data_output_spectra) 859 DO j = 1, ny/2 860 DO k = 1, n 861 spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k) 862 ENDDO 863 ENDDO 864 865 ENDIF 866 867 ! 868 !-- n_sp_y is needed by data_output_spectra_y 869 n_sp_y = n 870 871 END SUBROUTINE calc_spectra_y 918 919 ENDIF 920 921 ! 922 !-- n_sp_y is needed by data_output_spectra_y 923 n_sp_y = n 924 925 END SUBROUTINE calc_spectra_y 872 926 #endif 873 927
Note: See TracChangeset
for help on using the changeset viewer.
