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

Last change on this file since 1026 was 965, checked in by raasch, 12 years ago

last commit documented

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