#include "options.inc"

c#define convert_restart
c
c
#ifdef convert_restart
c
c     convert binary restart file of averaged blue correction
c     (annual mean version) to NetCDF format
c
#else
c
c     convert binary file containing averaged blue correction
c     (annual mean version) written during model run to NetCDF format
c
#endif
c
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


      subroutine barrier
      end

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


      program cv_blue_mean
      use kind_mod
      implicit none

      real (kind=r_knd), allocatable, dimension(:,:,:) :: corr
      real (kind=r_knd), allocatable, dimension(:)   :: xt,yt,zt
      real (kind=r_knd) :: dd
      real (kind=r_knd) :: spval=-9.9e12

      integer :: imt,jmt,km
      integer :: counter
      integer, allocatable, dimension(:,:)     :: kmt
      integer :: i,j,k,io
      character (len=80) :: infile,outfile

#include "netcdf.inc"
      character (len=80) :: name,unit
      integer ncid,iret,dims(4),start(4),count(4)
      integer :: lon_tdim,lon_tid
      integer :: lat_tdim,lat_tid
      integer :: depth_tdim,depth_tid
      integer :: itimedim,itimeid
      integer :: blueid,n,len
      integer :: year,month,day,domain

      character (len=15) :: c15
      character (len=1) :: c1
      character (len=2) :: c2
      character (len=4) :: c4


#ifdef convert_restart
      print*,' Converting blue_mean restart written by SPFLAME ',
     &       ' to NetCDF file '
      call get_commandline_arg(1,infile,iret)
      if (iret/=0) goto 99
      print*,' reading from file ',infile(1:len_trim(infile))
      call get_commandline_arg(2,outfile,iret)
      if (iret/=0) goto 99
      print*,' writing to file ',outfile(1:len_trim(outfile))
#else
      print*,' Converting blue_mean correction written by SPFLAME ',
     &       ' to NetCDF file '
      call get_commandline_arg(1,infile,iret)
      if (iret/=0) goto 99
      print*,' reading from file ',infile(1:len_trim(infile))

      read(infile,'(a15,i3,a2,i4,a1,i2,a1,i2,a4)') 
     &            c15,domain,c2,year,c1,month,c1,day,c4

      write(outfile,
     &   '("blue_mean_corr_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            domain,year,month,day
      do i=1,len_trim(outfile);
        if (outfile(i:i)==' ')outfile(i:i)='0';
      enddo
      print*,' sub domain = ',domain
      print*,' year/month/day  = ',year,month,day
      print*,' writing to file ',outfile

#endif

#ifdef convert_restart
c     read in the grid from a binary file
      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

      print*,' found dimensions: '
      print*,' imt=',imt,' jmt=',jmt,' km=',km

      allocate( xt(imt), yt(jmt), zt(km) )
      
      read (io) (dd,i=1,imt)
     &,         (dd,j=1,jmt)
     &,         (dd,i=1,imt)
     &,         (dd,j=1,jmt)
     &,         (dd,k=1,km)
     &,         (dd,k=0,km)
     &,         (xt(i),i=1,imt)
     &,         (dd,i=1,imt)
     &,         (yt(j),j=1,jmt)
     &,         (dd,j=1,jmt)
     &,         (zt(k),k=1,km)
     &,         (dd,k=1,km)

      close (io)

      allocate( kmt(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
c      check i,j,k
      if (i/=imt .or. j /= jmt .or. k/= km) then
         print*,' Error: '
         print*,' read imt,jmt,km=',i,j,k
         print*,' but previously ',imt,jmt,km
        call halt_stop(' reading kmt file ')
       endif
      read (io) 
      read (io) kmt
      close(io)

      allocate( corr( imt,jmt,km) )

      call getunit (io, infile, 'usr ieee')
      read (io) counter
      do k=1,km
        read(io) corr(:,:,k),corr(:,:,k)
      enddo
      close(io)
#else
c
c     all contained in one file
c
      call getunit (io, infile, 'usr ieee')
      read(io) imt,jmt,km
      print*,' imt=',imt,' jmt=',jmt,' km=',km
      allocate( xt(imt),yt(jmt),zt(km))
      allocate( kmt(imt,jmt) )
      allocate( corr(imt,jmt,km) )
      read(io) xt,yt,zt
      read(io) counter
      print*,' counter for averages=',counter
      read(io) kmt
      do k=1,km
         read(io) corr(:,:,k)
      enddo
      close(io)

#endif


      dd=365./2. 

      print*,' writing to file ',outfile

      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c       dimensions
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', 1, 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(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c     attributes of the grid
      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 = '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 = '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,
     &'31-DEC-1899 00:00:00', iret)
c     variables
      dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
      blueid=ncvdef (ncid,'blue_corr', NCFLOAT,4,dims,iret)
c     attributes of variables
      name = 'Correction to momentum equation'; unit = 'kg/m^3'
      call dvcdf_here(ncid,blueid,name,24,unit,16,spval)
      call ncendf(ncid, iret)

      zt=zt/100.

#ifdef real_8
      iret= nf_put_vara_double(ncid,lon_tid,1,imt,xt)
      iret= nf_put_vara_double(ncid,lat_tid,1,jmt,yt)
      iret= nf_put_vara_double(ncid,depth_tid,1,km,zt)
      iret= nf_put_vara_double(ncid,itimeid,1,1,dd)
#else
      call ncvpt(ncid, Lon_tid, 1, imt,xt, iret)
      call ncvpt(ncid, Lat_tid, 1, jmt,yt, iret)
      call ncvpt(ncid, depth_tid, 1, km,zt, iret)
      call ncvpt(ncid, itimeid, 1, 1,dd, iret)
#endif

c     convert to SI units 
      corr=corr*1000.   ! g/cm^3 -> kg/m^3

c     apply land mask and write to file
      do k=1,km
       where( kmt <k ) corr(:,:,k) = spval
       start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
       iret= nf_put_vara_double(ncid,blueid,start,count,corr(:,:,k))
#else
       iret= nf_put_vara_real(ncid,blueid,start,count,corr(:,:,k))
#endif
      enddo
      call ncclos (ncid, iret)

      stop

 99   print*,' program needs command line parameter'
#ifdef convert_restart
      print*,' 1. input binary file (blue_mean restart file)'
      print*,' 1. output NetCDF file '
#else
      print*,' 1. input binary file'
#endif

      end program cv_blue_mean



      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
      real (kind=4):: 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




