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

Last change on this file since 203 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

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