#include "options.inc"

c
c     convert unformatted written data containing
c     float positions to NetCDF format
c

      subroutine barrier
      end

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


      program cv_float
      implicit none
      character (len=80) :: infile,outfile
      integer :: ierr
      character (len=32) :: stamp
      integer :: nfl,itm
      integer :: io,k

#include "netcdf.inc"
      integer :: iret,ncid,timedim,nrdim,xid,yid,zid,timeid,dims(4)
      integer :: corner(2), edges(2)
      character (len=80) :: name,unit
      real :: tt

      real (kind=4) , allocatable :: pxyz(:,:)
      real (kind=4) :: time,spval = -1.e20

      print*,' Converting float data written by SPFLAME ',
     &       ' to NetCDF file '

      call get_commandline_arg(1,infile,ierr)
      if (ierr/=0) goto 99
      print*,' reading from file ',infile(1:len_trim(infile))
      call get_commandline_arg(2,outfile,ierr)
      if (ierr/=0) goto 99
      print*,' writing to file ',outfile(1:len_trim(outfile))

      call getunit (io, infile,'u s r ieee')

      itm=0
 10   continue
      do k=1,3; read (io,end=20); enddo
      itm=itm+1
      goto 10
 20   continue

      print*,' found ',itm,' time steps in file ',infile
      rewind(io)
      read(io) stamp
      read(io) nfl
      rewind(io)
      print*,' found ',nfl,' floats in file ',infile

      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c     dimensions
      Timedim   = ncddef(ncid,'Time', nf_unlimited, iret)
      nrdim     = ncddef(ncid,'Nr  ', nfl, iret)
c     variables
      Timeid   =  ncvdef(ncid,'Time', NCFLOAT,1,timedim,iret)
      name = 'Time'; unit = 'days'
      call ncaptc(ncid,Timeid, 'long_name', NCCHAR,4, name, iret) 
      call ncaptc(ncid,Timeid, 'units',    NCCHAR,4, unit, iret) 
      call ncaptc(ncid,Timeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
      dims=(/nrdim,timedim,1,1/)
      xid=ncvdef (ncid,'lon', NCFLOAT,2,dims,iret)
      yid=ncvdef (ncid,'lat', NCFLOAT,2,dims,iret)
      zid=ncvdef (ncid,'depth', NCFLOAT,2,dims,iret)
      name = 'Longitude on T grid'; unit = 'degrees_W'
      call dvcdf_here(ncid,xid,name,len_trim(name),unit,
     &            len_trim(unit),spval)
      name = 'Latitude on T grid'; unit = 'degrees_N'
      call dvcdf_here(ncid,yid,name,len_trim(name),unit,
     &            len_trim(unit),spval)
      name = 'Depth on T grid'; unit = 'm'
      call dvcdf_here(ncid,zid,name,len_trim(name),unit,
     &            len_trim(unit),spval)
      call ncendf(ncid, iret)

      allocate( pxyz(3,nfl) )

      do k=1,itm
       read (io) stamp
       call read_stamp(stamp,tt); time=tt
       print*,' at stamp=',stamp,
     &          ' (days since origin : ',time,')'
       read (io)
       read (io) pxyz

       iret= nf_put_vara_real (ncid,timeid, k,1,time)
       corner=(/1,k/); edges = (/nfl,1/)
       pxyz(1,:)=pxyz(1,:)-360.
       iret= nf_put_vara_real (ncid,xid,corner,edges,pxyz(1,:))
       iret= nf_put_vara_real (ncid,yid,corner,edges,pxyz(2,:))
       pxyz(3,:)=pxyz(3,:)/(-100.)
       iret= nf_put_vara_real (ncid,zid,corner,edges,pxyz(3,:))
      enddo

      call ncclos (ncid, iret)
      print*,' done'
      stop

 99   print*,' program needs command line parameter'
      print*,' 1. input binary file name (written by spflame)'
      print*,' 2. output NetCDF file name '

      end program cv_float



      subroutine dvcdf_here(ncid,ivarid,name,iname,unit,iunit,spval)
c-----------------------------------------------------------------------
c     define some standard attributes of variable ivarid in NetCDF file ncid 
c-----------------------------------------------------------------------
      implicit none
      integer ncid,ivarid,iname,iunit,iret
      character*(*) name, unit
      real (kind=4):: spval,vv
#include "netcdf.inc"
      vv=spval
      call ncaptc(ncid,ivarid, 'long_name', NCCHAR,iname , name, iret) 
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncaptc(ncid,ivarid, 'units',     NCCHAR,iunit, unit, iret) 
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncapt (ncid,ivarid, 'missing_value',NCFLOAT,1,vv,iret)
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncapt (ncid,ivarid, '_FillValue', NCFLOAT, 1,vv, iret)
        if (iret.ne.0) print*,nf_strerror(iret)
      end subroutine dvcdf_here





