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

Last change on this file since 1311 was 1310, checked in by raasch, 10 years ago

update of GPL copyright

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