Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |