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

Last change on this file since 1037 was 1037, checked in by raasch, 9 years ago

last commit documented

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