Home

Context Navigation

calc_spectra.f90 @ 667

Last change on this file since 667 was 667, checked in by suehring, 13 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.

Context Navigation

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

Download in other formats: