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

Last change on this file since 799 was 708, checked in by raasch, 13 years ago

last commit documented

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