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

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

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:


lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:


Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

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