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

Last change on this file since 691 was 691, checked in by maronga, 10 years ago

Bugfix for precursor atmosphere/ocean runs, re-adjustments for lcxt4

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