Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |