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

Last change on this file since 1267 was 1245, checked in by raasch, 11 years ago

last commit documented

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