Home

Context Navigation

spectra_mod.f90 @ 4883

Last change on this file since 4883 was 4843, checked in by raasch, 4 years 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 hellstea $
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.

Context Navigation

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

Download in other formats: