!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! this program reads several input netcdf files, presumably containing "compressed
! by gathering" data, and combines them into a single output file
#define __NF_ASRT__(ierr) call nfu_check_err(ierr,__FILE__,__LINE__)
program decompress

  use nfu_mod
  use nfu_compress_mod
  implicit none
  include 'netcdf.inc'

  integer, parameter :: PATH_MAX = 1024 ! max len of the file name; 

  character(PATH_MAX), allocatable :: files(:) ! names of all files on the command line
  character(PATH_MAX)              :: outfile  ! name of the output file
  integer :: nfiles    ! number of files on command line
  integer :: debug = 0 ! debug output verbosity level
  integer, allocatable :: input(:)             ! netcdf IDs of input files
  integer :: i,iret,ncid,dimid,varid,varid1,xtype,dimlen,vsize,ndims,nvars,natts
  integer :: dimids(NF_MAX_VAR_DIMS), dimlens(NF_MAX_VAR_DIMS)
  character(NF_MAX_NAME) :: dimname,varname,attname
  logical :: is_dim, is_compressed, has_records, add_missing_value = .FALSE.
  real   , allocatable :: buffer(:)
  logical, allocatable :: mask(:)
  real :: missing

  ! get command line options and list of files
  call parse_command_line() ! modigies global data!

  call assert(nfiles>0,'at least one input file must be specified')
  if(debug>0) then
     do i = 1,nfiles
        write(*,'("input file",i3,":",a)')i, '"'//trim(files(i))//'"'
     enddo
     write(*,'("output file:",a)')'"'//trim(outfile)//'"'
  endif

  ! open all input files
  allocate(input(nfiles))
  do i = 1,nfiles
     __NF_ASRT__(nf_open(files(i),NF_NOWRITE,input(i)))
  enddo

  ! create output file
  __NF_ASRT__(nf_create(outfile,NF_CLOBBER,ncid))

  ! Create netcdf structure in the output NetCDF file, using last input file
  ! as a template.

  ! clone all dimensions except compressed ones; compressed are just skipped
  __NF_ASRT__(nf_inq_ndims(input(nfiles),ndims))
  do dimid = 1,ndims
     __NF_ASRT__(nfu_inq_dim(input(nfiles),dimid,dimname=dimname,dimlen=dimlen,is_unlim=has_records))
     if(nfu_inq_att(input(nfiles),dimname,'compress')==NF_NOERR) cycle
     if(has_records)&
          dimlen=NF_UNLIMITED
     if(debug>0)&
          write(*,*)'defining dimension "'//trim(dimname)//'" with length',dimlen
     __NF_ASRT__(nf_def_dim(ncid,dimname,dimlen,i)) ! i is just a space for id
  enddo

  ! clone all variable definitions, replacing compressed dimensions with sets
  ! of uncompressed ones
  __NF_ASRT__(nf_inq_nvars(input(nfiles),nvars))
  do varid = 1,nvars
     iret = nfu_inq_compressed_var(input(nfiles),varid,varname,xtype,ndims,dimids,&
          natts=natts,is_dim=is_dim, is_compressed=is_compressed)
     if(debug>0)&
          write(*,*)'defining variable "'//trim(varname)//'"'
     __NF_ASRT__(iret) ! just because the line is going to be too long
     if(is_dim.and.is_compressed) cycle
     do i = 1,ndims
        __NF_ASRT__(nf_inq_dimname(input(nfiles),dimids(i),dimname))
        __NF_ASRT__(nf_inq_dimid(ncid,dimname,dimids(i)))
     enddo
     __NF_ASRT__(nf_def_var(ncid,varname,xtype,ndims,dimids,varid1))
     do i=1,natts
        __NF_ASRT__(nf_inq_attname(input(nfiles),varid,i,attname))
        __NF_ASRT__(nf_copy_att(input(nfiles),varid,attname,ncid,varid1))
     enddo
     if(add_missing_value.and.is_compressed) then
        ! check if missing_value or _FillValue attributes are present 
        if ( nf_inq_atttype(input(nfiles),varid,'missing_value',iret)/=NF_NOERR .and. &
             nf_inq_atttype(input(nfiles),varid,'_FillValue',iret)/=NF_NOERR ) then
           ! if not, define the missing value attribute
           select case(xtype)
           case(NF_DOUBLE) 
              missing = NF_FILL_DOUBLE
           case(NF_FLOAT)
              missing = NF_FILL_FLOAT
           case(NF_INT)
              missing = NF_FILL_INT
           end select
	   ! and add it to the output variable
	   __NF_ASRT__(nf_put_att_double(ncid,varid1,'missing_value',xtype,1,missing))
	endif
     endif
  enddo

  ! clone all global attributes
  __NF_ASRT__(nf_inq_natts(input(nfiles),natts))
  do i = 1,natts
     __NF_ASRT__(nf_inq_attname(input(nfiles),NF_GLOBAL,i,attname))
     __NF_ASRT__(nf_copy_att(input(nfiles),NF_GLOBAL,attname,ncid,NF_GLOBAL))
  enddo

  ! ---- end of definition stage
  __NF_ASRT__(nf_enddef(ncid))


  ! grow unlimited dimension, if necessary
  do varid=1,nvars
     __NF_ASRT__(nfu_inq_compressed_var(input(nfiles),varid,name=varname,dimlens=dimlens,has_records=has_records))
     if(has_records) then
        ! just write an integer at the very end of the variable -- that extends the 
        ! record dimensions as well
        __NF_ASRT__(nf_inq_varid(ncid,varname,varid))
        __NF_ASRT__(nf_put_var1_int(ncid,varid,dimlens,0))
        exit ! this loop
     endif
  enddo

  ! gather and copy data
  __NF_ASRT__(nf_inq_nvars(ncid,nvars))
  do varid = 1,nvars
!     __NF_ASRT__(nfu_inq_compressed_var(ncid,varid,name=varname,varsize=vsize))
     __NF_ASRT__(nfu_inq_var(ncid,varid,varname,xtype,varsize=vsize))
     if(debug>0) &
          write(*,*)'processing var "'//trim(varname)//'"'
     allocate(buffer(vsize),mask(vsize))
     mask(:) = .false.

     ! obtain the missing value 
     if(nfu_get_att(ncid,varname,'missing_value',missing)==NF_NOERR) then
        ! do nothing, the value is already in the "missing" variable
     else if(nfu_get_att(ncid,varname,'_FillValue',missing)==NF_NOERR) then
        ! do nothing, the value is already in the "missing" variable
     else
        ! get fill value for the type instead of the missing value
        select case(xtype)
        case(NF_DOUBLE) 
           missing = NF_FILL_DOUBLE
        case(NF_FLOAT)
           missing = NF_FILL_FLOAT
        case(NF_INT)
           missing = NF_FILL_INT
        end select
     endif
     ! fill the buffer with the missing value
     buffer=missing

     ! read the variable
     do i=1,nfiles
        __NF_ASRT__(nfu_get_compressed_var_r8n(input(i),varname,buffer,mask))
     enddo
     ! write the variable
     __NF_ASRT__(nfu_put_var_r8(ncid,varname,buffer))
     deallocate(buffer,mask)
  enddo
     
  __NF_ASRT__(nf_close(ncid))

contains ! -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

! ---- parses command line arguments, getting options and gathering list of
! file names
! NOTE: updates global variables.
subroutine parse_command_line()
  character(PATH_MAX) :: arg, param

  integer :: nargs     ! number of command-line arguments
  logical :: do_interpret_arguments
  integer :: i, iostat

  integer, external :: iargc

  nargs = iargc()
  if(nargs==0) then
     call usage()
     call exit(1)
  endif
  
  allocate(files(nargs))  ! allocate storage for all file names
  do_interpret_arguments = .true.
  add_missing_value = .false.
  i=1        ! counter of all command-line arguments
  nfiles = 0 ! counter of input files
  do while (i<=nargs)
     call getarg(i,arg)
     if(debug>1) write(*,*)'argument ',i, trim(arg)
     if(arg(1:1)=='-'.and.do_interpret_arguments) then
        select case(trim(arg))
        case('--')
           do_interpret_arguments = .false.

        case('-D','--debug-level')
           call assert(i<nargs,trim(arg)//' flag must be followed by integer verbosity level')
           call getarg(i+1,param)
           read(param,*,iostat=iostat) debug
           call assert(iostat==0,trim(arg)//' flag must be followed by integer verbosity level')
           i=i+1

        case('-m','--add-missing-value')
	   add_missing_value = .TRUE.

        case ('-h','-?','--help')
           call usage()
           call exit(1)

        case default
           call usage()
           call assert(.false.,'argument "'//trim(arg)//'" is illegal')
        end select
     else
        ! argument is either input or output file name, add it to the list
        nfiles = nfiles+1
        files(nfiles) = arg
     endif
     i = i+1
  enddo
  if (nfiles>0) then
     outfile = files(nfiles)
     nfiles  = nfiles-1
  endif
end subroutine


! ---- prints usage information
subroutine usage()
! this program reads several input netcdf files, presumably containing "compressed
! by gathering" data, and combines them into a single output file
  character(len=PATH_MAX) :: name
  call getarg(0,name)
  write(*,'(a)')'Converts one or several compressed-by-gathering netcdf file into'
  write(*,'(a)')'one regular netcdf. Normally used to convert lm3 restarts into a'
  write(*,'(a)')'form suitable for visualization applications.'
  write(*,'(a)')
  write(*,'(a)')'Usage:'
  write(*,'(a)')'  '//trim(name)//' [-D debug-level] [-m] in.nc [...] out.nc'
  write(*,'(a)')
  write(*,'(a)')'-D debug-level   Specifies level of debug output verbosity'
  write(*,'(a)')'-m               Forces adding a missing_value attribute to the variables'
  write(*,'(a)')'                 that do not have it'
  write(*,'(a)')'in.nc            Input file name(s)'
  write(*,'(a)')'out.nc           Output file name'
  write(*,'(a)')
  write(*,'(a)')'WARNING: output file is overwritten.'
end subroutine


! ===========================================================================
! ---- prints error message an exits if condition is not satisfied
subroutine assert(cond,message)
  logical     , intent(in) :: cond    ! condition to check
  character(*), intent(in) :: message ! error message to print if condition is not satisfied
  
  if(.not.cond) then
     write(*,*) 'ERROR :: ',trim(message)
     call exit(1)
  endif
end subroutine

end program decompress