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

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