source: palm/trunk/SOURCE/spectra_mod.f90 @ 4725

Last change on this file since 4725 was 4629, checked in by raasch, 4 years ago

support for MPI Fortran77 interface (mpif.h) removed

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