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

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

changes for coupling with independent precursor runs; z_i calculation with Sullivan criterion

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