#include "options.inc"

c
c     convert binary output of time series written by SPFLAME
c     to NetCDF format
c
c     real variables are 8 or 4 bytes long set by
c#define real_8



      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



      module cv_anna_module
      use kind_mod
      integer :: imt,jmt,km,nr_var,itm,nt,nr_var_given,itt
      real (kind=r_knd),allocatable::xt(:),yt(:),xu(:),yu(:),zt(:),zw(:)
      real (kind=r_knd) :: spval
      integer :: year,month,day,domain,min,sec,hour
      character (len=80) :: infile,outfile
      integer :: block_size=0
      end module cv_anna_module

      program cv_anna
      use cv_anna_module
      implicit none
#include "netcdf.inc"

      character (len=80) :: ident
      character  :: c9*9, c1*1, c2*2, c4*4,c8*8,c6*6
      integer :: i,j,k,n,nn,ierr,io,fsize
      character (len=32) :: stamp
      real (kind=r_knd) :: time,tt
      real (kind=r_knd) ,allocatable :: buf(:,:)

      integer :: ncid,iret,varid
      integer :: dims(4), corner(4), edges(4)
      character (len=80) :: name, unit, text

      print*,' Converting averaged quantities written by SPFLAME ',
     &       ' to NetCDF file '

      call get_commandline_arg(1,infile,ierr)
      if (ierr/=0) goto 99
      nr_var_given = -999
      call  get_commandline_arg_int(2,nr_var_given,iret)

      print*,' reading from file ',infile(1:len_trim(infile))
      read(infile,'(a5,i3,a2,i4,a1,i2,a1,i2,a1,i8)') 
     &            c9,domain,c2,year,c1,month,c1,day,c1,itt
      write(outfile,'("anna_",i3,"_y",i4,"m",i2,"d",i2,"_",i8,".cdf")') 
     &            domain,year,month,day,itt
      call replace_space_zero(outfile)


      print*,' writing to file ',outfile(1:len_trim(outfile))

      call getunit(io,infile,'usr ieee')

      read(io) imt,jmt,km,nt,nr_var
c      read(io) imt,jmt,km,nr_var; nt=2
      print*,' imt=',imt,' jmt=',jmt,' km=',km,' nt=',nt
      if (nr_var_given /= nr_var .and. nr_var_given >0) then
       print*,' file says that ',nr_var,' variables in file '
       print*,' however, that was overuled by command line input'
       nr_var = nr_var_given
       print*,' -> searching for ',nr_var,' variables in file '
      else
       print*,' found ',nr_var,' variables in file '
      endif

      allocate( xt(imt),yt(jmt),xu(imt),yu(jmt),zt(km),zw(km))
      read(io) xt,yt,xu,yu,zt,zw,spval
      zt=zt/100.
      zw=zw/100.

      call cv_anna_define(io)

      fsize=(imt*2+jmt*2+km*2)*4+1024

      itm=0
      iret=nf_open(outfile,NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)

 20   continue
      read(io,end=30) stamp,itt

      print*,' reading time step ', stamp
      print*,' estimated file size: ',fsize/1e6,' Mbytes'


c    check is file is too large 
      fsize=fsize+block_size*4
      if (fsize> 1500e6) then
        print*,' NetCDF file is getting too large'
        print*,' (estimated file size is ',fsize/1e6, 'Mbytes)'
        print*,' closing file ',outfile(1:len_trim(outfile))
        call ncclos (ncid, iret)
         read  (stamp, '(a6,i2,a1,i2,a1,i4,a8,i2,a1,i2,a1,i2)')
     &       c6, month, c1, day, c1, year, c8, hour,
     &       c1, min, c1, sec
        write(outfile,
     &        '("anna_",i3,"_y",i4,"m",i2,"d",i2,"_",i8,".cdf")') 
     &            domain,year,month,day,itt
        call replace_space_zero(outfile)
        print*,' opening file ',outfile(1:len_trim(outfile))
        backspace (io)
        call cv_anna_define(io)
        read(io)
        fsize=(imt*2+jmt*2+km*2)*4+1024
        itm=0
        iret=nf_open(outfile,NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
      endif



      itm=itm+1

      call read_stamp(stamp,tt); time=tt
      iret=nf_inq_varid(ncid,'Time',varid)
      call nf_put_var(ncid,varid, itm,1,time)

      do n=1,nr_var

        read(io) ident

        if (ident(1:3) == 'dml' ) then
         allocate( buf(imt,jmt) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'dml',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'uvts') then
         allocate( buf(imt,jmt) )
         do k=1,km
          corner=(/1,1,k,itm/); edges=(/imt,jmt,1,1/)
          read(io) buf
          iret=nf_inq_varid(ncid,'u',varid)
          call  nf_put_var(ncid,varid ,corner, edges,buf)
          read(io) buf
          iret=nf_inq_varid(ncid,'v',varid)
          call  nf_put_var(ncid,varid ,corner, edges,buf)
          do nn=1,nt
           read(io) buf
           if (nn==1) then 
              name='temp'
           elseif (nn==2) then  
               name = 'salt'
           else
            write(name,'("tracer_",i2)') nn
            call replace_space_zero(name)
           endif
           iret=nf_inq_varid(ncid,name,varid)
           call  nf_put_var(ncid,varid ,corner, edges,buf)
          enddo
         enddo
         deallocate( buf)
        elseif (ident(1:4) == 'potv') then
         allocate( buf(imt,jmt) )
         do k=1,km
          read(io) buf
          corner=(/1,1,k,itm/); edges=(/imt,jmt,1,1/)
          iret=nf_inq_varid(ncid,'potv',varid)
          call  nf_put_var(ncid,varid ,corner, edges,buf)
         enddo
         deallocate( buf)
        elseif (ident(1:4) == 'bvor') then
         allocate( buf(imt,jmt) )
         do k=1,km
          read(io) buf
          corner=(/1,1,k,itm/); edges=(/imt,jmt,1,1/)
          iret=nf_inq_varid(ncid,'bvort',varid)
          call  nf_put_var(ncid,varid ,corner, edges,buf)
         enddo
         deallocate( buf)
        elseif (ident(1:4) == 'pv35') then
         allocate( buf(jmt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/jmt,km,1,1/)
         iret=nf_inq_varid(ncid,'potv35W',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'bv35') then
         allocate( buf(jmt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/jmt,km,1,1/)
         iret=nf_inq_varid(ncid,'bvort35W',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'pv0N') then
         allocate( buf(imt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/imt,km,1,1/)
         iret=nf_inq_varid(ncid,'potv0N',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'bv0N') then
         allocate( buf(imt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/imt,km,1,1/)
         iret=nf_inq_varid(ncid,'bpot0N',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 't100') then
         allocate( buf(imt,jmt) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'t100',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'u100') then
         allocate( buf(imt,jmt) )
         read(io) buf ; where (buf /= spval) buf=buf/100.0
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'u100',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'v100') then
         allocate( buf(imt,jmt) )
         read(io) buf; where (buf /= spval) buf=buf/100.0
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'v100',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:3) == 's50') then
         allocate( buf(imt,jmt) )
         read(io) buf; where ( buf /=spval) buf=buf*1000.0+35.0
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'s50',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:3) == 't50') then
         allocate( buf(imt,jmt) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'t50',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:3) == 'u50') then
         allocate( buf(imt,jmt) )
         read(io) buf; where (buf /= spval) buf=buf/100.0
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'u50',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:3) == 'v50') then
         allocate( buf(imt,jmt) )
         read(io) buf; where (buf /= spval) buf=buf/100.0
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'v50',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:2) == 's0') then
         allocate( buf(imt,jmt) )
         read(io) buf; where ( buf /=spval) buf=buf*1000.0+35.0
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'s0',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:2) == 't0') then
         allocate( buf(imt,jmt) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'t0',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:2) == 'u0' .and. ident(1:3) /= 'u0N' ) then
         allocate( buf(imt,jmt) )
         read(io) buf; where (buf /= spval) buf=buf/100.0
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'u0',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:2) == 'v0') then
         allocate( buf(imt,jmt) )
         read(io) buf; where (buf /= spval) buf=buf/100.0
         corner=(/1,1,itm,1/); edges=(/imt,jmt,1,1/)
         iret=nf_inq_varid(ncid,'v0',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:3) == 'u0N') then
         allocate( buf(imt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/imt,km,1,1/)
         iret=nf_inq_varid(ncid,'u0N',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'Av0N') then
         allocate( buf(imt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/imt,km,1,1/)
         iret=nf_inq_varid(ncid,'Av0N',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'T57N') then
         allocate( buf(imt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/imt,km,1,1/)
         iret=nf_inq_varid(ncid,'T57N',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'u35W') then
         allocate( buf(jmt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/jmt,km,1,1/)
         iret=nf_inq_varid(ncid,'u35W',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'v35W') then
         allocate( buf(jmt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/jmt,km,1,1/)
         iret=nf_inq_varid(ncid,'v35W',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'w35W') then
         allocate( buf(jmt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/jmt,km,1,1/)
         iret=nf_inq_varid(ncid,'w35W',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:3) == 'u5W') then
         allocate( buf(jmt,km) )
         read(io) buf
         corner=(/1,1,itm,1/); edges=(/jmt,km,1,1/)
         iret=nf_inq_varid(ncid,'u5W',varid)
         call  nf_put_var(ncid,varid ,corner, edges,buf)
         deallocate( buf)
        elseif (ident(1:4) == 'bats') then
         allocate( buf(km,nt) )
         read(io) buf
         corner=(/1,itm,1,1/); edges=(/km,1,1,1/)
         do k=1,nt
          if (k==1) then
           name='t_bats'
          elseif (k==2) then
           name='s_bats'
          else
           write(name,'("tr",i2,"_bats")') k
           call replace_space_zero(name)
          endif
          iret= nf_inq_varid(ncid,name,varid)
          call  nf_put_var(ncid,varid ,corner, edges,buf(:,k))
         enddo
         deallocate( buf)

        else ! ident == 
         print*,' problem reading variable ',ident,' from file'
        endif

      enddo



      goto 20

 30   close(io)
      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. (optional) number of variables in file'
      end program cv_anna






      subroutine cv_anna_define(io)
      use cv_anna_module
      implicit none
#include "netcdf.inc"

      character (len=80) :: ident
      integer :: ncid,iret,io
      integer :: n,k,nn
      integer :: lat_tdim,lat_udim,lon_tdim,lon_udim
      integer :: lat_tid,lat_uid,lon_tid,lon_uid
      integer :: depth_tdim,depth_wdim,depth_tid,depth_wid
      integer :: itimedim,itimeid,varid
      integer :: dims(4)
      character (len=80) :: name, unit
      integer :: blk 

      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c     dimensions
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, 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)  = Lon_udim
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
      dims(1)  = depth_wdim
      depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,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 = 'Latitude on U grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_uid, '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 = 'Longitude on U grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of T Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of W Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_wid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_wid, '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)


      block_size= 2048

      read(io) ! stamp
      blk=1 ! number of blocks to rewind
      do n=1,nr_var
        read(io) ident
        read(io)
        blk=blk+2
        if     (ident(1:3) == 'dml' ) then
         print*,' found mixed layer depth '
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         varid=ncvdef (ncid,'dml', NCFLOAT,3,dims,iret)
         name = 'Mixed layer depth'; unit = 'm'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
         block_size=block_size+imt*jmt
        elseif (ident(1:4) == 'uvts') then
         print*,' found complete u,v,t,s fields '
         dims=(/lon_udim,lat_udim,depth_tdim,itimedim/)
         varid=ncvdef (ncid,'u', NCFLOAT,4,dims,iret)
         name = 'Zonal velocity'; unit = 'm/s'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
         varid=ncvdef (ncid,'v', NCFLOAT,4,dims,iret)
         name = 'Meridional velocity'; unit = 'm/s'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
         dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
         varid=ncvdef (ncid,'temp', NCFLOAT,4,dims,iret)
         name = 'Temperature'; unit = 'deg C'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
         varid=ncvdef (ncid,'salt', NCFLOAT,4,dims,iret)
         name = 'Salinity'; unit = 'psu'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
         if (nt>2) then
          do nn=3,nt
           write(name,'("tracer_",i2)') nn
           call replace_space_zero(name)
           varid=ncvdef (ncid,name, NCFLOAT,4,dims,iret)
           name = 'Passive tracer'; unit = ' '
           call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
          enddo
         endif
         block_size=block_size+imt*jmt*km*(2+nt)

         do k=1,km
          read(io); 
          if (k>1) read(io)
          do nn=1,nt
           read(io)
          enddo
          blk=blk+nt+1
          if (k>1) blk=blk+1
         enddo

        elseif (ident(1:4) == 'potv') then
         print*,' found potential vorticity'
         block_size=block_size+imt*jmt*km
         do k=2,km
          blk=blk+1
          read(io)
         enddo
         dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
         varid=ncvdef (ncid,'potv', NCFLOAT,4,dims,iret)
         name = 'Potential vorticity'; unit = '1/m/s^2'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)

        elseif (ident(1:4) == 'bvor') then
         print*,' found barotropic vorticity'
         block_size=block_size+imt*jmt*km
         do k=2,km
          blk=blk+1
          read(io)
         enddo
         dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
         varid=ncvdef (ncid,'bvort', NCFLOAT,4,dims,iret)
         name = 'Barotropic vorticity'; unit = '1/m/s^2'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'pv35') then
         block_size=block_size+jmt*km
         print*,' found potential vorticity at 35W'
         dims=(/lat_udim,depth_tdim,itimedim,1/)
         varid=ncvdef (ncid,'potv35W', NCFLOAT,3,dims,iret)
         name = 'Potential Vorticity at 35W'; unit = '1/m/s^2'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'bv35') then
         block_size=block_size+jmt*km
         print*,' found barotropic vorticity at 35W'
         dims=(/lat_udim,depth_tdim,itimedim,1/)
         varid=ncvdef (ncid,'bvort35W', NCFLOAT,3,dims,iret)
         name = 'Potential Vorticity at 35W'; unit = '1/m/s^2'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'pv0N') then
         block_size=block_size+imt*km
         print*,' found potential vorticity at 0N'
         dims=(/lon_tdim,depth_tdim,itimedim,1/)
         varid=ncvdef (ncid,'potv0N', NCFLOAT,3,dims,iret)
         name = 'Potential vorticity at Eq.'; unit = '1/m/s^2'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'bv0N') then
         block_size=block_size+imt*km
         print*,' found barotropic vorticity at 0N'
         dims=(/lon_tdim,depth_tdim,itimedim,1/)
         varid=ncvdef (ncid,'bpot0N', NCFLOAT,3,dims,iret)
         name = 'Barotropic vorticity at Eq.'; unit = '1/m/s^2'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 't100') then
         block_size=block_size+imt*jmt
         print*,' found temperature in 100m depth '
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         varid=ncvdef (ncid,'t100', NCFLOAT,3,dims,iret)
         name = 'Temperature in 100m'; unit = 'deg C'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'u100') then
         block_size=block_size+imt*jmt
         print*,' found zonal velocity in 100m depth '
         dims=(/lon_udim,lat_udim,itimedim,1/)
         varid =ncvdef (ncid,'u100', NCFLOAT,3,dims,iret)
         name = 'Zonal velocity in 100m'; unit = 'm/s'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &                  unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'v100') then
         print*,' found meridional velocity in 100m depth '
         block_size=block_size+imt*jmt
         dims=(/lon_udim,lat_udim,itimedim,1/)
         varid =ncvdef (ncid,'v100', NCFLOAT,3,dims,iret)
         name = 'Meridional velocity in 100m'; unit = 'm/s'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &                  unit,len_trim(unit),spval)
        elseif (ident(1:3) == 's50') then
         print*,' found salinity in 50m depth '
         block_size=block_size+imt*jmt
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         varid =ncvdef (ncid,'s50', NCFLOAT,3,dims,iret)
         name = 'Salinity in 50m'; unit = 'psu'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &                 unit,len_trim(unit),spval)
        elseif (ident(1:3) == 't50') then
         print*,' found temperature in 50m depth '
         block_size=block_size+imt*jmt
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         varid =ncvdef (ncid,'t50', NCFLOAT,3,dims,iret)
         name = 'Temperature in 50m'; unit = 'deg C'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:3) == 'u50') then
         print*,' found zonal velocity in 50m depth '
         block_size=block_size+imt*jmt
         dims=(/lon_udim,lat_udim,itimedim,1/)
         varid =ncvdef (ncid,'u50', NCFLOAT,3,dims,iret)
         name = 'Zonal velocity in 50m'; unit = 'm/s'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:3) == 'v50') then
         print*,' found meridional velocity in 50m depth '
         block_size=block_size+imt*jmt
         dims=(/lon_udim,lat_udim,itimedim,1/)
         varid =ncvdef (ncid,'v50', NCFLOAT,3,dims,iret)
         name = 'Meridional velocity in 50m'; unit = 'm/s'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:2) == 's0') then
         print*,' found salinity in 0m depth '
         block_size=block_size+imt*jmt
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         varid =ncvdef (ncid,'s0', NCFLOAT,3,dims,iret)
         name = 'Salinity in 0m'; unit = 'psu'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:2) == 't0') then
         print*,' found temperature in 0m depth '
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         varid =ncvdef (ncid,'t0', NCFLOAT,3,dims,iret)
         name = 'Temperature in 0m'; unit = 'deg C'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:2) == 'u0' .and. ident(1:3) /= 'u0N' ) then
         print*,' found zonal velocity in 0m depth '
         block_size=block_size+imt*jmt
         dims=(/lon_udim,lat_udim,itimedim,1/)
         varid =ncvdef (ncid,'u0', NCFLOAT,3,dims,iret)
         name = 'Zonal velocity in 0m'; unit = 'm/s'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:2) == 'v0') then
         print*,' found meridional velocity in 0m depth '
         block_size=block_size+imt*jmt
         dims=(/lon_udim,lat_udim,itimedim,1/)
         varid =ncvdef (ncid,'v0', NCFLOAT,3,dims,iret)
         name = 'Meridional velocity in 0m'; unit = 'm/s'
         call dvcdf_here(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:3) == 'u0N')  then
         print*,' found zonal vel. at eq.'
         block_size=block_size+imt*km
         dims=(/lon_udim,depth_tdim,itimedim,1/)
         varid=ncvdef (ncid,'u0N', NCFLOAT,3,dims,iret)
         name = 'Zonal velocity at Eq.'; unit = 'm/s'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'Av0N')  then
         print*,' found vertical viscosity at eq.'
         block_size=block_size+imt*km
         dims=(/lon_udim,depth_wdim,itimedim,1/)
         varid=ncvdef (ncid,'Av0N', NCFLOAT,3,dims,iret)
         name = 'Vertical viscosity at Eq.'; unit = 'm/s'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'T57N')  then
         print*,' found temperature at 57N'
         block_size=block_size+imt*km
         dims=(/lon_tdim,depth_tdim,itimedim,1/)
         varid=ncvdef (ncid,'T57N', NCFLOAT,3,dims,iret)
         name = 'Temperature at 57N'; unit = 'deg C'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'u35W') then
         print*,' found zonal vel. at 35W'
         block_size=block_size+jmt*km
         dims=(/lat_udim,depth_tdim,itimedim,1/)
         varid=ncvdef (ncid,'u35W', NCFLOAT,3,dims,iret)
         name = 'Zonal velocity at 35W'; unit = 'm/s'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'v35W') then
         print*,' found merid vel. at 35W'
         block_size=block_size+jmt*km
         dims=(/lat_udim,depth_tdim,itimedim,1/)
         varid=ncvdef (ncid,'v35W', NCFLOAT,3,dims,iret)
         name = 'Merid. velocity at 35W'; unit = 'm/s'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'w35W') then
         print*,' found vert. vel. at 35W'
         block_size=block_size+jmt*km
         dims=(/lat_udim,depth_wdim,itimedim,1/)
         varid=ncvdef (ncid,'w35W', NCFLOAT,3,dims,iret)
         name = 'Vertical velocity at 35W'; unit = 'm/s'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:3) == 'u5W') then
         print*,' found zonal vel. at 5W'
         block_size=block_size+jmt*km
         dims=(/lat_udim,depth_tdim,itimedim,1/)
         varid=ncvdef (ncid,'u5W', NCFLOAT,3,dims,iret)
         name = 'Zonal velocity at 5W'; unit = 'm/s'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
        elseif (ident(1:4) == 'bats') then
         print*,' found station BATS'
         block_size=block_size+km*nt
         dims=(/depth_tdim,itimedim,1,1/)
         varid=ncvdef (ncid,'t_bats', NCFLOAT,2,dims,iret)
         name = 'Temperature at 65W, 32N'; unit = 'deg C'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &                    unit,len_trim(unit),spval)
         varid=ncvdef (ncid,'s_bats', NCFLOAT,2,dims,iret)
         name = 'Salinity at 65W, 32N'; unit = 'psu'
         call dvcdf_here(ncid,varid,name,len_trim(name),
     &                    unit,len_trim(unit),spval)
         do k=3,nt
          write(name,'("tr",i2,"_bats")') k
          call replace_space_zero(name)
          varid=ncvdef (ncid,name, NCFLOAT,2,dims,iret)
          write(name,'("Passive tracer ",i2)') k; unit = ' '
          call dvcdf_here(ncid,varid,name,len_trim(name),
     &                    unit,len_trim(unit),spval)
         enddo
        else
          print*,' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          print*,' Warning: do not know how to read ',ident,' from file'
          print*,' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        endif

      enddo

      call ncendf(ncid, iret)

      call  nf_put_var(ncid,lat_tid ,1, jmt,yt)
      call  nf_put_var(ncid,lat_uid ,1, jmt,yu)
      call  nf_put_var(ncid,lon_tid ,1, imt,xt)
      call  nf_put_var(ncid,lon_uid ,1, imt,xu)
      call  nf_put_var(ncid,depth_tid ,1, km,zt)
      call  nf_put_var(ncid,depth_wid ,1, km,zw)

      call ncclos (ncid, iret)


      do k=1,blk
         backspace(io) 
      enddo

      end subroutine cv_anna_define








      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


      subroutine nf_put_var(ncid,varid,corner,edges,dd)
      use kind_mod
      implicit none
#include "netcdf.inc"
      integer :: ncid,varid,corner,edges,iret
      real (kind=r_knd) :: dd
#ifdef real_8
      iret= nf_put_vara_double(ncid,varid,corner,edges,dd)
#else
      iret= nf_put_vara_real(ncid,varid,corner,edges,dd)
#endif
      end subroutine nf_put_var


      subroutine get_file_size(size)
c-----------------------------------------------------------------------
c     returns the file size in bytes
c-----------------------------------------------------------------------
      implicit none
      integer :: size
      
      end subroutine get_file_size

