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

Last change on this file since 230 was 226, checked in by raasch, 16 years ago

preparations for the next release

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