Home

Context Navigation

spectrum.f90 @ 1825

Last change on this file since 1825 was 1818, checked in by maronga, 9 years ago
last commit documented / copyright update
Property svn:keywords set to `Id`
File size: 20.9 KB

Line
1	!> @file spectrum.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: spectrum.f90 1818 2016-04-06 15:53:27Z gronemeier $
26	!
27	! 1815 2016-04-06 13:49:59Z raasch
28	! bugfix: preprocessor directives included for the non-parallel case
29	!
30	! 1808 2016-04-05 19:44:00Z raasch
31	! MPI module used by default on all machines
32	!
33	! 1786 2016-03-08 05:49:27Z raasch
34	! routine is modularized, filename renamed from calc_spectra to spectrum,
35	! privious data module spectrum moved from modules.f90 to here,
36	! cpp-direktives for spectra removed, immediate return if no spectra levels are
37	! given
38	!
39	! 1682 2015-10-07 23:56:08Z knoop
40	! Code annotations made doxygen readable
41	!
42	! 1575 2015-03-27 09:56:27Z raasch
43	! adjustments for psolver-queries
44	!
45	! 1511 2014-12-16 15:54:16Z suehring
46	! Bugfix concerning spectra normalization
47	!
48	! 1431 2014-07-15 14:47:17Z suehring
49	! Wavenumber-integrated spectra coincide with respective variance.
50	!
51	! 1342 2014-03-26 17:04:47Z kanani
52	! REAL constants defined as wp-kinds
53	!
54	! 1324 2014-03-21 09:13:16Z suehring
55	! Bugfix: nzb_x, nzb_yd, nyn_x, nyn_x, nzt_x, nzt_yd belong to transpose_indices
56	!
57	! 1320 2014-03-20 08:40:49Z raasch
58	! ONLY-attribute added to USE-statements,
59	! kind-parameters added to all INTEGER and REAL declaration statements,
60	! kinds are defined in new module kinds,
61	! revision history before 2012 removed,
62	! comment fields (!:) to be used for variable explanations added to
63	! all variable declaration statements
64	!
65	! 1318 2014-03-17 13:35:16Z raasch
66	! module interfaces removed
67	!
68	! 1216 2013-08-26 09:31:42Z raasch
69	! resorting of array moved to separate routine resort_for_zx,
70	! one argument removed from the transpose_..d routines
71	!
72	! 1120 2013-04-05 15:11:35Z raasch
73	! bugfix: calls of fft_x\|y replaced by fft_x\|y_1d
74	!
75	! 1036 2012-10-22 13:43:42Z raasch
76	! code put under GPL (PALM 3.9)
77	!
78	! 1003 2012-09-14 14:35:53Z raasch
79	! adjustment of array tend for cases with unequal subdomain sizes removed
80	!
81	! Revision 1.1 2001/01/05 15:08:07 raasch
82	! Initial revision
83	!
84	!
85	! Description:
86	! ------------
87	!> Calculate horizontal spectra along x and y.
88	!> ATTENTION: 1d-decomposition along y still needs improvement, because in that
89	!> case the gridpoint number along z still depends on the PE number
90	!> because transpose_xz has to be used (and possibly also
91	!> transpose_zyd needs modification).
92	!------------------------------------------------------------------------------!
93	MODULE spectrum
94
95	USE kinds
96
97	PRIVATE
98
99	CHARACTER (LEN=6), DIMENSION(1:5) :: header_char = (/ 'PS(u) ', 'PS(v) ',&
100	'PS(w) ', 'PS(pt)', 'PS(q) ' /)
101	CHARACTER (LEN=2), DIMENSION(10) :: spectra_direction = 'x'
102	CHARACTER (LEN=10), DIMENSION(10) :: data_output_sp = ' '
103	CHARACTER (LEN=25), DIMENSION(1:5) :: utext_char = &
104	(/ '-power spectrum of u ', &
105	'-power spectrum of v ', &
106	'-power spectrum of w ', &
107	'-power spectrum of ^1185 ', &
108	'-power spectrum of q ' /)
109	CHARACTER (LEN=39), DIMENSION(1:5) :: ytext_char = &
110	(/ 'k ^2236 ^2566^2569<u(k) in m>2s>->2 ', &
111	'k ^2236 ^2566^2569<v(k) in m>2s>->2 ', &
112	'k ^2236 ^2566^2569<w(k) in m>2s>->2 ', &
113	'k ^2236 ^2566^2569<^1185(k) in m>2s>->2', &
114	'k ^2236 ^2566^2569<q(k) in m>2s>->2 ' /)
115
116	INTEGER(iwp) :: klist_x = 0, klist_y = 0, n_sp_x = 0, n_sp_y = 0
117
118	INTEGER(iwp) :: comp_spectra_level(100) = 999999, &
119	lstyles(100) = (/ 0, 7, 3, 10, 1, 4, 9, 2, 6, 8, &
120	0, 7, 3, 10, 1, 4, 9, 2, 6, 8, &
121	0, 7, 3, 10, 1, 4, 9, 2, 6, 8, &
122	0, 7, 3, 10, 1, 4, 9, 2, 6, 8, &
123	0, 7, 3, 10, 1, 4, 9, 2, 6, 8, &
124	0, 7, 3, 10, 1, 4, 9, 2, 6, 8, &
125	0, 7, 3, 10, 1, 4, 9, 2, 6, 8, &
126	0, 7, 3, 10, 1, 4, 9, 2, 6, 8, &
127	0, 7, 3, 10, 1, 4, 9, 2, 6, 8, &
128	0, 7, 3, 10, 1, 4, 9, 2, 6, 8 /), &
129	plot_spectra_level(100) = 999999
130
131	REAL(wp) :: time_to_start_sp = 0.0_wp
132
133	PUBLIC comp_spectra_level, data_output_sp, header_char, klist_x, klist_y, &
134	lstyles, n_sp_x, n_sp_y, plot_spectra_level, spectra_direction, &
135	utext_char, ytext_char
136
137	SAVE
138
139	INTERFACE calc_spectra
140	MODULE PROCEDURE calc_spectra
141	END INTERFACE calc_spectra
142
143	INTERFACE preprocess_spectra
144	MODULE PROCEDURE preprocess_spectra
145	END INTERFACE preprocess_spectra
146
147	INTERFACE calc_spectra_x
148	MODULE PROCEDURE calc_spectra_x
149	END INTERFACE calc_spectra_x
150
151	INTERFACE calc_spectra_y
152	MODULE PROCEDURE calc_spectra_y
153	END INTERFACE calc_spectra_y
154
155	PUBLIC calc_spectra
156
157
158	CONTAINS
159
160	SUBROUTINE calc_spectra
161
162	USE arrays_3d, &
163	ONLY: d, tend
164
165	USE control_parameters, &
166	ONLY: average_count_sp, bc_lr_cyc, bc_ns_cyc, message_string, psolver
167
168	USE cpulog, &
169	ONLY: cpu_log, log_point
170
171	USE fft_xy, &
172	ONLY: fft_init
173
174	USE indices, &
175	ONLY: nxl, nxr, nyn, nys, nzb, nzt
176
177	USE kinds
178
179	USE pegrid, &
180	ONLY: myid, pdims
181
182	IMPLICIT NONE
183
184	INTEGER(iwp) :: m !<
185	INTEGER(iwp) :: pr !<
186
187
188	!
189	!-- Check if user gave any levels for spectra to be calculated
190	IF ( comp_spectra_level(1) == 999999 ) RETURN
191
192	CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
193
194	!
195	!-- Initialize ffts
196	CALL fft_init
197
198	!
199	!-- Reallocate array d in required size
200	IF ( psolver(1:9) == 'multigrid' ) THEN
201	DEALLOCATE( d )
202	ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
203	ENDIF
204
205	m = 1
206	DO WHILE ( data_output_sp(m) /= ' ' .AND. m <= 10 )
207	!
208	!-- Transposition from z --> x ( y --> x in case of a 1d-decomposition
209	!-- along x)
210	IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN
211
212	!
213	!-- Calculation of spectra works for cyclic boundary conditions only
214	IF ( .NOT. bc_lr_cyc ) THEN
215
216	message_string = 'non-cyclic lateral boundaries along x do'// &
217	' not & allow calculation of spectra along x'
218	CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
219	ENDIF
220
221	CALL preprocess_spectra( m, pr )
222
223	#if defined( __parallel )
224	IF ( pdims(2) /= 1 ) THEN
225	CALL resort_for_zx( d, tend )
226	CALL transpose_zx( tend, d )
227	ELSE
228	CALL transpose_yxd( d, d )
229	ENDIF
230	CALL calc_spectra_x( d, pr, m )
231	#else
232	message_string = 'sorry, calculation of spectra in non paral' // &
233	'lel mode& is still not realized'
234	CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
235	#endif
236
237	ENDIF
238
239	!
240	!-- Transposition from z --> y (d is rearranged only in case of a
241	!-- 1d-decomposition along x)
242	IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN
243
244	!
245	!-- Calculation of spectra works for cyclic boundary conditions only
246	IF ( .NOT. bc_ns_cyc ) THEN
247	IF ( myid == 0 ) THEN
248	message_string = 'non-cyclic lateral boundaries along y' // &
249	' do not & allow calculation of spectr' // &
250	'a along y'
251	CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
252	ENDIF
253	CALL local_stop
254	ENDIF
255
256	CALL preprocess_spectra( m, pr )
257
258	#if defined( __parallel )
259	CALL transpose_zyd( d, d )
260	CALL calc_spectra_y( d, pr, m )
261	#else
262	message_string = 'sorry, calculation of spectra in non paral' // &
263	'lel mode& is still not realized'
264	CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
265	#endif
266
267	ENDIF
268
269	!
270	!-- Increase counter for next spectrum
271	m = m + 1
272
273	ENDDO
274
275	!
276	!-- Increase counter for averaging process in routine plot_spectra
277	average_count_sp = average_count_sp + 1
278
279	CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
280
281	END SUBROUTINE calc_spectra
282
283
284	!------------------------------------------------------------------------------!
285	! Description:
286	! ------------
287	!> @todo Missing subroutine description.
288	!------------------------------------------------------------------------------!
289	SUBROUTINE preprocess_spectra( m, pr )
290
291	USE arrays_3d, &
292	ONLY: d, pt, q, u, v, w
293
294	USE indices, &
295	ONLY: ngp_2dh, nxl, nxr, nyn, nys, nzb, nzt
296
297	USE kinds
298
299	#if defined( __parallel )
300	#if defined( __mpifh )
301	INCLUDE "mpif.h"
302	#else
303	USE MPI
304	#endif
305	#endif
306	USE pegrid, &
307	ONLY: collective_wait, comm2d, ierr
308
309	USE statistics, &
310	ONLY: hom, var_d
311
312
313	IMPLICIT NONE
314
315	INTEGER(iwp) :: i !<
316	INTEGER(iwp) :: j !<
317	INTEGER(iwp) :: k !<
318	INTEGER(iwp) :: m !<
319	INTEGER(iwp) :: pr !<
320
321	REAL(wp), DIMENSION(nzb:nzt+1) :: var_d_l
322
323	SELECT CASE ( TRIM( data_output_sp(m) ) )
324
325	CASE ( 'u' )
326	pr = 1
327	d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
328
329	CASE ( 'v' )
330	pr = 2
331	d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
332
333	CASE ( 'w' )
334	pr = 3
335	d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
336
337	CASE ( 'pt' )
338	pr = 4
339	d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
340
341	CASE ( 'q' )
342	pr = 41
343	d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
344
345	CASE DEFAULT
346	!
347	!-- The DEFAULT case is reached either if the parameter data_output_sp(m)
348	!-- contains a wrong character string or if the user has coded a special
349	!-- case in the user interface. There, the subroutine user_spectra
350	!-- checks which of these two conditions applies.
351	CALL user_spectra( 'preprocess', m, pr )
352
353	END SELECT
354
355	!
356	!-- Subtract horizontal mean from the array, for which spectra have to be
357	!-- calculated
358	var_d_l(:) = 0.0_wp
359	DO i = nxl, nxr
360	DO j = nys, nyn
361	DO k = nzb+1, nzt
362	d(k,j,i) = d(k,j,i) - hom(k,1,pr,0)
363	var_d_l(k) = var_d_l(k) + d(k,j,i) * d(k,j,i)
364	ENDDO
365	ENDDO
366	ENDDO
367	!
368	!-- Compute total variance from local variances
369	var_d(:) = 0.0_wp
370	#if defined( __parallel )
371	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
372	CALL MPI_ALLREDUCE( var_d_l(0), var_d(0), nzt+1-nzb, MPI_REAL, MPI_SUM, &
373	comm2d, ierr )
374	#else
375	var_d(:) = var_d_l(:)
376	#endif
377	var_d(:) = var_d(:) / ngp_2dh(0)
378
379	END SUBROUTINE preprocess_spectra
380
381
382	!------------------------------------------------------------------------------!
383	! Description:
384	! ------------
385	!> @todo Missing subroutine description.
386	!------------------------------------------------------------------------------!
387	SUBROUTINE calc_spectra_x( ddd, pr, m )
388
389	USE control_parameters, &
390	ONLY: fft_method
391
392	USE fft_xy, &
393	ONLY: fft_x_1d
394
395	USE grid_variables, &
396	ONLY: dx
397
398	USE indices, &
399	ONLY: nx, ny
400
401	USE kinds
402
403	#if defined( __parallel )
404	#if defined( __mpifh )
405	INCLUDE "mpif.h"
406	#else
407	USE MPI
408	#endif
409	#endif
410	USE pegrid, &
411	ONLY: comm2d, ierr, myid
412
413	USE statistics, &
414	ONLY: spectrum_x, var_d
415
416	USE transpose_indices, &
417	ONLY: nyn_x, nys_x, nzb_x, nzt_x
418
419
420	IMPLICIT NONE
421
422	INTEGER(iwp) :: i !<
423	INTEGER(iwp) :: ishape(1) !<
424	INTEGER(iwp) :: j !<
425	INTEGER(iwp) :: k !<
426	INTEGER(iwp) :: m !<
427	INTEGER(iwp) :: n !<
428	INTEGER(iwp) :: pr !<
429
430	REAL(wp) :: exponent !<
431	REAL(wp) :: sum_spec_dum !< wavenumber-integrated spectrum
432
433	REAL(wp), DIMENSION(0:nx) :: work !<
434
435	REAL(wp), DIMENSION(0:nx/2) :: sums_spectra_l !<
436
437	REAL(wp), DIMENSION(0:nx/2,100) :: sums_spectra !<
438
439	REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ddd !<
440
441	!
442	!-- Exponent for geometric average
443	exponent = 1.0_wp / ( ny + 1.0_wp )
444
445	!
446	!-- Loop over all levels defined by the user
447	n = 1
448	DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 )
449
450	k = comp_spectra_level(n)
451
452	!
453	!-- Calculate FFT only if the corresponding level is situated on this PE
454	IF ( k >= nzb_x .AND. k <= nzt_x ) THEN
455
456	DO j = nys_x, nyn_x
457
458	work = ddd(0:nx,j,k)
459	CALL fft_x_1d( work, 'forward' )
460
461	ddd(0,j,k) = dx * work(0)**2
462	DO i = 1, nx/2
463	ddd(i,j,k) = dx * ( work(i)2 + work(nx+1-i)2 )
464	ENDDO
465
466	ENDDO
467
468	!
469	!-- Local sum and geometric average of these spectra
470	!-- (WARNING: no global sum should be performed, because floating
471	!-- point overflow may occur)
472	DO i = 0, nx/2
473
474	sums_spectra_l(i) = 1.0_wp
475	DO j = nys_x, nyn_x
476	sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
477	ENDDO
478
479	ENDDO
480
481	ELSE
482
483	sums_spectra_l = 1.0_wp
484
485	ENDIF
486
487	!
488	!-- Global sum of spectra on PE0 (from where they are written on file)
489	sums_spectra(:,n) = 0.0_wp
490	#if defined( __parallel )
491	CALL MPI_BARRIER( comm2d, ierr ) ! Necessary?
492	CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1, &
493	MPI_REAL, MPI_PROD, 0, comm2d, ierr )
494	#else
495	sums_spectra(:,n) = sums_spectra_l
496	#endif
497	!
498	!-- Normalize spectra by variance
499	sum_spec_dum = SUM( sums_spectra(:,n) )
500	IF ( sum_spec_dum /= 0.0_wp ) THEN
501	sums_spectra(:,n) = sums_spectra(:,n) * var_d(k) / sum_spec_dum
502	ENDIF
503	n = n + 1
504
505	ENDDO
506	n = n - 1
507
508	IF ( myid == 0 ) THEN
509	!
510	!-- Sum of spectra for later averaging (see routine data_output_spectra)
511	DO i = 1, nx/2
512	DO k = 1, n
513	spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k)
514	ENDDO
515	ENDDO
516
517	ENDIF
518	!
519	!-- n_sp_x is needed by data_output_spectra_x
520	n_sp_x = n
521
522	END SUBROUTINE calc_spectra_x
523
524
525	!------------------------------------------------------------------------------!
526	! Description:
527	! ------------
528	!> @todo Missing subroutine description.
529	!------------------------------------------------------------------------------!
530	SUBROUTINE calc_spectra_y( ddd, pr, m )
531
532	USE control_parameters, &
533	ONLY: fft_method
534
535	USE fft_xy, &
536	ONLY: fft_y_1d
537
538	USE grid_variables, &
539	ONLY: dy
540
541	USE indices, &
542	ONLY: nx, ny
543
544	USE kinds
545
546	#if defined( __parallel )
547	#if defined( __mpifh )
548	INCLUDE "mpif.h"
549	#else
550	USE MPI
551	#endif
552	#endif
553	USE pegrid, &
554	ONLY: comm2d, ierr, myid
555
556	USE statistics, &
557	ONLY: spectrum_y, var_d
558
559	USE transpose_indices, &
560	ONLY: nxl_yd, nxr_yd, nzb_yd, nzt_yd
561
562
563	IMPLICIT NONE
564
565	INTEGER(iwp) :: i !<
566	INTEGER(iwp) :: j !<
567	INTEGER(iwp) :: jshape(1) !<
568	INTEGER(iwp) :: k !<
569	INTEGER(iwp) :: m !<
570	INTEGER(iwp) :: n !<
571	INTEGER(iwp) :: pr !<
572
573	REAL(wp) :: exponent !<
574	REAL(wp) :: sum_spec_dum !< wavenumber-integrated spectrum
575
576	REAL(wp), DIMENSION(0:ny) :: work !<
577
578	REAL(wp), DIMENSION(0:ny/2) :: sums_spectra_l !<
579
580	REAL(wp), DIMENSION(0:ny/2,100) :: sums_spectra !<
581
582	REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !<
583
584
585	!
586	!-- Exponent for geometric average
587	exponent = 1.0_wp / ( nx + 1.0_wp )
588
589	!
590	!-- Loop over all levels defined by the user
591	n = 1
592	DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 )
593
594	k = comp_spectra_level(n)
595
596	!
597	!-- Calculate FFT only if the corresponding level is situated on this PE
598	IF ( k >= nzb_yd .AND. k <= nzt_yd ) THEN
599
600	DO i = nxl_yd, nxr_yd
601
602	work = ddd(0:ny,i,k)
603	CALL fft_y_1d( work, 'forward' )
604
605	ddd(0,i,k) = dy * work(0)**2
606	DO j = 1, ny/2
607	ddd(j,i,k) = dy * ( work(j)2 + work(ny+1-j)2 )
608	ENDDO
609
610	ENDDO
611
612	!
613	!-- Local sum and geometric average of these spectra
614	!-- (WARNING: no global sum should be performed, because floating
615	!-- point overflow may occur)
616	DO j = 0, ny/2
617
618	sums_spectra_l(j) = 1.0_wp
619	DO i = nxl_yd, nxr_yd
620	sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
621	ENDDO
622
623	ENDDO
624
625	ELSE
626
627	sums_spectra_l = 1.0_wp
628
629	ENDIF
630
631	!
632	!-- Global sum of spectra on PE0 (from where they are written on file)
633	sums_spectra(:,n) = 0.0_wp
634	#if defined( __parallel )
635	CALL MPI_BARRIER( comm2d, ierr ) ! Necessary?
636	CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1, &
637	MPI_REAL, MPI_PROD, 0, comm2d, ierr )
638	#else
639	sums_spectra(:,n) = sums_spectra_l
640	#endif
641	!
642	!-- Normalize spectra by variance
643	sum_spec_dum = SUM( sums_spectra(:,n) )
644	IF ( SUM(sums_spectra(:,n)) /= 0.0_wp ) THEN
645	sums_spectra(:,n) = sums_spectra(:,n) * &
646	var_d(k) / SUM(sums_spectra(:,n))
647	ENDIF
648	n = n + 1
649
650	ENDDO
651	n = n - 1
652
653
654	IF ( myid == 0 ) THEN
655	!
656	!-- Sum of spectra for later averaging (see routine data_output_spectra)
657	DO j = 1, ny/2
658	DO k = 1, n
659	spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k)
660	ENDDO
661	ENDDO
662
663	ENDIF
664
665	!
666	!-- n_sp_y is needed by data_output_spectra_y
667	n_sp_y = n
668
669	END SUBROUTINE calc_spectra_y
670
671	END MODULE spectrum

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

Context Navigation

source: palm/trunk/SOURCE/spectrum.f90 @ 1825

Download in other formats: