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

Last change on this file since 1026 was 1008, checked in by franke, 12 years ago

last commit documented

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