source: palm/trunk/SOURCE/spectra_mod.f90

Last change on this file was 4843, checked in by raasch, 8 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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