#include "options.inc"

c
c     convert unformatted written data containing
c     surface and hydrostatic pressure to NetCDF format
c

c     real variables are 8 or 4 bytes long set by
c#define real_8
c
c
      module kind_mod
#ifdef real_8
      integer, parameter::  r_knd  = selected_real_kind(13)
#else
      integer, parameter::  r_knd  = selected_real_kind(6)
#endif
      end module kind_mod


      program cv_press
      use kind_mod
      implicit none
      character (len=80) :: infile,outfile
      integer :: ierr
      character (len=32) :: stamp
      integer :: imt,jmt,km,itm
      real (kind=r_knd),allocatable :: xt(:), yt(:), time(:)
      real (kind=r_knd), dimension(:,:,:), allocatable :: sfp, hsp
      real (kind=r_knd) :: spval
      real (kind=r_knd) :: tt
      integer :: io,n
#include "netcdf.inc"
      integer ncid,iret,i,j,k
      integer lat_tdim,itimedim, lon_tdim
      integer lat_tid,itimeid,lon_tid,sfpid,hspid
      integer dims(4), corner(4), edges(4)
      character name*24, unit*16, text*80
      print*,' Converting surface pressure 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,6; 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) imt,jmt
      read(io) spval
      allocate( xt(imt), yt(jmt) )
      read(io) xt,yt
      rewind(io)
      allocate( time(itm), sfp(imt,jmt,itm), hsp(imt,jmt,itm) )
      do n=1,itm
       read(io) stamp
       call read_stamp(stamp,tt)
       time(n)=tt
       print*,' at stamp=',stamp,
     &          ' (days since origin : ',time(n),')'

       do k=1,3; read(io); enddo
       read(io) sfp(:,:,n)
       read(io) hsp(:,:,n)
      enddo
      close(io)


      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c     dimensions
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      iTimedim  = ncddef(ncid, 'Time', nf_unlimited, iret)
c     grid variables
      dims(1)  = Lon_tdim
      Lon_tid  = ncvdef (ncid,'Longitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)

      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      sfpid=ncvdef (ncid,'sfp', NCFLOAT,3,dims,iret)
      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      hspid=ncvdef (ncid,'hsp', NCFLOAT,3,dims,iret)
c     attributes of the grid
      name = 'Latitude on T grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Longitude on T grid      '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Time                    '; unit = 'days            '
      call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
      call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c     attributes of variables
      name = 'Sea surface pressure'; unit = 'bar'
      call dvcdf_here(ncid,sfpid,name,24,unit,16,spval)
      name = 'Hydrostatic pressure'; unit = 'bar'
      call dvcdf_here(ncid,hspid,name,24,unit,16,spval)
      call ncendf(ncid, iret)

      corner=(/1,1,1,1/); 
      edges=(/imt,jmt,itm,1/)
#ifdef real_8
      iret= nf_put_vara_double (ncid,lat_tid,1,jmt,yt)
      iret= nf_put_vara_double (ncid,lon_tid,1,imt,xt)
      iret= nf_put_vara_double (ncid,itimeid,1,itm,time)
      iret= nf_put_vara_double (ncid,sfpid ,corner, edges,sfp)
      iret= nf_put_vara_double (ncid,hspid ,corner, edges,hsp)
#else
      iret= nf_put_vara_real (ncid,lat_tid,1,jmt,yt)
      iret= nf_put_vara_real (ncid,lon_tid,1,imt,xt)
      iret= nf_put_vara_real (ncid,itimeid,1,itm,time)
      iret= nf_put_vara_real (ncid,sfpid ,corner, edges,sfp)
      iret= nf_put_vara_real (ncid,hspid ,corner, edges,hsp)
#endif
      call ncclos (ncid, iret)

      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_press





      subroutine dvcdf_here(ncid,ivarid,name,iname,unit,iunit,spval)
c-----------------------------------------------------------------------
c define some standard attributes of variable ivarid in NetCDF file ncid 
c-----------------------------------------------------------------------
      use kind_mod
      implicit none
      integer ncid,ivarid,iname,iunit,iret
      character*(*) name, unit
      real (kind=r_knd):: spval
#include "netcdf.inc"
      call ncaptc(ncid,ivarid, 'long_name', NCCHAR,iname , name, iret) 
      call ncaptc(ncid,ivarid, 'units',     NCCHAR,iunit, unit, iret) 
#ifdef real_8
      iret=nf_put_att_double(NCID,iVARID, 'missing_value', NCFLOAT, 
     &                       1,spval)
      iret=nf_put_att_double(NCID,iVARID, '_FillValue', NCFLOAT, 
     &                       1,spval)
#else
      iret=nf_put_att_real(NCID,iVARID, 'missing_value', NCFLOAT, 
     &                       1,spval)
      iret=nf_put_att_real(NCID,iVARID, '_FillValue', NCFLOAT, 
     &                       1,spval)
#endif
c        if (iret.ne.0) print*,nf_strerror(iret)
      end subroutine dvcdf_here



      subroutine barrier
      end

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

