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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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