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

Last change on this file since 1036 was 1036, checked in by raasch, 11 years ago

code has been put under the GNU General Public License (v3)

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