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

Last change on this file since 776 was 772, checked in by heinze, 13 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 16.5 KB
Line 
1 SUBROUTINE data_output_3d( av )
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7! Former revisions:
8! -----------------
9! $Id: data_output_3d.f90 772 2011-10-27 11:08:46Z heinze $
10!
11! 771 2011-10-27 10:56:21Z heinze
12! +lpt
13!
14! 759 2011-09-15 13:58:31Z raasch
15! Splitting of parallel I/O
16!
17! 727 2011-04-20 20:05:25Z suehring
18! Exchange ghost layers also for p_av.
19!
20! 725 2011-04-11 09:37:01Z suehring
21! Exchange ghost layers for p regardless of used pressure solver (except SOR).
22!
23! 691 2011-03-04 08:45:30Z maronga
24! Replaced simulated_time by time_since_reference_point
25!
26! 673 2011-01-18 16:19:48Z suehring
27! When using Multigrid or SOR solver an additional CALL exchange_horiz is
28! is needed for pressure output.
29!
30! 667 2010-12-23 12:06:00Z suehring/gryschka
31! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
32! allocation of arrays.  Calls of exchange_horiz are modified.
33! Skip-value skip_do_avs changed to a dynamic adaption of ghost points.
34!
35! 646 2010-12-15 13:03:52Z raasch
36! bugfix: missing define statements for netcdf added
37!
38! 493 2010-03-01 08:30:24Z raasch
39! NetCDF4 support (parallel output)
40!
41! 355 2009-07-17 01:03:01Z letzel
42! simulated_time in NetCDF output replaced by time_since_reference_point.
43! Output of NetCDF messages with aid of message handling routine.
44! Output of messages replaced by message handling routine.
45! Bugfix: to_be_resorted => s_av for time-averaged scalars
46!
47! 96 2007-06-04 08:07:41Z raasch
48! Output of density and salinity
49!
50! 75 2007-03-22 09:54:05Z raasch
51! 2nd+3rd argument removed from exchange horiz
52!
53! RCS Log replace by Id keyword, revision history cleaned up
54!
55! Revision 1.3  2006/06/02 15:18:59  raasch
56! +argument "found", -argument grid in call of routine user_data_output_3d
57!
58! Revision 1.2  2006/02/23 10:23:07  raasch
59! Former subroutine plot_3d renamed data_output_3d, pl.. renamed do..,
60! .._anz renamed .._n,
61! output extended to (almost) all quantities, output of user-defined quantities
62!
63! Revision 1.1  1997/09/03 06:29:36  raasch
64! Initial revision
65!
66!
67! Description:
68! ------------
69! Output of the 3D-arrays in NetCDF and/or AVS format.
70!------------------------------------------------------------------------------!
71
72    USE array_kind
73    USE arrays_3d
74    USE averaging
75    USE cloud_parameters
76    USE control_parameters
77    USE cpulog
78    USE indices
79    USE interfaces
80    USE netcdf_control
81    USE particle_attributes
82    USE pegrid
83
84    IMPLICIT NONE
85
86    CHARACTER (LEN=9) ::  simulated_time_mod
87
88    INTEGER           ::  av, i, if, j, k, n, pos, prec, psi
89
90    LOGICAL           ::  found, resorted
91
92    REAL              ::  mean_r, s_r3, s_r4
93
94    REAL(spk), DIMENSION(:,:,:), ALLOCATABLE  ::  local_pf
95
96    REAL, DIMENSION(:,:,:), POINTER ::  to_be_resorted
97
98!
99!-- Return, if nothing to output
100    IF ( do3d_no(av) == 0 )  RETURN
101
102    CALL cpu_log (log_point(14),'data_output_3d','start')
103
104!
105!-- Open output file.
106!-- Also creates coordinate and fld-file for AVS.
107!-- For classic or 64bit NetCDF output or output of other (old) data formats,
108!-- for a run on more than one PE, each PE opens its own file and
109!-- writes the data of its subdomain in binary format (regardless of the format
110!-- the user has requested). After the run, these files are combined to one
111!-- file by combine_plot_fields in the format requested by the user (netcdf
112!-- and/or avs).
113!-- For NetCDF4/HDF5 output, data is written in parallel into one file.
114    IF ( netcdf_output )  THEN
115       IF ( netcdf_data_format < 3 )  THEN
116          CALL check_open( 30 )
117          IF ( myid == 0 )  CALL check_open( 106+av*10 )
118       ELSE
119          CALL check_open( 106+av*10 )
120       ENDIF
121    ELSE
122       IF ( avs_output  .OR.  ( numprocs > 1 ) )  CALL check_open( 30 )
123    ENDIF
124
125!
126!-- Allocate a temporary array with the desired output dimensions.
127    ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) )
128
129!
130!-- Update the NetCDF time axis
131#if defined( __netcdf )
132    IF ( myid == 0  .OR.  netcdf_data_format > 2 )  THEN
133       do3d_time_count(av) = do3d_time_count(av) + 1
134       IF ( netcdf_output )  THEN
135          nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av), &
136                                  (/ time_since_reference_point /),  &
137                                  start = (/ do3d_time_count(av) /), &
138                                  count = (/ 1 /) )
139          CALL handle_netcdf_error( 'data_output_3d', 376 )
140       ENDIF
141    ENDIF
142#endif
143
144!
145!-- Loop over all variables to be written.
146    if = 1
147
148    DO  WHILE ( do3d(av,if)(1:1) /= ' ' )
149!
150!--    Set the precision for data compression.
151       IF ( do3d_compress )  THEN
152          DO  i = 1, 100
153             IF ( plot_3d_precision(i)%variable == do3d(av,if) )  THEN
154                prec = plot_3d_precision(i)%precision
155                EXIT
156             ENDIF
157          ENDDO
158       ENDIF
159
160!
161!--    Store the array chosen on the temporary array.
162       resorted = .FALSE.
163       SELECT CASE ( TRIM( do3d(av,if) ) )
164
165          CASE ( 'e' )
166             IF ( av == 0 )  THEN
167                to_be_resorted => e
168             ELSE
169                to_be_resorted => e_av
170             ENDIF
171
172          CASE ( 'lpt' )
173             IF ( av == 0 )  THEN
174                to_be_resorted => pt
175             ELSE
176                to_be_resorted => lpt_av
177             ENDIF
178
179          CASE ( 'p' )
180             IF ( av == 0 )  THEN
181                IF ( psolver /= 'sor' )  CALL exchange_horiz( p, nbgp )
182                to_be_resorted => p
183             ELSE
184                IF ( psolver /= 'sor' )  CALL exchange_horiz( p_av, nbgp )
185                to_be_resorted => p_av
186             ENDIF
187
188          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
189             IF ( av == 0 )  THEN
190                tend = prt_count
191                CALL exchange_horiz( tend, nbgp )
192                DO  i = nxlg, nxrg
193                   DO  j = nysg, nyng
194                      DO  k = nzb, nz_do3d
195                         local_pf(i,j,k) = tend(k,j,i)
196                      ENDDO
197                   ENDDO
198                ENDDO
199                resorted = .TRUE.
200             ELSE
201                CALL exchange_horiz( pc_av, nbgp )
202                to_be_resorted => pc_av
203             ENDIF
204
205          CASE ( 'pr' )  ! mean particle radius
206             IF ( av == 0 )  THEN
207                DO  i = nxl, nxr
208                   DO  j = nys, nyn
209                      DO  k = nzb, nzt+1
210                         psi = prt_start_index(k,j,i)
211                         s_r3 = 0.0
212                         s_r4 = 0.0
213                         DO  n = psi, psi+prt_count(k,j,i)-1
214                            s_r3 = s_r3 + particles(n)%radius**3
215                            s_r4 = s_r4 + particles(n)%radius**4
216                         ENDDO
217                         IF ( s_r3 /= 0.0 )  THEN
218                            mean_r = s_r4 / s_r3
219                         ELSE
220                            mean_r = 0.0
221                         ENDIF
222                         tend(k,j,i) = mean_r
223                      ENDDO
224                   ENDDO
225                ENDDO
226                CALL exchange_horiz( tend, nbgp )
227                DO  i = nxlg, nxrg
228                   DO  j = nysg, nyng
229                      DO  k = nzb, nzt+1
230                         local_pf(i,j,k) = tend(k,j,i)
231                      ENDDO
232                   ENDDO
233                ENDDO
234                resorted = .TRUE.
235             ELSE
236                CALL exchange_horiz( pr_av, nbgp )
237                to_be_resorted => pr_av
238             ENDIF
239
240          CASE ( 'pt' )
241             IF ( av == 0 )  THEN
242                IF ( .NOT. cloud_physics ) THEN
243                   to_be_resorted => pt
244                ELSE
245                   DO  i = nxlg, nxrg
246                      DO  j = nysg, nyng
247                         DO  k = nzb, nz_do3d
248                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp *    &
249                                                          pt_d_t(k) * &
250                                                          ql(k,j,i)
251                         ENDDO
252                      ENDDO
253                   ENDDO
254                   resorted = .TRUE.
255                ENDIF
256             ELSE
257                to_be_resorted => pt_av
258             ENDIF
259
260          CASE ( 'q' )
261             IF ( av == 0 )  THEN
262                to_be_resorted => q
263             ELSE
264                to_be_resorted => q_av
265             ENDIF
266
267          CASE ( 'ql' )
268             IF ( av == 0 )  THEN
269                to_be_resorted => ql
270             ELSE
271                to_be_resorted => ql_av
272             ENDIF
273
274          CASE ( 'ql_c' )
275             IF ( av == 0 )  THEN
276                to_be_resorted => ql_c
277             ELSE
278                to_be_resorted => ql_c_av
279             ENDIF
280
281          CASE ( 'ql_v' )
282             IF ( av == 0 )  THEN
283                to_be_resorted => ql_v
284             ELSE
285                to_be_resorted => ql_v_av
286             ENDIF
287
288          CASE ( 'ql_vp' )
289             IF ( av == 0 )  THEN
290                to_be_resorted => ql_vp
291             ELSE
292                to_be_resorted => ql_vp_av
293             ENDIF
294
295          CASE ( 'qv' )
296             IF ( av == 0 )  THEN
297                DO  i = nxlg, nxrg
298                   DO  j = nysg, nyng
299                      DO  k = nzb, nz_do3d
300                         local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
301                      ENDDO
302                   ENDDO
303                ENDDO
304                resorted = .TRUE.
305             ELSE
306                to_be_resorted => qv_av
307             ENDIF
308
309          CASE ( 'rho' )
310             IF ( av == 0 )  THEN
311                to_be_resorted => rho
312             ELSE
313                to_be_resorted => rho_av
314             ENDIF
315
316          CASE ( 's' )
317             IF ( av == 0 )  THEN
318                to_be_resorted => q
319             ELSE
320                to_be_resorted => s_av
321             ENDIF
322
323          CASE ( 'sa' )
324             IF ( av == 0 )  THEN
325                to_be_resorted => sa
326             ELSE
327                to_be_resorted => sa_av
328             ENDIF
329
330          CASE ( 'u' )
331             IF ( av == 0 )  THEN
332                to_be_resorted => u
333             ELSE
334                to_be_resorted => u_av
335             ENDIF
336
337          CASE ( 'v' )
338             IF ( av == 0 )  THEN
339                to_be_resorted => v
340             ELSE
341                to_be_resorted => v_av
342             ENDIF
343
344          CASE ( 'vpt' )
345             IF ( av == 0 )  THEN
346                to_be_resorted => vpt
347             ELSE
348                to_be_resorted => vpt_av
349             ENDIF
350
351          CASE ( 'w' )
352             IF ( av == 0 )  THEN
353                to_be_resorted => w
354             ELSE
355                to_be_resorted => w_av
356             ENDIF
357
358          CASE DEFAULT
359!
360!--          User defined quantity
361             CALL user_data_output_3d( av, do3d(av,if), found, local_pf, &
362                                       nz_do3d )
363             resorted = .TRUE.
364
365             IF ( .NOT. found )  THEN
366                message_string =  'no output available for: ' //   &
367                                  TRIM( do3d(av,if) )
368                CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
369             ENDIF
370
371       END SELECT
372
373!
374!--    Resort the array to be output, if not done above
375       IF ( .NOT. resorted )  THEN
376          DO  i = nxlg, nxrg
377             DO  j = nysg, nyng
378                DO  k = nzb, nz_do3d
379                   local_pf(i,j,k) = to_be_resorted(k,j,i)
380                ENDDO
381             ENDDO
382          ENDDO
383       ENDIF
384
385!
386!--    Output of the volume data information for the AVS-FLD-file.
387       do3d_avs_n = do3d_avs_n + 1
388       IF ( myid == 0  .AND.  avs_output )  THEN
389!
390!--       AVS-labels must not contain any colons. Hence they must be removed
391!--       from the time character string.
392          simulated_time_mod = simulated_time_chr
393          DO  WHILE ( SCAN( simulated_time_mod, ':' ) /= 0 )
394             pos = SCAN( simulated_time_mod, ':' )
395             simulated_time_mod(pos:pos) = '/'
396          ENDDO
397
398          IF ( av == 0 )  THEN
399             WRITE ( 33, 3300 )  do3d_avs_n, TRIM( avs_data_file ), &
400                                 skip_do_avs, TRIM( do3d(av,if) ), &
401                                 TRIM( simulated_time_mod )
402          ELSE
403             WRITE ( 33, 3300 )  do3d_avs_n, TRIM( avs_data_file ), &
404                                 skip_do_avs, TRIM( do3d(av,if) ) // &
405                                 ' averaged', TRIM( simulated_time_mod )
406          ENDIF
407!
408!--       Determine the Skip-value for the next array. Record end and start
409!--       require 4 byte each.
410          skip_do_avs = skip_do_avs + ( ( ( nx+2*nbgp ) * ( ny+2*nbgp ) * &
411                                          ( nz_do3d+1 ) ) * 4 + 8 )
412       ENDIF
413
414!
415!--    Output of the 3D-array. (compressed/uncompressed)
416       IF ( do3d_compress )  THEN
417!
418!--       Compression, output of compression information on FLD-file and output
419!--       of compressed data.
420          CALL write_compressed( local_pf, 30, 33, myid, nxl, nxr, nyn, nys, &
421                                 nzb, nz_do3d, prec, nbgp )
422       ELSE
423!
424!--       Uncompressed output.
425#if defined( __parallel )
426          IF ( netcdf_output )  THEN
427             IF ( netcdf_data_format < 3 )  THEN
428!
429!--             Classic or 64bit format. Data is output in parallel in FORTRAN
430!--             binary format here, and later collected into one file by
431!--             combine_plot_fields
432                IF ( myid == 0 )  THEN
433                   WRITE ( 30 )  time_since_reference_point,                   &
434                                 do3d_time_count(av), av
435                ENDIF
436                DO  i = 0, io_blocks-1
437                   IF ( i == io_group )  THEN
438                      WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb, nz_do3d
439                      WRITE ( 30 )  local_pf
440                   ENDIF
441#if defined( __parallel )
442                   CALL MPI_BARRIER( comm2d, ierr )
443#endif
444                ENDDO
445
446             ELSE
447#if defined( __netcdf )
448!
449!--             Output in NetCDF4/HDF5 format.
450!--             Do not output redundant ghost point data except for the
451!--             boundaries of the total domain.
452                IF ( nxr == nx  .AND.  nyn /= ny )  THEN
453                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
454                                  local_pf(nxl:nxrg,nys:nyn,nzb:nz_do3d),    &
455                      start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
456                      count = (/ nxr-nxl+1+nbgp, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
457                ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
458                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
459                                  local_pf(nxl:nxr,nys:nyng,nzb:nz_do3d),    &
460                      start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
461                      count = (/ nxr-nxl+1, nyn-nys+1+nbgp, nz_do3d-nzb+1, 1 /) )
462                ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
463                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
464                                  local_pf(nxl:nxrg,nys:nyng,nzb:nz_do3d),  &
465                      start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
466                      count = (/ nxr-nxl+1+nbgp, nyn-nys+1+nbgp, nz_do3d-nzb+1, 1 /) )
467                ELSE
468                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
469                                  local_pf(nxl:nxr,nys:nyn,nzb:nz_do3d),      &
470                      start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
471                      count = (/ nxr-nxl+1, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
472                ENDIF
473                CALL handle_netcdf_error( 'data_output_3d', 386 )
474#endif
475             ENDIF
476          ENDIF
477#else
478          IF ( avs_output )  THEN
479             WRITE ( 30 )  local_pf(nxl:nxr+1,nys:nyn+1,:)
480          ENDIF
481#if defined( __netcdf )
482          IF ( netcdf_output )  THEN
483
484             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),    &
485                               local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),  &
486                               start = (/ 1, 1, 1, do3d_time_count(av) /), &
487                               count = (/ nx+2, ny+2, nz_do3d-nzb+1, 1 /) )
488             CALL handle_netcdf_error( 'data_output_3d', 446 )
489
490          ENDIF
491#endif
492#endif
493       ENDIF
494
495       if = if + 1
496
497    ENDDO
498
499!
500!-- Deallocate temporary array.
501    DEALLOCATE( local_pf )
502
503
504    CALL cpu_log (log_point(14),'data_output_3d','stop','nobarrier')
505
506!
507!-- Formats.
5083300 FORMAT ('variable ',I4,'  file=',A,'  filetype=unformatted  skip=',I12/ &
509             'label = ',A,A)
510
511 END SUBROUTINE data_output_3d
Note: See TracBrowser for help on using the repository browser.