Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |