source: palm/trunk/SOURCE/data_output_spectra.f90 @ 170

Last change on this file since 170 was 144, checked in by letzel, 16 years ago

User-defined spectra.

Bugfix: extra '*' removed in user_statistics sample code (user_interface).

  • Property svn:keywords set to Id
File size: 16.3 KB
Line 
1 SUBROUTINE data_output_spectra
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! user-defined spectra
7!
8! Former revisions:
9! -----------------
10! $Id: data_output_spectra.f90 144 2008-01-04 04:29:45Z raasch $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.7  2006/04/11 14:56:38  raasch
14! pl_spectra renamed data_output_sp
15!
16! Revision 1.1  2001/01/05 15:14:20  raasch
17! Initial revision
18!
19!
20! Description:
21! ------------
22! Writing spectra data on file, using a special format which allows
23! plotting of these data with PROFIL-graphic-software
24!------------------------------------------------------------------------------!
25#if defined( __spectra )
26
27    USE arrays_3d
28    USE control_parameters
29    USE cpulog
30    USE interfaces
31    USE netcdf_control
32    USE pegrid
33    USE spectrum
34    USE statistics
35
36
37    IMPLICIT NONE
38
39    INTEGER :: m, pr, cranz_x, cranz_y
40    LOGICAL :: frame_x, frame_y
41
42    CALL cpu_log( log_point(31), 'data_output_spectra', 'start' )
43
44!
45!-- Output is only performed on PE0
46    IF ( myid == 0 )  THEN
47
48!
49!--    Open file for spectra output in NetCDF format
50       IF ( netcdf_output )  CALL check_open( 107 )
51
52!
53!--    Increment the counter for number of output times
54       dosp_time_count = dosp_time_count + 1
55
56#if defined( __netcdf )
57!
58!--    Update the spectra time axis
59       nc_stat = NF90_PUT_VAR( id_set_sp, id_var_time_sp, (/ simulated_time /),&
60                               start = (/ dosp_time_count /), count = (/ 1 /) )
61       IF (nc_stat /= NF90_NOERR)  CALL handle_netcdf_error( 47 )
62#endif
63
64       IF ( profil_output )  THEN
65!
66!--       Compute RAHMEN-Parameter CRANZ for x- and y-spectra separately
67          cranz_x = 0; cranz_y = 0; frame_x = .FALSE.; frame_y = .FALSE.
68
69          m = 1
70          DO WHILE ( data_output_sp(m) /= ' ' .AND. m <= 10 )
71
72             IF ( INDEX( spectra_direction(m), 'x' ) /= 0 )  THEN
73                cranz_x = cranz_x + 1
74             ENDIF
75             IF ( INDEX( spectra_direction(m), 'y' ) /= 0 )  THEN
76                cranz_y = cranz_y + 1
77             ENDIF
78
79             m = m + 1
80
81          ENDDO
82
83       ENDIF
84
85!
86!--    If necessary, calculate time average and reset average counter
87       IF ( average_count_sp == 0 )  THEN
88          PRINT*, '+++ data_output_spectra: no spectra data available'
89       ENDIF
90       IF ( average_count_sp /= 1 )  THEN
91          spectrum_x = spectrum_x / REAL( average_count_sp )
92          spectrum_y = spectrum_y / REAL( average_count_sp )
93          average_count_sp = 0
94       ENDIF
95
96!
97!--    Loop over all spectra defined by the user
98       m = 1
99       DO WHILE ( data_output_sp(m) /= ' '  .AND.  m <= 10 )
100
101          SELECT CASE ( TRIM( data_output_sp(m) ) )
102
103             CASE ( 'u' )
104                pr = 1
105
106             CASE ( 'v' )
107                pr = 2
108
109             CASE ( 'w' )
110                pr = 3
111
112             CASE ( 'pt' )
113                pr = 4
114
115             CASE ( 'q' )
116                pr = 5
117
118             CASE DEFAULT
119!
120!--             The DEFAULT case is reached either if the parameter
121!--             data_output_sp(m) contains a wrong character string or if the
122!--             user has coded a special case in the user interface. There, the
123!--             subroutine user_spectra checks which of these two conditions
124!--             applies.
125                CALL user_spectra( 'data_output', m, pr )
126
127          END SELECT
128
129!
130!--       Output of spectra in NetCDF format
131          IF ( netcdf_output )  THEN
132!
133!--          Output of x-spectra
134             IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN
135                CALL output_spectra_netcdf( m, 'x' )
136             ENDIF
137!
138!--          Output of y-spectra
139             IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN
140                CALL output_spectra_netcdf( m, 'y' )
141             ENDIF
142          ENDIF
143
144!
145!--       Output of spectra in profil format
146          IF ( profil_output )  THEN
147!
148!--          Output of x-spectra
149             IF ( INDEX( spectra_direction(m), 'x' ) /= 0 ) THEN
150                CALL data_output_spectra_x( m, cranz_x, pr, frame_x )
151             ENDIF
152
153!
154!--          Output of y-spectra
155             IF ( INDEX( spectra_direction(m), 'y' ) /= 0 ) THEN
156                CALL data_output_spectra_y( m, cranz_y, pr, frame_y )
157             ENDIF
158          ENDIF
159
160!
161!--       Increase counter for next spectrum
162          m = m + 1
163
164       ENDDO
165
166!
167!--    Reset spectra values
168       spectrum_x = 0.0; spectrum_y = 0.0
169
170    ENDIF
171
172    CALL cpu_log( log_point(31), 'data_output_spectra', 'stop' )
173
174#if defined( __parallel )
175!    CALL MPI_BARRIER( comm2d, ierr )  ! really necessary
176#endif
177
178#endif
179 END SUBROUTINE data_output_spectra
180
181
182 SUBROUTINE output_spectra_netcdf( nsp, direction )
183#if defined( __netcdf )
184
185    USE constants
186    USE control_parameters
187    USE grid_variables
188    USE indices
189    USE netcdf_control
190    USE spectrum
191    USE statistics
192
193    IMPLICIT NONE
194
195    CHARACTER (LEN=1), INTENT(IN) ::  direction
196
197    INTEGER, INTENT(IN) ::  nsp
198
199    INTEGER ::  i, k
200
201    REAL ::  frequency
202
203    REAL, DIMENSION(nx/2) ::  netcdf_data_x
204    REAL, DIMENSION(ny/2) ::  netcdf_data_y
205
206
207    IF ( direction == 'x' )  THEN
208
209       DO  k = 1, n_sp_x
210
211          DO  i = 1, nx/2
212             frequency = 2.0 * pi * i / ( dx * ( nx + 1 ) )
213             netcdf_data_x(i) = frequency * spectrum_x(i,k,nsp)
214          ENDDO
215
216          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_dospx(nsp), netcdf_data_x, &
217                                  start = (/ 1, k, dosp_time_count /), &
218                                  count = (/ nx/2, 1, 1 /) )
219          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 348 )
220
221       ENDDO
222
223    ENDIF
224
225    IF ( direction == 'y' )  THEN
226
227       DO  k = 1, n_sp_y
228
229          DO  i = 1, ny/2
230             frequency = 2.0 * pi * i / ( dy * ( ny + 1 ) )
231             netcdf_data_y(i) = frequency * spectrum_y(i,k,nsp)
232          ENDDO
233
234          nc_stat = NF90_PUT_VAR( id_set_sp, id_var_dospy(nsp), netcdf_data_y, &
235                                  start = (/ 1, k, dosp_time_count /), &
236                                  count = (/ ny/2, 1, 1 /) )
237          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 349 )
238
239       ENDDO
240
241    ENDIF
242
243#endif
244 END SUBROUTINE output_spectra_netcdf
245
246
247#if defined( __spectra )
248 SUBROUTINE data_output_spectra_x( m, cranz, pr, frame_written )
249
250    USE arrays_3d
251    USE constants
252    USE control_parameters
253    USE grid_variables
254    USE indices
255    USE pegrid
256    USE singleton
257    USE spectrum
258    USE statistics
259    USE transpose_indices
260
261    IMPLICIT NONE
262
263    CHARACTER (LEN=30) ::  atext
264    INTEGER            ::  i, j, k, m, pr
265    LOGICAL            ::  frame_written
266    REAL               ::  frequency = 0.0
267
268!
269!-- Variables needed for PROFIL-namelist
270    INTEGER                  :: cranz, labforx = 3, labfory = 3, legpos = 3, &
271                                timodex = 1
272    INTEGER, DIMENSION(1:10) :: cucol = 1, klist = 999999, lstyle = 0
273    LOGICAL                  :: datleg = .TRUE., grid = .TRUE., &
274                                lclose = .TRUE., rand = .TRUE., &
275                                swap = .TRUE., twoxa = .TRUE.,  &
276                                xlog = .TRUE., ylog = .TRUE.
277    CHARACTER (LEN=80)       :: rtext, utext, xtext = 'k in m>->1', ytext
278    REAL                     :: gwid = 0.1, rlegfak = 0.7, uxmin, uxmax, &
279                                uymin, uymax
280    REAL, DIMENSION(1:10)    :: lwid = 0.6
281    REAL, DIMENSION(10)      :: uyma, uymi
282
283    NAMELIST /RAHMEN/  cranz, datleg, rtext, swap
284    NAMELIST /CROSS/   rand, cucol, grid, gwid, klist, labforx, labfory,      &
285                       legpos, lclose, lstyle, lwid, rlegfak, timodex, utext, &
286                       uxmin, uxmax, uymin, uymax, twoxa, xlog, xtext, ylog,  &
287                       ytext
288
289
290    rtext = '\0.5 ' // run_description_header
291
292!
293!-- Open parameter- and data-file
294    CALL check_open( 81 )
295    CALL check_open( 82 )
296
297!
298!-- Write file header,
299!-- write RAHMEN-parameters (pr=3: w-array is on zw, other arrys on zu,
300!-- pr serves as an index for output of strings (axis-labels) of the
301!-- different quantities u, v, w, pt and q)
302    DO  k = 1, n_sp_x
303       IF ( k < 10 )  THEN
304          IF ( pr == 3 )  THEN
305             WRITE ( 82, 100 )  '#', k, header_char( pr ),        &
306                                INT( zw(comp_spectra_level(k)) ), &
307                                simulated_time_chr
308          ELSE
309             WRITE ( 82, 100 )  '#', k, header_char( pr ),        &
310                                INT( zu(comp_spectra_level(k)) ), &
311                                simulated_time_chr
312          ENDIF
313       ELSE
314          IF ( pr == 3 )  THEN
315             WRITE ( 82, 101 )  '#', k, header_char( pr ),        &
316                                INT( zw(comp_spectra_level(k)) ), &
317                                simulated_time_chr
318          ELSE
319             WRITE ( 82, 101 )  '#', k, header_char( pr ),        &
320                                INT( zu(comp_spectra_level(k)) ), &
321                                simulated_time_chr
322          ENDIF
323       ENDIF
324    ENDDO
325
326    IF ( .NOT. frame_written )  THEN
327       WRITE ( 81, RAHMEN )
328       frame_written = .TRUE.
329    ENDIF
330
331!
332!-- Write all data and calculate uymi and uyma. They serve to calculate
333!-- the CROSS-parameters uymin and uymax
334    uymi = 999.999; uyma = -999.999
335    DO  i = 1, nx/2
336       frequency = 2.0 * pi * i / ( dx * ( nx + 1 ) )
337       WRITE ( 82, 102 )  frequency, ( frequency * spectrum_x(i,k,m), k = 1, &
338                          n_sp_x )
339       DO  k = 1, n_sp_x
340          uymi(k) = MIN( uymi(k), frequency * spectrum_x(i,k,m) )
341          uyma(k) = MAX( uyma(k), frequency * spectrum_x(i,k,m) )
342       ENDDO
343    ENDDO
344
345!
346!-- Determine CROSS-parameters
347    cucol(1:n_sp_x)  = (/ ( k, k = 1, n_sp_x ) /)
348    lstyle(1:n_sp_x) = (/ ( lstyles(k), k = 1, n_sp_x ) /)
349
350!
351!-- Calculate klist-values from the available comp_spectra_level values
352    i = 1; k = 1
353    DO WHILE ( i <= 10  .AND.  plot_spectra_level(i) /= 999999 )
354       DO WHILE ( k <= n_sp_x  .AND. &
355                  plot_spectra_level(i) >= comp_spectra_level(k) )
356          IF ( plot_spectra_level(i) == comp_spectra_level(k) )  THEN
357             klist(i) = k + klist_x
358          ELSE
359             uymi(k) =  999.999
360             uyma(k) = -999.999
361          ENDIF
362          k = k + 1
363       ENDDO
364       i = i + 1
365    ENDDO
366    uymi(k:n_sp_x) =  999.999
367    uyma(k:n_sp_x) = -999.999
368    utext = 'x'//utext_char( pr )
369    IF ( averaging_interval_sp /= 0.0 ) THEN
370       WRITE ( atext, 104 )  averaging_interval_sp
371       utext = TRIM(utext) // ',  ' // TRIM( atext )
372    ENDIF
373    uxmin = 0.8 * 2.0 * pi        / ( dx * ( nx + 1 ) )
374    uxmax = 1.2 * 2.0 * pi * nx/2 / ( dx * ( nx + 1 ) )
375    uymin = 0.8 * MIN (  999.999, MINVAL ( uymi ) )
376    uymax = 1.2 * MAX ( -999.999, MAXVAL ( uyma ) )
377    ytext = ytext_char( pr )
378
379!
380!-- Output of CROSS-parameters
381    WRITE ( 81, CROSS )
382
383!
384!-- Increase counter by the number of profiles written in the actual block
385    klist_x = klist_x + n_sp_x
386
387!
388!-- Write end-mark
389    WRITE ( 82, 103 )
390
391!
392!-- Close parameter- and data-file
393    CALL close_file( 81 )
394    CALL close_file( 82 )
395
396!
397!-- Formats
398100 FORMAT (A,I1,1X,A,1X,I4,'m ',A)
399101 FORMAT (A,I2,1X,A,1X,I4,'m ',A)
400102 FORMAT (E15.7,10(1X,E15.7))
401103 FORMAT ('NEXT')
402104 FORMAT ('time averaged over',F7.1,' s')
403
404 END SUBROUTINE data_output_spectra_x
405
406
407 SUBROUTINE data_output_spectra_y( m, cranz, pr, frame_written )
408
409    USE arrays_3d
410    USE constants
411    USE control_parameters
412    USE grid_variables
413    USE indices
414    USE pegrid
415    USE singleton
416    USE spectrum
417    USE statistics
418    USE transpose_indices
419
420    IMPLICIT NONE
421
422    CHARACTER (LEN=30) ::  atext
423    INTEGER            :: i, j, k, m, pr
424    LOGICAL            :: frame_written
425    REAL               :: frequency = 0.0
426
427!
428!-- Variables needed for PROFIL-namelist
429    INTEGER                  :: cranz, labforx = 3, labfory = 3, legpos = 3, &
430                                timodex = 1
431    INTEGER, DIMENSION(1:10) :: cucol = 1, klist = 999999, lstyle = 0
432    LOGICAL                  :: datleg = .TRUE., grid = .TRUE., &
433                                lclose = .TRUE., rand = .TRUE., &
434                                swap = .TRUE., twoxa = .TRUE.,  &
435                                xlog = .TRUE., ylog = .TRUE.
436    CHARACTER (LEN=80)       :: rtext, utext, xtext = 'k in m>->1', ytext
437    REAL                     :: gwid = 0.1, rlegfak = 0.7, uxmin, uxmax, &
438                                uymin, uymax
439    REAL, DIMENSION(1:10)    :: lwid = 0.6
440    REAL, DIMENSION(10)      :: uyma, uymi
441
442    NAMELIST /RAHMEN/  cranz, datleg, rtext, swap
443    NAMELIST /CROSS/   rand, cucol, grid, gwid, klist, labforx, labfory,      &
444                       legpos, lclose, lstyle, lwid, rlegfak, timodex, utext, &
445                       uxmin, uxmax, uymin, uymax, twoxa, xlog, xtext, ylog,  &
446                       ytext
447
448
449    rtext = '\0.5 ' // run_description_header
450
451!
452!-- Open parameter- and data-file
453    CALL check_open( 83 )
454    CALL check_open( 84 )
455
456!
457!-- Write file header,
458!-- write RAHMEN-parameters (pr=3: w-array is on zw, other arrys on zu,
459!-- pr serves as an index for output of strings (axis-labels) of the
460!-- different quantities u, v, w, pt and q)
461    DO  k = 1, n_sp_y
462       IF ( k < 10 )  THEN
463          IF ( pr == 3 ) THEN
464             WRITE ( 84, 100 )  '#', k, header_char( pr ),        &
465                                INT( zw(comp_spectra_level(k)) ), &
466                                simulated_time_chr
467          ELSE
468             WRITE ( 84, 100 )  '#', k, header_char( pr ),        &
469                                INT( zu(comp_spectra_level(k)) ), &
470                                simulated_time_chr
471          ENDIF
472       ELSE
473          IF ( pr == 3 )  THEN
474             WRITE ( 84, 101 )  '#', k, header_char( pr ),        &
475                                INT( zw(comp_spectra_level(k)) ), &
476                                simulated_time_chr
477          ELSE
478             WRITE ( 84, 101 )  '#', k, header_char( pr ),        &
479                                INT( zu(comp_spectra_level(k)) ), &
480                                simulated_time_chr
481          ENDIF
482       ENDIF
483    ENDDO
484
485    IF ( .NOT. frame_written )  THEN
486       WRITE ( 83, RAHMEN )
487       frame_written = .TRUE.
488    ENDIF
489
490!
491!-- Write all data and calculate uymi and uyma. They serve to calculate
492!-- the CROSS-parameters uymin and uymax
493    uymi = 999.999; uyma = -999.999
494    DO  j = 1, ny/2
495       frequency = 2.0 * pi * j / ( dy * ( ny + 1 ) )
496       WRITE ( 84, 102 ) frequency, ( frequency * spectrum_y(j,k,m), &
497                                      k = 1, n_sp_y ) 
498       DO k = 1, n_sp_y
499          uymi(k) = MIN( uymi(k), frequency * spectrum_y(j,k,m) )
500          uyma(k) = MAX( uyma(k), frequency * spectrum_y(j,k,m) )
501       ENDDO
502    ENDDO
503
504!
505!-- Determine CROSS-parameters
506    cucol(1:n_sp_y)  = (/ ( k, k = 1, n_sp_y ) /)
507    lstyle(1:n_sp_y) = (/ ( lstyles(k), k = 1, n_sp_y ) /)
508
509!
510!-- Calculate klist-values from the available comp_spectra_level values
511    j = 1; k = 1
512    DO WHILE ( j <= 10  .AND.  plot_spectra_level(j) /= 999999 )
513       DO WHILE ( k <= n_sp_y  .AND. &
514                  plot_spectra_level(j) >= comp_spectra_level(k) )
515          IF ( plot_spectra_level(j) == comp_spectra_level(k) )  THEN
516             klist(j) = k + klist_y
517          ELSE
518             uymi(k) =  999.999
519             uyma(k) = -999.999
520          ENDIF
521          k = k + 1
522       ENDDO
523       j = j + 1
524    ENDDO
525    uymi(k:n_sp_y) =  999.999
526    uyma(k:n_sp_y) = -999.999
527    utext = 'y'//utext_char( pr )
528    IF ( averaging_interval_sp /= 0.0 )  THEN
529       WRITE ( atext, 104 )  averaging_interval_sp
530       utext = TRIM(utext) // ',  ' // TRIM( atext )
531    ENDIF
532    uxmin = 0.8 * 2.0 * pi        / ( dy * ( ny + 1 ) )
533    uxmax = 1.2 * 2.0 * pi * ny/2 / ( dy * ( ny + 1 ) )
534    uymin = 0.8 * MIN (  999.999, MINVAL ( uymi ) )
535    uymax = 1.2 * MAX ( -999.999, MAXVAL ( uyma ) )
536    ytext = ytext_char( pr )
537
538!
539!-- Output CROSS-parameters
540    WRITE ( 83, CROSS )
541
542!
543!-- Increase counter by the number of profiles written in the actual block
544    klist_y = klist_y + n_sp_y
545
546!
547!-- Write end-mark
548    WRITE ( 84, 103 ) 
549
550!
551!-- Close parameter- and data-file
552    CALL close_file( 83 )
553    CALL close_file( 84 )
554
555!
556!-- Formats
557100 FORMAT (A,I1,1X,A,1X,I4,'m ',A)
558101 FORMAT (A,I2,1X,A,1X,I4,'m ',A)
559102 FORMAT (E15.7,10(1X,E15.7))
560103 FORMAT ('NEXT')
561104 FORMAT ('time averaged over',F7.1,' s')
562
563 END SUBROUTINE data_output_spectra_y
564#endif
Note: See TracBrowser for help on using the repository browser.