source: palm/trunk/UTIL/combine_plot_fields.f90 @ 3222

Last change on this file since 3222 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

  • Property svn:keywords set to Id
File size: 30.1 KB
Line 
1 PROGRAM combine_plot_fields
2
3!------------------------------------------------------------------------------!
4! This file is part of the PALM model system.
5!
6! PALM is free software: you can redistribute it and/or modify it under the
7! terms of the GNU General Public License as published by the Free Software
8! Foundation, either version 3 of the License, or (at your option) any later
9! version.
10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
18! Copyright 1997-2018  Leibniz Universitaet Hannover
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: combine_plot_fields.f90 2718 2018-01-02 08:49:38Z suehring $
28! Corrected "Former revisions" section
29!
30! 2696 2017-12-14 17:12:51Z kanani
31! Change in file header (GPL part)
32!
33! 2669 2017-12-06 16:03:27Z raasch
34! data of 3d-nest runs are completely processed now
35!
36! 2523 2017-10-05 14:42:47Z kanani
37! Increased LEN for CHARACTER variable var_name, equal to the value in PALM
38!
39! 2512 2017-10-04 08:26:59Z raasch
40! PALM output does not contain ghost layer data any more
41! avs- and iso2d-related parts removed, handling of compressed data removed
42!
43! 2365 2017-08-21 14:59:59Z kanani
44! Vertical grid nesting implemented (SadiqHuq)
45!
46! 1809 2016-04-05 20:13:28Z raasch
47!
48! 1808 2016-04-05 19:44:00Z raasch
49! cpu measurements are done with standard FORTRAN routine on every machine
50!
51! 1551 2015-03-03 14:18:16Z maronga
52! Adjustments for data output of soil model quantities
53!
54! 1468 2014-09-24 14:06:57Z maronga
55! Adapted for use on up to 6-digit processor cores (not tested)
56!
57! 1394 2014-05-06 10:17:31Z keck
58! KIND-attribute added to all INTEGER and REAL declaration statements
59!
60! 1046 2012-11-09 14:38:45Z maronga
61! code put under GPL (PALM 3.9)
62!
63! 691 2011-03-04 08:45:30Z maronga
64! Bugfix: check for precursor ocean runs, removed typo
65!
66! 493 2010-03-01 08:30:24Z raasch
67! Exit in case of already complete NetCDF data (due to parallel output in PALM)
68! cpu measurements included
69!
70! 210 2008-11-06 08:54:02Z raasch
71! Size of pf3d adjusted to the required output size (1 gridpoint less, along
72! all three dimensions), because output of a subset of the data
73! (pf3d(nxa:nxe...) in the NF90_PUT_VAR statement caused segmentation fault
74! with the INTEL compiler.
75! Subdomain data are read into temporary arrays pf_tmp/pf3d_tmp in order to
76! avoid INTEL compiler warnings about (automatic) creation of temporary arrays
77! Bugfix: three misplaced #endif directives
78!
79! 114 2007-10-10 00:03:15Z raasch
80! Bugfix: model_string needed a default value
81!
82! Aug 07    Loop for processing of output by coupled runs, id_string does not
83!           contain modus any longer
84!
85! 18/01/06  Output of time-averaged data
86!
87! 25/05/05  Errors removed
88!
89! 26/04/05  Output in NetCDF format, iso2d and avs output only if parameter
90!           file exists
91!
92! 31/10/01  All comments and messages translated into English
93!
94! 23/02/99  Keine Bearbeitung komprimierter 3D-Daten
95! Ursprungsversion vom 28/07/97
96!
97!
98! Description:
99! ------------
100! This routine combines data of the PALM-subdomains into one file. In PALM
101! every processor element opens its own file and writes 2D- or 3D-binary-data
102! into it (different files are opened for xy-, xz-, yz-cross-sections and
103! 3D-data). For plotting or analyzing these PE-data have to be collected and
104! to be put into single files, which is done by this routine.
105! Output format is NetCDF. Additionally, a data are output in a binary format
106! readable by ISO2D-software (cross-sections) and by AVS (3D-data).
107!------------------------------------------------------------------------------!
108
109#if defined( __netcdf )
110    USE netcdf
111#endif
112
113    IMPLICIT NONE
114
115!
116!-- Local variables
117    CHARACTER (LEN=2)    ::  modus
118    CHARACTER (LEN=4)    ::  model_string
119    CHARACTER (LEN=6)    ::  id_string
120    CHARACTER (LEN=30)   ::  dimname, var_name
121    CHARACTER (LEN=40)   ::  filename
122
123    CHARACTER (LEN=2000), DIMENSION(0:1) ::  var_list
124
125    INTEGER, PARAMETER ::  iwp = 4, spk = SELECTED_REAL_KIND( 6 ), wp = 8
126
127    INTEGER(iwp) ::  av, danz, i, id, j, k, model, models, nc_stat,            &
128                     nxa, nxag, nxe, nxeg, nya, nyag, nye, nyeg,               &
129                     nza, nze, pos, time_step, xa, xe, xxa, xxe, ya, ya_do,    &
130                     ya_tot, ye, ye_do, ye_tot, yya, yye, za, ze, zza, zze
131
132    INTEGER(8)                        ::  count, count_rate
133
134    INTEGER(iwp), DIMENSION(0:1)      ::  current_level, current_var,          &
135                                          fanz, id_set, id_var_time, num_var
136    INTEGER(iwp), DIMENSION(4)        ::  id_dims_loc
137    INTEGER(iwp), DIMENSION(0:1,4)    ::  id_dims
138    INTEGER(iwp), DIMENSION(0:1,1000) ::  id_var, levels
139
140    LOGICAL  ::  found, nest3d, netcdf_output, netcdf_parallel, netcdf_0,      &
141                 netcdf_1, vnest
142
143    REAL(wp) ::  cpu_start_time, cpu_end_time, dx, simulated_time
144    REAL(wp),  DIMENSION(:,:), ALLOCATABLE   ::  pf, pf_tmp
145    REAL(spk), DIMENSION(:,:,:), ALLOCATABLE ::  pf3d, pf3d_tmp
146
147
148
149    PRINT*, ''
150    PRINT*, ''
151    PRINT*, '*** combine_plot_fields ***'
152
153!
154!-- Find out if a coupled or nested run has been carried out
155    INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
156    INQUIRE( FILE='VNESTING_PORT_OPENED', EXIST=vnest )
157    INQUIRE( FILE='3DNESTING', EXIST=nest3d )
158    IF ( found )  THEN
159       models = 2
160       PRINT*, '    coupled run'
161    ELSEIF ( vnest )  THEN
162       models = 2
163       PRINT*, '    Vertically nested grid coupling'
164    ELSEIF ( nest3d )  THEN
165       OPEN( 90, FILE='3DNESTING', FORM='FORMATTED' )
166       READ ( 90, '(I2)' )  models
167       CLOSE ( 90 )
168       PRINT*, '    3d-nest run'
169       PRINT*, '    number of nest domains = ', models
170    ELSE
171       models = 1
172       PRINT*, '    uncoupled run'
173    ENDIF
174
175!
176!-- Find out if a precursor ocean run has been carried out
177    INQUIRE( FILE='PRECURSOR_OCEAN', EXIST=found )
178    IF ( found )  THEN
179       model_string = '_O'
180       PRINT*, '    precursor ocean run'
181    ELSE
182       model_string = ''
183    ENDIF
184
185!
186!-- Do everything for each model
187    DO model = 1, models
188!
189!--    Set the model string used to identify the filenames
190       IF ( found  .OR.  vnest )  THEN
191          PRINT*, ''
192          IF ( model == 2 )  THEN
193             IF ( vnest )  THEN
194                model_string = '_N'
195                PRINT*, '    now combining FINE data'
196                PRINT*, '    ========================'
197             ELSE
198                model_string = '_O'
199                PRINT*, '    now combining ocean data'
200                PRINT*, '    ========================'
201             ENDIF
202          ELSE
203             IF ( vnest )  THEN
204                PRINT*, '    now combining COARSE data'
205                PRINT*, '    ============================='
206             ELSE
207                PRINT*, '    now combining atmosphere data'
208                PRINT*, '    ============================='
209             ENDIF
210          ENDIF
211       ELSEIF ( nest3d )  THEN
212          PRINT*, ''
213          PRINT*, '--> processing nest id = ', model
214          IF ( model == 1 )  THEN
215             model_string = ''
216          ELSE
217             WRITE( model_string, '(A2,I2.2)' )  '_N', model
218          ENDIF
219       ENDIF
220!
221!--    2D-arrays for ISO2D
222!--    Main loop for the three different cross-sections, starting with
223!--    xy-section
224       modus = 'XY'
225       PRINT*, ''
226       DO  WHILE ( modus == 'XY'  .OR.  modus == 'XZ'  .OR.  modus == 'YZ' )
227!
228!--       Take current time
229          CALL SYSTEM_CLOCK( count, count_rate )
230          cpu_start_time = REAL( count ) / REAL( count_rate )
231
232          netcdf_parallel = .FALSE.
233!
234!--       Check, if file from PE0 exists. If it does not exist, PALM did not
235!--       create any output for this cross-section.
236          danz = 0
237          WRITE (id_string,'(I6.6)')  danz
238          INQUIRE ( &
239               FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
240               EXIST=found )
241!
242!--       Find out the number of files (equal to the number of PEs which
243!--       have been used in PALM) and open them
244          DO  WHILE ( found )
245
246             OPEN ( danz+110, &
247                  FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
248                  FORM='UNFORMATTED' )
249             danz = danz + 1
250             WRITE (id_string,'(I6.6)')  danz
251             INQUIRE ( &
252                  FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
253                  EXIST=found )
254
255          ENDDO
256
257!
258!--       Inquire whether a NetCDF file exists
259          INQUIRE( FILE='DATA_2D_'//modus//'_NETCDF'//TRIM( model_string ), &
260               EXIST=netcdf_0 )
261
262!
263!--       Inquire whether a NetCDF file for time-averaged data exists
264          INQUIRE( FILE='DATA_2D_'//modus//'_AV_NETCDF'//TRIM( model_string ),&
265               EXIST=netcdf_1 )
266
267          IF ( netcdf_0  .OR.  netcdf_1 )  THEN
268             netcdf_output = .TRUE.
269!
270!--          Inquire whether the NetCDF file is already complete (parallel
271!--          output)
272             INQUIRE( FILE='NO_COMBINE_PLOT_FIELDS_'//modus, &
273                      EXIST=netcdf_parallel )
274             IF ( netcdf_parallel )  THEN
275                netcdf_parallel = .TRUE.
276             ELSE
277                netcdf_parallel = .FALSE.
278             ENDIF
279          ELSE
280             netcdf_output = .FALSE.
281          ENDIF
282
283!
284!--       Info-output
285          PRINT*, ''
286#if defined( __netcdf )
287          IF ( netcdf_output )  THEN
288             IF ( netcdf_parallel )  THEN
289             PRINT*, '    NetCDF ' // modus // '-data are in one file ', &
290                          '(NetCDF4-format) - merging not neccessary'
291             ELSE
292                PRINT*, '    NetCDF output enabled'
293             ENDIF
294          ENDIF
295#else
296          IF ( netcdf_output )  THEN
297             PRINT*, '--- Sorry, no NetCDF support on this host'
298             netcdf_output = .FALSE.
299          ENDIF
300#endif
301          IF ( .NOT. netcdf_parallel )  THEN
302             IF ( danz /= 0 )  THEN
303                PRINT*, '    ',modus,'-section:  ', danz, ' file(s) found'
304             ELSE
305                PRINT*, '    no ', modus, '-section data available'
306             ENDIF
307          ENDIF
308
309          IF ( netcdf_output  .AND. .NOT. netcdf_parallel  .AND.  danz /= 0 ) &
310          THEN
311#if defined( __netcdf )
312             DO  av = 0, 1
313
314                IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
315                IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
316
317!
318!--             Open NetCDF dataset
319                IF ( av == 0 )  THEN
320                   filename = 'DATA_2D_'//modus//'_NETCDF' &
321                        //TRIM( model_string )
322                ELSE
323                   filename = 'DATA_2D_'//modus//'_AV_NETCDF' &
324                        //TRIM( model_string )
325                ENDIF
326                nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
327                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 1 )
328
329!
330!--             Get the list of variables (order of variables corresponds with
331!--             the order of data on the binary file)
332                var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
333                nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
334                     var_list(av) )
335                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 2 )
336
337!
338!--             Inquire id of the time coordinate variable
339                nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
340                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 3 )
341
342!
343!--             Count number of variables; there is one more semicolon in the
344!--             string than variable names
345                num_var(av) = -1
346                DO  i = 1, LEN( var_list(av) )
347                   IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) +1
348                ENDDO
349
350!
351!--             Extract the variable names from the list and inquire their
352!--             NetCDF IDs
353                pos = INDEX( var_list(av), ';' )
354!
355!--             Loop over all variables
356                DO  i = 1, num_var(av)
357
358!
359!--                Extract variable name from list
360                   var_list(av) = var_list(av)(pos+1:)
361                   pos = INDEX( var_list(av), ';' )
362                   var_name = var_list(av)(1:pos-1)
363
364!
365!--                Get variable ID from name
366                   nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
367                        id_var(av,i) )
368                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 4 )
369
370!
371!--                Get number of x/y/z levels for that variable
372                   nc_stat = NF90_INQUIRE_VARIABLE( id_set(av), id_var(av,i), &
373                        dimids = id_dims_loc )
374                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 5 )
375                   id_dims(av,:) = id_dims_loc
376
377!
378!--                Inquire dimension ID
379                   DO  j = 1, 4
380                      nc_stat = NF90_INQUIRE_DIMENSION( id_set(av), &
381                           id_dims(av,j), dimname, levels(av,i) )
382                      IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 6 )
383
384                      IF ( modus == 'XY' .AND. INDEX(dimname, 'z') /= 0 )  EXIT
385                      IF ( modus == 'XZ' .AND. INDEX(dimname, 'y') /= 0 )  EXIT
386                      IF ( modus == 'YZ' .AND. INDEX(dimname, 'x') /= 0 )  EXIT
387                   ENDDO
388
389                ENDDO
390
391             ENDDO   ! av = 0, 1
392
393#endif
394          ENDIF
395
396!
397!--       Read the arrays, as long as the end of the file is reached
398          IF ( .NOT. netcdf_parallel )  THEN
399
400             fanz          =         0
401             current_level =         1
402             current_var   = 999999999
403
404             DO  WHILE ( danz /= 0 )
405
406!
407!--             Loop over all files (reading data of the subdomains)
408                DO  id = 0, danz-1
409!
410!--                File from PE0 contains special information at the beginning,
411!--                concerning the lower and upper indices of the total-domain
412!--                used in PALM (nxa, nxe, nya, nye).
413!--                Allocate necessary arrays, open the output file and write
414!--                the coordinate informations needed by ISO2D.
415                   IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 ) THEN
416
417                      READ ( id+110 )  nxa, nxe, nya, nye
418                      ALLOCATE ( pf(nxa:nxe,nya:nye) )
419!
420!--                   Set actual domain bounds to total domain
421                      ya_do = nya
422                      ye_do = nye
423
424                   ENDIF
425!
426!--                Read output time
427                   IF ( netcdf_output  .AND.  id == 0 )  THEN
428                      READ ( id+110, END=998 )  simulated_time, time_step, av
429                   ENDIF
430!
431!--                Read subdomain indices
432                   READ ( id+110, END=998 )  xa, xe, ya, ye, ya_tot, ye_tot
433
434!
435!--                IF the PE made no output (in case that no part of the
436!--                cross-section is situated on this PE), indices have the
437!--                value -1
438                   IF ( .NOT. ( xa == -1  .AND.  xe == -1  .AND. &
439                                ya == -1  .AND.  ye == -1 ) )  THEN
440
441
442!
443!--                   Read the subdomain grid-point values
444                      ALLOCATE( pf_tmp(xa:xe,ya:ye) )
445                      READ ( id+110 )  pf_tmp
446!
447!--                   Calculate indices on atmospheric grid (e.g. for soil model
448!--                   quantities)
449                      IF ( ya /= ya_tot .OR. ye /= ye_tot )  THEN
450                         ye_do = ye - ya
451                         ya_do = ya
452                         pf(xa:xe,ya_do:ye_do) = pf_tmp
453                      ELSE
454                         ye_do = nye
455                         ya_do = nya
456                         pf(xa:xe,ya:ye) = pf_tmp
457                      ENDIF
458
459                      DEALLOCATE( pf_tmp )
460                   ENDIF
461                   IF ( id == 0 )  fanz(av) = fanz(av) + 1
462
463                ENDDO
464!
465!--             Write data in NetCDF format
466                IF ( netcdf_output )  THEN
467#if defined( __netcdf )
468!
469!--                Check if a new time step has begun; if yes write data to
470!--                time axis
471                   IF ( current_var(av) > num_var(av) )  THEN
472                      current_var(av) = 1
473                      nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
474                                              (/ simulated_time /),        &
475                                              start = (/ time_step /),     &
476                                              count = (/ 1 /) )
477                      IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 7 )
478                   ENDIF
479
480!
481!--                Now write the data; this is mode dependent
482                   SELECT CASE ( modus )
483
484                      CASE ( 'XY' )
485                         nc_stat = NF90_PUT_VAR( id_set(av),                   &
486                                           id_var(av,current_var(av)),         &
487                                           pf(nxa:nxe,nya:nye),                &
488                             start = (/ 1, 1, current_level(av), time_step /), &
489                                      count = (/ nxe-nxa+1, nye-nya+1, 1, 1 /) )
490                         IF ( nc_stat /= NF90_NOERR )  THEN
491                            CALL handle_netcdf_error( 8 )
492                         ENDIF
493                 
494                      CASE ( 'XZ' )
495                         nc_stat = NF90_PUT_VAR( id_set(av),                   &
496                                           id_var(av,current_var(av)),         &
497                                           pf(nxa:nxe,ya_do:ye_do),            &
498                             start = (/ 1, current_level(av), 1, time_step /), &
499                                      count = (/ nxe-nxa+1, 1, ye_do-ya_do+1, 1 /) )
500                         IF ( nc_stat /= NF90_NOERR )  THEN
501                            CALL handle_netcdf_error( 9 )
502                         ENDIF
503
504                      CASE ( 'YZ' )
505                         nc_stat = NF90_PUT_VAR( id_set(av),                   &
506                                           id_var(av,current_var(av)),         &
507                                           pf(nxa:nxe,ya_do:ye_do),            &
508                             start = (/ current_level(av), 1, 1, time_step /), &
509                                      count = (/ 1, nxe-nxa+1, ye_do-ya_do+1, 1 /) )
510                         IF ( nc_stat /= NF90_NOERR )  THEN
511                            CALL handle_netcdf_error( 10 )
512                         ENDIF
513
514                   END SELECT
515
516!
517!--                Data is written, check if max level is reached
518                   current_level(av) = current_level(av) + 1
519                   IF ( current_level(av) > levels(av,current_var(av)) )  THEN
520                      current_level(av) = 1
521                      current_var(av)   = current_var(av) + 1
522                   ENDIF
523
524#endif
525                ENDIF
526
527             ENDDO
528
529          ENDIF
530
531998       IF ( danz /= 0  .AND.  .NOT. netcdf_parallel )  THEN
532!
533!--          Print the number of the arrays processed
534             WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
535             IF ( fanz(1) /= 0 )  THEN
536                WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
537             ENDIF
538
539!
540!--          Close all files and deallocate arrays
541             DO  id = 0, danz-1
542                CLOSE ( id+110 )
543             ENDDO
544             CLOSE ( 2 )
545             DEALLOCATE ( pf )
546
547!
548!--          Close the NetCDF file
549             IF ( netcdf_output )  THEN
550#if defined( __netcdf )
551                IF ( netcdf_0 )  THEN
552                   nc_stat = NF90_CLOSE( id_set(0) )
553                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 11 )
554                ENDIF
555                IF ( netcdf_1 )  THEN
556                   nc_stat = NF90_CLOSE( id_set(1) )
557                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 12 )
558                ENDIF
559#endif
560             ENDIF
561          ENDIF
562
563!
564!--       Output required cpu time
565          IF ( danz /= 0  .AND.  .NOT. netcdf_parallel )  THEN
566             CALL SYSTEM_CLOCK( count, count_rate )
567             cpu_end_time = REAL( count ) / REAL( count_rate )
568             WRITE (*,'(5X,A,F9.3,A)')  'Required cpu-time: ', &
569                                        cpu_end_time-cpu_start_time, ' sec'
570          ENDIF
571
572!
573!--       Choose the next cross-section
574          SELECT CASE ( modus )
575             CASE ( 'XY' )
576                modus = 'XZ'
577             CASE ( 'XZ' )
578                modus = 'YZ'
579             CASE ( 'YZ' )
580                modus = 'no'
581          END SELECT
582
583       ENDDO
584
585
586!
587!--    Combine the 3D-arrays
588       netcdf_parallel = .FALSE.
589
590!
591!--    Info-output
592       PRINT*, ' '
593
594!
595!--    Take current time
596       CALL SYSTEM_CLOCK( count, count_rate )
597       cpu_start_time = REAL( count ) / REAL( count_rate )
598
599!
600!--    Inquire whether a NetCDF file exists
601       INQUIRE( FILE='DATA_3D_NETCDF'//TRIM( model_string ), EXIST=netcdf_0 )
602
603!
604!--    Inquire whether a NetCDF file for time-averaged data exists
605       INQUIRE( FILE='DATA_3D_AV_NETCDF'//TRIM( model_string ), EXIST=netcdf_1 )
606
607       IF ( netcdf_0  .OR.  netcdf_1 )  THEN
608          netcdf_output = .TRUE.
609!
610!--       Inquire whether the NetCDF file is already complete (parallel output)
611          INQUIRE( FILE='NO_COMBINE_PLOT_FIELDS_3D', EXIST=netcdf_parallel )
612          IF ( netcdf_parallel )  THEN
613             netcdf_parallel = .TRUE.
614          ELSE
615             netcdf_parallel = .FALSE.
616          ENDIF
617       ELSE
618          netcdf_output = .FALSE.
619       ENDIF
620
621!
622!--    Check, if file from PE0 exists; not neccessary in case of parallel
623!--    PALM output
624       IF ( .NOT. netcdf_parallel )  THEN
625          danz = 0
626          WRITE (id_string,'(I6.6)')  danz
627          INQUIRE ( &
628               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM( id_string ),  &
629               EXIST=found )
630       ELSE
631          found = .FALSE.
632       ENDIF
633
634!
635!--    Find out the number of files and open them
636       DO  WHILE ( found )
637
638          OPEN ( danz+110, &
639               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM(id_string), &
640               FORM='UNFORMATTED')
641          danz = danz + 1
642          WRITE (id_string,'(I6.6)')  danz
643          INQUIRE ( &
644               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM(id_string), &
645               EXIST=found )
646
647       ENDDO
648
649#if defined( __netcdf )
650       IF ( netcdf_output )  THEN
651          IF ( netcdf_parallel )  THEN
652             PRINT*, '    NetCDF data are in one file (NetCDF4-format)', &
653                          ' - merging not neccessary'
654          ELSE
655             PRINT*, '    NetCDF output enabled'
656          ENDIF
657       ENDIF
658#else
659       IF ( netcdf_output )  THEN
660          PRINT*, '--- Sorry, no NetCDF support on this host'
661          netcdf_output = .FALSE.
662       ENDIF
663#endif
664       IF ( .NOT. netcdf_parallel )  THEN
665          IF ( danz /= 0 )  THEN
666             PRINT*, '    3D-data:     ', danz, ' file(s) found'
667          ELSE
668             PRINT*, '    no 3D-data file available'
669          ENDIF
670       ENDIF
671
672       IF ( netcdf_output  .AND. .NOT. netcdf_parallel  .AND.  danz /= 0 )  THEN
673#if defined( __netcdf )
674          DO  av = 0, 1
675
676             IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
677             IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
678
679!
680!--          Open NetCDF dataset
681             IF ( av == 0 )  THEN
682                filename = 'DATA_3D_NETCDF'//TRIM( model_string )
683             ELSE
684                filename = 'DATA_3D_AV_NETCDF'//TRIM( model_string )
685             ENDIF
686             nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
687             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 13 )
688
689
690!
691!--          Get the list of variables (order of variables corresponds with the
692!--          order of data on the binary file)
693             var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
694             nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
695                  var_list(av) )
696             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 14 )
697
698!
699!--          Inquire id of the time coordinate variable
700             nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
701             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 15 )
702
703!
704!--          Count number of variables; there is one more semicolon in the
705!--          string than variable names
706             num_var(av) = -1
707             DO  i = 1, LEN( var_list(av) )
708                IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) + 1
709             ENDDO
710
711!
712!--          Extract the variable names from the list and inquire their NetCDF
713!--          IDs
714             pos = INDEX( var_list(av), ';' )
715!
716!--          Loop over all variables
717             DO  i = 1, num_var(av)
718
719!
720!--             Extract variable name from list
721                var_list(av) = var_list(av)(pos+1:)
722                pos = INDEX( var_list(av), ';' )
723                var_name = var_list(av)(1:pos-1)
724
725!
726!--             Get variable ID from name
727                nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
728                     id_var(av,i) )
729                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 16 )
730
731             ENDDO
732
733          ENDDO    ! av=0,1
734
735#endif
736       ENDIF
737
738!
739!--    Read arrays, until the end of the file is reached
740       IF ( .NOT. netcdf_parallel )  THEN
741
742          current_var = 999999999
743          fanz = 0
744          DO  WHILE ( danz /= 0 )
745
746!
747!--          Loop over all files
748             DO  id = 0, danz-1
749!
750!--             File from PE0 contains at the beginning the index bounds
751!--             of PALM's total domain.
752!--             Allocate the array for storing the total domain data
753                IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
754!                   READ ( id+110 )  nxag, nxeg, nyag, nyeg
755                   READ ( id+110 )  nxa, nxe, nya, nye, nza, nze
756                   ALLOCATE ( pf3d(nxa:nxe,nya:nye,nza:nze) )
757                ENDIF
758
759!
760!--             Read output time
761                IF ( netcdf_output  .AND.  id == 0 )  THEN
762                   IF ( netcdf_1 )  THEN
763                      READ ( id+110, END=999 )  simulated_time, time_step, av
764                   ELSE
765!
766!--                   For compatibility with earlier PALM versions
767                      READ ( id+110, END=999 )  simulated_time, time_step
768                      av = 0
769                   ENDIF
770                ENDIF
771
772!
773!--             Read subdomain indices and grid point values
774                READ ( id+110, END=999 )  xa, xe, ya, ye, za, ze
775                ALLOCATE( pf3d_tmp(xa:xe,ya:ye,za:ze) )
776                READ ( id+110 )  pf3d_tmp
777
778                xxa = MAX( nxa, xa )
779                xxe = MIN( nxe, xe )
780                yya = MAX( nya, ya )
781                yye = MIN( nye, ye )
782                DO  k = za, ze
783                   DO  j = yya, yye
784                      DO  i = xxa, xxe
785                         pf3d(i,j,k) = pf3d_tmp(i,j,k)
786                      ENDDO
787                   ENDDO
788                ENDDO
789
790                DEALLOCATE( pf3d_tmp )
791                IF ( id == 0 )  fanz(av) = fanz(av) + 1
792
793             ENDDO
794
795!
796!--          Write data of the total domain in NetCDF format
797             IF ( netcdf_output )  THEN
798#if defined( __netcdf )
799!
800!--             Check if a new time step has begun; if yes write data to time
801!--             axis
802                IF ( current_var(av) > num_var(av) )  THEN
803                   current_var(av) = 1
804                   nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
805                                      (/ simulated_time /),&
806                                      start = (/ time_step /), count = (/ 1 /) )
807                   IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 17 )
808                ENDIF
809
810!
811!--             Now write the data
812                nc_stat = NF90_PUT_VAR( id_set(av), id_var(av,current_var(av)),&
813                                        pf3d(nxa:nxe,nya:nye,za:ze), start = (/ 1, 1, 1, time_step /),&
814                              count = (/ nxe-nxa+1, nye-nya+1, ze-za+1, 1 /) )
815                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 18 )
816
817                current_var(av) = current_var(av) + 1
818
819#endif
820             ENDIF
821
822          ENDDO
823
824       ENDIF
825
826999    IF ( danz /= 0  .AND.  .NOT. netcdf_parallel )  THEN
827!
828!--       Print the number of arrays processed
829          WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
830          IF ( fanz(1) /= 0 )  THEN
831             WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
832          ENDIF
833!
834!--       Close all files and deallocate array
835          DO  id = 0, danz-1
836             CLOSE ( id+110 )
837          ENDDO
838          CLOSE ( 2 )
839          DEALLOCATE ( pf3d )
840!
841!--       Close the NetCDF file
842          IF ( netcdf_output )  THEN
843#if defined( __netcdf )
844             IF ( netcdf_0 )  THEN
845                nc_stat = NF90_CLOSE( id_set(0) )
846                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 19 )
847             ENDIF
848             IF ( netcdf_1 )  THEN
849                nc_stat = NF90_CLOSE( id_set(1) )
850                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 20 )
851             ENDIF
852#endif
853          ENDIF
854
855!
856!--       Output required cpu time
857          CALL SYSTEM_CLOCK( count, count_rate )
858          cpu_end_time = REAL( count ) / REAL( count_rate )
859          WRITE (*,'(5X,A,F9.3,A)')  'Required cpu-time: ', &
860                                     cpu_end_time-cpu_start_time, ' sec'
861
862       ENDIF
863
864    ENDDO  ! models
865
866
867 CONTAINS
868
869
870    SUBROUTINE handle_netcdf_error( errno )
871!
872!--    Prints out a text message corresponding to the current NetCDF status
873
874       IMPLICIT NONE
875
876       INTEGER, INTENT(IN) ::  errno
877
878#if defined( __netcdf )
879       IF ( nc_stat /= NF90_NOERR )  THEN
880          PRINT*, '+++ combine_plot_fields  netcdf: ', av, errno, &
881                  TRIM( nf90_strerror( nc_stat ) )
882       ENDIF
883#endif
884
885    END SUBROUTINE handle_netcdf_error
886
887
888 END PROGRAM combine_plot_fields
889
890
891
Note: See TracBrowser for help on using the repository browser.