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

Last change on this file since 1379 was 1343, checked in by kanani, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 14.9 KB
Line 
1 SUBROUTINE calc_spectra
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: calc_spectra.f90 1343 2014-03-26 17:07:58Z raasch $
27!
28! 1342 2014-03-26 17:04:47Z kanani
29! REAL constants defined as wp-kinds
30!
31! 1324 2014-03-21 09:13:16Z suehring
32! Bugfix: nzb_x, nzb_yd, nyn_x, nyn_x, nzt_x, nzt_yd belong to transpose_indices
33!
34! 1320 2014-03-20 08:40:49Z raasch
35! ONLY-attribute added to USE-statements,
36! kind-parameters added to all INTEGER and REAL declaration statements,
37! kinds are defined in new module kinds,
38! revision history before 2012 removed,
39! comment fields (!:) to be used for variable explanations added to
40! all variable declaration statements
41!
42! 1318 2014-03-17 13:35:16Z raasch
43! module interfaces removed
44!
45! 1216 2013-08-26 09:31:42Z raasch
46! resorting of array moved to separate routine resort_for_zx,
47! one argument removed from the transpose_..d routines
48!
49! 1120 2013-04-05 15:11:35Z raasch
50! bugfix: calls of fft_x|y replaced by fft_x|y_1d
51!
52! 1036 2012-10-22 13:43:42Z raasch
53! code put under GPL (PALM 3.9)
54!
55! 1003 2012-09-14 14:35:53Z raasch
56! adjustment of array tend for cases with unequal subdomain sizes removed
57!
58! Revision 1.1  2001/01/05 15:08:07  raasch
59! Initial revision
60!
61!
62! Description:
63! ------------
64! Calculate horizontal spectra along x and y.
65! ATTENTION: 1d-decomposition along y still needs improvement, because in that
66!            case the gridpoint number along z still depends on the PE number
67!            because transpose_xz has to be used (and possibly also
68!            transpose_zyd needs modification).
69!------------------------------------------------------------------------------!
70
71#if defined( __spectra )
72    USE arrays_3d,                                                             &
73        ONLY:  d, tend
74
75    USE control_parameters,                                                    &
76        ONLY:  average_count_sp, bc_lr_cyc, bc_ns_cyc, message_string, psolver
77
78    USE cpulog,                                                                &
79        ONLY:  cpu_log, log_point
80
81    USE fft_xy,                                                                &
82        ONLY:  fft_init
83
84    USE indices,                                                               &
85        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
86
87    USE kinds
88
89    USE pegrid
90
91    USE spectrum,                                                              &
92        ONLY:  data_output_sp, spectra_direction
93
94
95    IMPLICIT NONE
96
97    INTEGER(iwp) ::  m  !:
98    INTEGER(iwp) ::  pr !:
99
100
101    CALL cpu_log( log_point(30), 'calc_spectra', 'start' )
102
103!
104!-- Initialize ffts
105    CALL fft_init
106
107!
108!-- Reallocate array d in required size
109    IF ( psolver == 'multigrid' )  THEN
110       DEALLOCATE( d )
111       ALLOCATE( d(nzb+1:nzt,nys:nyn,nxl:nxr) )
112    ENDIF
113
114    m = 1
115    DO WHILE ( data_output_sp(m) /= ' '  .AND.  m <= 10 )
116!
117!--    Transposition from z --> x  ( y --> x in case of a 1d-decomposition
118!--    along x)
119       IF ( INDEX( spectra_direction(m), 'x' ) /= 0 )  THEN
120
121!
122!--       Calculation of spectra works for cyclic boundary conditions only
123          IF ( .NOT. bc_lr_cyc )  THEN
124
125             message_string = 'non-cyclic lateral boundaries along x do not'// &
126                              '& allow calculation of spectra along x'
127             CALL message( 'calc_spectra', 'PA0160', 1, 2, 0, 6, 0 )
128          ENDIF
129
130          CALL preprocess_spectra( m, pr )
131
132#if defined( __parallel )
133          IF ( pdims(2) /= 1 )  THEN
134             CALL resort_for_zx( d, tend )
135             CALL transpose_zx( tend, d )
136          ELSE
137             CALL transpose_yxd( d, d )
138          ENDIF
139          CALL calc_spectra_x( d, pr, m )
140#else
141          message_string = 'sorry, calculation of spectra in non parallel ' // &
142                           'mode& is still not realized'
143          CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )     
144#endif
145
146       ENDIF
147
148!
149!--    Transposition from z --> y (d is rearranged only in case of a
150!--    1d-decomposition along x)
151       IF ( INDEX( spectra_direction(m), 'y' ) /= 0 )  THEN
152
153!
154!--       Calculation of spectra works for cyclic boundary conditions only
155          IF ( .NOT. bc_ns_cyc )  THEN
156             IF ( myid == 0 )  THEN
157                message_string = 'non-cyclic lateral boundaries along y do' // &
158                                 ' not & allow calculation of spectra along y' 
159                CALL message( 'calc_spectra', 'PA0162', 1, 2, 0, 6, 0 )
160             ENDIF
161             CALL local_stop
162          ENDIF
163
164          CALL preprocess_spectra( m, pr )
165
166#if defined( __parallel )
167          CALL transpose_zyd( d, d )
168          CALL calc_spectra_y( d, pr, m )
169#else
170          message_string = 'sorry, calculation of spectra in non parallel' //  &
171                           'mode& is still not realized'
172          CALL message( 'calc_spectra', 'PA0161', 1, 2, 0, 6, 0 )
173#endif
174
175       ENDIF
176
177!
178!--    Increase counter for next spectrum
179       m = m + 1
180         
181    ENDDO
182
183!
184!-- Increase counter for averaging process in routine plot_spectra
185    average_count_sp = average_count_sp + 1
186
187    CALL cpu_log( log_point(30), 'calc_spectra', 'stop' )
188
189#endif
190 END SUBROUTINE calc_spectra
191
192
193#if defined( __spectra )
194 SUBROUTINE preprocess_spectra( m, pr )
195
196    USE arrays_3d,                                                             &
197        ONLY:  d, pt, q, u, v, w
198
199    USE indices,                                                               &
200        ONLY:  nxl, nxr, nyn, nys, nzb, nzt
201
202    USE kinds
203
204    USE pegrid
205
206    USE spectrum,                                                              &
207        ONLY:  data_output_sp
208
209    USE statistics,                                                            &
210        ONLY:  sums
211
212
213    IMPLICIT NONE
214
215    INTEGER(iwp) :: i  !:
216    INTEGER(iwp) :: j  !:
217    INTEGER(iwp) :: k  !:
218    INTEGER(iwp) :: m  !:
219    INTEGER(iwp) :: pr !:
220
221    SELECT CASE ( TRIM( data_output_sp(m) ) )
222         
223    CASE ( 'u' )
224       pr = 1
225       d(nzb+1:nzt,nys:nyn,nxl:nxr) = u(nzb+1:nzt,nys:nyn,nxl:nxr)
226       
227    CASE ( 'v' )
228       pr = 2
229       d(nzb+1:nzt,nys:nyn,nxl:nxr) = v(nzb+1:nzt,nys:nyn,nxl:nxr)
230       
231    CASE ( 'w' )
232       pr = 3
233       d(nzb+1:nzt,nys:nyn,nxl:nxr) = w(nzb+1:nzt,nys:nyn,nxl:nxr)
234       
235    CASE ( 'pt' )
236       pr = 4
237       d(nzb+1:nzt,nys:nyn,nxl:nxr) = pt(nzb+1:nzt,nys:nyn,nxl:nxr)
238       
239    CASE ( 'q' )
240       pr = 41
241       d(nzb+1:nzt,nys:nyn,nxl:nxr) = q(nzb+1:nzt,nys:nyn,nxl:nxr)
242       
243    CASE DEFAULT
244!
245!--    The DEFAULT case is reached either if the parameter data_output_sp(m)
246!--    contains a wrong character string or if the user has coded a special
247!--    case in the user interface. There, the subroutine user_spectra
248!--    checks which of these two conditions applies.
249       CALL user_spectra( 'preprocess', m, pr )
250         
251    END SELECT
252
253!
254!-- Subtract horizontal mean from the array, for which spectra have to be
255!-- calculated
256    DO  i = nxl, nxr
257       DO  j = nys, nyn
258          DO  k = nzb+1, nzt
259             d(k,j,i) = d(k,j,i) - sums(k,pr)
260          ENDDO
261       ENDDO
262    ENDDO
263
264 END SUBROUTINE preprocess_spectra
265
266
267 SUBROUTINE calc_spectra_x( ddd, pr, m )
268
269    USE arrays_3d,                                                             &
270        ONLY: 
271
272    USE control_parameters,                                                    &
273        ONLY:  fft_method
274
275    USE fft_xy,                                                                &
276        ONLY:  fft_x_1d
277
278    USE grid_variables,                                                        &
279        ONLY:  dx
280
281    USE indices,                                                               &
282        ONLY:  nx, ny
283
284    USE kinds
285
286    USE pegrid
287
288    USE spectrum,                                                              &
289        ONLY:  comp_spectra_level, n_sp_x
290
291    USE statistics,                                                            &
292        ONLY:  spectrum_x
293
294    USE transpose_indices,                                                     &
295        ONLY:  nyn_x, nys_x, nzb_x, nzt_x
296
297
298    IMPLICIT NONE
299
300    INTEGER(iwp) ::  i         !:
301    INTEGER(iwp) ::  ishape(1) !:
302    INTEGER(iwp) ::  j         !:
303    INTEGER(iwp) ::  k         !:
304    INTEGER(iwp) ::  m         !:
305    INTEGER(iwp) ::  n         !:
306    INTEGER(iwp) ::  pr        !:
307
308    REAL(wp) ::  fac      !:
309    REAL(wp) ::  exponent !:
310   
311    REAL(wp), DIMENSION(0:nx) ::  work !:
312   
313    REAL(wp), DIMENSION(0:nx/2) ::  sums_spectra_l !:
314   
315    REAL(wp), DIMENSION(0:nx/2,100) ::  sums_spectra !:
316   
317    REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ddd !:
318
319!
320!-- Exponent for geometric average
321    exponent = 1.0_wp / ( ny + 1.0_wp )
322
323!
324!-- Loop over all levels defined by the user
325    n = 1
326    DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
327
328       k = comp_spectra_level(n)
329
330!
331!--    Calculate FFT only if the corresponding level is situated on this PE
332       IF ( k >= nzb_x  .AND.  k <= nzt_x )  THEN
333         
334          DO  j = nys_x, nyn_x
335
336             work = ddd(0:nx,j,k)
337             CALL fft_x_1d( work, 'forward' )
338
339             ddd(0,j,k) = dx * work(0)**2
340             DO  i = 1, nx/2
341                ddd(i,j,k) = dx * ( work(i)**2 + work(nx+1-i)**2 )
342             ENDDO
343
344          ENDDO
345
346!
347!--       Local sum and geometric average of these spectra
348!--       (WARNING: no global sum should be performed, because floating
349!--       point overflow may occur)
350          DO  i = 0, nx/2
351
352             sums_spectra_l(i) = 1.0_wp
353             DO  j = nys_x, nyn_x
354                sums_spectra_l(i) = sums_spectra_l(i) * ddd(i,j,k)**exponent
355             ENDDO
356
357          ENDDO
358         
359       ELSE
360
361          sums_spectra_l = 1.0_wp
362
363       ENDIF
364
365!
366!--    Global sum of spectra on PE0 (from where they are written on file)
367       sums_spectra(:,n) = 0.0_wp
368#if defined( __parallel )   
369       CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
370       CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), nx/2+1,          &
371                        MPI_REAL, MPI_PROD, 0, comm2d, ierr )
372#else
373       sums_spectra(:,n) = sums_spectra_l
374#endif
375
376       n = n + 1
377
378    ENDDO
379    n = n - 1
380
381    IF ( myid == 0 )  THEN
382!
383!--    Sum of spectra for later averaging (see routine data_output_spectra)
384!--    Temperton fft results need to be normalized
385       IF ( fft_method == 'temperton-algorithm' )  THEN
386          fac = nx + 1.0_wp
387       ELSE
388          fac = 1.0_wp
389       ENDIF
390       DO  i = 1, nx/2
391          DO k = 1, n
392             spectrum_x(i,k,m) = spectrum_x(i,k,m) + sums_spectra(i,k) * fac
393          ENDDO
394       ENDDO
395
396    ENDIF
397
398!
399!-- n_sp_x is needed by data_output_spectra_x
400    n_sp_x = n
401
402 END SUBROUTINE calc_spectra_x
403
404
405 SUBROUTINE calc_spectra_y( ddd, pr, m )
406
407    USE arrays_3d,                                                             &
408        ONLY: 
409
410    USE control_parameters,                                                    &
411        ONLY:  fft_method
412
413    USE fft_xy,                                                                &
414        ONLY:  fft_y_1d
415
416    USE grid_variables,                                                        &
417        ONLY:  dy
418
419    USE indices,                                                               &
420        ONLY:  nx, ny
421
422    USE kinds
423
424    USE pegrid
425
426    USE spectrum,                                                              &
427        ONLY:  comp_spectra_level, n_sp_y
428
429    USE statistics,                                                            &
430        ONLY:  spectrum_y
431
432    USE transpose_indices,                                                     &
433        ONLY:  nxl_yd, nxr_yd, nzb_yd, nzt_yd
434
435
436    IMPLICIT NONE
437
438    INTEGER(iwp) ::  i         !:
439    INTEGER(iwp) ::  j         !:
440    INTEGER(iwp) ::  jshape(1) !:
441    INTEGER(iwp) ::  k         !:
442    INTEGER(iwp) ::  m         !:
443    INTEGER(iwp) ::  n         !:
444    INTEGER(iwp) ::  pr        !:
445
446    REAL(wp) ::  fac      !:
447    REAL(wp) ::  exponent !:
448   
449    REAL(wp), DIMENSION(0:ny) ::  work !:
450   
451    REAL(wp), DIMENSION(0:ny/2) ::  sums_spectra_l !:
452   
453    REAL(wp), DIMENSION(0:ny/2,100) ::  sums_spectra !:
454   
455    REAL(wp), DIMENSION(0:ny,nxl_yd:nxr_yd,nzb_yd:nzt_yd) :: ddd !:
456
457
458!
459!-- Exponent for geometric average
460    exponent = 1.0_wp / ( nx + 1.0_wp )
461
462!
463!-- Loop over all levels defined by the user
464    n = 1
465    DO WHILE ( comp_spectra_level(n) /= 999999  .AND.  n <= 100 )
466
467       k = comp_spectra_level(n)
468
469!
470!--    Calculate FFT only if the corresponding level is situated on this PE
471       IF ( k >= nzb_yd  .AND.  k <= nzt_yd )  THEN
472         
473          DO  i = nxl_yd, nxr_yd
474
475             work = ddd(0:ny,i,k)
476             CALL fft_y_1d( work, 'forward' )
477
478             ddd(0,i,k) = dy * work(0)**2
479             DO  j = 1, ny/2
480                ddd(j,i,k) = dy * ( work(j)**2 + work(ny+1-j)**2 )
481             ENDDO
482
483          ENDDO
484
485!
486!--       Local sum and geometric average of these spectra
487!--       (WARNING: no global sum should be performed, because floating
488!--       point overflow may occur)
489          DO  j = 0, ny/2
490
491             sums_spectra_l(j) = 1.0_wp
492             DO  i = nxl_yd, nxr_yd
493                sums_spectra_l(j) = sums_spectra_l(j) * ddd(j,i,k)**exponent
494             ENDDO
495
496          ENDDO
497         
498       ELSE
499
500          sums_spectra_l = 1.0_wp
501
502       ENDIF
503
504!
505!--    Global sum of spectra on PE0 (from where they are written on file)
506       sums_spectra(:,n) = 0.0_wp
507#if defined( __parallel )   
508       CALL MPI_BARRIER( comm2d, ierr )  ! Necessary?
509       CALL MPI_REDUCE( sums_spectra_l(0), sums_spectra(0,n), ny/2+1,          &
510                        MPI_REAL, MPI_PROD, 0, comm2d, ierr )
511#else
512       sums_spectra(:,n) = sums_spectra_l
513#endif
514
515       n = n + 1
516
517    ENDDO
518    n = n - 1
519
520
521    IF ( myid == 0 )  THEN
522!
523!--    Sum of spectra for later averaging (see routine data_output_spectra)
524!--    Temperton fft results need to be normalized
525       IF ( fft_method == 'temperton-algorithm' )  THEN
526          fac = ny + 1.0_wp
527       ELSE
528          fac = 1.0_wp
529       ENDIF
530       DO  j = 1, ny/2
531          DO k = 1, n
532             spectrum_y(j,k,m) = spectrum_y(j,k,m) + sums_spectra(j,k) * fac
533          ENDDO
534       ENDDO
535
536    ENDIF
537
538!
539!-- n_sp_y is needed by data_output_spectra_y
540    n_sp_y = n
541
542 END SUBROUTINE calc_spectra_y
543#endif
Note: See TracBrowser for help on using the repository browser.