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

Last change on this file since 667 was 667, checked in by suehring, 11 years ago

summary:


Gryschka:

  • Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
  • Exchange of u and v from ocean surface to atmosphere surface
  • Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
  • Inflow turbulence is now defined by flucuations around spanwise mean
  • Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

  • New advection added ( Wicker and Skamarock 5th order ), therefore:
    • New module advec_ws.f90
    • Modified exchange of ghost boundaries.
    • Modified evaluation of turbulent fluxes
    • New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90


Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90


Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90


Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90


When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90


Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90


Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90


Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90


Steering parameter dissipation_control added in inipar.

Makefile


Module advec_ws added.

Modules


Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90


New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90


Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90


Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90


Output of total array size was adapted to nbgp.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90


Calls of exchange_horiz are modified.

diffusion_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90


Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90


Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90


determination of target_id's moved to init_pegrid

init_pt_anomaly.f90


Call of exchange_horiz are modified.

init_rankine.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90


Output of advection scheme.

poismg.f90


Calls of exchange_horiz are modified.

prandtl_fluxes.f90


Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90


+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90


Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90


Allocation with nbgp.

wall_fluxes.f90


nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90


Array bounds and nx, ny adapted with nbgp.

sor.f90


bug removed in declaration of ddzw(), nz replaced by nzt+1

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