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

Last change on this file since 1078 was 1037, checked in by raasch, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 11.8 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: calc_spectra.f90 1037 2012-10-22 14:10:22Z maronga $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 1003 2012-09-14 14:35:53Z raasch
32! adjustment of array tend for cases with unequal subdomain sizes removed
33!
34! 707 2011-03-29 11:39:40Z raasch
35! bc_lr/ns replaced by bc_lr/ns_cyc
36!
37! 667 2010-12-23 12:06:00Z suehring/gryschka
38! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation
39! of tend
40!
41! 274 2009-03-26 15:11:21Z heinze
42! Output of messages replaced by message handling routine
43!
44! 225 2009-01-26 14:44:20Z raasch
45! Bugfix: array d is reallocated in case that multigrid is used
46!
47! 192 2008-08-27 16:51:49Z letzel
48! bugfix in calc_spectra_x: exponent = 1.0 / ( ny + 1.0 )
49! allow 100 spectra levels instead of 10 for consistency with
50! define_netcdf_header
51! user-defined spectra, arguments removed from transpose routines
52!
53! February 2007
54! RCS Log replace by Id keyword, revision history cleaned up
55!
56! Revision 1.9  2006/04/11 14:56:00  raasch
57! pl_spectra renamed data_output_sp
58!
59! Revision 1.1  2001/01/05 15:08:07  raasch
60! Initial revision
61!
62!
63! Description:
64! ------------
65! Calculate horizontal spectra along x and y.
66! ATTENTION: 1d-decomposition along y still needs improvement, because in that
67!            case the gridpoint number along z still depends on the PE number
68!            because transpose_xz has to be used (and possibly also
69!            transpose_zyd needs modification).
70!------------------------------------------------------------------------------!
71
72#if defined( __spectra )
73    USE arrays_3d
74    USE control_parameters
75    USE cpulog
76    USE fft_xy
77    USE indices
78    USE interfaces
79    USE pegrid
80    USE spectrum
81
82    IMPLICIT NONE
83
84    INTEGER ::  m, pr
85
86
87    CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
88
89!
90!-- Initialize ffts
91    CALL fft_init
92
93!
94!-- Reallocate array d in required size
95    IF ( psolver == 'multigrid' )  THEN
96       DEALLOCATE( d )
97       ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
98    ENDIF
99
100    m = 1
101    DO WHILE ( data_output_sp(m) /= ' '  .AND.  m <= 10 )
102!
103!--    Transposition from z --> x  ( y --> x in case of a 1d-decomposition
104!--    along x)
105       IF ( INDEX( spectra_direction(m), 'x' ) /= 0 )  THEN
106
107!
108!--       Calculation of spectra works for cyclic boundary conditions only
109          IF ( .NOT. bc_lr_cyc )  THEN
110
111             message_string = 'non-cyclic lateral boundaries along x do not'// &
112                              '& allow calculation of spectra along x'
113             CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
114          ENDIF
115
116          CALL preprocess_spectra( m, pr )
117
118#if defined( __parallel )
119          IF ( pdims(2) /= 1 )  THEN
120             CALL transpose_zx( d, tend, d )
121          ELSE
122             CALL transpose_yxd( d, tend, d )
123          ENDIF
124          CALL calc_spectra_x( d, pr, m )
125#else
126          message_string = 'sorry, calculation of spectra in non parallel ' // &
127                           'mode& is still not realized'
128          CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )     
129#endif
130
131       ENDIF
132
133!
134!--    Transposition from z --> y (d is rearranged only in case of a
135!--    1d-decomposition along x)
136       IF ( INDEX( spectra_direction(m), 'y' ) /= 0 )  THEN
137
138!
139!--       Calculation of spectra works for cyclic boundary conditions only
140          IF ( .NOT. bc_ns_cyc )  THEN
141             IF ( myid == 0 )  THEN
142                message_string = 'non-cyclic lateral boundaries along y do' // &
143                                 ' not & allow calculation of spectra along y' 
144                CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
145             ENDIF
146             CALL local_stop
147          ENDIF
148
149          CALL preprocess_spectra( m, pr )
150
151#if defined( __parallel )
152          CALL transpose_zyd( d, tend, d )
153          CALL calc_spectra_y( d, pr, m )
154#else
155          message_string = 'sorry, calculation of spectra in non parallel' // &
156                           'mode& is still not realized'
157          CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
158#endif
159
160       ENDIF
161
162!
163!--    Increase counter for next spectrum
164       m = m + 1
165         
166    ENDDO
167
168!
169!-- Increase counter for averaging process in routine plot_spectra
170    average_count_sp = average_count_sp + 1
171
172    CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
173
174#endif
175 END SUBROUTINE calc_spectra
176
177
178#if defined( __spectra )
179 SUBROUTINE preprocess_spectra( m, pr )
180
181    USE arrays_3d
182    USE indices
183    USE pegrid
184    USE spectrum
185    USE statistics
186
187    IMPLICIT NONE
188
189    INTEGER :: i, j, k, m, pr
190
191    SELECT CASE ( TRIM( data_output_sp(m) ) )
192         
193    CASE ( 'u' )
194       pr = 1
195       d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
196       
197    CASE ( 'v' )
198       pr = 2
199       d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
200       
201    CASE ( 'w' )
202       pr = 3
203       d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
204       
205    CASE ( 'pt' )
206       pr = 4
207       d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
208       
209    CASE ( 'q' )
210       pr = 41
211       d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
212       
213    CASE DEFAULT
214!
215!--    The DEFAULT case is reached either if the parameter data_output_sp(m)
216!--    contains a wrong character string or if the user has coded a special
217!--    case in the user interface. There, the subroutine user_spectra
218!--    checks which of these two conditions applies.
219       CALL user_spectra( 'preprocess', m, pr )
220         
221    END SELECT
222
223!
224!-- Subtract horizontal mean from the array, for which spectra have to be
225!-- calculated
226    DO  i = nxl, nxr
227       DO  j = nys, nyn
228          DO  k = nzb+1, nzt
229             d(k,j,i) = d(k,j,i) - sums(k,pr)
230          ENDDO
231       ENDDO
232    ENDDO
233
234 END SUBROUTINE preprocess_spectra
235
236
237 SUBROUTINE calc_spectra_x( ddd, pr, m )
238
239    USE arrays_3d
240    USE constants
241    USE control_parameters
242    USE fft_xy
243    USE grid_variables
244    USE indices
245    USE pegrid
246    USE spectrum
247    USE statistics
248    USE transpose_indices
249
250    IMPLICIT NONE
251
252    INTEGER                    ::  i, ishape(1), j, k, m, n, pr
253
254    REAL                       ::  fac, exponent
255    REAL, DIMENSION(0:nx)      ::  work
256    REAL, DIMENSION(0:nx/2)    ::  sums_spectra_l
257    REAL, DIMENSION(0:nx/2,100)::  sums_spectra
258
259    REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ddd
260
261!
262!-- Exponent for geometric average
263    exponent = 1.0 / ( ny + 1.0 )
264
265!
266!-- Loop over all levels defined by the user
267    n = 1
268    DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
269
270       k = comp_spectra_level(n)
271
272!
273!--    Calculate FFT only if the corresponding level is situated on this PE
274       IF ( k >= nzb_x  .AND.  k <= nzt_x )  THEN
275         
276          DO  j = nys_x, nyn_x
277
278             work = ddd(0:nx,j,k)
279             CALL fft_x( work, 'forward' )
280
281             ddd(0,j,k) = dx * work(0)**2
282             DO  i = 1, nx/2
283                ddd(i,j,k) = dx * ( work(i)**2 + work(nx+1-i)**2 )
284             ENDDO
285
286          ENDDO
287
288!
289!--       Local sum and geometric average of these spectra
290!--       (WARNING: no global sum should be performed, because floating
291!--       point overflow may occur)
292          DO  i = 0, nx/2
293
294             sums_spectra_l(i) = 1.0
295             DO  j = nys_x, nyn_x
296                sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
297             ENDDO
298
299          ENDDO
300         
301       ELSE
302
303          sums_spectra_l = 1.0
304
305       ENDIF
306
307!
308!--    Global sum of spectra on PE0 (from where they are written on file)
309       sums_spectra(:,n) = 0.0
310#if defined( __parallel )   
311       CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
312       CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1, &
313                        MPI_REAL, MPI_PROD, 0, comm2d, ierr )
314#else
315       sums_spectra(:,n) = sums_spectra_l
316#endif
317
318       n = n + 1
319
320    ENDDO
321    n = n - 1
322
323    IF ( myid == 0 )  THEN
324!
325!--    Sum of spectra for later averaging (see routine data_output_spectra)
326!--    Temperton fft results need to be normalized
327       IF ( fft_method == 'temperton-algorithm' )  THEN
328          fac = nx + 1.0
329       ELSE
330          fac = 1.0
331       ENDIF
332       DO  i = 1, nx/2
333          DO k = 1, n
334             spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k) * fac
335          ENDDO
336       ENDDO
337
338    ENDIF
339
340!
341!-- n_sp_x is needed by data_output_spectra_x
342    n_sp_x = n
343
344 END SUBROUTINE calc_spectra_x
345
346
347 SUBROUTINE calc_spectra_y( ddd, pr, m )
348
349    USE arrays_3d
350    USE constants
351    USE control_parameters
352    USE fft_xy
353    USE grid_variables
354    USE indices
355    USE pegrid
356    USE spectrum
357    USE statistics
358    USE transpose_indices
359
360    IMPLICIT NONE
361
362    INTEGER :: i, j, jshape(1), k, m, n, pr
363
364    REAL                       ::  fac, exponent
365    REAL, DIMENSION(0:ny)      ::  work
366    REAL, DIMENSION(0:ny/2)    ::  sums_spectra_l
367    REAL, DIMENSION(0:ny/2,100)::  sums_spectra
368
369    REAL, DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd
370
371
372!
373!-- Exponent for geometric average
374    exponent = 1.0 / ( nx + 1.0 )
375
376!
377!-- Loop over all levels defined by the user
378    n = 1
379    DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
380
381       k = comp_spectra_level(n)
382
383!
384!--    Calculate FFT only if the corresponding level is situated on this PE
385       IF ( k >= nzb_yd  .AND.  k <= nzt_yd )  THEN
386         
387          DO  i = nxl_yd, nxr_yd
388
389             work = ddd(0:ny,i,k)
390             CALL fft_y( work, 'forward' )
391
392             ddd(0,i,k) = dy * work(0)**2
393             DO  j = 1, ny/2
394                ddd(j,i,k) = dy * ( work(j)**2 + work(ny+1-j)**2 )
395             ENDDO
396
397          ENDDO
398
399!
400!--       Local sum and geometric average of these spectra
401!--       (WARNING: no global sum should be performed, because floating
402!--       point overflow may occur)
403          DO  j = 0, ny/2
404
405             sums_spectra_l(j) = 1.0
406             DO  i = nxl_yd, nxr_yd
407                sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
408             ENDDO
409
410          ENDDO
411         
412       ELSE
413
414          sums_spectra_l = 1.0
415
416       ENDIF
417
418!
419!--    Global sum of spectra on PE0 (from where they are written on file)
420       sums_spectra(:,n) = 0.0
421#if defined( __parallel )   
422       CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
423       CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1, &
424                        MPI_REAL, MPI_PROD, 0, comm2d, ierr )
425#else
426       sums_spectra(:,n) = sums_spectra_l
427#endif
428
429       n = n + 1
430
431    ENDDO
432    n = n - 1
433
434
435    IF ( myid == 0 )  THEN
436!
437!--    Sum of spectra for later averaging (see routine data_output_spectra)
438!--    Temperton fft results need to be normalized
439       IF ( fft_method == 'temperton-algorithm' )  THEN
440          fac = ny + 1.0
441       ELSE
442          fac = 1.0
443       ENDIF
444       DO  j = 1, ny/2
445          DO k = 1, n
446             spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k) * fac
447          ENDDO
448       ENDDO
449
450    ENDIF
451
452!
453!-- n_sp_y is needed by data_output_spectra_y
454    n_sp_y = n
455
456 END SUBROUTINE calc_spectra_y
457#endif
Note: See TracBrowser for help on using the repository browser.