#include "options.inc"


      program nc_ufo
c=======================================================================
c     unfold netcdf variable 
c=======================================================================
      implicit none
      real,parameter :: version = 0.1
      integer :: ncin,varin,dimid(10),count(10),start(10)
      integer :: ncout,varout,imtf,imt,jmt,jmtf,km,itm
      integer :: nr_com,n,ierr,iok,i,j,ndims
      integer, allocatable :: kmt(:,:),kmu(:,:)
      character (len=80) :: fname,varname,infile,outfile
#include "netcdf.inc"
      integer :: iret,xtype,k,xid
      real, allocatable :: bufin(:,:,:,:),bufout(:,:,:,:)
      real, allocatable :: xt(:),xu(:)
      logical :: ugrid = .false.
      real :: spval

      print'(" This is nc_ufo.x :: version ",f5.3)',version
c  read command line input
      call get_commandline_argnr(nr_com)
      if (nr_com<3) goto 100

      call get_commandline_arg(1,infile,ierr)
      print*,' reading from file ',infile(1:len_trim(infile))
      call get_commandline_arg(2,outfile,ierr)
      print*,' writing to file ',outfile(1:len_trim(outfile))
      call get_commandline_arg(3,varname,ierr)
      print*,' variable ',varname(1:len_trim(varname))
      if (nr_com>3) then
        call get_commandline_arg(4,fname,ierr)
        if (fname == 'ugrid' ) ugrid=.true.
      endif
      if (ugrid) print*,' variable is on ugrid '
       
c   open KMT mask file
      print*,' reading from file kmt.ufl '
      iok=10
      open(iok,file='kmt.ufl',form='unformatted',status='old')
      read (iok); read (iok); read (iok) imt,jmt,km
      print*,' found in kmt.ufl imt=',imt,' jmt=',jmt,' km=',km
      allocate( kmt(imt,jmt), kmu(imt,jmt) )
      read (iok) ; read (iok) kmt
      close(iok)
      print*,' reading from file grid.ufl '
      open(iok,file='grid.ufl',form='unformatted',status='old')
      read (iok)
      read (iok) i,j,k
      if (i/=imt.or.j/=jmt.or.k/=km) then
        print*,' grid.ufl does not match kmt.ufl '
        stop
      endif
      allocate(xt(imt),xu(imt) )
      read (iok) (spval,i=1,(imt+jmt)*2+km*2+1),xt, xu
      close(iok)
      kmu=0.
      do j=1,jmt-1
        do i=1,imt-1
           kmu(i,j) = min (kmt(i,j), kmt(i+1,j), 
     &                     kmt(i,j+1), kmt(i+1,j+1))
        enddo
      enddo
      if (ugrid) then 
        kmt=kmu
        xt=xu
      endif

c   open input files and search for variable
      iret=nf_open(infile,NF_NOWRITE,ncin)
      if (iret /=0 ) goto 110
      iret=nf_inq_varid(ncin,varname,varin)
      if (iret /=0 ) goto 120

c   create output file
      ncout = nccre (outfile, NCNOCLOB, iret)
      if (iret /=0 ) goto 110
      iret=nf_set_fill(ncout, NF_NOFILL, iret)
      call copy_file_attributes(ncin,ncout)
c   define unfolded Longitude
      iret=nf_inq_vardimid(ncin,varin,dimid)
      iret=nf_inq_dimlen(ncin,dimid(1),imtf)
      iret=nf_inq_dimlen(ncin,dimid(2),jmtf)
      print*,' found imtf = ',imtf
      if (imt <= imtf ) then
       print*,' variable already unfolded or wrong kmt mask'
       stop
      endif
      if (jmt /= jmtf) then
       print*,' wrong kmt mask'
       stop
      endif
      iret=NF_INQ_DIMNAME (NCin, dimid(1), fname)
      call copy_dimension(ncin,ncout,fname,.true.,imt)
c   copy rest of dimensions and variables 
      call copy_variable_novalues(ncin,ncout,varname,.false.,0)
      iret = nf_enddef(ncout)

      print*,' writing new x--coord ',fname(1:len_trim(fname))
      iret=nf_inq_varid(ncout,fname,xid)
        if (iret.ne.0) print*,nf_strerror(iret)
      iret= nf_put_vara_real(ncout,xid,1,imt,xt)
        if (iret.ne.0) print*,nf_strerror(iret)

      print*,' writing new variable ',varname(1:len_trim(varname))
      start=1
      iret=nf_inq_varid(ncout,varname,varout)
      iret=nf_inq_vardimid(ncin,varin,dimid)
      iret = NF_INQ_VARTYPE (ncin, VARin,xtype)
      iret= nf_get_att_real (ncin,varin ,'missing_value',spval)
      iret=nf_inq_varndims(ncin,varin,ndims)
      if (ndims == 3) then
       iret=nf_inq_dimlen(ncin,dimid(3),itm)
       allocate( bufin(imtf,jmt,itm,1) )
       allocate( bufout(imt,jmt,itm,1) )
       iret= nf_get_vara_real(ncin,varin,start,(/imtf,jmt,itm/),bufin)
       bufout(1:imtf,:,:,1)= bufin(1:imtf,:,:,1)
c       bufout(imtf-1:imt,:,:,1)= bufout(1:imt-imtf+2,:,:,1)
       bufout(imtf:imt,:,:,1)= bufout(2:imt-imtf+2,:,:,1)
       do n=1,itm
          where( kmt ==0 ) bufout(:,:,n,1) = spval
       enddo
       iret= nf_put_vara_real(ncout,varout,start,(/imt,jmt,itm/),bufout)
      elseif (ndims==4) then
       iret=nf_inq_dimlen(ncin,dimid(4),itm)
       allocate( bufin(imtf,jmt,km,itm) )
       allocate( bufout(imt,jmt,km,itm) )
       iret= nf_get_vara_real(
     &             ncin,varin,start,(/imtf,jmt,km,itm/),bufin)
       bufout(1:imtf,:,:,:)= bufin(1:imtf,:,:,:)
c       bufout(imtf-1:imt,:,:,:)= bufout(1:imt-imtf+2,:,:,:)
       bufout(imtf:imt,:,:,:)= bufout(2:imt-imtf+2,:,:,:)
       do n=1,itm
        do k=1,km
         where( kmt <k ) bufout(:,:,k,n) = spval
        enddo
       enddo
       iret= nf_put_vara_real(
     &             ncout,varout,start,(/imt,jmt,km,itm/),bufout)
      endif
      print*,' ok'

      call ncclos (ncout, iret); if (iret /=0 ) goto 130
      call ncclos (ncin, iret); if (iret /=0 ) goto 130
      print*,'all done, bye'
      stop

 100  print*,' program needs command line parameter'
      print*,' input netcdf file names, output netcdf file name,'
      print*,' variable name to be unfolded [and switch ugrid]'
      stop
 110  print*,' problem opening NetCDF file ',fname(1:len_trim(fname))
      print*,' netcdf lib says: ',nf_strerror(iret)
      stop
 120  print*,' problem finding variable ',varname(1:len_trim(varname)),
     &   ' in file ',fname(1:len_trim(fname))
      print*,' netcdf lib says: ',nf_strerror(iret)
      stop
 130  print*,' problem closing NetCDF file '
      print*,' netcdf lib says: ',nf_strerror(iret)
      stop
      end program nc_ufo





c   dummy routines

      subroutine barrier
      end

      subroutine halt_stop(s)
      character*(*) s
      print*,s
      stop
      end



