!
! File: sift_create.f
! Author: Christopher Moore
!
! Creates netCDF output file for Desktop SIFT.
! This routine includes timeseries at user-
! selected locations.
! See sift_create_max.f for max_amp and max_speed
! for each grid.
!
      subroutine siftcreate(ncfn,
     &             h1a,h2a,ma,na,nxa1,nya1,
     &             h1b,h2b,mb1,nb1,nxb1,nyb1,
     &             h1c,h2c,mc1,nc1,nxc1,nyc1,
     &             numts,tsgrid,tsi,tsj,ncid,tstimeid,tsid,
     &             maxhaid,minhaid,maxspaid,
     &             maxhbid,minhbid,maxspbid,
     &             maxhcid,minhcid,maxspcid,
     &             floodid,inunthreshold,
     &             latimgA,latimgB,latimgC)
      implicit none
      include 'netcdf.inc'

      character*(*) ncfn
      integer*4 iret, ncid
      integer*4 i, ii, j, jj, npics
      integer*4 ma, na, nxa1, nya1
      integer*4 mb1, nb1, nxb1, nyb1
      integer*4 mc1, nc1, nxc1, nyc1
      real*8 h1a(nxa1), h2a(nya1)
      real*8 h1b(nxb1), h2b(nyb1)
      real*8 h1c(nxc1), h2c(nyc1)
      integer*4 numts, tsgrid(numts), tsi(numts), tsj(numts),intval(1)
      integer*4 maxhaid,minhaid,maxspaid
      integer*4 maxhbid,minhbid,maxspbid
      integer*4 maxhcid,minhcid,maxspcid,floodid,inunthreshold
      integer*4 tsid, tstimeid, tslatid, tslonid,tsiid,tsjid,tsgridid
      integer*4 lona_id,lata_id,lonb_id,latb_id,lonc_id,latc_id
      integer*4 LON_dim, LAT_dim, TS_dim, TIME_dim
      integer*4 LON_dims(1),LAT_dims(1),TIME_dims(1)
      integer*4 MAX_dims(2), TS_dims(2)
      real*8 tslon(numts), tslat(numts)
      real*4 realval(1)
      real*8 dval(1)
      logical latimgA,latimgB,latimgC

!.....enter define mode
      iret = nf_create(ncfn, NF_CLOBBER, ncid)
      call check_err(iret)

!.....Timeseries dims and vars
!.....define dimensions
      if (numts.gt.0) then
        iret = nf_def_dim(ncid, 'NUMTS', numts, TS_dim)
        call check_err(iret)
      end if
      iret = nf_def_dim(ncid, 'TIME', NF_UNLIMITED, TIME_dim)
      call check_err(iret)
!.....define variables
      TIME_dims(1) = TIME_dim
      iret = nf_def_var(ncid, 'TIME', NF_DOUBLE, 1, TIME_dims,
     &  tstimeid)
      iret = nf_put_att_text(ncid, tstimeid, 'long_name', 4, 'time')
      iret = nf_put_att_text(ncid, tstimeid, 'units', 7, 
     &     'seconds')
      iret = nf_put_att_text(ncid, tstimeid, 'calendar', 9, 'gregorian')
      call check_err(iret)

!.....Create flag for timestep of first inundation
      iret = nf_def_var(ncid, 'inundation',NF_DOUBLE,1,TIME_dims,
     &  floodid)
      call check_err(iret)
      iret = nf_put_att_text(ncid, floodid, 'units', 4, 'km^2')
      intval(1) = -1
      iret = nf_put_att_int(ncid, floodid, 'TimestepOfFirstInundation',
     &  NF_INT, 1, intval)
      intval(1) = inunthreshold
      iret = nf_put_att_int(ncid, floodid, 'InundationCellsThreshold',
     &  NF_INT, 1, intval)
      dval(1) = 0
      iret = nf_put_att_double(ncid,floodid,'MaxInundationArea',
     &  NF_DOUBLE, 1, dval(1))
      call check_err(iret)

      if (numts.gt.0) then
        TIME_dims(1) = TS_dim
        iret=nf_def_var(ncid, 'tslon', NF_DOUBLE, 1,TIME_dims,tslonid)
        call check_err(iret)
        iret=nf_def_var(ncid, 'tsi', NF_INT, 1,TIME_dims,tsiid)
        call check_err(iret)
        iret=nf_put_att_text(ncid,tsiid,'long_name',
     &    28,'1-based index into longitude')
        call check_err(iret)
        iret=nf_def_var(ncid, 'tslat', NF_DOUBLE, 1,TIME_dims,tslatid)
        call check_err(iret)
        iret=nf_def_var(ncid, 'tsj', NF_INT, 1,TIME_dims,tsjid)
        call check_err(iret)
        iret=nf_put_att_text(ncid,tsjid,'long_name',
     &    27,'1-based index into latitude')
        call check_err(iret)
        iret=nf_def_var(ncid, 'tsgrid', NF_INT, 1,TIME_dims,tsgridid)
        call check_err(iret)
        iret=nf_put_att_text(ncid,tsgridid,'long_name',
     &    11,'Grid Number')
        call check_err(iret)
        iret=nf_put_att_text(ncid,tsgridid,'units',
     &    9,'(123=ABC)')
        call check_err(iret)
        TIME_dims(1) = TIME_dim
        TS_dims(1) = TS_dim
        TS_dims(2) = TIME_dim
        iret = nf_def_var(ncid, 'ts', NF_REAL, 2, TS_dims, tsid)
        call check_err(iret)
        iret = nf_put_att_text(ncid, tsid, 'long_name', 14,
     &   'Wave Amplitude')
        call check_err(iret)
        iret = nf_put_att_text(ncid, tsid, 'units', 11, 'CENTIMETERS')
        call check_err(iret)
      endif

!.....Create max_amp and max_speed vars and axes
      call SiftCreateMax(ncid,'lona','lata',ma,na,lona_id,lata_id,
     & 'max_amp_a','min_amp_a','max_speed_a',maxhaid,minhaid,maxspaid)
      call SiftCreateMax(ncid,'lonb','latb',mb1,nb1,lonb_id,latb_id,
     & 'max_amp_b','min_amp_b','max_speed_b',maxhbid,minhbid,maxspbid)
      call SiftCreateMax(ncid,'lonc','latc',mc1,nc1,lonc_id,latc_id,
     & 'max_amp_c','min_amp_c','max_speed_c',maxhcid,minhcid,maxspcid)


!.....leave define mode
      iret = nf_enddef(ncid)
      call check_err(iret)

!.....store LON and LAT
      call SiftWriteAxes(ncid,ma,na,nxa1,nya1,h1a,h2a,lona_id,lata_id,
     & latimgA)
      call SiftWriteAxes(ncid,mb1,nb1,nxb1,nyb1,h1b,h2b,lonb_id,latb_id,
     & latimgB)
      call SiftWriteAxes(ncid,mc1,nc1,nxc1,nyc1,h1c,h2c,lonc_id,latc_id,
     & latimgC)

!.....store timeseries lon, lat
      if (numts.gt.0) then
        do i=1,numts
          if (tsgrid(i).eq.1) then
            if (latimgA) then
              tslon(i) = h1a(tsi(i))
              tslat(i) = h2a(na-tsj(i)+1)
            else
              tslon(i) = h1a(tsi(i))
              tslat(i) = h2a(tsj(i))
            end if
          else if (tsgrid(i).eq.2) then
            if (latimgB) then
              tslon(i) = h1b(tsi(i))
              tslat(i) = h2b(nb1-tsj(i)+1)
            else
              tslon(i) = h1b(tsi(i))
              tslat(i) = h2b(tsj(i))
            end if
          else
            if (latimgC) then
              tslon(i) = h1c(tsi(i))
              tslat(i) = h2c(nc1-tsj(i)+1)
            else
              tslon(i) = h1c(tsi(i))
              tslat(i) = h2c(tsj(i))
            end if
          end if
        enddo
        iret = nf_put_var_double(ncid, tslonid, tslon)
        call check_err(iret)
        iret = nf_put_var_double(ncid, tslatid, tslat)
        call check_err(iret)
        iret = nf_put_var_int(ncid, tsiid, tsi)
        call check_err(iret)
        iret = nf_put_var_int(ncid, tsjid, tsj)
        call check_err(iret)
        iret = nf_put_var_int(ncid, tsgridid, tsgrid)
        call check_err(iret)
      endif

      iret = nf_sync(ncid)
      call check_err(iret)

      return
      end

