Home

Context Navigation

source: palm/trunk/SOURCE/calc_spectra.f90 @ 1320

Last change on this file since 1320 was 1320, checked in by raasch, 11 years ago
ONLY-attribute added to USE-statements, kind-parameters added to all INTEGER and REAL declaration statements, kinds are defined in new module kinds, old module precision_kind is removed, revision history before 2012 removed, comment fields (!:) to be used for variable explanations added to all variable declaration statements
Property svn:keywords set to `Id`
File size: 14.5 KB

Line
1	SUBROUTINE calc_spectra
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	! ONLY-attribute added to USE-statements,
23	! kind-parameters added to all INTEGER and REAL declaration statements,
24	! kinds are defined in new module kinds,
25	! revision history before 2012 removed,
26	! comment fields (!:) to be used for variable explanations added to
27	! all variable declaration statements
28	!
29	! Former revisions:
30	! -----------------
31	! $Id: calc_spectra.f90 1320 2014-03-20 08:40:49Z raasch $
32	!
33	! 1318 2014-03-17 13:35:16Z raasch
34	! module interfaces removed
35	!
36	! 1216 2013-08-26 09:31:42Z raasch
37	! resorting of array moved to separate routine resort_for_zx,
38	! one argument removed from the transpose_..d routines
39	!
40	! 1120 2013-04-05 15:11:35Z raasch
41	! bugfix: calls of fft_x\|y replaced by fft_x\|y_1d
42	!
43	! 1036 2012-10-22 13:43:42Z raasch
44	! code put under GPL (PALM 3.9)
45	!
46	! 1003 2012-09-14 14:35:53Z raasch
47	! adjustment of array tend for cases with unequal subdomain sizes removed
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
57	! case the gridpoint number along z still depends on the PE number
58	! because transpose_xz has to be used (and possibly also
59	! transpose_zyd needs modification).
60	!------------------------------------------------------------------------------!
61
62	#if defined( __spectra )
63	USE arrays_3d, &
64	ONLY: d, tend
65
66	USE control_parameters, &
67	ONLY: average_count_sp, bc_lr_cyc, bc_ns_cyc, message_string, psolver
68
69	USE cpulog, &
70	ONLY: cpu_log, log_point
71
72	USE fft_xy, &
73	ONLY: fft_init
74
75	USE indices, &
76	ONLY: nxl, nxr, nyn, nys, nzb, nzt, nzt_x, nzt_yd
77
78	USE kinds
79
80	USE pegrid
81
82	USE spectrum, &
83	ONLY: data_output_sp, spectra_direction
84
85
86	IMPLICIT NONE
87
88	INTEGER(iwp) :: m !:
89	INTEGER(iwp) :: pr !:
90
91
92	CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
93
94	!
95	!-- Initialize ffts
96	CALL fft_init
97
98	!
99	!-- Reallocate array d in required size
100	IF ( psolver == 'multigrid' ) THEN
101	DEALLOCATE( d )
102	ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
103	ENDIF
104
105	m = 1
106	DO WHILE ( data_output_sp(m) /= ' ' .AND. m <= 10 )
107	!
108	!-- Transposition from z --> x ( y --> x in case of a 1d-decomposition
109	!-- along x)
110	IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN
111
112	!
113	!-- Calculation of spectra works for cyclic boundary conditions only
114	IF ( .NOT. bc_lr_cyc ) THEN
115
116	message_string = 'non-cyclic lateral boundaries along x do not'// &
117	'& allow calculation of spectra along x'
118	CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
119	ENDIF
120
121	CALL preprocess_spectra( m, pr )
122
123	#if defined( __parallel )
124	IF ( pdims(2) /= 1 ) THEN
125	CALL resort_for_zx( d, tend )
126	CALL transpose_zx( tend, d )
127	ELSE
128	CALL transpose_yxd( d, d )
129	ENDIF
130	CALL calc_spectra_x( d, pr, m )
131	#else
132	message_string = 'sorry, calculation of spectra in non parallel ' // &
133	'mode& is still not realized'
134	CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
135	#endif
136
137	ENDIF
138
139	!
140	!-- Transposition from z --> y (d is rearranged only in case of a
141	!-- 1d-decomposition along x)
142	IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN
143
144	!
145	!-- Calculation of spectra works for cyclic boundary conditions only
146	IF ( .NOT. bc_ns_cyc ) THEN
147	IF ( myid == 0 ) THEN
148	message_string = 'non-cyclic lateral boundaries along y do' // &
149	' not & allow calculation of spectra along y'
150	CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
151	ENDIF
152	CALL local_stop
153	ENDIF
154
155	CALL preprocess_spectra( m, pr )
156
157	#if defined( __parallel )
158	CALL transpose_zyd( d, d )
159	CALL calc_spectra_y( d, pr, m )
160	#else
161	message_string = 'sorry, calculation of spectra in non parallel' // &
162	'mode& is still not realized'
163	CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
164	#endif
165
166	ENDIF
167
168	!
169	!-- Increase counter for next spectrum
170	m = m + 1
171
172	ENDDO
173
174	!
175	!-- Increase counter for averaging process in routine plot_spectra
176	average_count_sp = average_count_sp + 1
177
178	CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
179
180	#endif
181	END SUBROUTINE calc_spectra
182
183
184	#if defined( __spectra )
185	SUBROUTINE preprocess_spectra( m, pr )
186
187	USE arrays_3d, &
188	ONLY: d, pt, q, u, v, w
189
190	USE indices, &
191	ONLY: nxl, nxr, nyn, nys, nzb, nzt
192
193	USE kinds
194
195	USE pegrid
196
197	USE spectrum, &
198	ONLY: data_output_sp
199
200	USE statistics, &
201	ONLY: sums
202
203
204	IMPLICIT NONE
205
206	INTEGER(iwp) :: i !:
207	INTEGER(iwp) :: j !:
208	INTEGER(iwp) :: k !:
209	INTEGER(iwp) :: m !:
210	INTEGER(iwp) :: pr !:
211
212	SELECT CASE ( TRIM( data_output_sp(m) ) )
213
214	CASE ( 'u' )
215	pr = 1
216	d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
217
218	CASE ( 'v' )
219	pr = 2
220	d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
221
222	CASE ( 'w' )
223	pr = 3
224	d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
225
226	CASE ( 'pt' )
227	pr = 4
228	d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
229
230	CASE ( 'q' )
231	pr = 41
232	d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
233
234	CASE DEFAULT
235	!
236	!-- The DEFAULT case is reached either if the parameter data_output_sp(m)
237	!-- contains a wrong character string or if the user has coded a special
238	!-- case in the user interface. There, the subroutine user_spectra
239	!-- checks which of these two conditions applies.
240	CALL user_spectra( 'preprocess', m, pr )
241
242	END SELECT
243
244	!
245	!-- Subtract horizontal mean from the array, for which spectra have to be
246	!-- calculated
247	DO i = nxl, nxr
248	DO j = nys, nyn
249	DO k = nzb+1, nzt
250	d(k,j,i) = d(k,j,i) - sums(k,pr)
251	ENDDO
252	ENDDO
253	ENDDO
254
255	END SUBROUTINE preprocess_spectra
256
257
258	SUBROUTINE calc_spectra_x( ddd, pr, m )
259
260	USE arrays_3d, &
261	ONLY:
262
263	USE control_parameters, &
264	ONLY: fft_method
265
266	USE fft_xy, &
267	ONLY: fft_x_1d
268
269	USE grid_variables, &
270	ONLY: dx
271
272	USE indices, &
273	ONLY: nx, ny, nyn_x, nys_x, nzb_x, nzt_x
274
275	USE kinds
276
277	USE pegrid
278
279	USE spectrum, &
280	ONLY: comp_spectra_level, n_sp_x
281
282	USE statistics, &
283	ONLY: spectrum_x
284
285	USE transpose_indices
286
287
288	IMPLICIT NONE
289
290	INTEGER(iwp) :: i !:
291	INTEGER(iwp) :: ishape(1) !:
292	INTEGER(iwp) :: j !:
293	INTEGER(iwp) :: k !:
294	INTEGER(iwp) :: m !:
295	INTEGER(iwp) :: n !:
296	INTEGER(iwp) :: pr !:
297
298	REAL(wp) :: fac !:
299	REAL(wp) :: exponent !:
300
301	REAL(wp), DIMENSION(0:nx) :: work !:
302
303	REAL(wp), DIMENSION(0:nx/2) :: sums_spectra_l !:
304
305	REAL(wp), DIMENSION(0:nx/2,100) :: sums_spectra !:
306
307	REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ddd !:
308
309	!
310	!-- Exponent for geometric average
311	exponent = 1.0 / ( ny + 1.0 )
312
313	!
314	!-- Loop over all levels defined by the user
315	n = 1
316	DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 )
317
318	k = comp_spectra_level(n)
319
320	!
321	!-- Calculate FFT only if the corresponding level is situated on this PE
322	IF ( k >= nzb_x .AND. k <= nzt_x ) THEN
323
324	DO j = nys_x, nyn_x
325
326	work = ddd(0:nx,j,k)
327	CALL fft_x_1d( work, 'forward' )
328
329	ddd(0,j,k) = dx * work(0)**2
330	DO i = 1, nx/2
331	ddd(i,j,k) = dx * ( work(i)2 + work(nx+1-i)2 )
332	ENDDO
333
334	ENDDO
335
336	!
337	!-- Local sum and geometric average of these spectra
338	!-- (WARNING: no global sum should be performed, because floating
339	!-- point overflow may occur)
340	DO i = 0, nx/2
341
342	sums_spectra_l(i) = 1.0
343	DO j = nys_x, nyn_x
344	sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
345	ENDDO
346
347	ENDDO
348
349	ELSE
350
351	sums_spectra_l = 1.0
352
353	ENDIF
354
355	!
356	!-- Global sum of spectra on PE0 (from where they are written on file)
357	sums_spectra(:,n) = 0.0
358	#if defined( __parallel )
359	CALL MPI_BARRIER( comm2d, ierr ) ! Necessary?
360	CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1, &
361	MPI_REAL, MPI_PROD, 0, comm2d, ierr )
362	#else
363	sums_spectra(:,n) = sums_spectra_l
364	#endif
365
366	n = n + 1
367
368	ENDDO
369	n = n - 1
370
371	IF ( myid == 0 ) THEN
372	!
373	!-- Sum of spectra for later averaging (see routine data_output_spectra)
374	!-- Temperton fft results need to be normalized
375	IF ( fft_method == 'temperton-algorithm' ) THEN
376	fac = nx + 1.0
377	ELSE
378	fac = 1.0
379	ENDIF
380	DO i = 1, nx/2
381	DO k = 1, n
382	spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k) * fac
383	ENDDO
384	ENDDO
385
386	ENDIF
387
388	!
389	!-- n_sp_x is needed by data_output_spectra_x
390	n_sp_x = n
391
392	END SUBROUTINE calc_spectra_x
393
394
395	SUBROUTINE calc_spectra_y( ddd, pr, m )
396
397	USE arrays_3d, &
398	ONLY:
399
400	USE control_parameters, &
401	ONLY: fft_method
402
403	USE fft_xy, &
404	ONLY: fft_y_1d
405
406	USE grid_variables, &
407	ONLY: dy
408
409	USE indices, &
410	ONLY: nx, ny, nxl_yd, nxr_yd, nzb_yd, nzt_yd
411
412	USE kinds
413
414	USE pegrid
415
416	USE spectrum, &
417	ONLY: comp_spectra_level, n_sp_y
418
419	USE statistics, &
420	ONLY: spectrum_y
421
422	USE transpose_indices
423
424
425	IMPLICIT NONE
426
427	INTEGER(iwp) :: i !:
428	INTEGER(iwp) :: j !:
429	INTEGER(iwp) :: jshape(1) !:
430	INTEGER(iwp) :: k !:
431	INTEGER(iwp) :: m !:
432	INTEGER(iwp) :: n !:
433	INTEGER(iwp) :: pr !:
434
435	REAL(wp) :: fac !:
436	REAL(wp) :: exponent !:
437
438	REAL(wp), DIMENSION(0:ny) :: work !:
439
440	REAL(wp), DIMENSION(0:ny/2) :: sums_spectra_l !:
441
442	REAL(wp), DIMENSION(0:ny/2,100) :: sums_spectra !:
443
444	REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !:
445
446
447	!
448	!-- Exponent for geometric average
449	exponent = 1.0 / ( nx + 1.0 )
450
451	!
452	!-- Loop over all levels defined by the user
453	n = 1
454	DO WHILE ( comp_spectra_level(n) /= 999999 .AND. n <= 100 )
455
456	k = comp_spectra_level(n)
457
458	!
459	!-- Calculate FFT only if the corresponding level is situated on this PE
460	IF ( k >= nzb_yd .AND. k <= nzt_yd ) THEN
461
462	DO i = nxl_yd, nxr_yd
463
464	work = ddd(0:ny,i,k)
465	CALL fft_y_1d( work, 'forward' )
466
467	ddd(0,i,k) = dy * work(0)**2
468	DO j = 1, ny/2
469	ddd(j,i,k) = dy * ( work(j)2 + work(ny+1-j)2 )
470	ENDDO
471
472	ENDDO
473
474	!
475	!-- Local sum and geometric average of these spectra
476	!-- (WARNING: no global sum should be performed, because floating
477	!-- point overflow may occur)
478	DO j = 0, ny/2
479
480	sums_spectra_l(j) = 1.0
481	DO i = nxl_yd, nxr_yd
482	sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
483	ENDDO
484
485	ENDDO
486
487	ELSE
488
489	sums_spectra_l = 1.0
490
491	ENDIF
492
493	!
494	!-- Global sum of spectra on PE0 (from where they are written on file)
495	sums_spectra(:,n) = 0.0
496	#if defined( __parallel )
497	CALL MPI_BARRIER( comm2d, ierr ) ! Necessary?
498	CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1, &
499	MPI_REAL, MPI_PROD, 0, comm2d, ierr )
500	#else
501	sums_spectra(:,n) = sums_spectra_l
502	#endif
503
504	n = n + 1
505
506	ENDDO
507	n = n - 1
508
509
510	IF ( myid == 0 ) THEN
511	!
512	!-- Sum of spectra for later averaging (see routine data_output_spectra)
513	!-- Temperton fft results need to be normalized
514	IF ( fft_method == 'temperton-algorithm' ) THEN
515	fac = ny + 1.0
516	ELSE
517	fac = 1.0
518	ENDIF
519	DO j = 1, ny/2
520	DO k = 1, n
521	spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k) * fac
522	ENDDO
523	ENDDO
524
525	ENDIF
526
527	!
528	!-- n_sp_y is needed by data_output_spectra_y
529	n_sp_y = n
530
531	END SUBROUTINE calc_spectra_y
532	#endif

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |