!
! File: sift_save.f
! Author: Christopher Moore
!
! At each output timestep, saves timeseries
! at user-selected locations to netCDF output file for Desktop SIFT
!
      subroutine siftsave(ncid,TIME_id,T_stp,t,
     & qa,da,nxa,nya,ma,na,
     & qb,db,nxb,nyb,mb1,nb1,
     & qc,dc,nxc,nyc,mc1,nc1,
     & tsgrid,tsi,tsj,numts,tsid,
     & floodid,inun,h1,h2,
     & latimgA,latimgB,latimgC,dry)

      implicit none
      include 'netcdf.inc'

      integer*4 i,j,k

c.....netCDF id and error status return
      integer*4  ncid, iret

c.....variable ids
      integer*4  TIME_id, floodid

c.....netCDF dimension sizes for dimensions used with record variables
      integer*4  nxa,nya,nxb,nyb,nxc,nyc
      integer*4  ma,na,mb1,nb1,mc1,nc1
c.....starts and counts for array sections of record variables
      integer*4  TIME_start(1), TIME_count(1)

c..... timeseries
      integer*4 numts, tsid
      integer*4 tsgrid(numts),tsi(numts), tsj(numts)
      real*4 ts(1),tst(numts,1)
      integer*4 TS_start(2),TS_count(2)
      
c.....flooding inundation area
      real*8 dval(1), area, maxarea
      real*8 h1(nxc),h2(nyc)
      real*8 PI,RE,GA
      integer*4 inun(nxc,nyc)

c.....data variables
      integer*4  T_stp
      real*8 t, dry
      real*8  TIME(1)
      real*8 qa(nxa, nya),qb(nxb,nyb),qc(nxc,nyc)
      real*8 da(nxa, nya),db(nxb,nyb),dc(nxc,nyc)
c..... image-ordered indices
      logical latimgA, latimgB, latimgC
      save maxarea
      data maxarea /0/
      common /const/ PI,RE,GA

!      write(3,*) ' sift_save'
!      call flush(3)

c.....store TIME
      TIME_start(1) = T_stp
      TIME_count(1) = 1
      TIME(1)=t
      iret = nf_put_vara_double(ncid,TIME_id,TIME_start,TIME_count,TIME)
      call check_err(iret)
c.....store inundation area
      area = 0
      do i=1,mc1-1
        do j= 1,nc1-1
          if (inun(i,j) .gt. 0) then
            area = area +
     &      (h1(i+1)-h1(i))*(h2(j+1)-h2(j))*cos(PI/180*h2(j))*111.320**2
          endif
        end do
      end do
      i=mc1
      do j=1,nc1-1
        if (inun(i,j) .gt. 0) then
          area = area +
     &    (h1(i)-h1(i-1))*(h2(j+1)-h2(j))*cos(PI/180*h2(j))*111.320**2
        endif
      end do
      j=nc1
      do i=1,mc1-1
        if (inun(i,j) .gt. 0) then
          area = area + 
     &    (h1(i+1)-h1(i))*(h2(j)-h2(j-1))*cos(PI/180*h2(j))*111.320**2
        endif
      end do
      if (inun(mc1,nc1) .gt. 0) then
        area = area + 
     &  (h1(mc1)-h1(mc1-1))*(h2(nc1)-h2(nc1-1))*
     &   cos(PI/180*h2(nc1))*111.320**2
      endif
      dval(1)=area
      iret=nf_put_vara_double(ncid,floodid,TIME_start,TIME_count,
     &  dval)
      call check_err(iret)
      if (area.gt.maxarea) then
        maxarea = area
        iret = nf_put_att_double(ncid,floodid,'MaxInundationArea',
     &    NF_REAL, 1, dval)
      endif
      
c.....store timeseries
      TS_start(1) = 1
      TS_start(2) = T_stp
      TS_count(1) = numts
      TS_count(2) = 1
      if (numts.gt.0) then
        do k=1,numts
          tst(k,1)= -9.9999998e+33
          if (tsgrid(k).eq.1) then
            if (latimgA) then
              if (qa(tsi(k),na-tsj(k)+1).GE.dry) then
                tst(k,1)=(qa(tsi(k),na-tsj(k)+1)
     &           -da(tsi(k),na-tsj(k)+1))*100.0
              endif
            else
              if (qa(tsi(k),tsj(k)).GE.dry) then
                tst(k,1)=(qa(tsi(k),tsj(k))-da(tsi(k),tsj(k)))*100.0
              endif
            endif
          else if (tsgrid(k).eq.2) then
            if (latimgB) then
              if(qb(tsi(k),nb1-tsj(k)+1).GE.dry) then
                tst(k,1)=(qb(tsi(k),nb1-tsj(k)+1)
     &           -db(tsi(k),nb1-tsj(k)+1))*100.0
              endif
            else
              if (qb(tsi(k),tsj(k)).GE.dry) then
                tst(k,1)=(qb(tsi(k),tsj(k))-db(tsi(k),tsj(k)))*100.0
              endif
            endif
          else
            if (latimgC) then
              if (qc(tsi(k),nc1-tsj(k)+1).GE.dry) then
                tst(k,1)=(qc(tsi(k),nc1-tsj(k)+1)
     &           -dc(tsi(k),nc1-tsj(k)+1))*100.0
              endif
            else
              if (qc(tsi(k),tsj(k)).GE.dry) then
                tst(k,1)=(qc(tsi(k),tsj(k))-dc(tsi(k),tsj(k)))*100.0
              endif
            endif
          end if
        end do
        iret = nf_put_vara_real(ncid,tsid,TS_start,TS_count,tst)
        call check_err(iret)
      endif


      end


