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

Last change on this file since 1035 was 1035, checked in by raasch, 9 years ago

revisions r1031 and r1034 documented

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