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

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

Bugfixes


Missing calculation of mean particle weighting factor for output added. (data_output_2d, data_output_3d, data_output_mask, sum_up_3d_data)
Calculation of mean particle radius for output now considers the weighting factor. (data_output_mask)
Calculation of sugrid-scale buoyancy flux for humidity and cloud droplets corrected. (flow_statistics)
Factor in calculation of enhancement factor for collision efficencies corrected. (lpm_collision_kernels)
Calculation of buoyancy production now considers the liquid water mixing ratio in case of cloud droplets. (production_e)

Changes


Calculation of buoyancy flux for humidity in case of WS-scheme is now using turbulent fluxes of WS-scheme. (flow_statistics)
Calculation of the collision kernels now in SI units. (lpm_collision_kernels)

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