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

Last change on this file since 1319 was 1319, checked in by raasch, 10 years ago

last commit documented

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