#include "options.inc"

c
c     convert unformatted written data containing
c     overturning to NetCDF format
c
c
c      or:
c
c     read meridional velocity from snapshot or
c     averaged file to compute overturning
#ifdef partial_cell
c     from an experiment with partial cells
#endif
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_vsf_module
      use kind_mod
      use time_manager_module
      integer :: imt,jmt,km,itm
      real (kind=r_knd),allocatable :: xt(:), yt(:), xu(:), yu(:),
     & dxtdeg(:), dytdeg(:), dxudeg(:), dyudeg(:),
     & dxt(:), dyt(:), dxu(:), dyu(:),
     & zt(:),zw(:),dzt(:),dzw(:),time(:)
      integer, dimension(:,:), allocatable :: kmt, kmu
      real (kind=r_knd), dimension(:,:), allocatable :: htp, hup
      real (kind=r_knd), dimension(:,:,:), allocatable :: dht, dhu
      real (kind=r_knd), dimension(:,:,:), allocatable :: vsf
      real (kind=r_knd), parameter :: radius   = 6370.0e5
      real (kind=r_knd), parameter :: pi       = 
     &        3.14159265358979323846264338327950588
      real (kind=r_knd), parameter :: radian   = 360./(2.*pi)
      real (kind=r_knd), parameter :: degtcm   = radius/radian
      real (kind=r_knd) :: spval = -9.9e15
      character (len=32) :: istamp ='m/d/y= 1/ 1/1900, h:m:s= 0: 0: 0' 
      end module cv_vsf_module




      program cv_vsf
      use cv_vsf_module
      implicit none
      character (len=80) :: infile,outfile
      integer :: ierr
      print*,' Converting overturning written by SPFLAME ',
     &       ' to NetCDF file '
      call set_calendar_type(no_leap)
      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 read_grid
c      call read_vel(infile)
      call read_vsf(infile)
      call write_vsf(outfile)
      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_vsf


      subroutine read_vsf(infile)
c
c-----------------------------------------------------------------------
c     reads meridional overturning vsf
c-----------------------------------------------------------------------
c
      use cv_vsf_module
      implicit none
      character (len=*) :: infile
      character (len=32) :: stamp
      type( time_type ) :: t
      integer :: io,k,n,j,i
      real  (kind=r_knd):: x

      call getunit (io, infile,'u s r ieee')
      itm=0
 10   continue
      do k=1,8; read (io,end=20); enddo
      itm=itm+1
      goto 10
 20   continue
      print*,' found ',itm,' time steps in file ',infile
      rewind(io)
      allocate( vsf(jmt,km,itm), time(itm) )
      vsf=spval
      do n=1,itm
       read (io) stamp ! iotext, expnam
       t =  get_stamp (stamp) - get_stamp(istamp)
       time(n)=t%days + t%seconds/60./60./24.
       print*,' n=',n,' days since origin=',time(n),' date ',stamp
       read (io) j, k
       if ((j/=jmt) .or. (k/=km)) then
        print*,' found jmt=',j,' km=',k,' in file ',infile
        print*,' previous definition was jmt=',jmt,' km=',km
        stop
       endif
       read (io) ! stamp ! iotext, expnam
       read (io) ! zw
       read (io) !current_stamp, iotext, expnam
       read (io) !yu
       read (io) !current_stamp, iotext, expnam
       read (io) vsf(:,:,n)

       do j=1,jmt
        do k=1,km
         x=0
         do i=2,imt-1
          if (kmu(i,j)>=k) x=1
         enddo
         vsf(j,k,n)=vsf(j,k,n)/1.e12*x+spval*(1-x)
        enddo
       enddo
      enddo
      close(io)

      end subroutine read_vsf


      subroutine read_vel(infile)
c
c-----------------------------------------------------------------------
c     reads meridional velocity and computes vsf
c-----------------------------------------------------------------------
c
      use cv_vsf_module
      implicit none
      character (len=*) :: infile
#include "netcdf.inc"
      integer :: iret,ncid,vid,dim(4),timeid,sshid
      integer corner(4), edges(4)
      character (len=80) :: s
      real (kind=r_knd), allocatable :: vbarx(:),ssh(:,:,:)
      real (kind=r_knd), dimension(:,:,:,:), allocatable :: v
      real (kind=r_knd) :: x
      integer :: k,n,i,j
      logical :: freesurf=.false.

      print*,' reading meridional velocity from file ',
     &     infile(1:len_trim(infile))
      iret=nf_open(infile,NF_NOWRITE,ncid)
      iret=nf_inq_varid(ncid,'v',vid)
      iret=NF_INQ_VARDIMID(ncid,vid,dim)
      iret= NF_INQ_DIMLEN (ncid, DIM(4), itm)
      s='';iret=NF_INQ_DIMNAME(ncid,dim(4), s) 
      iret=nf_inq_varid(ncid,s,timeid)
      print*,' time steps in file: ',itm
      allocate( v(imt,jmt,km,itm), time(itm) );v=0.;time=0.
      iret=nf_inq_varid(ncid,'ssh',sshid)
      if (iret == 0) then
        freesurf=.true.
        print*,' found free surface in snapshot'
        print*,' assuming that the exp. used free surface '
        allocate(ssh(imt,jmt,itm))
        corner=1; edges=(/imt,jmt,itm,1/)
#ifdef real_8
        iret= nf_get_vara_double(ncid,sshid, corner,edges,ssh)
        iret = nf_get_att_double(ncid, sshid, 'missing_value', spval)
#else
        iret= nf_get_vara_real (ncid,sshid, corner,edges,ssh)
        iret = nf_get_att_real(ncid, sshid, 'missing_value', spval)
#endif
        where( ssh/=spval) ssh=ssh*100. ! convert to cm
        where( ssh==spval) ssh=0.
      endif

      corner=1; edges=(/itm,1,1,1/)
#ifdef real_8
      iret= nf_get_vara_double(ncid,timeid, corner,edges,time)
#else
      iret= nf_get_vara_real (ncid,timeid, corner,edges,time)
#endif
      corner=1; edges=(/imt,jmt,km,itm/)
#ifdef real_8
      iret= nf_get_vara_double(ncid,vid, corner,edges,v)
      iret = nf_get_att_double(ncid, vid, 'missing_value', spval)
#else
      iret= nf_get_vara_real (ncid,vid, corner,edges,v)
      iret = nf_get_att_real(ncid, vid, 'missing_value', spval)
#endif
      where( v/=spval) v=v*100 ! convert to cm/s
      where( v==spval) v=0.
      call ncclos (ncid, iret)

      allocate( vbarx(km) );vbarx=0.
      allocate( vsf(jmt,km,itm) )
      vsf=spval

      do n=1,itm
      do j=2,jmt-1
       do k=1,km
        vbarx=0.
        do i=2,imt-1
         vbarx(k) = vbarx(k) + v(i,j,k,n)*dxu(i)
#ifdef partial_cell
     &                                          *dhu(i,k,j)
#else
     &                                          *dzt(k)
#endif
     &        *cos(yu(j)/radian)
        enddo
        if (k == 1) then
         if (freesurf) then
          do i=2,imt-1
           vbarx(k) = vbarx(k) + v(i,j,k,n)*dxu(i)*
     &                   0.5*((ssh(i,j  ,n)+ssh(i-1,j  ,n)*0.5)
     &                       +(ssh(i,j-1,n)+ssh(i-1,j-1,n)*0.5)) 
          enddo
         endif
         vsf(j,k,n) = vbarx(k)
	else
         vsf(j,k,n) = vsf(j,k-1,n) + vbarx(k)
	endif
       enddo
       do k=1,km
        x=0
        do i=2,imt-1
         if (kmu(i,j)>=k) x=1
        enddo
        vsf(j,k,n)=vsf(j,k,n)/1.e12*x+spval*(1-x)
       enddo
      enddo
      enddo
      deallocate(vbarx,v)
      end subroutine read_vel


      subroutine write_vsf(outfile)
      use cv_vsf_module
      implicit none
      character (len=*) :: outfile
      real  (kind=r_knd):: x
      integer :: i,j,k,n
#include "netcdf.inc"
      integer :: iret,ncid,lat_udim,lat_uid
      integer :: depth_tdim,depth_tid,itimedim,itimeid
      integer :: dims(4),vsfid
      character name*24, unit*16, text*80
      integer corner(4), edges(4)

      print*,' writing meridional overturning to file ',
     &     outfile(1:len_trim(outfile))
      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', nf_unlimited, iret)
c     grid variables
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
      dims=(/lat_udim,depth_tdim,itimedim,1/)
      vsfid=ncvdef (ncid,'vsf', NCFLOAT,3,dims,iret)
c     attributes of the grid
      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 = '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 = '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,
     &  '01-JAN-1900 00:00:00', iret)
c     attributes of variables
      name = 'Meridional overturning'; unit = 'Sv'
      call dvcdf_here(ncid,vsfid,name,24,unit,16,spval)
      call ncendf(ncid, iret)
      zt=zt/100.
      corner=(/1,1,1,1/); edges=(/jmt,km,itm,1/)
#ifdef real_8
      iret= nf_put_vara_double(ncid,lat_uid ,1, jmt,yu)
      iret= nf_put_vara_double(ncid,depth_tid ,1, km,zt)
      iret= nf_put_vara_double(ncid,itimeid ,1, itm,time)
      iret= nf_put_vara_double(ncid,vsfid ,corner, edges,vsf)
#else
      iret= nf_put_vara_real(ncid,lat_uid ,1, jmt,yu)
      iret= nf_put_vara_real(ncid,depth_tid ,1, km,zt)
      iret= nf_put_vara_real(ncid,itimeid ,1, itm,time)
      iret= nf_put_vara_real(ncid,vsfid ,corner, edges,vsf)
#endif
      call ncclos (ncid, iret)

      end subroutine write_vsf
 


      subroutine read_grid ()
c
c-----------------------------------------------------------------------
c     read in the grid from a binary file
c-----------------------------------------------------------------------
c
      use cv_vsf_module
      implicit none
      integer i,j,k,io,kb

      print*,
     & ' =>Reading the grid definition from file grid.dta'
      call getunit (io, 'grid.dta', 'u s r ieee')
      read (io) ! no checks for awi id
      read (io) imt, jmt, km
      allocate( xt(imt), yt(jmt), xu(imt), yu(jmt))
      allocate( dxtdeg(imt), dytdeg(jmt), dxudeg(imt), dyudeg(jmt))
      allocate( dxt(imt), dyt(jmt), dxu(imt), dyu(jmt))
      allocate( zt(km),zw(km),dzt(km),dzw(0:km) )
      
      read (io) (dxtdeg(i),i=1,imt)
     &,         (dytdeg(j),j=1,jmt)
     &,         (dxudeg(i),i=1,imt)
     &,         (dyudeg(j),j=1,jmt)
     &,         (dzt(k),k=1,km)
     &,         (dzw(k),k=0,km)
     &,         (xt(i),i=1,imt)
     &,         (xu(i),i=1,imt)
     &,         (yt(j),j=1,jmt)
     &,         (yu(j),j=1,jmt)
     &,         (zt(k),k=1,km)
     &,         (zw(k),k=1,km)

      do j=1,jmt
        dyt(j) = dytdeg(j)*degtcm
        dyu(j) = dyudeg(j)*degtcm
      enddo
      do i=1,imt
        dxt(i) = dxtdeg(i)*degtcm
        dxu(i) = dxudeg(i)*degtcm
      enddo

      close (io)

      allocate( kmt(imt,jmt), kmu(imt,jmt) )
      print*,
     & ' =>Reading a "kmt" field from binary file kmt.dta'
      call getunit (io, 'kmt.dta', 'usr ieee')
      read (io) ! no checks
c       skip header record and read dimensional info
      read (io) 
      read (io) i, j, k
      read (io) 
      read (io) kmt
      close(io)

      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

#ifdef partial_cell
      allocate( htp(imt,jmt), hup(imt,jmt) )
      allocate( dht(imt,km,jmt), dhu(imt,km,jmt) )

      print*,
     & ' =>Reading a "htp" field from binary file htp.dta'
      call getunit (io, 'htp.dta', 'usr ieee')
      read (io) ! no checks
c       skip header record and read dimensional info
      read (io) 
      read (io) i, j, k
      read (io) 
      read (io) htp
      close(io)

      hup=0.
      do j=1,jmt-1
        do i=1,imt-1
          hup(i,j) = min( htp(i,j),htp(i+1,j),htp(i,j+1),htp(i+1,j+1))
	enddo
      enddo

      do k=1,km
       dht(:,k,:)=dzt(k)
       dhu(:,k,:)=dzt(k)
      enddo

      do j=1,jmt
	do i=1,imt
	  kb = kmt(i,j)
	  if (kb .gt. 1) then
	    dht(i,kb,j) = htp(i,j) - zw(kb-1)
	  endif
	enddo
      enddo
c

      do j=1,jmt-1
        do k=1,km
          do i=1,imt-1
	    dhu(i,k,j) = min(dht(i,k,j), dht(i+1,k,j), dht(i,k,j+1)
     &,                      dht(i+1,k,j+1))
	  enddo
          dhu(1,k,j)   = dht(1,k,j)
          dhu(imt,k,j) = dht(imt,k,j)
	enddo
      enddo
#endif

      end subroutine read_grid



      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
