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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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