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

Last change on this file since 790 was 790, checked in by raasch, 12 years ago

Bugfix for output of mean particle radius + preliminary works for implementing the Wang collision kernel

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