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

 PROGRAM combine_plot_fields

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! 18/01/06  Output of time-averaged data
!
! 25/05/05  Errors removed
!
! 26/04/05  Output in NetCDF format, iso2d and avs output only if parameter
!           file exists
!
! 31/10/01  All comments and messages translated into English
!
! 23/02/99  Keine Bearbeitung komprimierter 3D-Daten
! Ursprungsversion vom 28/07/97
!
!
! Description:
! ------------
! This routine combines data of the PALM-subdomains into one file. In PALM
! every processor element opens its own file and writes 2D- or 3D-binary-data
! into it (different files are opened for xy-, xz-, yz-cross-sections and
! 3D-data). For plotting or analyzing these PE-data have to be collected and
! to be put into single files, which is done by this routine.
! Output format is NetCDF. Additionally, a data are output in a binary format
! readable by ISO2D-software (cross-sections) and by AVS (3D-data).
!
! This routine must be compiled with:
! decalpha:
!    f95 -cpp -D__netcdf -fast -r8 -I/usr/local/netcdf-3.5.1/include
!    -L/usr/local/netcdf-3.5.1/lib -lnetcdf
! IBM-Regatta:
!    xlf95 -qsuffix=cpp=f90 -WF,-D__netcdf -qrealsize=8 -q64 -qmaxmem=-1 -Q
!    -I /aws/dataformats/netcdf-3.6.0-p1/64-32/include-64
!    -L/aws/dataformats/netcdf-3.6.0-p1/64-32/lib -lnetcdf -O3
! IBM-Regatta KISTI:
!    xlf95 -qsuffix=cpp=f90 -WF,-D__netcdf -qrealsize=8 -q64 -qmaxmem=-1 -Q
!    -I /applic/netcdf64/src/f90
!    -L/applic/lib/NETCDF64 -lnetcdf -O3
! IBM-Regatta Yonsei (gfdl5):
!    xlf95 -qsuffix=cpp=f90 -WF,-D__netcdf -qrealsize=8 -q64 -qmaxmem=-1 -Q
!    -I /usr1/users/raasch/pub/netcdf-3.6.0-p1/include
!    -L/usr1/users/raasch/pub/netcdf-3.6.0-p1/lib -lnetcdf -O3
! IMUK:
!    ifort combine...f90 -o combine...x
!    -cpp -D__netcdf -I /muksoft/packages/netcdf/linux/include -axW -r8 -nbs
!    -Vaxlib -L /muksoft/packages/netcdf/linux/lib -lnetcdf
! NEC-SX6:
!    sxf90 combine...f90 -o combine...x 
!    -I /pool/SX-6/netcdf/netcdf-3.6.0-p1/include  -C hopt -Wf '-A idbl4'
!    -D__netcdf -L/pool/SX-6/netcdf/netcdf-3.6.0-p1/lib -lnetcdf
! Sun Fire X4600
!    pgf95 combine...f90 -o combine...x
!    -Mpreprocess -D__netcdf -I /home/usr5/mkanda/netcdf-3.6.1/src/f90 -r8
!    -fast -L/home/usr5/mkanda/netcdf-3.6.1/src/libsrc -lnetcdf
!------------------------------------------------------------------------------!

#if defined( __netcdf )
    USE netcdf
#endif

    IMPLICIT NONE

!
!-- Local variables
    CHARACTER (LEN=2)    ::  modus
    CHARACTER (LEN=7)    ::  id_char
    CHARACTER (LEN=10)   ::  dimname, var_name
    CHARACTER (LEN=40)   ::  filename

    CHARACTER (LEN=2000), DIMENSION(0:1) ::  var_list

    INTEGER, PARAMETER ::  spk = SELECTED_REAL_KIND( 6 )

    INTEGER ::  av, danz, i, id,             &
                j, nc_stat, nxa, nxag, nxe, nxeg, nya,   &
                nyag, nye, nyeg, nza, nzag, nze, nzeg, pos, time_step, xa, xe, &
                ya, ye, za, ze

    INTEGER, DIMENSION(0:1) ::  current_level, current_var, fanz, id_set, id_var_time, num_var

    INTEGER, DIMENSION(4) ::  id_dims_loc

    INTEGER, DIMENSION(0:1,4) ::  id_dims

    INTEGER, DIMENSION(0:1,1000) ::  id_var, levels

    LOGICAL ::  avs_output, compressed, found, iso2d_output, netcdf_output, &
                netcdf_0, netcdf_1

    REAL ::  dx, simulated_time
    REAL, DIMENSION(:),   ALLOCATABLE   ::  eta, ho, hu
    REAL, DIMENSION(:,:), ALLOCATABLE   ::  pf
    REAL(spk), DIMENSION(:,:,:), ALLOCATABLE ::  pf3d


!
!-- 2D-arrays for ISO2D
!-- Main loop for the three different cross-sections, starting with xy-section
    modus = 'XY'
    PRINT*, ''
    DO  WHILE ( modus == 'XY'  .OR.  modus == 'XZ'  .OR.  modus == 'YZ' )
!
!--    Check, if file from PE0 exists. If it does not exist, PALM did not
!--    create any output for this cross-section.
       danz = 0
       WRITE (id_char,'(A2,''_'',I4.4)')  modus, danz
       INQUIRE ( FILE='PLOT2D_'//id_char, EXIST=found )
!
!--    Find out the number of files (equal to the number of PEs which
!--    have been used in PALM) and open them
       DO  WHILE ( found )

          OPEN ( danz+110, FILE='PLOT2D_'//id_char, FORM='UNFORMATTED' )
          danz = danz + 1
          WRITE (id_char,'(A2,''_'',I4.4)')  modus, danz
          INQUIRE ( FILE='PLOT2D_'//id_char, EXIST=found )

       ENDDO

!
!--    Inquire whether an iso2d parameter file exists
       INQUIRE( FILE='PLOT2D_' // modus // '_GLOBAL', EXIST=iso2d_output )

!
!--    Inquire whether a NetCDF file exists
       INQUIRE( FILE='DATA_2D_' // modus // '_NETCDF', EXIST=netcdf_0 )

!
!--    Inquire whether a NetCDF file for time-averaged data exists
       INQUIRE( FILE='DATA_2D_' // modus // '_AV_NETCDF', EXIST=netcdf_1 )

       IF ( netcdf_0  .OR.  netcdf_1 )  THEN
          netcdf_output = .TRUE.
       ELSE
          netcdf_output = .FALSE.
       ENDIF

!
!--    Info-output
       PRINT*, ''
       PRINT*, '*** combine_plot_fields ***'
#if defined( __netcdf )
       IF ( netcdf_output )  PRINT*, '    NetCDF output enabled'
#else
       IF ( netcdf_output )  THEN
          PRINT*, '--- Sorry, no NetCDF support on this host'
          netcdf_output = .FALSE.
       ENDIF
#endif
       IF ( danz /= 0 )  THEN
          PRINT*, '    ',modus,'-section:  ', danz, ' file(s) found'
       ELSE
          PRINT*, '    no ', modus, '-section data available'
       ENDIF

       IF ( netcdf_output  .AND.  danz /= 0 )  THEN
#if defined( __netcdf )
          DO  av = 0, 1

             IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
             IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE

!
!--          Open NetCDF dataset
             IF ( av == 0 )  THEN
                filename = 'DATA_2D_' // modus // '_NETCDF'
             ELSE
                filename = 'DATA_2D_' // modus // '_AV_NETCDF'
             ENDIF
             nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 1 )

!
!--          Get the list of variables (order of variables corresponds with the
!--          order of data on the binary file)
             var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
             nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
                                     var_list(av) )
             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 2 )

!
!--          Inquire id of the time coordinate variable
             nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 3 )

!
!--          Count number of variables; there is one more semicolon in the
!--          string than variable names
             num_var(av) = -1
             DO  i = 1, LEN( var_list(av) )
                IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) + 1
             ENDDO

!
!--          Extract the variable names from the list and inquire their
!--          NetCDF IDs
             pos = INDEX( var_list(av), ';' )
!
!--          Loop over all variables
             DO  i = 1, num_var(av)

!
!--             Extract variable name from list
                var_list(av) = var_list(av)(pos+1:)
                pos = INDEX( var_list(av), ';' )
                var_name = var_list(av)(1:pos-1)

!
!--             Get variable ID from name
                nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
                                          id_var(av,i) )
                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 4 )

!
!--             Get number of x/y/z levels for that variable
                nc_stat = NF90_INQUIRE_VARIABLE( id_set(av), id_var(av,i), &
                                                 dimids = id_dims_loc )
                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 5 )
                id_dims(av,:) = id_dims_loc

!
!--             Inquire dimension ID
                DO  j = 1, 4
                   nc_stat = NF90_INQUIRE_DIMENSION( id_set(av), id_dims(av,j),&
                                                     dimname, levels(av,i) )
                   IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 6 )

                   IF ( modus == 'XY' .AND. INDEX(dimname, 'z') /= 0 )  EXIT
                   IF ( modus == 'XZ' .AND. INDEX(dimname, 'y') /= 0 )  EXIT
                   IF ( modus == 'YZ' .AND. INDEX(dimname, 'x') /= 0 )  EXIT
                ENDDO

             ENDDO

          ENDDO   ! av = 0, 1

       ENDIF
#endif

!
!--    Read the arrays, as long as the end of the file is reached
       fanz          =         0
       current_level =         1
       current_var   = 999999999

       DO  WHILE ( danz /= 0 )

!
!--       Loop over all files (reading data of the subdomains)
          DO  id = 0, danz-1
!
!--          File from PE0 contains special information at the beginning,
!--          concerning the lower and upper indices of the total-domain used
!--          in PALM (nxag, nxeg, nyag, nyeg) and the lower and upper
!--          indices of the array to be writte by this routine (nxa, nxe, nya,
!--          nye). Usually in the horizontal directions nxag=-1 and nxa=0
!--          while all other variables have the same value (i.e. nxeg=nxe).
!--          Allocate necessary arrays, open the output file and write
!--          the coordinate informations needed by ISO2D.
             IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
                READ ( id+110 )  nxag, nxeg, nyag, nyeg
                READ ( id+110 )  nxa, nxe, nya, nye
                ALLOCATE ( eta(nya:nye), ho(nxa:nxe), hu(nxa:nxe), &
                           pf(nxag:nxeg,nyag:nyeg) )
                READ ( id+110 )  dx, eta, hu, ho

                IF ( iso2d_output )  THEN
                   OPEN ( 2, FILE='PLOT2D_'//modus, FORM='UNFORMATTED' )
                   WRITE ( 2 )  dx, eta, hu, ho
                ENDIF
             ENDIF
!
!--          Read output time
             IF ( netcdf_output  .AND.  id == 0 )  THEN
                IF ( netcdf_1 )  THEN
                   READ ( id+110, END=998 )  simulated_time, time_step, av
                ELSE
!
!--                For compatibility with earlier PALM versions
                   READ ( id+110, END=998 )  simulated_time, time_step
                   av = 0
                ENDIF
             ENDIF
!
!--          Read subdomain indices
             READ ( id+110, END=998 )  xa, xe, ya, ye
!
!--          IF the PE made no output (in case that no part of the
!--          cross-section is situated on this PE), indices have the
!--          value -1
             IF ( .NOT. ( xa == -1  .AND.  xe == -1  .AND. &
                          ya == -1  .AND.  ye == -1 ) )  THEN
!
!--              Read the subdomain grid-point values
                 READ ( id+110 )  pf(xa:xe,ya:ye)
             ENDIF
             IF ( id == 0 )  fanz(av) = fanz(av) + 1

          ENDDO
!
!--       Write the data of the total domain cross-section
          IF ( iso2d_output )  WRITE ( 2 )  pf(nxa:nxe,nya:nye)
       
!
!--       Write same data in NetCDF format
          IF ( netcdf_output )  THEN
#if defined( __netcdf )
!
!--          Check if a new time step has begun; if yes write data to time axis
             IF ( current_var(av) > num_var(av) )  THEN
                current_var(av) = 1
                nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
                                        (/ simulated_time /),        &
                                        start = (/ time_step /),     &
                                        count = (/ 1 /) )
                IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 7 )
             ENDIF

!
!--          Now write the data; this is mode dependent
             SELECT CASE ( modus )

                CASE ( 'XY' )
                   nc_stat = NF90_PUT_VAR( id_set(av),                         &
                                           id_var(av,current_var(av)),         &
                                           pf(nxa:nxe,nya:nye),                &
                             start = (/ 1, 1, current_level(av), time_step /), &
                                      count = (/ nxe-nxa+1, nye-nya+1, 1, 1 /) )
                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 8 )
                  
                CASE ( 'XZ' )
                   nc_stat = NF90_PUT_VAR( id_set(av),                         &
                                           id_var(av,current_var(av)),         &
                                           pf(nxa:nxe,nya:nye),                &
                             start = (/ 1, current_level(av), 1, time_step /), &
                                      count = (/ nxe-nxa+1, 1, nye-nya+1, 1 /) )
                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 9 )

                CASE ( 'YZ' )
                   nc_stat = NF90_PUT_VAR( id_set(av),                         &
                                           id_var(av,current_var(av)),         &
                                           pf(nxa:nxe,nya:nye),                &
                             start = (/ current_level(av), 1, 1, time_step /), &
                                      count = (/ 1, nxe-nxa+1, nye-nya+1, 1 /) )
                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 10 )

             END SELECT

!
!--          Data is written, check if max level is reached
             current_level(av) = current_level(av) + 1
             IF ( current_level(av) > levels(av,current_var(av)) )  THEN
                current_level(av) = 1
                current_var(av)   = current_var(av) + 1
             ENDIF

          ENDIF
#endif

       ENDDO

998    IF ( danz /= 0 )  THEN
!
!--       Print the number of the arrays processed
          WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
          IF ( fanz(1) /= 0 )  THEN
             WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
          ENDIF

!
!--       Close all files and deallocate arrays
          DO  id = 0, danz-1
             CLOSE ( id+110 )
          ENDDO
          CLOSE ( 2 )
          DEALLOCATE ( eta, ho, hu, pf )
       ENDIF

!
!--    Close the NetCDF file
       IF ( netcdf_output  .AND.  danz /= 0 )  THEN
#if defined( __netcdf )
          IF ( netcdf_0 )  THEN
             nc_stat = NF90_CLOSE( id_set(0) )
             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 11 )
          ENDIF
          IF ( netcdf_1 )  THEN
             nc_stat = NF90_CLOSE( id_set(1) )
             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 12 )
          ENDIF
#endif
       ENDIF

!
!--    Choose the next cross-section
       SELECT CASE ( modus )
          CASE ( 'XY' )
             modus = 'XZ'
          CASE ( 'XZ' )
             modus = 'YZ'
          CASE ( 'YZ' )
             modus = 'no'
       END SELECT

    ENDDO


!
!-- Combine the 3D-arrays

!
!-- Inquire whether an avs fld file exists
    INQUIRE( FILE='PLOT3D_FLD', EXIST=avs_output )

!
!-- Inquire whether a NetCDF file exists
    INQUIRE( FILE='DATA_3D_NETCDF', EXIST=netcdf_0 )

!
!-- Inquire whether a NetCDF file for time-averaged data exists
    INQUIRE( FILE='DATA_3D_AV_NETCDF', EXIST=netcdf_1 )

    IF ( netcdf_0  .OR.  netcdf_1 )  THEN
       netcdf_output = .TRUE.
    ELSE
       netcdf_output = .FALSE.
    ENDIF

!
!-- Check, if file from PE0 exists
    danz = 0
    WRITE (id_char,'(I4.4)')  danz
    INQUIRE ( FILE='PLOT3D_DATA_'//TRIM( id_char ), EXIST=found )

!
!-- Combination only works, if data are not compressed. In that case, 
!-- PALM created a flag file (PLOT3D_COMPRESSED)
    INQUIRE ( FILE='PLOT3D_COMPRESSED', EXIST=compressed )

!
!-- Find out the number of files and open them
    DO  WHILE ( found  .AND.  .NOT. compressed )

       OPEN ( danz+110, FILE='PLOT3D_DATA_'//TRIM( id_char ), &
              FORM='UNFORMATTED')
       danz = danz + 1
       WRITE (id_char,'(I4.4)')  danz
       INQUIRE ( FILE='PLOT3D_DATA_'//TRIM( id_char ), EXIST=found )

    ENDDO

!
!-- Info-output
    PRINT*, ' '
    PRINT*, '*** combine_plot_fields ***'
#if defined( __netcdf )
       IF ( netcdf_output )  PRINT*, '    NetCDF output enabled'
#else
       IF ( netcdf_output )  THEN
          PRINT*, '--- Sorry, no NetCDF support on this host'
          netcdf_output = .FALSE.
       ENDIF
#endif
    IF ( danz /= 0 )  THEN
       PRINT*, '    3D-data:     ', danz, ' file(s) found'
    ELSE
       IF ( found .AND. compressed )  THEN
          PRINT*, '+++ no 3D-data processing, since data are compressed'
       ELSE
          PRINT*, '    no 3D-data file available'
       ENDIF
    ENDIF

    IF ( netcdf_output  .AND.  danz /= 0 )  THEN
#if defined( __netcdf )
       DO  av = 0, 1

          IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
          IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE

!
!--       Open NetCDF dataset
          IF ( av == 0 )  THEN
             filename = 'DATA_3D_NETCDF'
          ELSE
             filename = 'DATA_3D_AV_NETCDF'
          ENDIF
          nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 13 )


!
!--       Get the list of variables (order of variables corresponds with the
!--       order of data on the binary file)
          var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
          nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
                                  var_list(av) )
          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 14 )

!
!--       Inquire id of the time coordinate variable
          nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 15 )

!
!--       Count number of variables; there is one more semicolon in the string
!--       than variable names
          num_var(av) = -1
          DO  i = 1, LEN( var_list(av) )
             IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) + 1
          ENDDO

!
!--       Extract the variable names from the list and inquire their NetCDF IDs
          pos = INDEX( var_list(av), ';' )
!
!--       Loop over all variables
          DO  i = 1, num_var(av)

!
!--          Extract variable name from list
             var_list(av) = var_list(av)(pos+1:)
             pos = INDEX( var_list(av), ';' )
             var_name = var_list(av)(1:pos-1)

!
!--          Get variable ID from name
!             print*, '*** find id for "',TRIM( var_name ),'" begin'
             nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
                                       id_var(av,i) )
             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 16 )
!             print*, '*** find id for "',TRIM( var_name ),'" end'

          ENDDO

       ENDDO    ! av=0,1

    ENDIF
#endif

!
!-- Read arrays, until the end of the file is reached
    current_var = 999999999
    fanz = 0
    DO  WHILE ( danz /= 0 )

!
!--    Loop over all files
       DO  id = 0, danz-1
!
!--       File from PE0 contains special information at the beginning,
!--       concerning the lower and upper indices of the total-domain used
!--       in PALM (nxag, nxeg, nyag, nyeg, nzag, nzeg) and the lower and upper
!--       indices of the array to be written by this routine (nxa, nxe, nya,
!--       nye, nza, nze). Usually nxag=-1 and nxa=0, nyag=-1 and nya=0,
!--       nzeg=nz and nze=nz_plot3d.
!--       Allocate necessary array and open the output file.
          IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
             READ ( id+110 )  nxag, nxeg, nyag, nyeg, nzag, nzeg
             READ ( id+110 )  nxa, nxe, nya, nye, nza, nze
             ALLOCATE ( pf3d(nxag:nxeg,nyag:nyeg,nzag:nzeg) )
             IF ( avs_output )  THEN
                OPEN ( 2, FILE='PLOT3D_DATA', FORM='UNFORMATTED' )
             ENDIF
          ENDIF

!
!--       Read output time
          IF ( netcdf_output  .AND.  id == 0 )  THEN
             IF ( netcdf_1 )  THEN
                READ ( id+110, END=999 )  simulated_time, time_step, av
             ELSE
!
!--             For compatibility with earlier PALM versions
                READ ( id+110, END=999 )  simulated_time, time_step
                av = 0
             ENDIF
          ENDIF

!
!--       Read subdomain indices and grid point values
          READ ( id+110, END=999 )  xa, xe, ya, ye, za, ze
          READ ( id+110 )  pf3d(xa:xe,ya:ye,za:ze)
          IF ( id == 0 )  fanz(av) = fanz(av) + 1

       ENDDO

!
!--    Write data of the total domain
       IF ( avs_output )  WRITE ( 2 )  pf3d(nxa:nxe,nya:nye,nza:nze)
       
!
!--    Write same data in NetCDF format
       IF ( netcdf_output )  THEN
#if defined( __netcdf )
!
!--       Check if a new time step has begun; if yes write data to time axis
          IF ( current_var(av) > num_var(av) )  THEN
             current_var(av) = 1
             nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
                                     (/ simulated_time /),&
                                     start = (/ time_step /), count = (/ 1 /) )
             IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 17 )
          ENDIF

!
!--       Now write the data
          nc_stat = NF90_PUT_VAR( id_set(av), id_var(av,current_var(av)), &
                                  pf3d(nxa:nxe,nya:nye,nza:nze),      &
                                  start = (/ 1, 1, 1, time_step /),   &
                              count = (/ nxe-nxa+1, nye-nya+1, nze-nza+1, 1 /) )
          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 18 )

          current_var(av) = current_var(av) + 1

#endif
       ENDIF

    ENDDO

999 IF ( danz /= 0 )  THEN
!
!--    Print the number of arrays processed
       WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
       IF ( fanz(1) /= 0 )  THEN
          WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
       ENDIF
!
!--    Close all files and deallocate array
       DO  id = 0, danz-1
          CLOSE ( id+110 )
       ENDDO
       CLOSE ( 2 )
       DEALLOCATE ( pf3d )
!
!--    Close the NetCDF file
       IF ( netcdf_output )  THEN
#if defined( __netcdf )
          IF ( netcdf_0 )  THEN
             nc_stat = NF90_CLOSE( id_set(0) )
             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 19 )
          ENDIF
          IF ( netcdf_1 )  THEN
             nc_stat = NF90_CLOSE( id_set(1) )
             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 20 )
          ENDIF
#endif
       ENDIF
    ENDIF


 CONTAINS


    SUBROUTINE handle_netcdf_error( errno )
!
!--    Prints out a text message corresponding to the current NetCDF status

       IMPLICIT NONE

       INTEGER, INTENT(IN) ::  errno

       IF ( nc_stat /= NF90_NOERR )  THEN
          PRINT*, '+++ combine_plot_fields  netcdf: ', av, errno, &
                  TRIM( nf90_strerror( nc_stat ) )
       ENDIF

    END SUBROUTINE handle_netcdf_error


 END PROGRAM combine_plot_fields