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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

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