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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

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