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

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

bugfixes concerning cpu time measurements and calculation of spectra in case of multigrid method switched on

  • Property svn:keywords set to Id
File size: 10.9 KB
Line 
1 SUBROUTINE calc_spectra
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Bugfix: array d is reallocated in case that multigrid is used
7!
8! Former revisions:
9! -----------------
10! $Id: calc_spectra.f90 225 2009-01-26 14:44:20Z raasch $
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!-- Reallocate array d in required size
60    IF ( psolver == 'multigrid' )  THEN
61       DEALLOCATE( d )
62       ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra) )
63    ENDIF
64
65!
66!-- Enlarge the size of tend, used as a working array for the transpositions
67    IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
68       DEALLOCATE( tend )
69       ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
70    ENDIF
71
72    m = 1
73    DO WHILE ( data_output_sp(m) /= ' '  .AND.  m <= 10 )
74!
75!--    Transposition from z --> x  ( y --> x in case of a 1d-decomposition
76!--    along x)
77       IF ( INDEX( spectra_direction(m), 'x' ) /= 0 )  THEN
78
79!
80!--       Calculation of spectra works for cyclic boundary conditions only
81          IF ( bc_lr /= 'cyclic' )  THEN
82             IF ( myid == 0 )  THEN
83                PRINT*, '+++ calc_spectra:'
84                PRINT*, '    non-cyclic lateral boundaries along x do not ', &
85                             'allow calculation of spectra along x' 
86             ENDIF
87             CALL local_stop
88          ENDIF
89
90          CALL preprocess_spectra( m, pr )
91
92#if defined( __parallel )
93          IF ( pdims(2) /= 1 )  THEN
94             CALL transpose_zx( d, tend, d )
95          ELSE
96             CALL transpose_yxd( d, tend, d )
97          ENDIF
98          CALL calc_spectra_x( d, pr, m )
99#else
100          PRINT*, '+++ calc_spectra: sorry, calculation of spectra ', &
101                                    'in non parallel mode'
102          PRINT*, '                  is still not realized'
103          CALL local_stop
104#endif
105
106       ENDIF
107
108!
109!--    Transposition from z --> y (d is rearranged only in case of a
110!--    1d-decomposition along x)
111       IF ( INDEX( spectra_direction(m), 'y' ) /= 0 )  THEN
112
113!
114!--       Calculation of spectra works for cyclic boundary conditions only
115          IF ( bc_ns /= 'cyclic' )  THEN
116             IF ( myid == 0 )  THEN
117                PRINT*, '+++ calc_spectra:'
118                PRINT*, '    non-cyclic lateral boundaries along y do not ', &
119                             'allow calculation of spectra along y' 
120             ENDIF
121             CALL local_stop
122          ENDIF
123
124          CALL preprocess_spectra( m, pr )
125
126#if defined( __parallel )
127          CALL transpose_zyd( d, tend, d )
128          CALL calc_spectra_y( d, pr, m )
129#else
130          PRINT*, '+++ calc_spectra: sorry, calculation of spectra', &
131                                    'in non parallel mode'
132          PRINT*, '                  still not realized'
133          CALL local_stop
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.