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

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