source: palm/trunk/SOURCE/data_output_3d.f90 @ 375

Last change on this file since 375 was 355, checked in by letzel, 15 years ago
  • Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)
  • Property svn:keywords set to Id
File size: 12.5 KB
RevLine 
[1]1 SUBROUTINE data_output_3d( av )
2
3!------------------------------------------------------------------------------!
[254]4! Current revisions:
[1]5! -----------------
[291]6! simulated_time in NetCDF output replaced by time_since_reference_point.
[263]7! Output of NetCDF messages with aid of message handling routine.
[254]8! Output of messages replaced by message handling routine.
[355]9! Bugfix: to_be_resorted => s_av for time-averaged scalars
[1]10!
[254]11!
[1]12! Former revisions:
13! -----------------
[3]14! $Id: data_output_3d.f90 355 2009-07-17 01:03:01Z heinze $
[77]15!
[98]16! 96 2007-06-04 08:07:41Z raasch
17! Output of density and salinity
18!
[77]19! 75 2007-03-22 09:54:05Z raasch
20! 2nd+3rd argument removed from exchange horiz
21!
[3]22! RCS Log replace by Id keyword, revision history cleaned up
23!
[1]24! Revision 1.3  2006/06/02 15:18:59  raasch
25! +argument "found", -argument grid in call of routine user_data_output_3d
26!
27! Revision 1.2  2006/02/23 10:23:07  raasch
28! Former subroutine plot_3d renamed data_output_3d, pl.. renamed do..,
29! .._anz renamed .._n,
30! output extended to (almost) all quantities, output of user-defined quantities
31!
32! Revision 1.1  1997/09/03 06:29:36  raasch
33! Initial revision
34!
35!
36! Description:
37! ------------
38! Output of the 3D-arrays in NetCDF and/or AVS format.
39!------------------------------------------------------------------------------!
40
41    USE array_kind
42    USE arrays_3d
43    USE averaging
44    USE cloud_parameters
45    USE control_parameters
46    USE cpulog
47    USE indices
48    USE interfaces
49    USE netcdf_control
50    USE particle_attributes
51    USE pegrid
52
53    IMPLICIT NONE
54
55    CHARACTER (LEN=9) ::  simulated_time_mod
56
57    INTEGER           ::  av, i, if, j, k, n, pos, prec, psi
58
59    LOGICAL           ::  found, resorted
60
61    REAL              ::  mean_r, s_r3, s_r4
62
63    REAL(spk), DIMENSION(:,:,:), ALLOCATABLE  ::  local_pf
64
65    REAL, DIMENSION(:,:,:), POINTER ::  to_be_resorted
66
67!
68!-- Return, if nothing to output
69    IF ( do3d_no(av) == 0 )  RETURN
70
71    CALL cpu_log (log_point(14),'data_output_3d','start')
72
73!
74!-- Open output file.
75!-- Also creates coordinate and fld-file for AVS.
76!-- In case of a run on more than one PE, each PE opens its own file and
77!-- writes the data of its subdomain in binary format (regardless of the format
78!-- the user has requested). After the run, these files are combined to one
79!-- file by combine_plot_fields in the format requested by the user (netcdf
80!-- and/or avs).
81    IF ( avs_output  .OR.  ( numprocs > 1 ) )  CALL check_open( 30 )
82
83#if defined( __netcdf )
84    IF ( myid == 0  .AND.  netcdf_output )  CALL check_open( 106+av*10 )
85#endif
86
87!
88!-- Allocate a temporary array with the desired output dimensions.
89    ALLOCATE( local_pf(nxl-1:nxr+1,nys-1:nyn+1,nzb:nz_do3d) )
90
91!
92!-- Update the NetCDF time axis
93#if defined( __netcdf )
94    do3d_time_count(av) = do3d_time_count(av) + 1
95    IF ( myid == 0  .AND.  netcdf_output )  THEN
96       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av), &
[291]97                               (/ time_since_reference_point /),  &
[1]98                               start = (/ do3d_time_count(av) /), &
99                               count = (/ 1 /) )
[314]100       CALL handle_netcdf_error( 'data_output_3d', 376 )
[1]101    ENDIF
102#endif
103
104!
105!-- Loop over all variables to be written.
106    if = 1
107
108    DO  WHILE ( do3d(av,if)(1:1) /= ' ' )
109!
110!--    Set the precision for data compression.
111       IF ( do3d_compress )  THEN
112          DO  i = 1, 100
113             IF ( plot_3d_precision(i)%variable == do3d(av,if) )  THEN
114                prec = plot_3d_precision(i)%precision
115                EXIT
116             ENDIF
117          ENDDO
118       ENDIF
119
120!
121!--    Store the array chosen on the temporary array.
122       resorted = .FALSE.
123       SELECT CASE ( TRIM( do3d(av,if) ) )
124
125          CASE ( 'e' )
126             IF ( av == 0 )  THEN
127                to_be_resorted => e
128             ELSE
129                to_be_resorted => e_av
130             ENDIF
131
132          CASE ( 'p' )
133             IF ( av == 0 )  THEN
134                to_be_resorted => p
135             ELSE
136                to_be_resorted => p_av
137             ENDIF
138
139          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
140             IF ( av == 0 )  THEN
141                tend = prt_count
[75]142                CALL exchange_horiz( tend )
[1]143                DO  i = nxl-1, nxr+1
144                   DO  j = nys-1, nyn+1
145                      DO  k = nzb, nz_do3d
146                         local_pf(i,j,k) = tend(k,j,i)
147                      ENDDO
148                   ENDDO
149                ENDDO
150                resorted = .TRUE.
151             ELSE
[75]152                CALL exchange_horiz( pc_av )
[1]153                to_be_resorted => pc_av
154             ENDIF
155
156          CASE ( 'pr' )  ! mean particle radius
157             IF ( av == 0 )  THEN
158                DO  i = nxl, nxr
159                   DO  j = nys, nyn
160                      DO  k = nzb, nzt+1
161                         psi = prt_start_index(k,j,i)
162                         s_r3 = 0.0
163                         s_r4 = 0.0
164                         DO  n = psi, psi+prt_count(k,j,i)-1
165                            s_r3 = s_r3 + particles(n)%radius**3
166                            s_r4 = s_r4 + particles(n)%radius**4
167                         ENDDO
168                         IF ( s_r3 /= 0.0 )  THEN
169                            mean_r = s_r4 / s_r3
170                         ELSE
171                            mean_r = 0.0
172                         ENDIF
173                         tend(k,j,i) = mean_r
174                      ENDDO
175                   ENDDO
176                ENDDO
[75]177                CALL exchange_horiz( tend )
[1]178                DO  i = nxl-1, nxr+1
179                   DO  j = nys-1, nyn+1
180                      DO  k = nzb, nzt+1
181                         local_pf(i,j,k) = tend(k,j,i)
182                      ENDDO
183                   ENDDO
184                ENDDO
185                resorted = .TRUE.
186             ELSE
[75]187                CALL exchange_horiz( pr_av )
[1]188                to_be_resorted => pr_av
189             ENDIF
190
191          CASE ( 'pt' )
192             IF ( av == 0 )  THEN
193                IF ( .NOT. cloud_physics ) THEN
194                   to_be_resorted => pt
195                ELSE
196                   DO  i = nxl-1, nxr+1
197                      DO  j = nys-1, nyn+1
198                         DO  k = nzb, nz_do3d
199                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp *    &
200                                                          pt_d_t(k) * &
201                                                          ql(k,j,i)
202                         ENDDO
203                      ENDDO
204                   ENDDO
205                   resorted = .TRUE.
206                ENDIF
207             ELSE
208                to_be_resorted => pt_av
209             ENDIF
210
211          CASE ( 'q' )
212             IF ( av == 0 )  THEN
213                to_be_resorted => q
214             ELSE
215                to_be_resorted => q_av
216             ENDIF
217             
218          CASE ( 'ql' )
219             IF ( av == 0 )  THEN
220                to_be_resorted => ql
221             ELSE
222                to_be_resorted => ql_av
223             ENDIF
224
225          CASE ( 'ql_c' )
226             IF ( av == 0 )  THEN
227                to_be_resorted => ql_c
228             ELSE
229                to_be_resorted => ql_c_av
230             ENDIF
231
232          CASE ( 'ql_v' )
233             IF ( av == 0 )  THEN
234                to_be_resorted => ql_v
235             ELSE
236                to_be_resorted => ql_v_av
237             ENDIF
238
239          CASE ( 'ql_vp' )
240             IF ( av == 0 )  THEN
241                to_be_resorted => ql_vp
242             ELSE
243                to_be_resorted => ql_vp_av
244             ENDIF
245
246          CASE ( 'qv' )
247             IF ( av == 0 )  THEN
248                DO  i = nxl-1, nxr+1
249                   DO  j = nys-1, nyn+1
250                      DO  k = nzb, nz_do3d
251                         local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
252                      ENDDO
253                   ENDDO
254                ENDDO
255                resorted = .TRUE.
256             ELSE
257                to_be_resorted => qv_av
258             ENDIF
259
[96]260          CASE ( 'rho' )
261             IF ( av == 0 )  THEN
262                to_be_resorted => rho
263             ELSE
264                to_be_resorted => rho_av
265             ENDIF
266             
[1]267          CASE ( 's' )
268             IF ( av == 0 )  THEN
269                to_be_resorted => q
270             ELSE
[355]271                to_be_resorted => s_av
[1]272             ENDIF
273             
[96]274          CASE ( 'sa' )
275             IF ( av == 0 )  THEN
276                to_be_resorted => sa
277             ELSE
278                to_be_resorted => sa_av
279             ENDIF
280             
[1]281          CASE ( 'u' )
282             IF ( av == 0 )  THEN
283                to_be_resorted => u
284             ELSE
285                to_be_resorted => u_av
286             ENDIF
287
288          CASE ( 'v' )
289             IF ( av == 0 )  THEN
290                to_be_resorted => v
291             ELSE
292                to_be_resorted => v_av
293             ENDIF
294
295          CASE ( 'vpt' )
296             IF ( av == 0 )  THEN
297                to_be_resorted => vpt
298             ELSE
299                to_be_resorted => vpt_av
300             ENDIF
301
302          CASE ( 'w' )
303             IF ( av == 0 )  THEN
304                to_be_resorted => w
305             ELSE
306                to_be_resorted => w_av
307             ENDIF
308
309          CASE DEFAULT
310!
311!--          User defined quantity
312             CALL user_data_output_3d( av, do3d(av,if), found, local_pf, &
313                                       nz_do3d )
314             resorted = .TRUE.
315
[254]316             IF ( .NOT. found )  THEN
[274]317                message_string =  'no output available for: ' //   &
318                                  TRIM( do3d(av,if) )
[254]319                CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
[1]320             ENDIF
321
322       END SELECT
323
324!
325!--    Resort the array to be output, if not done above
326       IF ( .NOT. resorted )  THEN
327          DO  i = nxl-1, nxr+1
328             DO  j = nys-1, nyn+1
329                DO  k = nzb, nz_do3d
330                   local_pf(i,j,k) = to_be_resorted(k,j,i)
331                ENDDO
332             ENDDO
333          ENDDO
334       ENDIF
335
336!
337!--    Output of the volume data information for the AVS-FLD-file.
338       do3d_avs_n = do3d_avs_n + 1
339       IF ( myid == 0  .AND.  avs_output )  THEN
340!
341!--       AVS-labels must not contain any colons. Hence they must be removed
342!--       from the time character string.
343          simulated_time_mod = simulated_time_chr
344          DO  WHILE ( SCAN( simulated_time_mod, ':' ) /= 0 )
345             pos = SCAN( simulated_time_mod, ':' )
346             simulated_time_mod(pos:pos) = '/'
347          ENDDO
348
349          IF ( av == 0 )  THEN
350             WRITE ( 33, 3300 )  do3d_avs_n, TRIM( avs_data_file ), &
351                                 skip_do_avs, TRIM( do3d(av,if) ), &
352                                 TRIM( simulated_time_mod )
353          ELSE
354             WRITE ( 33, 3300 )  do3d_avs_n, TRIM( avs_data_file ), &
355                                 skip_do_avs, TRIM( do3d(av,if) ) // &
356                                 ' averaged', TRIM( simulated_time_mod )
357          ENDIF
358!
359!--       Determine the Skip-value for the next array. Record end and start
360!--       require 4 byte each.
361          skip_do_avs = skip_do_avs + ( ((nx+2)*(ny+2)*(nz_do3d+1)) * 4 + 8 )
362       ENDIF
363
364!
365!--    Output of the 3D-array. (compressed/uncompressed)
366       IF ( do3d_compress )  THEN
367!
368!--       Compression, output of compression information on FLD-file and output
369!--       of compressed data.
370          CALL write_compressed( local_pf, 30, 33, myid, nxl, nxr, nyn, nys, &
371                                 nzb, nz_do3d, prec )
372       ELSE
373!
374!--       Uncompressed output.
375#if defined( __parallel )
376          IF ( netcdf_output  .AND.  myid == 0 )  THEN
377             WRITE ( 30 )  simulated_time, do3d_time_count(av), av
378          ENDIF
379          WRITE ( 30 )  nxl-1, nxr+1, nys-1, nyn+1, nzb, nz_do3d
380          WRITE ( 30 )  local_pf
381#else
382          IF ( avs_output )  THEN
383             WRITE ( 30 )  local_pf(nxl:nxr+1,nys:nyn+1,:)
384          ENDIF
385#if defined( __netcdf )
386          IF ( netcdf_output )  THEN
387             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),       &
388                                  local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),  &
389                                  start = (/ 1, 1, 1, do3d_time_count(av) /), &
390                                  count = (/ nx+2, ny+2, nz_do3d-nzb+1, 1 /) )
[314]391             CALL handle_netcdf_error( 'data_output_3d', 386 )
[1]392          ENDIF
393#endif
394#endif
395       ENDIF
396
397       if = if + 1
398
399    ENDDO
400
401!
402!-- Deallocate temporary array.
403    DEALLOCATE( local_pf )
404
405
406    CALL cpu_log (log_point(14),'data_output_3d','stop','nobarrier')
407
408!
409!-- Formats.
4103300 FORMAT ('variable ',I4,'  file=',A,'  filetype=unformatted  skip=',I12/ &
411             'label = ',A,A)
412
413 END SUBROUTINE data_output_3d
Note: See TracBrowser for help on using the repository browser.