#include "options.inc"


c
c     convert binary output of averaged blue correction
c     to NetCDF format
c
      module cv_blue_module
      integer :: imt,jmt,km

      real (kind=4), allocatable, dimension(:,:,:) :: mean_blue
      integer :: counter
      real (kind=4), allocatable, dimension(:)   :: xt,yt,zt
      integer, allocatable, dimension(:,:)     :: kmt
      end module cv_blue_module

      subroutine barrier
      end

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


      program cv_blue
      use cv_blue_module
      implicit none
      character (len=80) :: infile,outfile
      integer :: ierr,i
      character (len=14) :: c14
      character (len=1) :: c1
      character (len=2) :: c2
      character (len=4) :: c4
      integer :: year,month,day,domain

      print*,' Converting averaged quantities 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))

      read(infile,'(a14,i3,a2,i4,a1,i2,a1,i2,a4)') 
     &            c14,domain,c2,year,c1,month,c1,day,c4
      write(outfile,
     &   '("blue_averages_",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
      call read_blue_averages(infile)
      call write_blue_averages(outfile,year,month,day,domain)
      stop

 99   print*,' program needs command line parameter'
      print*,' 1. input binary file name (written by spflame)'
      end program cv_blue


      subroutine read_blue_averages(infile)
      use cv_blue_module
      implicit none
      character (len=*) :: infile
      real (kind=4), allocatable :: buf(:,:)
      integer :: i,j,k,io,n

      allocate( buf(imt,jmt) ); buf=0.
      io=10
      open(io,file=infile,form='unformatted')
      read(io) imt,jmt,km
      print*,' imt=',imt,' jmt=',jmt,' km=',km

      allocate( mean_blue(imt,km,jmt) )
      allocate( xt(imt),yt(jmt),zt(km))
      allocate( kmt(imt,jmt) )

      read(io) xt,yt,zt
      read(io) counter
      print*,' counter for averages=',counter
      read(io) kmt
      do k=1,km
         read(io) mean_blue(:,k,:)
      enddo
      close(io)
      deallocate(buf)
      end subroutine read_blue_averages


      subroutine write_blue_averages(outfile,year,month,day,domain)
      use cv_blue_module
      implicit none
      character (len=80) :: outfile,name,unit
      integer :: year,month,day,domain
#include "netcdf.inc"
      integer ncid,iret,i,j,k,dims(4),start(4),count(4)
      real (kind=4) :: spval=-9.9e12, tt
      integer :: lon_tdim,lon_tid
      integer :: lat_tdim,lat_tid
      integer :: depth_tdim,depth_tid
      integer :: itimedim,itimeid
      integer :: blueid,n,len
      integer :: 
     &days_in_month(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

      tt=(year-1900)*365 + day
      do n=1,month-1
       tt=tt+days_in_month(n)
      enddo
      tt=tt-days_in_month(month)/2. ! always monthly averages

      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)

      iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        if (iret.ne.0) print*,nf_strerror(iret)

      write( name,'("SFLAME sub domain # ",i3)') domain
      name=name(1:len_trim(name))//' monthly means'
      len=len_trim(name)
      iret= nf_put_att_text(ncid,nf_global,'history',len,name)

      call ncendf(ncid, iret)

      zt=zt/100.

      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,tt, iret)

c     convert to SI units and divide by N
      mean_blue=mean_blue*1000.  /counter  ! g/cm^3 -> kg/m^3

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

      end subroutine write_blue_averages




      subroutine dvcdf_here(ncid,ivarid,name,iname,unit,iunit,spval)
c-----------------------------------------------------------------------
c     define some standard attributes of variable ivarid in NetCDF file ncid 
c-----------------------------------------------------------------------
      implicit none
      integer ncid,ivarid,iname,iunit,iret
      character*(*) name, unit
      real (kind=4):: spval,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


