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

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