Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |