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

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

files re-formatted to follow the PALM coding standard

  • Property svn:keywords set to Id
File size: 32.5 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 4591 2020-07-06 15:56:08Z raasch $
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!
32! 4360 2020-01-07 11:25:50Z suehring
33! Corrected "Former revisions" section
34!
35! 3805 2019-03-20 15:26:35Z raasch
36! Unused variable removed
37!
38! 3655 2019-01-07 16:51:22Z knoop
39! Renamed output variables
40!
41! Revision 1.1  2001/01/05 15:08:07  raasch
42! Initial revision
43!
44!--------------------------------------------------------------------------------------------------!
45! Description:
46! ------------
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!--------------------------------------------------------------------------------------------------!
52 MODULE spectra_mod
53
54    USE kinds
55
56    PRIVATE
57
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            !<
65
66    LOGICAL ::  calculate_spectra   = .FALSE.  !< internal switch that spectra are calculated
67    LOGICAL ::  spectra_initialized = .FALSE.  !< internal switch that spectra related quantities are initialized
68
69    REAL(wp) ::  averaging_interval_sp = 9999999.9_wp  !< averaging interval for spectra output
70    REAL(wp) ::  dt_dosp = 9999999.9_wp                !< time interval for spectra output
71    REAL(wp) ::  skip_time_dosp = 9999999.9_wp         !< no output of spectra data before this interval has passed
72
73    REAL(wp), DIMENSION(:), ALLOCATABLE ::  var_d  !<
74
75    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  spectrum_x  !<
76    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  spectrum_y  !<
77
78    SAVE
79
80    INTERFACE calc_spectra
81       MODULE PROCEDURE calc_spectra
82    END INTERFACE calc_spectra
83
84    INTERFACE preprocess_spectra
85       MODULE PROCEDURE preprocess_spectra
86    END INTERFACE preprocess_spectra
87
88#if defined( __parallel )
89    INTERFACE calc_spectra_x
90       MODULE PROCEDURE calc_spectra_x
91    END INTERFACE calc_spectra_x
92
93    INTERFACE calc_spectra_y
94       MODULE PROCEDURE calc_spectra_y
95    END INTERFACE calc_spectra_y
96#endif
97
98    INTERFACE spectra_check_parameters
99       MODULE PROCEDURE spectra_check_parameters
100    END INTERFACE spectra_check_parameters
101
102    INTERFACE spectra_header
103       MODULE PROCEDURE spectra_header
104    END INTERFACE spectra_header
105
106    INTERFACE spectra_init
107       MODULE PROCEDURE spectra_init
108    END INTERFACE spectra_init
109
110    INTERFACE spectra_parin
111       MODULE PROCEDURE spectra_parin
112    END INTERFACE spectra_parin
113
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
133
134
135 CONTAINS
136
137!--------------------------------------------------------------------------------------------------!
138! Description:
139! ------------
140!> Parin for &spectra_par for calculating spectra
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!--------------------------------------------------------------------------------------------------!
238! Description:
239! ------------
240!> Initialization of spectra related variables
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
259       ENDIF
260
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!--------------------------------------------------------------------------------------------------!
277! Description:
278! ------------
279!> Check spectra related quantities
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!--------------------------------------------------------------------------------------------------!
311! Description:
312! ------------
313!> Header output for spectra
314!>
315!> @todo Output of netcdf data format and compression level
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!--------------------------------------------------------------------------------------------------!
373! Description:
374! ------------
375!> @todo Missing subroutine description.
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 )
469#else
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!--------------------------------------------------------------------------------------------------!
521! Description:
522! ------------
523!> @todo Missing subroutine description.
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 )
548#if !defined( __mpifh )
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 )
565#if defined( __mpifh )
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 )
631#else
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!--------------------------------------------------------------------------------------------------!
640! Description:
641! ------------
642!> @todo Missing subroutine description.
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
658
659#if !defined( __mpifh )
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
676
677#if defined( __mpifh )
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
876             DO  i = nxl_yd, nxr_yd
877                sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
878             ENDDO
879
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 )
895#else
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
917       ENDDO
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
926#endif
927
928 END MODULE spectra_mod
Note: See TracBrowser for help on using the repository browser.