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

Last change on this file since 343 was 274, checked in by heinze, 15 years ago

Indentation of the message calls corrected

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