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

Last change on this file since 707 was 707, checked in by raasch, 11 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:


Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:


localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

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