#include "options.inc"

c
c     convert unformatted written data containing
c     meridional heat transport 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_gyre_module
      use kind_mod
      use time_manager_module
      integer :: imt,jmt,km,itm,nt
      real (kind=r_knd) ,allocatable :: xt(:), yt(:), xu(:), yu(:),
     & dxtdeg(:), dytdeg(:), dxudeg(:), dyudeg(:),
     & dxt(:), dyt(:), dxu(:), dyu(:),
     & zt(:),zw(:),dzt(:),dzw(:)
      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
      real (kind=r_knd), allocatable :: time(:)
      real (kind=r_knd),dimension(:,:,:),allocatable::over_tr,gyre_tr
      real (kind=r_knd),dimension(:,:,:),allocatable::adv_tr
      real (kind=r_knd),dimension(:,:,:),allocatable::diff_tr,barotr_tr
      real (kind=r_knd),dimension(:,:,:),allocatable::total_tr,ekm_tr
      real (kind=r_knd),dimension(:,:,:),allocatable::barocl_tr
      character (len=32) :: istamp ='m/d/y= 1/ 1/1900, h:m:s= 0: 0: 0' 
      end module cv_gyre_module





      program cv_gyre
      use cv_gyre_module
      implicit none
      character (len=80) :: infile,outfile
      integer :: ierr
      print*,' Converting heat transport 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))
      spval=-9.9e12
      call read_grid
      call read_gyre(infile)
      call write_gyre(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_gyre



      subroutine read_gyre(infile)
      use cv_gyre_module
      implicit none
      character (len=*) :: infile
      character (len=32) :: stamp
      type( time_type ) :: t
      integer :: io,k,n,j,i
      real (kind=r_knd), allocatable :: buf(:,:,:)
      call getunit (io, infile,'u s r ieee')
      read (io)
      read (io) j,nt
      if (j/=jmt)  then
        print*,' found jmt=',j,' km=',k,' in file ',infile
        print*,' previous definition was jmt=',jmt,' km=',km
        stop
      endif
      rewind(io)
      itm=0
 10   continue
      do k=1,6; read (io,end=20); enddo
      itm=itm+1
      goto 10
 20   continue
      rewind(io)
      print*,' found ',itm,' time steps in file ',infile
      print*,' number of tracers  ',nt
      allocate(buf(8,jmt,nt), time(itm) )
      allocate(over_tr(jmt,nt,itm) )
      allocate(gyre_tr(jmt,nt,itm) )
      allocate(adv_tr(jmt,nt,itm) )
      allocate(diff_tr(jmt,nt,itm) )
      allocate(barotr_tr(jmt,nt,itm) )
      allocate(total_tr(jmt,nt,itm) )
      allocate(ekm_tr(jmt,nt,itm) )
      allocate(barocl_tr(jmt,nt,itm) )
      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) ! jmt,nt,retim
       read (io) ! stamp ! iotext, expnam
       read (io) buf
       read (io) !current_stamp, iotext, expnam
       read (io) ! buf, dummy
       do j=1,jmt
       over_tr(j,:,n)  =buf(1,j,:); gyre_tr(j,:,n)  =buf(2,j,:);
       barotr_tr(j,:,n)=buf(3,j,:); barocl_tr(j,:,n)=buf(4,j,:);
       ekm_tr(j,:,n)   =buf(5,j,:); adv_tr(j,:,n)   =buf(6,j,:);
       diff_tr(j,:,n)  =buf(7,j,:); total_tr(j,:,n) =buf(8,j,:);
       enddo
      enddo
      close(io)
      deallocate( buf)
      end subroutine read_gyre



      subroutine write_gyre(outfile)
      use cv_gyre_module
      implicit none
      character (len=*) :: outfile
#include "netcdf.inc"
      integer ncid,iret,i,j,k
      integer lat_udim,itimedim,lat_uid,itimeid
      integer n,id, advid, diffid, overid, gyreid, ekmid
      integer btrid,bclid
      integer dims(4), corner(4), edges(4)
      character name*24, unit*16,lname*24
      real (kind=r_knd)::factor, small=1.e-10
      real (kind=r_knd)::pwatts = 4.186e-15, csalt  = 1.e-10

      print*,' writing to NetCDF file ',outfile
      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c     dimensions
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, 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=(/lat_udim,itimedim,1,1/)
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 = '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
      do n=1,nt
      name='heat';unit='PW'
      if (n==2) name='salt'
      if (n==2) unit='10^10cm^2/s'
      if (n>2)  write(name,'("tracer_",i2)') n
      if (n>2) unit=' '
      do i=1,len_trim(name); if (name(i:i)==' ')name(i:i)='0'
      enddo
      k=len_trim(name)
      id=ncvdef (ncid,name(1:k)//'_tr', NCFLOAT,2,dims,iret)
      advid =ncvdef (ncid,name(1:k)//'_tr_adv', NCFLOAT,2,dims,iret)
      diffid=ncvdef (ncid,name(1:k)//'_tr_diff', NCFLOAT,2,dims,iret)
      overid=ncvdef (ncid,name(1:k)//'_tr_over', NCFLOAT,2,dims,iret)
      gyreid=ncvdef (ncid,name(1:k)//'_tr_gyre', NCFLOAT,2,dims,iret)
      ekmid =ncvdef (ncid,name(1:k)//'_tr_ekm', NCFLOAT,2,dims,iret)
      btrid =ncvdef (ncid,name(1:k)//'_tr_btr', NCFLOAT,2,dims,iret)
      bclid =ncvdef (ncid,name(1:k)//'_tr_bcl', NCFLOAT,2,dims,iret)
      lname = 'Northward '//name(1:k)//' transport'; 
      call dvcdf_here(ncid,id,lname,24,unit,16,spval)
      lname = 'Advective '//name(1:k)//' transport'; 
      call dvcdf_here(ncid,advid,lname,24,unit,16,spval)
      lname = 'Diffusive '//name(1:k)//' transport'; 
      call dvcdf_here(ncid,diffid,lname,24,unit,16,spval)
      lname = 'Overturning '//name(1:k)//' transport'; 
      call dvcdf_here(ncid,overid,lname,24,unit,16,spval)
      lname = 'Gyre '//name(1:k)//' transport';
      call dvcdf_here(ncid,gyreid,lname,24,unit,16,spval)
      lname = 'Ekman '//name(1:k)//' transport'; 
      call dvcdf_here(ncid,ekmid,lname,24,unit,16,spval)
      lname = 'Barotropic '//name(1:k)//' transport'; 
      call dvcdf_here(ncid,btrid,lname,24,unit,16,spval)
      lname = 'Baroclinic '//name(1:k)//' transport'; 
      call dvcdf_here(ncid,bclid,lname,24,unit,16,spval)
      enddo ! n
      call ncendf(ncid, iret)
#ifdef real_8
      iret= nf_put_vara_double(ncid,lat_uid ,1, jmt,yu)
      iret= nf_put_vara_double(ncid,itimeid ,1, itm,time)
#else
      iret= nf_put_vara_real(ncid,lat_uid ,1, jmt,yu)
      iret= nf_put_vara_real(ncid,itimeid ,1, itm,time)
#endif

      do n=1,nt
        name='heat'
        if (n==2) name='salt'
        if (n>2)  write(name,'("tracer_",i2)') n
        do i=1,len_trim(name); if (name(i:i)==' ')name(i:i)='0'
        enddo
        k=len_trim(name)
        print*,' Tracer ',name(1:k)
        corner=(/1,1,1,1/); edges=(/jmt,itm,1,1/)


#ifdef real_8
        iret=nf_inq_varid(ncid,name(1:k)//'_tr',id)
        iret= nf_put_vara_double (ncid,id ,corner, edges,
     &                          total_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_adv',id)
        iret= nf_put_vara_double (ncid,id ,corner, edges,
     &                          adv_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_diff',id)
        iret= nf_put_vara_double (ncid,id ,corner, edges,
     &                          diff_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_over',id)
        iret= nf_put_vara_double (ncid,id ,corner, edges,
     &                          over_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_gyre',id)
        iret= nf_put_vara_double (ncid,id ,corner, edges,
     &                          gyre_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_ekm',id)
        iret= nf_put_vara_double (ncid,id ,corner, edges,
     &                          ekm_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_btr',id)
        iret= nf_put_vara_double (ncid,id ,corner, edges,
     &                          barotr_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_bcl',id)
        iret= nf_put_vara_double (ncid,id ,corner, edges,
     &                          barocl_tr(:,n,:))
#else
        iret=nf_inq_varid(ncid,name(1:k)//'_tr',id)
        iret= nf_put_vara_real (ncid,id ,corner, edges,
     &                          total_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_adv',id)
        iret= nf_put_vara_real (ncid,id ,corner, edges,
     &                          adv_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_diff',id)
        iret= nf_put_vara_real (ncid,id ,corner, edges,
     &                          diff_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_over',id)
        iret= nf_put_vara_real (ncid,id ,corner, edges,
     &                          over_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_gyre',id)
        iret= nf_put_vara_real (ncid,id ,corner, edges,
     &                          gyre_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_ekm',id)
        iret= nf_put_vara_real (ncid,id ,corner, edges,
     &                          ekm_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_btr',id)
        iret= nf_put_vara_real (ncid,id ,corner, edges,
     &                          barotr_tr(:,n,:))
        iret=nf_inq_varid(ncid,name(1:k)//'_tr_bcl',id)
        iret= nf_put_vara_real (ncid,id ,corner, edges,
     &                          barocl_tr(:,n,:))
#endif
      enddo
      call ncclos (ncid, iret)
      end subroutine write_gyre






      subroutine read_grid ()
c
c-----------------------------------------------------------------------
c     read in the grid from a binary file
c-----------------------------------------------------------------------
c
      use cv_gyre_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)
      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
