program interpolate_adeck

      !----------------------------------------------------------------!
      !---                                                          ---!
      !--- Fortran program to interpolate a TC intensity forecast   ---!
      !--- given the intial intensity forecast at 6 hours... which  ---!
      !--- matches the current synoptic time.                       ---!
      !---                                                          ---!
      !--- Compile statement:                                       ---!
      !---      gfortran -o interpolate_adeck.x interpolate_adeck.f ---!
      !---                                                          ---!
      !--- Usage:                                                   ---!
      !---   ./interpolate_adeck.x int.file model.file H1 H2        ---!
      !--- where:                                                   ---!
      !---        1) int.file   = input file containing TC observed ---!
      !---                        intensity on line 1 and model     ---!
      !---                        intensity on line 2               ---!
      !---        2) model.file = input model intensity file name   ---!
      !---        3) H1 (real)  = beginning hour for interpolation  ---!
      !---        4) H2 (real)  = ending hour for interpolation     ---!
      !---                                                          ---!
      !--- VERSION: 1.0.0                                           ---!
      !---                                                          ---!
      !--- Edit history:                                            ---!
      !---        Matt Onderlinde - June 2018 - wrote original code ---!
      !---                                                          ---!
      !----------------------------------------------------------------!

      implicit none

      !--- Declare variables
      integer, parameter :: nsmooth = 10
      real               :: inth1, inth2            !--- beginning and end hours for linearly decreasing interpolation
      character*999      :: ifile, mfile, ofile, ininth1, ininth2
      character*1        :: dum
      integer            :: n, nf, t, reason, numint
      real               :: tcint, mint, ioffset, sum, curh
      real, allocatable  :: fhour(:), inten(:)
      real, allocatable  :: h3int(:), h3fhr(:), smth3int(:)

      !--- Get input file names from command line
      call get_command_argument(1,ifile)
      call get_command_argument(2,mfile)
      call get_command_argument(3,ininth1)
      call get_command_argument(4,ininth2)
      read(ininth1,*) inth1
      read(ininth2,*) inth2

      !--- Read in initial intensity data
      open(10, file=ifile, status='old')
      read(10,*) dum, tcint
      read(10,*) dum, mint
      close(10)

      !--- Read in 6-hour-old model forecast
      open(20, file=mfile, status='old')
      nf = 0
      do n = 1, 500
         read(20,*,IOSTAT=reason) dum
         if ( reason .lt. 0 ) then
             exit
         end if
         nf = nf + 1
      end do
      rewind(20)
      allocate(fhour(nf))
      allocate(inten(nf))
      do n = 1, nf
          read(20,*) fhour(n), inten(n)
      end do
      close(20)

      !--- Interpolate to 3-hourly
      numint = int(maxval(fhour)/3.0) + 1
      allocate(h3int(numint))
      allocate(h3fhr(numint))
      do n = 1, numint
          curh    = real(n*3) - 3.0
          h3fhr(n) = curh
          do t = 1, nf-1
              if ( fhour(t) .le. curh .and. fhour(t+1) .ge. curh ) then
                  h3int(n) = inten(t) + ((curh-fhour(t))/(fhour(t+1) -
     &                               fhour(t)))*(inten(t+1)-inten(t))
                  exit
              end if
          end do
      end do      

      !--- Perform 3 point smoother on interpolated intensity
      allocate(smth3int(numint))
      do n = 1, nsmooth
          smth3int = h3int
          do t = 2, numint-1
              if ( h3int(t-1) .gt. 0.0 .and. h3int(t+1) .gt. 0.0 ) then
                  smth3int(t) = h3int(t-1)*0.25 + h3int(t)*0.5 +
     &                          h3int(t+1)*0.25
              else
                  smth3int(t) = h3int(t)
              end if
          end do
          h3int = smth3int
      end do

      !--- Shift time 6 hours and set new t=0 values to HSU numbers
      h3fhr = h3fhr - 6.0

      !--- Compute initial intensity offset
      do n = 1, numint
          if ( (abs(h3fhr(n)-0.0)) .lt. 0.001 ) then
              ioffset = tcint - h3int(n)
              exit
          end if
      end do

      !--- Apply relaxation to deterministic forecast over time based on input time range
      do n = 1, numint
          if      ( h3fhr(n) .le. inth1 ) then
              h3int(n) = h3int(n) + ioffset
          else if ( h3fhr(n) .gt. inth1 .and. h3fhr(n) .lt. inth2 ) then
              h3int(n) = h3int(n)+ioffset*
     &                   (1.0-((h3fhr(n)-inth1)/(inth2-inth1)))
          end if
      end do

      !--- Set t=0 value to the HSU-prescribed value
      where(h3fhr .eq. 0.0 ) h3int = tcint

      !--- Write interpolated model intensity data
      ofile = trim(mfile)//".interp"
      open(100, file=trim(ofile), status='unknown')
      do n = 1, numint
          write(100,*) nint(h3fhr(n)), ",", nint(h3int(n))
      end do
      close(100)

      !--- Done
      end program