Home

Context Navigation

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

Last change on this file since 2039 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 26.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |