source: palm/tags/release-3.5/SOURCE/calc_spectra.f90 @ 3999

Last change on this file since 3999 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

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