!  2011 Copyleft by Patrick Corbitt Kerr
!  University of Notre Dame - CHL

!  Created as part of SURA-IOOS Project

program Toolkit

      implicit none
      integer                      :: iopt
      character(10)                :: os="unix"
      character(120)               :: filestationu=' '
      character(120)               :: fort61u=''
      character(120)               :: fort63u=' '
      character(120)               :: fort14u=' '     
      character(120)               :: RunID=' '
      character(120)               :: dfolder=' ' 
      integer                      :: cntrl=0,dtype=1
      integer                      :: dateu(5)

      dateu = (/ 0, 0, 0, 0, 0 /)      
      do
998   write(6,*)
      write(6,*)    '*************************************************'
      write(6,*)    '         NotreDame-CHL Toolkit by PCK  v.43      '    
      write(6,*)    '*************************************************'
      write(6,*)
      write(6,*)    '  1. Create Station List from IMEDS file '
      write(6,*)    '  2. Create IMEDS hwm from IMEDS watlev '
      write(6,*)    '  3. Convert a fort.61/62 to IMEDS format'
      write(6,*)    '  4. *.F, Filter/Remove Stations from an IMEDS file'
      write(6,*)    '  5. Number IMEDS Stations'
      write(6,*)    '  6. Merge IMEDS files'
      write(6,*)    '  7. Adjust Time or Elevation in an IMEDS file'
      write(6,*)    '  8. *.C, Clean an IMEDS file'      
      write(6,*)    '  ------------------------------------------------'
      write(6,*)    '  10. Create a fort.61/62 from a fort.63/64'
      write(6,*)    '  11. Create a wet.14 file'
      write(6,*)    '  12. Find Min-Max'
      write(6,*)    '  13. Create Maximum Envelope of Water (MEOW)'
      write(6,*)    '  14. Generate a max file from a TimeSeries File'
!      write(6,*)    '  15. Project a file onto a different grid'
      write(6,*)    '  16. Generate a Difference File'
      write(6,*)    '  17. *.P, Populate Stations file with Preferred Nodes'      
      write(6,*)    '  18. *.L, Popolate Stations file with new Locations'      
      write(6,*)    '  ------------------------------------------------'
      write(6,*)    '  21. Convert CRMS raw data to IMEDS format'
      write(6,*)    '  22. Convert NOAA raw data to IMEDS format'
      write(6,*)    '  23. Convert TCOON raw data to IMEDS format'
      write(6,*)    '  24. Convert USACE raw data to IMEDS format'
      write(6,*)    '  25. Convert USGS-DEPL raw data to IMEDS format'
      write(6,*)    '  26. Convert USGS-PERM raw data to IMEDS format'
      write(6,*)    '  27. Convert USACE-CHL raw data to IMEDS format'
      write(6,*)    '  28. Convert SLOSH raw data to IMEDS format'
      write(6,*)    '  29. Convert CSI raw data to IMEDS format'
      write(6,*)    '  30. Convert NDBC raw data to IMEDS format'      

      write(6,*)    '  ------------------------------------------------'
      write(6,*)    '  31. Create 61.IMEDS with IMEDS.ctrl'
      write(6,*)    '  ------------------------------------------------'      
      write(6,*)    '  41. Export Open Ocean Boundary Coordinates'
      write(6,*)    '  42. Extract OOB forcing from OSU'      
      write(6,*)    '  43. Extract OOB forcing from FES2004'            
      write(6,*)    '  44. Extract OOB forcing from EC2001'                  
      write(6,*)    '  45. Decompose a fort.53 into *.14 Format files'      
      write(6,*)    '  ------------------------------------------------'
      write(6,*)    '  51. Add Wave-Refraction to Fort.13'
      write(6,*)    '  52. River Flux Boundary Extraction'      
      write(6,*)    '  ------------------------------------------------'      
      write(6,*)    '  98. Change Operating System'
      write(6,*)    '  99. Terminate Program'

997     write(6,'(A)',ADVANCE="NO")    ' Select Option: '
        read(5,'(i10)',ERR=997) iopt
        write(6,*)
        select case(iopt)
          case(1)
            Call IMEDSlist(os)
          case(2)
            Call wl2hwm(os)
          case(3)
            Call IMEDS61(os, cntrl, dateu, filestationu, fort61u, dtype,dfolder)
          case(4)
            CALL IMEDSfilter(os)
          case(5)
            CALL IMEDSrename(os)
          case(6)
            Call IMEDSmerge(os)
          case(7)
            Call IMEDSadjust(os)
          case(8)
            Call IMEDSclean(os)
          case(10)
            Call TSextract(os,cntrl,filestationu,fort14u,fort63u,RunID,dfolder)            
          case(11)
            Call wetdry(os)
          case(12)
            Call findminmax(os)
          case(13)
            Call meow(os)
          case(14)
            Call generatemax(os)
          case(15)
            Call projector(os)
          case(16)
            call generatediff(os)
          case(17)
            call listnodes(os)
          case(18)
            call relocator(os)            
          case(21)
            Call CRMS2IMEDS
          case(22)
            Call NOAA2IMEDS
          case(23)
            CALL TCOON2IMEDS
          case(24)
            CALL USACE2IMEDS
          case(25)
            Call USGSDEPL2IMEDS
          case(26)
            Call USGSPERM2IMEDS
          case(27)
            Call USACECHL
          case(28)
            Call SLOSH2IMEDS
          case(29)
            Call CSI2IMEDS
          case(30)
            Call NDBC2IMEDS            
          case(31)
            Call IMEDScntrl(os)
          case(41)
            Call ExportOOB(os)            
          case(42)
            Call OSUextract
          case(43)
            Call FES2004extract
          case(44)
            Call EC2001extract
          case(45)
            Call Decompose53            
          case(51)
            Call AddWave13
          case(52)
            Call RiverExtract
          case(98)
            if (trim(os).eq."unix") then
              write(6,*) "Operating System changed to DOS"
              os ="dos"
            else
              write(6,*) "Operating System changed to UNIX"
              os ="unix"
            endif
          case(99)
            goto 999
          case default
            write(6,'(A)') 'Option Unacceptable, Try Again'
            goto 998
        end select
        write(6,*) 
        write(6,*) "Yatta!"
        write(6,*)
      enddo
      write(6,*)
999   write(6,*) 'Program Terminated'
   end program
!----------------------------------------------------------------------
!     1
!----------------------------------------------------------------------
   subroutine IMEDSlist(os)

      implicit none
      character(120)                :: filest,junkC,junkC2,junkC3
      character(1000)               :: junkCC
      character(120)                :: filewl
      character(10), intent(in)     :: os
      integer                       :: extension=1
      integer                       :: iflag, i, n, X1
      write(6,*) "Enter the name of the input watlev file"
      call choosefile(filewl,extension,os)
      X1 = len_trim(filewl)
      if (filewl(X1-4:X1).eq.'IMEDS') then
        filest=filewl(1:X1-5)//'Stations'
      else
        filest=trim(filewl)//'.Stations'
      endif
      open(unit=10,file=filewl,status='old',action='read')
      open(unit=20,file='.tempfile',status='replace',action='write')
      read(10,'(a)') junkC
      read(10,'(a)') junkC
      read(10,'(a)') junkC
      read(10,'(a)') junkC
      write(20,*) trim(junkC)        !Write station 
      n= 1
      iflag = 1
      do while (iflag.eq.1)
        do while (iflag.eq.1)
          read(10,'(a)',END=999) junkC2
          JunkC3 = adjustL(junkC2)
          if (JunkC3(1:1).ne."2") goto 998  !bad programming
        enddo
998     write(20,*) trim(junkC2)  !Write Station
        n = n + 1
      enddo
999   close(10)
      close(20)
      open(unit=20,file='.tempfile',status='old',action='read')
      open(unit=10,file=filest,status='replace',action='write')
      write(10,*) n
      do i = 1,n
        read(20,'(A)') junkCC
        write(10,'(A)') trim(junkCC)
      enddo
      close(10)
      close(20)

   end subroutine IMEDSLIST
!----------------------------------------------------------------------
!     2
!----------------------------------------------------------------------
   subroutine wl2hwm(os)

      implicit none
      character(120)                :: filewl, filehwm, templine,templine2
      character(10), intent(inout)  :: os
      character(120)                :: templine3
      character(120),allocatable    :: ct(:), mt(:)
      integer                       :: iflag, i, n, Defval=-999
      double precision              :: wse, mxe

      write(6,*) "Enter the name of the input watlev file"
      read(5,*) filewl
      write(6,*) "Enter the name of the output hwm file"
      read(5,*) filehwm
      write(6,*)

      open(unit=10,file=filewl,status='old',action='read')
      open(unit=20,file=filehwm,status='replace',action='write')
      read(10,*)
      write(20,*) '% IMEDS generic format - Water Level'
      read(10,'(a)') templine
      if (templine(28:30).eq.'sec') then
        n = 6
        write(20,*) "% ", "year ","month ","day ","hour ","min ","sec ", "hwm (m)"
      elseif (templine(24:26).eq.'min') then
        n = 5
        write(20,*) "% ", "year ","month ","day ","hour ","min ", "hwm (m)"
      elseif (templine(19:22).eq.'hour') then
        n = 4
        write(20,*) "% ", "year ","month ","day ","hour ", "hwm (m)"
      else
        n = 3
        write(20,*) "% ", "year ","month ","day ", "hwm (m)"
      endif
      read(10,'(a)') templine
      write(20,*) trim(templine)        !Write Header Line 3
      read(10,'(a)') templine
      write(20,*) trim(templine)        !Write station 
      allocate(ct(1:n), mt(1:n))
      iflag = 1
      do while (iflag.eq.1)
        mxe = -999.0d0
        do while (iflag.eq.1)
          read(10,'(a)',END=999) templine2
          templine3 = adjustL(templine2)
          if (templine3(1:1).ne."2") goto 998  !bad programming
          read(templine2(1:len_trim(templine2)),*) (ct(i),i=1,n) , wse
          if (wse.gt.mxe) then
            do i=1,n
              mt(i) = ct(i)
            enddo
            mxe=wse
          endif
        enddo
998     write(20,*) (trim(mt(i))//"    ",i=1,n) , mxe
        write(20,*) trim(templine2)  !Write Station
      enddo
999   write(20,*) (trim(mt(i))//"    ",i=1,n) , mxe
      close(10)
      close(20)
   end subroutine wl2hwm
!----------------------------------------------------------------------
!     3
!----------------------------------------------------------------------
   subroutine IMEDS61(os, cntrl, dateu, filestationu, fort61u, dtype,dfolder)

      implicit none
      character(120)                :: filestation, fort61
      character(10), intent(inout)  :: os
      character(120), intent(inout) :: filestationu
      character(120), intent(inout) :: fort61u, dfolder      
      integer, intent(inout)        :: cntrl, dtype
      integer, intent(inout)        :: dateu(5)
      character(120), allocatable   :: stnid(:)
      character(22)                 :: junkC
      integer                       :: Interval, nCol
      integer                       :: junki, Nsteps, Nsta, i, j, k
      integer                       :: yy, mo, dd, hh, mi, ss=0
      integer,allocatable           :: calen(:), dates(:,:), cd(:)
      double precision              :: junkd, day1, day2
      integer                       :: extension,iopt
      integer                       :: defval=-999
      double precision, allocatable :: lat(:), lon(:), time(:)
      double precision, allocatable :: val(:,:), val2(:,:)
  

      if (cntrl.eq.0) then
        write(6,*) "What is the date of the simulation start? (YYYY, MO, DD, HH, MI )"
        read(5,*) yy, mo, dd, hh, mi
        write(6,*)
        write(6,*) "Enter the name of the Stations ID file"
        extension=2
        call choosefile(filestation,extension,os)
        write(6,*)
        write(6,*) "Enter the name of the fort.61 file"
        extension=3
        call choosefile(fort61,extension,os)
        write(6,*)
      else
        filestation = trim(filestationu)
        fort61 = trim(fort61u)
        yy = dateu(1)
        mo = dateu(2)
        dd = dateu(3)
        hh = dateu(4)
        mi = dateu(5)        
      endif
!
      open(unit=61,file=fort61,status='old',action='read')
      read(61,*)
      read(61,*) Nsteps, Nsta, Interval, Interval, nCol
      allocate(stnid(1:Nsta),lat(1:Nsta),lon(1:Nsta))
      allocate(val(1:Nsteps,1:Nsta), time(1:Nsteps), dates(1:5,1:Nsteps))
      allocate(val2(1:Nsteps,1:Nsta))
!! Collect data from ADICRC output file
      i = 0
      do
         i = i + 1
         read(61,*,ERR=999,END=999) time(i)
         write(6,*) "Collecting data for timesnap ", i, "."
         do j = 1,Nsta
           if (nCol.eq.1) then
             read(61,*,ERR=999,END=999) junkI, val(i,j)
           else
             read(61,*,ERR=999,END=999) junkI, val(i,j), val2(i,j)
           endif
!           read(61,*,ERR=999,END=999) junkI, junkC
!           if (junkC.eq."NaN") then
!             val(i,j) = -99999.0d0
!           else
!             read(junkC,*) val(i,j)
!           endif
         enddo
      enddo
999   Nsteps = i-1
      close(61)
      
!! Break time into component parts
      allocate(calen(1:12), cd(1:13))
      calen = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
      cd = (/ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 /)
      if( mod(yy,4) == 0 ) then
        calen(2) = 29 !leap year
        cd = (/ 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366 /)
      endif
      day1 = cd(mo) + dd + hh/24.0d0 + mi/1440.0d0
      do i = 1,Nsteps
         day2 = time(i)/86400.0d0 + day1
         dates(1,i) = yy
         do j = 1,13
           if(cd(j).gt.day2-1) GoTo 500           
         enddo
  500    dates(2,i) = j-1
         dates(3,i) = floor(day2-cd(j-1))
         dates(4,i) = floor(24.0d0*((day2-cd(j-1))-dates(3,i)))
         dates(5,i) = floor(60.0d0*((24.0d0*((day2-cd(j-1))-dates(3,i))-dates(4,i))))
      enddo
      
!! Collect station information from station id file
      open(unit=11,file=trim(dfolder)//trim(filestation),status='old',action='read')
      read(11,*)
      do j = 1,Nsta
         read(11,*) stnid(j), lat(j), lon(j)
      enddo
      close(11)
      
!! Reformat data as IMEDS generic output
      open(unit=21,file=TRIM(fort61)//'.IMEDS',status='replace',action='write')
      if (cntrl.eq.0) then
        write(6,*) "IMEDS file types:"
        write(6,*) "  1. watlev - Water Level"
        write(6,*) "  2. hs     - Significant Wave Height"
        write(6,*) "  3. dir    - Mean Wave Direction"
        write(6,*) "  4. tp     - Peak Period"
        write(6,*) "  5. tm     - Mean Period"
        write(6,*) "  6. hwm    - High Water Mark"
        write(6,*) "  7. cur    - Current Velocity"
        write(6,*) "Select Option: "
        read(5,*) iopt
      else
        iopt=dtype
      endif
      select case(iopt)
        case(1)
          write(21,*) "% IMEDS generic format - Water Level"
          write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
        case(2)
          write(21,*) "% IMEDS generic format - wave-height"
         write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "hs (m)"
        case(3)
          write(21,*) "% IMEDS generic format - mean wave-direction"
          write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "dir (degrees)"
        case(4)
          write(21,*) "% IMEDS generic format - Peak Period"
          write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "tp (s)"
        case(5)
          write(21,*) "% IMEDS generic format - Mean Period"
          write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "tp (s)"
        case(6)
          write(21,*) "% IMEDS generic format - High Wate Mark"
          write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "hwm (m)"          
        case(7)
          write(21,*) "% IMEDS generic format - Current"
          write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "cur (m/s)"          
      end select
      write(21,*) "ADCIRC " , "UTC", " NAVD88"  !ADCIRC DEFAULT
      do j = 1,Nsta
         write(21,*) TRIM(stnid(j)), lat(j), lon(j)
         do i = 1,Nsteps
            select case(iopt)
            case(1,2,3,4,5,6)
              if (val(i,j).eq.-99999.0d0) then
                write(21,201) (dates(k,i),k=1,5),ss, defval
              else
                write(21,200) (dates(k,i),k=1,5),ss, val(i,j)
              endif
            case(7)
              if (val(i,j).eq.-99999.0d0) then
                write(21,203) (dates(k,i),k=1,5),ss, defval, defval
              else
                write(21,202) (dates(k,i),k=1,5),ss, val(i,j), val2(i,j)
              endif
            end select
         enddo
      enddo
200   Format(6i7,f14.4)
201   Format(6i7,i10)
202   Format(6i7,2f24.10)
203   Format(6i7,2i10)
      close(21)

   end subroutine IMEDS61
!----------------------------------------------------------------------
!     4
!----------------------------------------------------------------------
   subroutine IMEDSfilter(os)

      implicit none
      character(120)                :: fileIMEDS,filefilter,junkC
      character(120)                :: flagfile, tempC
      character(10), intent(inout)  :: os
      integer                       :: iflag, i, nsta
      integer                       :: extension, X1,ifilter
      integer, allocatable          :: filter(:)
      character(120),allocatable    :: StnID(:)
      real(8)                       :: junkR
      
      extension=1
      write(6,*) "Enter the name of the IMEDS file"
      call choosefile(fileIMEDS,extension,os)
      X1 = len_trim(fileIMEDS)-6
      filefilter = fileIMEDS(1:X1)//'.F.IMEDS'
      write(6,*) "Enter the name of the flag file"
      extension=4
      call choosefile(flagfile,extension,os)
      nsta=0
      open(10,file=trim(flagfile),status='old',action='read')
        do 
          read(10,*, END=997)
          nsta = nsta + 1
        enddo
997   close(10)
      nsta = nsta-4
      allocate(StnID(1:nsta),filter(1:nsta))
      open(10,file=trim(flagfile),status='old',action='read')
        read(10,*)
        read(10,*)
        read(10,*)
        read(10,*)
        do i = 1,nsta
          read(10,*) StnID(i), junkR, junkR, filter(i)
        enddo        
      close(10)
      open(10,file=trim(fileIMEDS),status='old',action='read')
      open(unit=20,file=trim(filefilter),status='replace',action='write')
      read(10,'(a)') junkC
      write(20,*) trim(junkC)
      read(10,'(a)') junkC
      write(20,*) trim(junkC)      
      read(10,'(a)') junkC
      write(20,*) trim(junkC)
      read(10,'(a)') junkC
      do
        read(junkC,*,ERR=999,END=999) tempC
        ifilter=1
        do i=1,nsta
          if (trim(tempC).eq.trim(StnID(i))) ifilter= filter(i)
        enddo     
        if (ifilter.eq.1) write(20,*) junkC
        do
          read(10,'(a)',END=999) junkC
          JunkC = adjustL(junkC)
          if (JunkC(1:1).ne."2") goto 998  !bad programming
          if (ifilter.eq.1) write(20,*) junkC
        enddo
998     continue
      enddo
999   close(10)
      close(20)

   end subroutine IMEDSfilter
!----------------------------------------------------------------------
!     5
!----------------------------------------------------------------------
   subroutine IMEDSrename(os)

      implicit none
      character(120)                :: fileIMEDS,filestations,junkC,Prefix,filenew
      character(10), intent(inout)  :: os
      integer                       :: iflag, i, nsta, X1
      integer, allocatable          :: filter(:)
      character(120)                :: StnID, NewID, junkC2,filenum
      integer                       :: extension=1
      double precision              :: Lat, Lon

      write(6,*) "Enter the name of the IMEDS file"
      call choosefile(fileIMEDS,extension,os)
      filestations = trim(fileIMEDS)//'.Stations'
      X1 = len_trim(fileIMEDS)      
      filenew = fileIMEDS(1:X1-6)//'.N.IMEDS'
      filenum = fileIMEDS(1:X1-6)//'.N'

      open(10,file=trim(fileIMEDS),status='old',action='read')
      open(unit=20,file=trim(filenew),status='replace',action='write')
      open(30,file=trim(filenum),status='replace',action='write')
      read(10,'(a)') junkC
      write(20,*) trim(junkC)
      read(10,'(a)') junkC
      write(20,*) trim(junkC)      
      read(10,'(a)') junkC
      write(20,*) trim(junkC)
      iflag = 1
      read(10,'(a)') junkC
      i = 10000
      do 
        i = i + 1
        read(junkC(1:len_trim(junkC)),*) StnID, Lat, Lon
        write(NewID,fmt='(I10)') i
        X1 = index(StnID,"_")
        Prefix = StnID(1:X1-1)
        write(20,*) trim(trim(Prefix)//adjustL(newID)), Lat, Lon
        write(30,*) trim(trim(Prefix)//adjustL(newID)), "   ", trim(StnID)
        do while (iflag.eq.1)
          read(10,'(a)',END=999) junkC
          JunkC = adjustL(junkC)
          if (JunkC(1:1).ne."2") goto 998  !bad programming
          write(20,*) trim(junkC)
        enddo
998     continue
      enddo
999   close(10)
      close(20)
      close(30)

   end subroutine IMEDSrename
!----------------------------------------------------------------------
!     6
!----------------------------------------------------------------------
   subroutine IMEDSmerge(os)

      implicit none
      character(120)                 :: junkC, Dsource, Dtz, Ddatum,fileIMEDS
      character(120)                 ::  Header1, Header2, imerge,filemerge
      integer                        :: i, iopt,extension
      character(10),intent(inout)      :: os
     
      write(6,*) "Note that this Merges files irrespective of Header Differences"
      write(6,*)
      write(6,*) "Select format option"
      write(6,*) " 1. watlev (m)"
      write(6,*) " 2. hs (m)"
      write(6,*) " 3. dir (deg)"
      write(6,*) " 4. tp (s)"
      write(6,*) " 5. tm (s)"
      write(6,*) " 6. hwm (m)"
      write(6,*) " 7. cur (m/s)"
      read(5,*)  iopt
      write(6,*)
      write(6,*) "Enter Source"
      read(5,*) Dsource
      write(6,*) "Enter Time Zone"
      Dtz="UTC"
!      read(5,*) Dtz
      write(6,*) "Enter Datum"
      Ddatum="NAVD88"
!      read(5,*) Ddatum
      write(6,*)
      write(6,*) "Enter the name of the mergefile"
      read(5,*) filemerge

      open(10,file=filemerge,status='replace',action='write')
      Select Case (iopt)
        case(1)
          write(10,'(A)') "% IMEDS generic format - Water Level"
          write(10,'(A)') "% year month day hour min sec watlev (m)"
        case(2)
          write(10,'(A)') "% IMEDS generic format - Significant Wave Height"
          write(10,'(A)') "% year month day hour min sec hs (m)"
        case(3)
          write(10,'(A)') "% IMEDS generic format - Mean Wave Direction"
          write(10,'(A)') "% year month day hour min sec dir (deg)"
        case(4)
          write(10,'(A)') "% IMEDS generic format - Peak Period"
          write(10,'(A)') "% year month day hour min sec tp (s)"
        case(5)
          write(10,'(A)') "% IMEDS generic format - Mean Period"
          write(10,'(A)') "% year month day hour min sec tp (s)"
        case(6)
          write(10,'(A)') "% IMEDS generic format - High Water Mark"
          write(10,'(A)') "% year month day hour min sec hwm (m)"
        case(7)
          write(10,'(A)') "% IMEDS generic format - Current Velocity"
          write(10,'(A)') "% year month day hour min sec cur (m/s)"
      end select
      write(10,101) trim(Dsource),trim(Dtz), trim(Ddatum)
101   Format(3A10)
      
      extension=1
      do
        write(6,*) "Enter the name of the IMEDS file" 
        call choosefile(fileIMEDS,extension,os)
        open(20,file=trim(fileIMEDS),status='old',action='read')
        do i=1,3
          read(20,'(A)') junkC
        enddo
        do
          read(20,'(A)',END=993) junkC
          write(10,'(A)') trim(junkC)
        enddo
993     close(20)
        write(6,*) "Do you Wish to Merge Another File? (Y or N)?"
        read(5,*) imerge
        if ((imerge(1:1).eq."n").or.(imerge(1:1).eq."N")) goto 994
      enddo
994   close(10)

   end subroutine IMEDSmerge
!----------------------------------------------------------------------
!     7
!----------------------------------------------------------------------
   subroutine IMEDSadjust(os)

      implicit none
      character(120)                 :: fileIMEDSold, fileIMEDSnew
      character(120)                 :: Header1, Header2, Model
      character(120)                 :: TimeZone, Datum
      character(120)                 :: StnID, JunkC
      real(8)                        :: Lat, Lon
      integer                        :: i, iopt,extension=1, X1
      real(8)                        :: MultR, DatAdj, TempR
      character(10),intent(inout)    :: os
      integer                        :: yy, mo, dd, hh, mi, ss, TZC
     
      write(6,*) "This program will adjust the time and/or Elevation"
      write(6,*)
      write(6,*) "Enter the name of the IMEDS file" 
      call choosefile(fileIMEDSold,extension,os)
      fileIMEDSnew = trim(fileIMEDSold)//'.A.IMEDS'
      open(10,file=trim(fileIMEDSold),status='old',action='read')
      read(10,'(A)') Header1
      read(10,'(A)') Header2
      read(10,*) Model, TimeZone, Datum
      X1 = len_trim(Header2)
      if (Header2(X1-1:X1-1).eq.'m') then
        Header2 = " % year month day hour min sec watlev (m)"
      elseif (Header2(X1-1:X1-1).eq.'t') then
        Header2 = " % year month day hour min sec watlev (ft)"
      endif      
      write(6,*)      
      write(6,*) ' Select Adjustment Option'
      write(6,*) ' 1. Change Time Zone'
      write(6,*) ' 2. Change Datum and Vertically Adjust Values'
      write(6,*) ' 3. Change Units from feet to meters'
      write(6,*) ' 4. Change Units from meters to feet'
      read(5,*) iopt
      MultR= 1.0d0
      Select Case (iopt)
        case(1) 
          write(6,*) 'Enter Time Zone'
          read(5,*) TimeZone
          write(6,*) 'Enter Time Zone Adjustment in Hours'
          read(5,*) TZC
        case(2)
          write(6,*) 'Enter Datum'
          read(5,*) Datum
          write(6,*) 'Enter Vertical Adjustment'
          read(5,*) DatAdj
        case(3)
          Header2 = " % year month day hour min sec watlev (m)"
          MultR = 0.3048d0
        case(4)
          Header2 = " % year month day hour min sec watlev (ft)"
          MultR = 1.0d0/0.3048d0
      end select

      open(20,file=trim(fileIMEDSnew),status='replace',action='write')
      write(20,'(A)') trim(Header1)
      write(20,'(A)') trim(Header2)
      write(20,*) trim(Model), " ", trim(TimeZone), " ", trim(Datum)
      read(10,'(A)') JunkC 
      do 
        read(junkC(1:len_trim(junkC)),*) StnID, Lat, Lon
        write(20,*) trim(StnID), Lat, Lon
        do
          read(10,'(a)',END=999) junkC
          JunkC = adjustL(junkC)
          mi=0
          ss=0
          if (JunkC(1:1).ne."2") goto 998  !bad programming
          read(junkC,*,ERR=995,END=995) yy, mo, dd, hh, mi, ss, tempR
          goto 997
995       read(junkC,*,ERR=996,END=996) yy, mo, dd, hh, mi, tempR
          goto 997
996       read(junkC,*) yy, mo, dd, hh, tempR
997       call  TZconvert(yy,mo,dd,hh,TZC)
          TempR = TempR*MultR
          TempR = TempR + DatAdj
          write(20,201) yy, mo, dd, hh, mi, ss, tempR
        enddo
998     continue
      enddo
999   close(10)
      close(20)
201   Format(6i7,f14.4)
202   Format(6i7,i14)

   end subroutine IMEDSadjust   
!----------------------------------------------------------------------
!     8
!----------------------------------------------------------------------
   subroutine IMEDSclean(os)

      implicit none
      character(120)                 :: Header1, Header2, fileIMEDSold
      character(120)                 :: fileIMEDSnew, fileIMEDSclean
      character(120)                 :: Model, TimeZone, Datum
      integer                        :: extension,tempI=-999, iexist
      character(10),intent(inout)    :: os
      character(120)                 :: JunkC, JunkC2, StnID, StnID2
      real(8)                        :: Lat, Lon, tempR
      integer                        :: yy, mo, dd, hh, mi, ss, X1
      integer                        :: yy2, mo2, dd2, hh2, mi2, ss2


     
      write(6,*) "This program will clean an IMEDSfile"
      write(6,*)
      write(6,*) "Enter the name of the IMEDS file to clean" 
      extension=1
      call choosefile(fileIMEDSold,extension,os)
      write(6,*) "Enter the name of the IMEDS flag file"
      extension=4
      call choosefile(fileIMEDSclean,extension,os)
      X1=len_trim(fileIMEDSold)-6
      fileIMEDSnew = fileIMEDSold(1:X1)//'.C.IMEDS'
      open(10,file=trim(fileIMEDSold),status='old',action='read')
      open(20,file=trim(fileIMEDSnew),status='replace',action='write')

      read(10,'(A)') Header1
      read(10,'(A)') Header2
      read(10,*) Model, TimeZone, Datum
      write(20,'(A)') trim(Header1)
      write(20,'(A)') trim(Header2)
      write(20,*) trim(Model), " ", trim(TimeZone), " ", trim(Datum)      

      read(10,'(A)') JunkC 
      do 
        read(junkC,*,ERR=999,END=999) StnID, Lat, Lon
        write(6,*) trim(StnID), Lat, Lon
        write(20,*) trim(StnID), Lat, Lon
        open(30,file=trim(fileIMEDSclean),status='old',action='read')      
          do
            read(30,*,END=901) StnID2
            if (StnID2.eq.StnID) goto 908
          enddo
901       do
            read(10,'(a)',END=999) junkC
            junkC = adjustL(JunkC)
            if (JunkC(1:1).ne."2") goto 900 !bad programming
            write(20,'(A)') JunkC
          enddo
908       do
907         read(30,'(A)',END=901) junkC2
            junkC2 = adjustL(junkC2)
            if (JunkC2(1:1).ne."2") goto 901
            read(JunkC2,*,END=900) yy2, mo2, dd2, hh2, mi2, ss2
            do
              read(10,'(a)',END=999) junkC
              junkC = adjustL(JunkC)
              if (JunkC(1:1).ne."2") goto 900  !bad programming
              read(junkC,*,END=999) yy, mo, dd, hh, mi, ss
              if ((yy2.eq.yy).and.(mo2.eq.mo).and.(dd2.eq.dd).and.(hh2.eq.hh).and.(mi2.eq.mi).and.(ss2.eq.ss)) then
                write(20,202) yy, mo, dd, hh, mi, ss, tempI
                goto 907
              else
                write(20,'(A)') JunkC
              endif
            enddo
          enddo
900     close(30)
      enddo
999   close(10)
      close(20)
202   Format(6i7,i14)

   end subroutine IMEDSclean     
!----------------------------------------------------------------------
!     10
!----------------------------------------------------------------------
   subroutine TSextract(os,cntrl,filestationu,fort14u,fort63u,RunID,dfolder)

      IMPLICIT NONE

      INTRINSIC           :: ABS
      INTRINSIC           :: TRIM
      character(10), intent(inout)     :: os
      character(120),intent(inout)     :: filestationu
      character(120),intent(inout)     :: fort14u
      character(120),intent(inout)     :: fort63u
      character(120),intent(inout)     :: RunID      
      character(120),intent(inout)     :: dfolder
      integer,intent(inout)            :: cntrl
      integer                          :: extension

      CHARACTER(LEN=1)    :: Ivector
      CHARACTER(LEN=1)    :: Isparse
      CHARACTER(LEN=100)  :: GridFile
      CHARACTER(LEN=100)  :: GlobalFile
      CHARACTER(LEN=1)    :: JunkC
      CHARACTER(LEN=500)  :: Line
      CHARACTER(LEN=100)  :: OutputFile
      CHARACTER(LEN=500)  :: StationsFile
      CHARACTER(LEN=50)   :: Comments
      CHARACTER(120)      :: JunkC2

      INTEGER,ALLOCATABLE :: Conn(:,:)
      INTEGER             :: E
      INTEGER             :: JunkI
      INTEGER             :: N
      INTEGER             :: Interval
      INTEGER             :: NumElems
      INTEGER             :: NumNodes
      INTEGER             :: NumNodesNonDefault
      INTEGER             :: NumSnaps
      INTEGER             :: NumStations
      INTEGER             :: S
      INTEGER             :: SS
      INTEGER,ALLOCATABLE :: StationsElem(:)
      INTEGER             :: UserSel

      REAL(8)             :: DefaultValue
      REAL(8),ALLOCATABLE :: Global1(:)
      REAL(8),ALLOCATABLE :: Global2(:)
      REAL(8),ALLOCATABLE :: GridLat(:)
      REAL(8),ALLOCATABLE :: GridLon(:)
      REAL(8)             :: JunkR
      REAL(8),ALLOCATABLE :: StationsLat(:)
      REAL(8),ALLOCATABLE :: StationsLon(:)
      REAL(8)             :: StationsTemp1
      REAL(8)             :: StationsTemp2
      REAL(8)             :: StationsTemp11
      REAL(8)             :: SubArea1
      REAL(8)             :: SubArea2
      REAL(8)             :: SubArea3
      REAL(8)             :: TotalArea
      REAL(8),ALLOCATABLE :: Weight(:,:)
      REAL(8)             :: X1
      REAL(8)             :: X2
      REAL(8)             :: X3
      REAL(8)             :: X4      
      REAL(8)             :: Y1
      REAL(8)             :: Y2
      REAL(8)             :: Y3
      REAL(8)             :: W1
      REAL(8)             :: W2
      REAL(8)             :: W3      
      REAL(8), allocatable         :: Bath(:), UseDepth(:)
      integer, allocatable         :: WN(:),Grid(:), SN(:), UseNode(:)
      character(120), allocatable  :: StnID(:)
      integer                      :: nCol
      REAL(8)                      :: SnapR
      Integer                      :: SnapI, XX1
      character(120)               :: filestationi

      ! Read Station File      
      if (cntrl.eq.0) then
        WRITE(*,'(A)',ADVANCE="YES") "Enter name of file with lat/lon or node numbers for stations: "
        extension = 2
        call choosefile(filestationi,extension,os)
        StationsFile = trim(filestationi)
        dfolder = ''
      else
        StationsFile = trim(filestationu)
      endif

      OPEN(UNIT=11,FILE=trim(dfolder)//TRIM(StationsFile),Status='OLD',ACTION="READ")
        READ(11,*) NumStations
        ALLOCATE(StationsLon(1:NumStations))
        ALLOCATE(StationsLat(1:NumStations))
        allocate(StnID(1:NumStations),UseNode(1:NumStations),UseDepth(1:NumStations))
        ALLOCATE(Weight(1:NumStations,1:3))
        allocate(grid(1:Numstations),SN(1:NumStations))
        allocate(WN(1:NumStations))
        ALLOCATE(StationsElem(1:NumStations))
        DO S=1,NumStations
           UseNode(S) = 0
           UseDepth(S) = -999.0d0
           READ(11,'(A)') Line
           READ(Line,*) StnID(S), StationsLat(S), StationsLon(S)
           READ(Line,*,ERR=110,END=110) StnID(S), StationsLat(S), StationsLon(S), UseNode(S)
           READ(Line,*,ERR=110,END=110) StnID(S), StationsLat(S), StationsLon(S), UseNode(S), UseDepth(S)
110        CONTINUE
        ENDDO
      CLOSE(UNIT=11,STATUS="KEEP")
 
      ! Read Grid File
      WRITE(*,'(A)',ADVANCE="YES") "Stations information was read successfully!"
      WRITE(*,'(A)',ADVANCE="YES") " "
      if (cntrl.eq.0) then
        WRITE(*,'(A)',ADVANCE="NO") "Enter name of ADCIRC grid file: "
        READ(*,'(A)') GridFile
      else
        GridFile = trim(fort14u)
      endif
      OPEN(UNIT=14,FILE=TRIM(GridFile),ACTION="READ")
        READ(14,'(A)') JunkC
        READ(14,*) NumElems, NumNodes
        ALLOCATE(GridLon(1:NumNodes),GridLat(1:NumNodes),Bath(1:NumNodes))
        allocate(Conn(1:NumElems,1:3))
        DO N=1,NumNodes
           READ(14,*) JunkI, GridLon(N), GridLat(N), Bath(N)
        ENDDO
        DO E=1,NumElems
           READ(14,*) JunkI, JunkI, Conn(E,1), Conn(E,2), Conn(E,3)
        ENDDO
      CLOSE(UNIT=14,STATUS="KEEP")

      ! Identify Interpolation Weights & Prepare Summary File      
      OPEN(unit=86,file=trim(StationsFile)//'.'//trim(GridFile)//'.Summary',status='replace',action='write')
      write(86,667) 'StationID','Inside Grid','Near Node','Near WetNode', 'Near Depth', &
      'Wet Depth','Station Lat','Station Lon','Node Lat','Node Lon','Wet Lat', 'Wet Lon'
667   Format(A50,11A15)
      DO S=1,NumStations
         IF (StationsLon(S).gt.-9999.0d0) THEN
         Grid(S) = 0
         SN(S) = 0
         WN(S) = 0
         inner: DO E=1,NumElems
            X1 = StationsLon(S)
            X2 = GridLon(Conn(E,2))
            X3 = GridLon(Conn(E,3))
            Y1 = StationsLat(S)
            Y2 = GridLat(Conn(E,2))
            Y3 = GridLat(Conn(E,3))
            SubArea1 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))
            X1 = GridLon(Conn(E,1))
            X2 = StationsLon(S)
            X3 = GridLon(Conn(E,3))
            Y1 = GridLat(Conn(E,1))
            Y2 = StationsLat(S)
            Y3 = GridLat(Conn(E,3))
            SubArea2 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))
            X1 = GridLon(Conn(E,1))
            X2 = GridLon(Conn(E,2))
            X3 = StationsLon(S)
            Y1 = GridLat(Conn(E,1))
            Y2 = GridLat(Conn(E,2))
            Y3 = StationsLat(S)
            SubArea3 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))
            X1 = GridLon(Conn(E,1))
            X2 = GridLon(Conn(E,2))
            X3 = GridLon(Conn(E,3))
            Y1 = GridLat(Conn(E,1))
            Y2 = GridLat(Conn(E,2))
            Y3 = GridLat(Conn(E,3))
            TotalArea = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))
            IF((SubArea1+SubArea2+SubArea3).LE.(1.01*TotalArea))THEN
               StationsElem(S) = E
               Weight(S,1) = ( (StationsLon(S)-X3)*(Y2-Y3)+(X2-X3)*(Y3-StationsLat(S)))/TotalArea
               Weight(S,2) = ( (StationsLon(S)-X1)*(Y3-Y1)-(StationsLat(S)-Y1)*(X3-X1))/TotalArea
               Weight(S,3) = (-(StationsLon(S)-X1)*(Y2-Y1)+(StationsLat(S)-Y1)*(X2-X1))/TotalArea
               PRINT*, Weight(S,1), Weight(S,2), Weight(S,3)
               Grid(S) = 1 !Station is Located within an Element
               EXIT inner
            ENDIF
         ENDDO inner
         ENDIF
         if ((Grid(S).eq.0)) then
           write(6,*) trim(StnID(S)), " is located outside of Grid"
         endif
         X3 = 1000000000.0d0
         X4 = 1000000000.0d0         
         do N=1,NumNodes
           X1 = GridLon(N)-StationsLon(S)
           Y1 = GridLat(N)-StationsLat(S)
           X2 = (X1**2.0d0+Y1**2.0d0)**0.5d0
           if (X2.lt.X3) then
             X3 = X2
             SN(S) = N !Identify Nearest Node
           endif
           if ((X2.lt.X4).and.(Bath(N).ge.UseDepth(S))) then
             X4 = X2
             WN(S) = N !Identify Nearest Node
           endif
         enddo
         write(86,666) adjustL(trim(StnID(S))),Grid(S),SN(S),WN(S),Bath(SN(S)),Bath(WN(S)), &
         StationsLat(S),StationsLon(S),GridLat(SN(S)),GridLon(SN(S)),GridLat(WN(S)),GridLon(WN(S))
666   Format(A50,3i15,8f15.10)
      ENDDO
      close(86)

      ! Read Global File
      WRITE(*,'(A)',ADVANCE="YES") "ADCIRC grid file was read successfully!"
      WRITE(*,'(A)',ADVANCE="YES") " "
      if (cntrl.eq.0) then
        WRITE(*,'(A)',ADVANCE="NO") "Enter name of file containing global values: "
        READ(*,'(A)',ADVANCE="YES") GlobalFile
        WRITE(*,'(A)',ADVANCE="NO") "Enter Run Tag: "
        READ(*,'(A)',ADVANCE="YES") RunID
      else
        GlobalFile=fort63u
      endif
      OPEN(UNIT=12,FILE=TRIM(GlobalFile),ACTION="READ")
      READ(12,'(A)') JunkC
      READ(12,'(A)') Line
      READ(Line,*) NumSnaps, JunkI, JunkR, Interval, nCol, Comments
      WRITE(*,'(A)',ADVANCE="YES") " "
      XX1 = len_trim(StationsFile)
      if (StationsFile(XX1-7:XX1).eq.'Stations') then
        if (nCol.eq.1) then
          Outputfile = StationsFile(1:XX1-9)//'.'//trim(RunID)//'.61'
        else
          Outputfile = StationsFile(1:XX1-9)//'.'//trim(RunID)//'.62'
        endif 
      else      
        if (nCol.eq.1) then
          Outputfile = TRIM(StationsFile)//'.'//trim(RunID)//'.61'
        else
          Outputfile = TRIM(StationsFile)//'.'//trim(RunID)//'.62'
        endif
      endif
      OPEN(UNIT=14,FILE=TRIM(OutputFile),ACTION="WRITE")

      WRITE(14,'(A)') "Interpolated values at stations"
      WRITE(14,'(5I11,1x,A50)') NumSnaps, NumStations, Interval, Interval, nCol, Comments
      WRITE(*,'(A)',ADVANCE="YES") " "

      DO SS=1,NumSnaps

         WRITE(*,'(A,I5,A,I5,A)',ADVANCE="YES") "Processing snap #",SS," of ",NumSnaps,"."

         ALLOCATE(Global1(1:NumNodes),Global2(1:NumNodes))

         READ(12,'(A)',ERR=999,END=999) Line
         
         Read(Line,*) SnapR, SnapI
         Read(Line,*,ERR=907,END=907) SnapR, SnapI, NumNodesNonDefault, DefaultValue
         Goto 908
907      NumNodesNonDefault = NumNodes
         DefaultValue = -99999.0d0       
908      WRITE(14,*) SnapR, SnapI
         DO N=1,NumNodes
            Global1(N) = DefaultValue
            Global2(N) = DefaultValue
         ENDDO
         if (nCol.eq.2) then
            DO N=1,NumNodesNonDefault
               READ(12,*,ERR=999,END=999) JunkI, StationsTemp1,StationsTemp2
               Global1(JunkI) = StationsTemp1
               Global2(JunkI) = StationsTemp2
            ENDDO
            DO S=1,NumStations
              IF ((UseNode(S).eq.0).and.(UseDepth(S).eq.-999.0d0)) THEN
                W1 = Weight(S,1)
                W2 = Weight(S,2)
                W3 = Weight(S,3)
                X1 = Global1(Conn(StationsElem(S),1))
                X2 = Global1(Conn(StationsElem(S),2))
                X3 = Global1(Conn(StationsElem(S),3))
                Y1 = Global2(Conn(StationsElem(S),1))
                Y2 = Global2(Conn(StationsElem(S),2))
                Y3 = Global2(Conn(StationsElem(S),3))
                if (X1.eq.DefaultValue) then
                  if (X2.eq.DefaultValue) then
                    if (X3.eq.DefaultValue) then
                      StationsTemp1 = DefaultValue                              
                      StationsTemp2 = DefaultValue                                                    
                    else
                      StationsTemp1 = X3
                      StationsTemp2 = Y3                      
                    endif
                  elseif (X3.eq.DefaultValue) then
                    StationsTemp1 = X2
                    StationsTemp2 = Y2                    
                  else
                    StationsTemp1 = X2*W2/(W2+W3) + X3*W3/(W2+W3)
                    StationsTemp2 = Y2*W2/(W2+W3) + Y3*W3/(W2+W3)                    
                  endif
                elseif (X2.eq.DefaultValue) then
                  if (X3.eq.DefaultValue) then
                    StationsTemp1 = X1
                    StationsTemp2 = Y1                    
                  else
                    StationsTemp1 = X1*W1/(W1+W3) + X3*W3/(W1+W3)
                    StationsTemp2 = Y1*W1/(W1+W3) + Y3*W3/(W1+W3)                    
                  endif
                elseif (X3.eq.DefaultValue) then
                   StationsTemp1 = X1*W1/(W1+W2) + X2*W2/(W1+W2)
                   StationsTemp2 = Y1*W1/(W1+W2) + Y2*W2/(W1+W2)                   
                else
                  StationsTemp1 = X1 * W1 + X2 * W2 + X3 * W3
                  StationsTemp2 = Y1 * W1 + Y2 * W2 + Y3 * W3
                endif
              elseif (UseDepth(S).eq.-999.0d0) then
                StationsTemp1 = Global1(UseNode(S))
                StationsTemp2 = Global2(UseNode(S))
              else !Use Preferred Node
                StationsTemp1 = Global1(WN(S))
                StationsTemp2 = Global2(WN(S))
              endif
              if (Grid(S).eq.0) then
                write(14,'(I10,5X,2E24.10)') S, DefaultValue, DefaultValue
              else
                write(14,'(I10,5X,2E24.10)') S, StationsTemp1, StationsTemp2
              endif
            ENDDO
         ELSE
            DO N=1,NumNodesNonDefault
               READ(12,*,ERR=999,END=999) JunkI, StationsTemp1
               Global1(JunkI) = StationsTemp1
            ENDDO
            DO S=1,NumStations
              IF ((UseNode(S).eq.0).and.(UseDepth(S).eq.-999.0d0)) THEN
                if (Grid(S).eq.1) then
                  X1 = Global1(Conn(StationsElem(S),1))
                  X2 = Global1(Conn(StationsElem(S),2))
                  X3 = Global1(Conn(StationsElem(S),3))
                  W1 = Weight(S,1)
                  W2 = Weight(S,2)
                  W3 = Weight(S,3)
                  if (X1.eq.DefaultValue) then
                    if (X2.eq.DefaultValue) then
                      if (X3.eq.DefaultValue) then
                        StationsTemp1 = DefaultValue                              
                      else
                        StationsTemp1 = X3
                      endif
                    elseif (X3.eq.DefaultValue) then
                      StationsTemp1 = X2
                    else
                      StationsTemp1 = X2*W2/(W2+W3) + X3*W3/(W2+W3)
                    endif
                  elseif (X2.eq.DefaultValue) then
                    if (X3.eq.DefaultValue) then
                      StationsTemp1 = X1
                    else
                      StationsTemp1 = X1*W1/(W1+W3) + X3*W3/(W1+W3)
                    endif
                  elseif (X3.eq.DefaultValue) then
                    StationsTemp1 = X1*W1/(W1+W2) + X2*W2/(W1+W2)
                  else
                    StationsTemp1 = X1 * W1 + X2 * W2 + X3 * W3
                  endif
                else
                  StationsTemp1 = Global1(SN(S))
                endif
              elseif (UseDepth(S).eq.-999.0d0) then
                StationsTemp1 = Global1(UseNode(S))
              else
                StationsTemp1 = Global1(WN(S))
              end if
              if (Grid(S).eq.0) then
                WRITE(14,'(I10,5X,1E24.10)') S, DefaultValue
              else
                WRITE(14,'(I10,5X,1E24.10)') S, StationsTemp1
              endif
            ENDDO
         ENDIF
         DEALLOCATE(Global1)
         DEALLOCATE(Global2)
      ENDDO

 999  CONTINUE

      CLOSE(UNIT=12,STATUS="KEEP")

      WRITE(*,'(A)',ADVANCE="YES") "Time series extracted successfully!"

      CLOSE(UNIT=14,STATUS="KEEP")

      DEALLOCATE(StationsElem)
      DEALLOCATE(StationsLat)
      DEALLOCATE(StationsLon)
      DEALLOCATE(Conn)
      DEALLOCATE(Weight)

   end subroutine TSextract
!----------------------------------------------------------------------
!     11
!----------------------------------------------------------------------
   Subroutine wetdry

       implicit none
       character(120)       :: file14, file63, filewet, junkC
       real(8)              :: junkR
       real(8), allocatable :: lat(:), Lon(:)
       integer, allocatable :: wetN(:), Conn(:,:)
       integer              :: i, j, junkI, NumNodes, NumElems, EC, wettag
       integer              :: N, E, NumSnaps, NumNodesNonDefault, SS

       write(6,*) "Enter the name of the fort.14 file"
       read(5,*) file14
       write(6,*) "Enter the name of the fort.63 file" 
       read(5,*) file63
       write(6,*) "Enter the wettag value"
       read(5,*) wettag      
       filewet = 'wet.14'

       open(10,file=file14,status='old',action='read')
         read(10,*)
         read(10,*) NumElems, NumNodes
         allocate(Lon(1:NumNodes),Lat(1:NumNodes))
         allocate(Conn(1:NumElems,1:3))
         allocate(WetN(1:NumNodes))
         do N = 1,NumNodes
           read(10,*) junkI, Lon(N), Lat(N)
         enddo
         do E = 1,NumElems
           read(10,*) junkI, junkI, Conn(E,1), Conn(E,2), Conn(E,3)
         enddo
       close(10)

       open(20,file=file63,status='old',action='read')
         read(20,*)
         read(20,*) NumSnaps
         do SS = 1,NumSnaps
           write(6,'(A,I5,A,I5,A)') "Processing snap #",SS," of ",NumSnaps,"."
           read(20,*,END=999) junkR, junkI, NumNodesNonDefault
           do j = 1,NumNodesNonDefault
             read(20,*) N
             wetN(N) = wettag
           enddo
         enddo
999      continue
       close(20)

       open(30,file=filewet,status='replace',action='write')
         write(30,'(A)') "Wetness Extent fort.14 style Grid"
         write(30,300) NumElems, NumNodes
         do N = 1,NumNodes
           write(30,301) N, Lon(N), Lat(N), WetN(N)
         enddo
         EC=3
         do E = 1,NumElems
           write(30,302) E, EC, Conn(E,1), Conn(E,2), Conn(E,3)
         enddo
300    Format(2i10) 
301    Format(i10,2f15.6,i5)
302    Format(5i10)
       close(30)

   end subroutine wetdry
!----------------------------------------------------------------------
!     12
!----------------------------------------------------------------------
   subroutine findminmax(os)
      implicit none
      character(10),intent(inout)      :: os
      character(120)       :: file14, GlobalFile, junkC, Line, Comments
      real(8)              :: junkR, SnapR, DefaultValue
      real(8), allocatable :: lat(:), Lon(:), maxR(:), minR(:), TempR(:)
      integer, allocatable :: wetN(:), Conn(:,:)
      integer              :: i, j, V, N, junkI, NumNodes, nCol
      integer              :: Interval, NumElems
      integer              :: NumSnaps, NumNodesNonDefault, SS, SnapI

      write(6,'(A)',ADVANCE="NO") "Enter the name of the fort.14 file: "
      read(5,'(A)',ADVANCE="YES") file14
      write(6,'(A)',ADVANCE="NO") "Enter the name of the global file: " 
      read(5,'(A)',ADVANCE="YES") GlobalFile      

      open(10,file=trim(file14),status='old',action='read')
        read(10,*)
        read(10,*) NumElems, NumNodes
        allocate(Lon(1:NumNodes),Lat(1:NumNodes))
        do N = 1,NumNodes
          read(10,*) junkI, Lon(N), Lat(N)
        enddo
      close(10)

      OPEN(UNIT=12,FILE=TRIM(GlobalFile),ACTION="READ")
      READ(12,'(A)') JunkC
      READ(12,'(A)') Line
      READ(Line,*) NumSnaps, JunkI, JunkR, Interval, nCol, Comments
      Allocate(maxR(1:nCol),minR(1:nCol),TempR(1:nCol))
      write(6,*)
      DO SS=1,NumSnaps
        WRITE(*,'(A,I5,A,I5,A)') "Processing snap #",SS," of ",NumSnaps,"."
        READ(12,'(A)',ERR=999,END=999) Line
        Read(Line,*) SnapR, SnapI
        Read(Line,*,ERR=907,END=907) SnapR, SnapI, NumNodesNonDefault
        Goto 908
907     NumNodesNonDefault = NumNodes
908     DO N=1,NumNodesNonDefault
          READ(12,*,ERR=999,END=999) JunkI, (TempR(V), V=1,nCol)
          DO V=1,nCol
            if ((TempR(V).gt.maxR(V)).and.(TempR(V).ne.-99999.0d0)) maxR(V) = TempR(V) 
            if ((TempR(V).lt.minR(V)).and.(TempR(V).ne.-99999.0d0)) minR(V) = TempR(V)            
          enddo
        ENDDO        
      enddo
      write(6,*)
      write(6,*) "Results:"
      write(6,219) "Vector", "Maximum", "Minimum"
999   DO V=1,nCol
        write(6,220) V, maxR(V), minR(V)
      enddo
219   Format(3a10)
220   Format(i10,2f10.4) 
   end subroutine findminmax
!----------------------------------------------------------------------
!     13
!----------------------------------------------------------------------
   subroutine meow(os)
      implicit none
      character(10),intent(inout)      :: os
      real*8           :: surge,r1,r2
      integer          :: nn,node,c1,c3,c4,c5
      integer          :: nf,tag,i,j
      character(60)    :: title, maxele
      character(35)    :: title35
      character(60)    :: maxele_filelist,output_file
      real*8,dimension(:),allocatable :: Z

      write(6,*)
      write(6,*) "Maximum Envelope of Water (MEOW) Compiler Program"
      write(6,*)
      write(6,*) "Enter the name of the MEOW Output File: "       
      read(5,*) output_file
      write(6,*) "Enter the name of the MEOW Input File List: "
      read(5,*)  maxele_filelist
      write(6,*) "Enter the number of Files in the Input File List: "
      read(5,*) nf
      open(1,file=output_file,action="write")           ! output file
      open(700,file=maxele_filelist,action="read")       ! list of all max elevation files

101   format(I10,E22.10)
      do i=1,nf
        read(700,*) maxele
        write(*,*) i,nf,maxele
        open(700+i,file=maxele)
          read(700+i,*) title
          read(700+i,*) c1,nn,r1,c3,c4
          read(700+i,*) r2,c5
          if (i==1) then
            allocate(Z(nn))
            do j=1,nn
              Z(j)=-99999d0
            end do
            write(1,*) title
            write(1,*) c1,nn,r1,c3,c4
            write(1,*) r2,c5
          end if
          write(*,*) "Reading the file:  ", title
          do j=1,nn
            read(700+i,*) node,surge
            if (surge>Z(j)) then
              Z(j)=surge
            end if
          end do
        close(700+i)
      end do
      do i=1,nn
        write(1,101) i,Z(i)
      end do
      close(1)
      close(700)
   end subroutine meow
!----------------------------------------------------------------------
!     14
!----------------------------------------------------------------------
   subroutine generatemax(os)
      implicit none
      character(10),intent(inout)      :: os
      character(120)       :: file14, GlobalFile, junkC, Line, Comments
      character(120)       :: Maxfile
      real(8)              :: junkR, SnapR, DV
      real(8), allocatable :: lat(:), Lon(:), maxR(:,:), TempR(:)
      integer, allocatable :: wetN(:), Conn(:,:)
      integer              :: i, j, V, N, junkI, NumNodes, nCol
      integer              :: Interval, NumElems
      integer              :: NumSnaps, NumNodesNonDefault, SS, SnapI

      write(6,'(A)',ADVANCE="NO") "Enter the name of the fort.14 file: "
      read(5,'(A)',ADVANCE="YES") file14
      write(6,'(A)',ADVANCE="NO") "Enter the name of the global file: " 
      read(5,'(A)',ADVANCE="YES") GlobalFile      
      write(6,'(A)',ADVANCE="NO") "Enter the name of the max file: " 
      read(5,'(A)',ADVANCE="YES") MaxFile            

      open(10,file=trim(file14),status='old',action='read')
        read(10,*)
        read(10,*) NumElems, NumNodes
        allocate(Lon(1:NumNodes),Lat(1:NumNodes))
        do N = 1,NumNodes
          read(10,*) junkI, Lon(N), Lat(N)
        enddo
      close(10)

      OPEN(UNIT=12,FILE=TRIM(GlobalFile),ACTION="READ")
      OPEN(UNIT=13,FILE=TRIM(MaxFile),Status='replace',ACTION="write")
      READ(12,'(A)') JunkC
      write(13,'(A)') JunkC
      READ(12,'(A)') Line
      READ(Line,*) NumSnaps, JunkI, JunkR, Interval, nCol, Comments
      write(13,131) 1, NumNodes, 1, 1, nCol, trim(Comments)
      write(13,132) 0, 0, NumNodes, -9999.0d0
131   Format(2i10, e10.6, 2i10, a100)
132   Format(e10.6,2i10,e10.6)
      Allocate(maxR(1:NumNodes,1:nCol),TempR(1:nCol))
      Do N=1,NumNodes
        do V=1,nCol
          maxR(N,V) = -99999.0d0
        enddo
      enddo
      write(6,*)
      DO SS=1,NumSnaps
        WRITE(*,'(A,I5,A,I5,A)') "Processing snap #",SS," of ",NumSnaps,"."
        READ(12,'(A)',ERR=999,END=999) Line
        Read(Line,*) SnapR, SnapI
        Read(Line,*,ERR=907,END=907) SnapR, SnapI, NumNodesNonDefault, DV
        Goto 908
907     NumNodesNonDefault = NumNodes
908     DO N=1,NumNodesNonDefault
          READ(12,*,ERR=999,END=999) i, (TempR(V), V=1,nCol)
          DO V=1,nCol
            if ((TempR(V).gt.maxR(i,V)).and.(TempR(V).ne.-99999.0d0)) then
              maxR(i,V) = TempR(V)
            endif 
          enddo
        ENDDO        
      enddo
      write(6,*)
999   Do N=1,NumNodes
          write(6,*) N, (maxR(N,V),V=1,nCol)
      enddo
219   Format(3a10)
   end subroutine generatemax
!----------------------------------------------------------------------
!     15
!----------------------------------------------------------------------
   subroutine projector(os)
      character(10),intent(inout)      :: os



   end subroutine projector
!----------------------------------------------------------------------
!     16
!----------------------------------------------------------------------
   subroutine generatediff(os)

       implicit none
       character(120)         :: file1,file2,file3,junkC
      character(10),intent(inout)      :: os
       real(8)                :: wse1,wse2,wse3
       integer                :: i

       write(6,*)
       write(6,*) "************************************************"
       write(6,*) "   This Program will Subtract File2 from File1"
       write(6,*) "************************************************"
       write(6,'(A)',ADVANCE="NO") "Enter Filename 1: "
       read(5,'(A)',ADVANCE="YES") file1
       write(6,'(A)',ADVANCE="NO") "Enter Filename 2: "
       read(5,'(A)',ADVANCE="YES") file2
       write(6,'(A)',ADVANCE="NO") "Enter Output File: "
       read(5,'(A)',ADVANCE="YES") file3
    
       open(10,file=file1,status='old',action='read')
       open(20,file=file2,status='old',action='read')
       open(30,file=file3,status='replace',action='write')
         do i = 1,3
           read(10,'(A)') junkC
           read(20,'(A)')
           write(30,'(A)') junkC
         enddo
         write(6,*) junkC
         do
           read(10,*,END=999) i, wse1
           read(20,*) i, wse2
           if ((wse1.eq.-99999.0d0).and.(wse2.eq.-99999.0d0)) then
             wse3  = -99999.0d0
           elseif(wse1.eq.-99999.0d0) then
             wse1 = 0.0d0
             wse3 = wse1 - wse2
           elseif(wse2.eq.-99999.0d0) then
             wse2 = 0.0d0
             wse3 = wse1 - wse2
           else
             wse3 = wse1 - wse2        
           endif
           write(30,100) i, wse3
         enddo
100    Format(i10,e22.10)
999    close(10)
       close(20)
       close(30)
   end subroutine generatediff
!----------------------------------------------------------------------
!     17
!----------------------------------------------------------------------     
   subroutine listnodes(os)

      implicit none
      character(120)                :: StationsFile,NodesFile,PFile
      character(120)                :: StnID
      character(120), allocatable   :: StationName(:)
      character(10), intent(in)     :: os
      integer                       :: extension, X1, NumNodes, i, S
      integer                       :: NumStations, UseNode
      integer, allocatable          :: NodeNumber(:)
      real (8)                      :: StationsLon, StationsLat

      WRITE(*,'(A)',ADVANCE="YES") "Enter name of stationfile: "
      extension = 2
      call choosefile(StationsFile,extension,os)
      write(6,*)

      WRITE(*,'(A)',ADVANCE="YES") "Enter name of preferred node file: "
      extension = 4
      call choosefile(NodesFile,extension,os)
      write(6,*)

      X1 = len_trim(StationsFile)      
      PFile = StationsFile(1:X1-8)//'P.Stations'


      OPEN(UNIT=21,FILE=TRIM(NodesFile),ACTION="READ")
        read(21,*) NumNodes
        allocate(StationName(1:NumNodes),NodeNumber(1:NumNodes))
        do i = 1,NumNodes
          read(21,*,END=999) StationName(i),NodeNumber(i)
        enddo
999   continue
      close(21)


      open(UNIT=11,FILE=TRIM(StationsFile),ACTION="READ")
      open(UNIT=31,FILE=TRIM(PFile),status='replace',ACTION="write")      
        read(11,*) NumStations
        write(31,*) NumStations
        do S=1,NumStations
           read(11,*) StnID, StationsLat, StationsLon
           UseNode = 0
           do i = 1,NumNodes
             if (trim(StnID).eq.trim(StationName(i))) UseNode = i
           enddo
           if (UseNode.gt.0) then
             write(31,*) trim(StnID),StationsLat,StationsLon,NodeNumber(UseNode)
           else
             write(31,*) trim(StnID),StationsLat,StationsLon
           endif             
        enddo
      close(UNIT=11,STATUS="KEEP")
      close(UNIT=31,STATUS="KEEP")      

   end subroutine listnodes
!----------------------------------------------------------------------
!     18
!----------------------------------------------------------------------     
   subroutine relocator(os)

      implicit none
      character(120)                :: StationsFile,RelocatorFile,LFile
      character(120)                :: StnID
      character(120), allocatable   :: StationName(:)
      character(10), intent(in)     :: os
      integer                       :: extension, X1, NumNodes, i, S
      integer                       :: NumStations, UseNode
      integer, allocatable          :: NodeNumber(:)
      real (8)                      :: StationsLon, StationsLat
      real (8), allocatable         :: NewLon(:), NewLat(:)      

      WRITE(*,'(A)',ADVANCE="YES") "Enter name of stationfile: "
      extension = 2
      call choosefile(StationsFile,extension,os)
      write(6,*)

      WRITE(*,'(A)',ADVANCE="YES") "Enter name of relocation file: "
      extension = 4
      call choosefile(RelocatorFile,extension,os)
      write(6,*)

      X1 = len_trim(StationsFile)      
      LFile = StationsFile(1:X1-8)//'L.Stations'


      OPEN(UNIT=21,FILE=TRIM(RelocatorFile),ACTION="READ")
        read(21,*) NumNodes
        allocate(StationName(1:NumNodes),NewLat(1:NumNodes),NewLon(1:NumNodes))
        do i = 1,NumNodes
          read(21,*,END=999) StationName(i),NewLat(i),NewLon(i)
        enddo
999   continue
      close(21)


      open(UNIT=11,FILE=TRIM(StationsFile),ACTION="READ")
      open(UNIT=31,FILE=TRIM(LFile),status='replace',ACTION="write")      
        read(11,*) NumStations
        write(31,*) NumStations
        do S=1,NumStations
           read(11,*) StnID, StationsLat, StationsLon
           UseNode = 0
           do i = 1,NumNodes
             if (trim(StnID).eq.trim(StationName(i))) UseNode = i
           enddo
           if (UseNode.gt.0) then
             write(31,*) trim(StnID),NewLat(UseNode),NewLon(UseNode)
           else
             write(31,*) trim(StnID),StationsLat,StationsLon
           endif             
        enddo
      close(UNIT=11,STATUS="KEEP")
      close(UNIT=31,STATUS="KEEP")      

   end subroutine relocator   
!----------------------------------------------------------------------
!     21
!----------------------------------------------------------------------   
   subroutine CRMS2IMEDS

      implicit none
      character(120)                :: filestations, fileCRMS, fileIMEDS, tempC
      character(1000)               :: templine
      character(120)                :: stnid, pstnid, sdate, stime, fstnid
      integer                       :: iflag, i , j, k, nsta 
      integer                       :: yyi, moi, ddi, hhi, mii, ssi
      integer                       :: yyf, mof, ddf, hhf, mif, ssf
      integer                       :: yyc, moc, ddc, hhc, mic, ssc
      integer                       :: defval=-999, nstid, TZC
      double precision              :: vadj, tadj
      integer, allocatable          :: X(:)
      integer,allocatable           :: calen(:), cd(:)
      double precision              :: IDATE, CDATE, FDATE, lat, lon, adjetd

      write(6,*) "What is the initial date? (Year, Month, Day, Hour, Min, Sec)"
      read(5,*) yyi, moi,ddi,hhi,mii,ssi
      write(6,*) "What is the final date? (Year, Month, Day, Hour, Min, Sec)"
      read(5,*) yyf, mof,ddf,hhf,mif,ssf
      write(6,*)
      write(6,*) "Enter the name of the CRMS input file"
      read(5,*) fileCRMS
      write(6,*) "Enter the name of the CRMS stations file"
      read(5,*) filestations
      write(6,*)
      Allocate(X(1:17))

      !Check for Cross Years
      if (yyi.ne.yyf) then
        write(6,*) "Cross-years currently not permitted with this release"
        stop
      endif
      allocate(calen(1:12), cd(1:13))
      calen = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
      cd = (/ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 /)
      if( mod(yyi,4) == 0 ) then
        calen(2) = 29 !leap year
        cd = (/ 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366 /)
      endif
      IDATE = cd(moi) + ddi + hhi/24.0d0 + mii/1440.0d0 + ssi/86400.0d0
      FDATE = cd(mof) + ddf + hhf/24.0d0 + mif/1440.0d0 + ssf/86400.0d0

      open(unit=23,file='notfound.txt',status='replace',action='write')

      open(unit=21,file='watlev_CRMS.IMEDS',status='replace',action='write')
      write(21,*) "% IMEDS generic format - Water Level"
      write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(21,*) "CRMS " , "UTC", " NAVD88"

      open(unit=20,file=fileCRMS,status='old',action='read')
      read(20,*)
      nsta = 0
      pstnid = "0"
      do 
        read(20,'(A1000)',ERR=999,END=999) templine
        X(1) = index(templine,",")
        do i = 2,17
          X(i) =X(i-1) + index(templine(X(i-1)+1:1000),",")
        enddo
        stnid = templine(1:X(1)-1)
        stnid = trim(AdjustL(stnid))
        sdate = templine(X(1)+1:X(2)-1)
        stime = templine(X(2)+1:X(3)-1)
        tempC = TRIM(templine(X(15)+1:X(16)-1))
        if (len_trim(tempC).eq.0) then
          adjetd = 0.0d0
        else 
          read(tempC(:),'(f10.5)') adjetd
        endif
          read(sdate(1:2),'(i10)') moc
          read(sdate(4:5),'(i10)') ddc
          read(sdate(7:10),'(i10)') yyc
          read(stime(1:2),'(i10)') hhc
          read(stime(4:5),'(i10)') mic
          read(stime(7:8),'(i10)') ssc
        CDATE = cd(moc) + ddc + hhc/24.0d0 + mic/1440.0d0 + ssc/86400.0d0
        if ((CDATE.ge.IDATE).and.(CDATE.le.FDATE).and.(yyc.eq.yyi)) then
          if (stnid.ne.pstnid) then
            write(6,*) "Processing Station Id: ", trim(stnid)
            pstnid=stnid
            !Search for Latitude,Longitude
            open(unit=22,file=filestations,status='old',action='read')             
              read(22,*) nstid
              read(22,*)
              lat = 0.0d0
              lon = 0.0d0
              vadj = -9999.0d0
              tadj = -9999.0d0              
              do i = 1,nstid
                read(22,*) fstnid, lon, lat, vadj, tadj
                fstnid = AdjustL(fstnid)
                if ((trim(fstnid)).eq.(trim(stnid))) goto 995
              enddo
              write(6,*) trim(stnid), " Could not be found"
              write(23,*) trim(stnid)
995         continue
            close(22)
            nsta = nsta + 1
            write(21,*) "CRMS_"//TRIM(stnid), lat, lon
          endif
          adjetd = 0.3048d0 * adjetd         ! Convert Feet to Meters
          if (vadj.ne.-9999.0d0) adjetd = adjetd + vadj
          if (tadj.ne.-9999.0d0) then
            TZC = tadj
          else
            TZC = 5   !CDT (Between 2nd SUN of March to 1st SUN of Nov) 
          endif 
          Call TZconvert(yyc,moc,ddc,hhc,TZC) !Adjust for Central Time to UTC
          write(21,220) yyc, moc, ddc, hhc, mic, ssc, adjetd
        endif
      enddo
220   Format(6i7,F14.4)
999   continue

      close(23)
      close(21)
      close(20)
      write(6,*) "Processed ", nsta, " Stations"

   end subroutine CRMS2IMEDS
!----------------------------------------------------------------------
!     22
!----------------------------------------------------------------------
   subroutine NOAA2IMEDS

      implicit none

      character(120)                :: filestations
      character(120)                :: sdate, stime, tempC, junkC
      character(200)                :: templine, templine2
      character(120),allocatable    :: stnid(:), datum(:), filename(:),state(:)
      integer                       :: iflag, i , j, k, nsta, iver 
      integer                       :: yy, mo, dd, hh, mi, ss
      double precision, allocatable :: lat(:), lon(:), conv(:)
      double precision              :: pred, ver
      integer                       :: defval=-999

      write(6,*) "Enter the name of the station list"
      read(5,*) filestations
!      filestations = "NOAA_Stations.txt"

! Read Station Data
      open(unit=20,file=filestations,status='old',action='read')
        read(20,*) nsta
        read(20,*)
        allocate(lat(1:nsta), lon(1:nsta),filename(1:nsta), conv(1:nsta))
        allocate(stnid(1:nsta), state(1:nsta))
        do i=1,nsta
          read(20,*) stnid(i), filename(i), state(i), lat(i), lon(i), conv(i)
        enddo
      close(20)

! Read & Write Raw Data to IMEDS format

      open(unit=21,file="watlev_NOAA_NAVD_VER.IMEDS",status='replace',action='write')
      write(21,*) "% IMEDS generic format - Water Level"
      write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(21,*) "NOAA " , "UTC", " NAVD88"

      open(unit=22,file="watlev_NOAA_MSL_VER.IMEDS",status='replace',action='write')
      write(22,*) "% IMEDS generic format - Water Level"
      write(22,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(22,*) "NOAA " , "UTC", " MSL"

      open(unit=23,file="watlev_NOAA_NAVD_PRE.IMEDS",status='replace',action='write')
      write(23,*) "% IMEDS generic format - Water Level"
      write(23,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(23,*) "NOAA " , "UTC", " NAVD88"

      open(unit=24,file="watlev_NOAA_MSL_PRE.IMEDS",status='replace',action='write')
      write(24,*) "% IMEDS generic format - Water Level"
      write(24,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(24,*) "NOAA " , "UTC", " MSL"

      iflag = 1
      ss = 0
      do i=1,nsta
        write(6,*) "Processing station Id: ", i, trim(stnid(i))
        open(unit=11,file=filename(i),status='old',action='read')
          do j=1,7
            read(11,'(a)') templine
          enddo
          iver = 1
          if (len_trim(templine).lt.35) iver = 0
          if (iver.eq.1) then
              if (conv(i).ne.-9999.0d0) then
              write(21,*) "NOAA_"//trim(stnid(i)), lat(i), lon(i)
            else
              write(22,*) "NOAA_"//trim(stnid(i)), lat(i), lon(i)
            endif
          else
              if (conv(i).ne.-9999.0d0) then
              write(23,*) "NOAA_"//trim(stnid(i)), lat(i), lon(i)
            else
              write(24,*) "NOAA_"//trim(stnid(i)), lat(i), lon(i)
            endif
          endif
          do while (iflag.eq.1)
            read(11,'(a)',END=996) templine
            if (len_trim(templine).eq.0) goto 996
            read(templine(9:12),'(i10)') yy
            read(templine(13:14),'(i10)') mo
            read(templine(15:16),'(i10)') dd
            read(templine(18:19),'(i10)') hh
            read(templine(21:22),'(i10)') mi
            templine2 = templine(24:len_trim(templine))
            if (iver.eq.1) then
              read(templine2,221,ERR=998) pred, ver
              if (conv(i).ne.-9999.0d0) then
                ver = ver-conv(i)
                write(21,220) yy, mo, dd, hh, mi, ss, ver
              else
                write(22,220) yy, mo, dd, hh, mi, ss, ver
              endif
              goto 997
998           if (conv(i).ne.-9999.0d0) then
                ver = ver-conv(i)        
                write(21,220) yy, mo, dd, hh, mi, ss, -999
              else
                write(22,220) yy, mo, dd, hh, mi, ss, -999
              endif
            else
              read(templine2,*) pred
              if (conv(i).ne.-9999.0d0) then
                pred = pred-conv(i)                      
                write(23,220) yy, mo, dd, hh, mi, ss, pred
              else
                write(24,220) yy, mo, dd, hh, mi, ss, pred
              endif              
            endif
997         continue
          enddo
996       continue
        close(11)
      enddo
      write(6,*) "FINISHED"
220   Format(6i7,F14.4)
221   Format(f7.3,1x,f7.3)
222   Format(6i6,i10)
      close(21)
      close(22)
      close(23)
      close(24)
      close(25)

   end subroutine NOAA2IMEDS
!----------------------------------------------------------------------
!     23
!----------------------------------------------------------------------
   subroutine TCOON2IMEDS
 
      implicit none
      character(120)                :: filestations
      character(120)                :: sdate, stime
      character(120),allocatable    :: stnid(:), filename(:), desc(:), datum(:)
      integer                       :: iflag, i , j, k, nsta
      integer                       :: yy, mo, dd, hh, mi, ss, dy
      double precision, allocatable :: lat(:), lon(:), conv(:)
      double precision              :: wse, junkD

      write(6,*) "Enter the name of the TCOON station list"
      read(5,*) filestations

! Read Station Data
      open(unit=20,file=filestations,status='old',action='read')
      read(20,*) nsta
      read(20,*)
      allocate(lat(1:nsta), lon(1:nsta),filename(1:nsta), conv(1:nsta))
      allocate(stnid(1:nsta),desc(1:nsta),datum(1:nsta))
      do i=1,nsta
        read(20,*) stnid(i), desc(i), filename(i), lat(i), lon(i), datum(i), conv(i)
      enddo
      close(20)

! Read & Write Raw Data to IMEDS format

      open(unit=21,file="watlev_TCOON_NAVD88.IMEDS",status='replace',action='write')
      write(21,*) "% IMEDS generic format - Water Level (surge)"
      write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(21,*) "TCOON " , "UTC", " NAVD88"

      iflag = 1
      ss = 0
      do i=1,nsta
        write(6,*) "Processing station Id: ", i, trim(stnid(i))
        open(unit=11,file=filename(i),status='old',action='read')
          do j=1,19
            read(11,*)
          enddo
          write(21,*) "TCOON_"//trim(stnid(i)), lat(i), lon(i)
          do while (iflag.eq.1)
            read(11,*,ERR=998, END=999) sdate, stime, wse

            read(sdate(1:4),'(i10)') yy
            read(sdate(5:7),'(i10)') dy
            read(stime(1:2),'(i10)') hh
            read(stime(3:4),'(i10)') mi
            if (conv(i).ne.-9999) wse = wse + conv(i)
            call getmonthday(dy,yy,mo,dd)
            write(21,200) yy, mo, dd, hh, mi, ss, wse
998         continue            
          enddo
999       continue
        close(11)
      enddo
200   Format(6i7,F14.4)      
      close(21)
      close(22)
   
   end subroutine TCOON2IMEDS
!----------------------------------------------------------------------
!     24
!----------------------------------------------------------------------
   subroutine USACE2IMEDS

      implicit none
      character(120)                :: filestations, tempC
      character(120)                :: sdate, stime
      character(120),allocatable    :: stnid(:), filename(:), datum(:)
      integer                       :: iflag, i , j, k, nsta
      integer                       :: X1, X2, X3,TZC=5
      integer                       :: yy, mo, dd, hh, mi, ss
      double precision, allocatable :: lat(:), lon(:), conv(:)
      double precision              :: wse, junkD

      write(6,*) "Enter the name of the USACE station list"
      read(5,*) filestations
!      filestations = "USACE_Stations.txt"

! Read Station Data
      open(unit=20,file=filestations,status='old',action='read')
      read(20,*) nsta
      read(20,*)
      allocate(lat(1:nsta), lon(1:nsta), filename(1:nsta))
      allocate(stnid(1:nsta), conv(1:nsta), datum(1:nsta))
      do i=1,nsta
        read(20,*) stnid(i), lat(i), lon(i), filename(i), datum(i), conv(i)
      enddo
      close(20)

! Read & Write Raw Data to IMEDS format

      open(unit=21,file="watlev_USACE.IMEDS",status='replace',action='write')
      write(21,*) "% IMEDS generic format - Water Level"
      write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(21,*) "USACE " , "UTC", " NAVD88"

      ss = 0
      do i=1,nsta
        write(6,*) "Processing station Id: ", i, trim(stnid(i))
        open(unit=11,file=filename(i),status='old',action='read')
          write(21,*) "USACE_"//trim(stnid(i)), lat(i), lon(i)
          do
            read(11,'(a)',END=999) tempC
            X1 = index(tempC,"/")
            X2 = index(tempC(X1+1:len_trim(tempC)),"/")+X1
            X3 = index(tempC,":")
            read(tempC(1:X1-1),'(i10)') mo
            read(tempC(X1+1:X2-1),'(i10)') dd
            read(tempC(X2+1:X2+4),'(i10)') yy
            read(tempC(X2+6:X3-1),'(i10)') hh
            read(tempC(X3+1:X3+2),'(i10)') mi
            read(tempC(X3+3:len_trim(tempC)),*,ERR=998,END=998) wse
            if (conv(i).ne.-9999) wse = wse + conv(i)
            wse = 0.3048d0*wse
            goto 997
998         wse = -999.0d0
997         Call TZconvert(yy,mo,dd,hh,TZC) !Adjust for Central Time to UTC
            write(21,200) yy, mo, dd, hh, mi, ss, wse
          enddo
999       continue
        close(11)
      enddo
200   Format(6i7,F14.4)     
      close(21)
      close(22)
      close(23)

   end subroutine 
!----------------------------------------------------------------------
!     25
!----------------------------------------------------------------------
   subroutine USGSDEPL2IMEDS

      implicit none
      character(120)                :: filestations
      character(120)                :: sdate, stime, junkC, TempC
      character(120),allocatable    :: stnid(:), datum(:), site(:),filename(:)
      integer                       :: iflag, i , j, k, nsta, iopt
      integer                       :: yy, mo, dd, hh, mi, ss, X1, X2, X3
      double precision, allocatable :: lat(:), lon(:), conv(:)
      double precision              :: wse
      integer                       :: TZC=5

      write(6,*) "Enter the name of the USGS-DEPL station list"
      read(5,*) filestations
!      filestations = "USGS-DEPL_Stations.txt"
      write(6,*)
      write(6,*) "Enter 1 for Ike or 2 for Rita"
      read(5,*) iopt
      write(6,*)

! Read Station Data
      open(unit=20,file=filestations,status='old',action='read')
      read(20,*) nsta
      read(20,*)
      allocate(lat(1:nsta), lon(1:nsta),filename(1:nsta))
      allocate(stnid(1:nsta),datum(1:nsta),site(1:nsta),conv(1:nsta))
      do i=1,nsta
        read(20,*) stnid(i), lat(i), lon(i), filename(i), datum(i), conv(i), junkC, site(i)
      enddo
      close(20)

! Read & Write Raw Data to IMEDS format

      open(unit=21,file="watlev_USGS-DEPL_surge.IMEDS",status='replace',action='write')
      write(21,*) "% IMEDS generic format - Water Level (surge)"
      write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(21,*) "USGS-DEPL " , "UTC", " NAVD88"

      open(unit=22,file="watlev_USGS-DEPL_beach.IMEDS",status='replace',action='write')
      write(22,*) "% IMEDS generic format - Water Level (beach)"
      write(22,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(22,*) "USGS-DEPL " , "UTC", " NAVD88"

      open(unit=23,file="watlev_USGS-DEPL_riverine.IMEDS",status='replace',action='write')
      write(23,*) "% IMEDS generic format - Water Level (riverine)"
      write(23,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(23,*) "USGS-DEPL " , "UTC", " NAVD88"
   
      iflag = 1
      ss = 0
      do i=1,nsta
        write(6,*) "Processing station Id: ", i, trim(stnid(i))
        open(unit=11,file=trim(filename(i)),status='old',action='read')
          select case (site(i))
            case ("surge")
              write(21,*) "USGS-DEPL_"//trim(stnid(i)), lat(i), lon(i)
            case("beach")
              write(22,*) "USGS-DEPL_"//trim(stnid(i)), lat(i), lon(i)                    
            case("riverine")
              write(23,*) "USGS-DEPL_"//trim(stnid(i)), lat(i), lon(i)
          end select  
          select case (iopt)
            case (1)
              do j=1,28
                read(11,*)
              enddo
              do while (iflag.eq.1)
                read(11,*,END=999) sdate, stime, wse
                read(sdate(7:10),'(i10)') yy
                read(sdate(1:2),'(i10)') mo
                read(sdate(4:5),'(i10)') dd
                read(stime(1:2),'(i10)') hh
                read(stime(4:5),'(i10)') mi
                Call TZconvert(yy,mo,dd,hh,TZC)
                wse = 0.3048d0*wse         !Convert Feet to Meters
                select case (site(i))
                  case ("surge")
                    write(21,200) yy, mo, dd, hh, mi, ss, wse
                  case("beach")
                    write(22,200) yy, mo, dd, hh, mi, ss, wse
                  case("riverine")
                    write(23,200) yy, mo, dd, hh, mi, ss, wse
                end select  
              enddo
            case (2)
              do j=1,23
                read(11,*)
              enddo
              do while (iflag.eq.1)
                read(11,'(a)',END=999) tempC
                X1 = index(tempC,"/")
                X2 = index(tempC(X1+1:len_trim(tempC)),"/")+X1
                X3 = index(tempC,":")
                read(tempC(1:X1-1),'(i10)') mo
                read(tempC(X1+1:X2-1),'(i10)') dd
                read(tempC(X2+1:X2+4),'(i10)') yy
                read(tempC(X3-2:X3-1),'(i10)') hh
                read(tempC(X3+1:X3+2),'(i10)') mi
                read(tempC(X3+4:X3+5),'(i10)') ss
                if (tempC(X3+7:X3+8).eq."PM") hh = hh + 12
                if ((tempC(X3+7:X3+8).eq."AM").and.(hh.eq.12)) hh = 0
                if (hh.eq.24) hh = 12
                read(tempC(X3+9:len_trim(tempC)),*,ERR=998,END=998) wse
                wse = 0.3048d0*wse
                goto 997
998             wse = -999.0d0
997             continue
                Call TZconvert(yy,mo,dd,hh,TZC)
                write(21,200) yy, mo, dd, hh, mi, ss, wse
              enddo
          end select
999       continue
        close(11)
      enddo
200   Format(6i7,F14.4)
      close(21)
      close(22)

   end subroutine USGSDEPL2IMEDS
!----------------------------------------------------------------------
!     26
!----------------------------------------------------------------------
   subroutine USGSPERM2IMEDS

      implicit none
      character(120)                :: filestations,file1,file2,file3
      character(120)                :: sdate, stime, tzone, Header1
      character(120)                :: junkC, tempC
      character(120),allocatable    :: stnid(:), filename(:), datum(:)
      integer                       :: iflag, i , j, k, nsta
      integer                       :: X1, X2, X3, TZC
      integer                       :: yy, mo, dd, hh, mi, ss
      double precision, allocatable :: lat(:), lon(:), conv(:)
      double precision              :: wse, junkD

      write(6,*) "Enter the name of the USGS-PERM station list"
      read(5,*) filestations
!      filestations = "USGS-PERM_Stations.txt"

! Read Station Data
      open(unit=20,file=filestations,status='old',action='read')
      read(20,*) nsta
      read(20,*)
      allocate(lat(1:nsta), lon(1:nsta), filename(1:nsta))
      allocate(stnid(1:nsta), conv(1:nsta), datum(1:nsta))
      do i=1,nsta
        read(20,*) stnid(i),lat(i),lon(i),filename(i),datum(i),conv(i)
      enddo
      close(20)

! Read & Write Raw Data to IMEDS format
      file1="watlev_USGS-PERM.IMEDS"
      Header1="% IMEDS generic format - Water Level"      

      open(unit=21,file=trim(file1),status='replace',action='write')
      write(21,*) trim(Header1)
      write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(21,*) "USGS-PERM " , "UTC", " NAVD88"

      ss = 0
      do i=1,nsta
        write(6,*) "Processing station Id: ", i, trim(stnid(i))
        open(unit=11,file=filename(i),status='old',action='read')
          write(21,*) "USGS-PERM_"//trim(stnid(i)), lat(i), lon(i)
          do  
            read(11,*) junkC
            if (JunkC(1:1).ne.('#')) goto 667
          enddo
667       read(11,*)          
          do
            read(11,'(A)',ERR=998,END=999) junkC
            junkC = adjustL(junkC)
            if ((JunkC(3:3).eq.'/').or.(JunkC(2:2).eq.'/')) then
              read(junkC(1:19),'(A)') tempC
              read(junkC(20:len_trim(junkC)),*) tzone, wse
              X1 = index(tempC,"/")
              X2 = index(tempC(X1+1:len_trim(tempC)),"/")+X1
              X3 = index(tempC,":")
              read(tempC(1:X1-1),'(i10)') mo
              read(tempC(X1+1:X2-1),'(i10)') dd
              read(tempC(X2+1:X2+4),'(i10)') yy
              read(tempC(X3-2:X3-1),'(i10)') hh
              read(tempC(X3+1:X3+2),'(i10)') mi
              read(tempC(X3+4:X3+5),'(i10)') ss
            else
              read(junkC,*) sdate, stime, tzone, wse
              read(sdate(1:4),'(i10)') yy
              read(sdate(5:6),'(i10)') mo
              read(sdate(7:8),'(i10)') dd
              read(stime(1:2),'(i10)') hh
              read(stime(3:4),'(i10)') mi
              read(stime(5:6),'(i10)') ss
            end if
            if (trim(AdjustL(tzone)).eq."CDT") then
                TZC=5
                Call TZconvert(yy,mo,dd,hh,TZC)
            elseif (trim(AdjustL(tzone)).eq."CST") then
                TZC=6
                Call TZconvert(yy,mo,dd,hh,TZC)
            else                
              write(6,*) tzone, trim(stnid(i))
            endif
            if (conv(i).ne.-9999) wse = wse + conv(i)
            wse  = 0.3048d0*wse    !Convert Feet to Meters
            write(21,200) yy, mo, dd, hh, mi, ss, wse
998         continue            
          enddo
999       continue
        close(11)
      enddo
200   Format(6i7,F14.4)      
      close(21)
      close(22)
      close(23)

   end subroutine USGSPERM2IMEDS
!----------------------------------------------------------------------
!     27
!----------------------------------------------------------------------
   subroutine USACECHL

      implicit none
      character(120)                :: filestations, line
      character(120),allocatable    :: stnid(:), filename1(:)
      character(120),allocatable    :: filename2(:)
      integer                       :: nsta, NumPts, i, j
      integer                       :: yy, mo, dd, hh, mi, ss
      double precision, allocatable :: lat(:), lon(:), conv(:)
      double precision              :: wse, hs, tm, tp

      write(6,*) "Enter the name of the USACE-CHL station list"
      read(5,*) filestations
!      filestations = "USACE-CHL_Stations.txt"

! Read Station Data
      open(unit=20,file=trim(filestations),status='old',action='read')
      read(20,*) nsta
      read(20,*)
      allocate(lat(1:nsta), lon(1:nsta), filename1(1:nsta))
      allocate(stnid(1:nsta), conv(1:nsta),filename2(1:nsta))
      do i=1,nsta
        read(20,*) stnid(i),filename1(i), filename2(i), lat(i),lon(i),conv(i)
      enddo
      close(20)

! Read & Write Raw Data to IMEDS format
      open(unit=21,file="hs_USACE-CHL.IMEDS",status='replace',action='write')
      write(21,*) "% IMEDS generic format - Wave Height"
      write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "hs (m)"
      write(21,*) "USACE-CHL " , "UTC", " NAVD88"

      open(unit=22,file="tp_USACE-CHL.IMEDS",status='replace',action='write')
      write(22,*) "% IMEDS generic format - Peak Period"
      write(22,*) "% ", "year ","month ","day ","hour ","min ","sec ", "tp (s)"
      write(22,*) "USACE-CHL " , "UTC", " NAVD88"

      open(unit=23,file="tm_USACE-CHL.IMEDS",status='replace',action='write')
      write(23,*) "% IMEDS generic format - Mean Period"
      write(23,*) "% ", "year ","month ","day ","hour ","min ","sec ", "tp (s)"
      write(23,*) "USACE-CHL " , "UTC", " NAVD88"

      open(unit=24,file="watlev_USACE-CHL.IMEDS",status='replace',action='write')
      write(24,*) "% IMEDS generic format - Water Level"
      write(24,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(24,*) "USACE-CHL " , "UTC", " NAVD88"      


      do i=1,nsta
        ss = 0
        mi = 0      
        write(6,*) "Processing station Id: ", i, trim(stnid(i))
        write(21,*) "USACE-CHL_"//trim(stnid(i)), lat(i), lon(i)
        write(22,*) "USACE-CHL_"//trim(stnid(i)), lat(i), lon(i)
        write(23,*) "USACE-CHL_"//trim(stnid(i)), lat(i), lon(i)
        write(24,*) "USACE-CHL_"//trim(stnid(i)), lat(i), lon(i)
        open(unit=11,file=filename1(i),status='old',action='read')
          read(11,*) NumPts
          do j = 1,8
            read(11,*)
          enddo
          do j=1,NumPts
            read(11,*) mo, dd, yy, hh, hs, tm, tp, wse
            yy = yy + 2000
            write(21,200) yy, mo, dd, hh, mi, ss, hs
            write(22,200) yy, mo, dd, hh, mi, ss, tp
            write(23,200) yy, mo, dd, hh, mi, ss, tm
          enddo
        close(11)
        open(unit=12,file=filename2(i),status='old',action='read')
          do j = 1,7
            read(12,*)
          enddo
          do j = 1,NumPts
            read(12,'(a)',err=999,end=999) line
            read(line(1:2),*) mo
            read(line(4:5),*) dd
            read(line(7:8),*) yy
            read(line(10:11),*) hh
            read(line(13:14),*) mi
            read(line(16:len_trim(line)),*,ERR=867) wse
            yy = yy + 2000
            wse = wse + conv(i)            
            goto 866
867         wse = -999
866         write(24,200) yy, mo, dd, hh, mi, ss, wse            
          enddo
999     close(12)
      enddo
200   Format(6i7,F14.4)      
      close(21)
      close(22)
      close(23)
      close(24)      

   end subroutine USACECHL
!----------------------------------------------------------------------
!     28
!----------------------------------------------------------------------
   subroutine SLOSH2IMEDS

      implicit none
      character(120)                :: filestations, Line
      character(120),allocatable    :: stnid(:), filename(:)
      integer                       :: nsta, i, j
      integer                       :: yy, mo, dd, hh, mi, ss
      double precision, allocatable :: lat(:), lon(:)
      double precision              :: wse

      write(6,*) "Enter the name of the SLOSH station list"
      read(5,*) filestations

! Read Station Data
      open(unit=20,file=trim(filestations),status='old',action='read')
      read(20,*) nsta
      read(20,*)
      allocate(lat(1:nsta), lon(1:nsta), filename(1:nsta))
      allocate(stnid(1:nsta))
      do i=1,nsta
        read(20,*) stnid(i),lat(i),lon(i),filename(i)
      enddo
      close(20)

! Read & Write Raw Data to IMEDS format
      open(unit=21,file="watlev_SLOSH.IMEDS",status='replace',action='write')
      write(21,*) "% IMEDS generic format - Water Level"
      write(21,*) "% ", "year ","month ","day ","hour ","min ","sec ", "watlev (m)"
      write(21,*) "SLOSH " , "UTC", " NAVD88"
      do i=1,nsta
        write(6,*) "Processing station Id: ", i, trim(stnid(i))
        open(unit=11,file=filename(i),status='old',action='read')
          write(21,*) trim(stnid(i)), lat(i), lon(i)
          do j = 1,4
            read(11,*)
          enddo
          do
            read(11,'(A)',END=999,ERR=999) Line
            read(Line(1:2),'(i10)') mo
            read(Line(4:5),'(i10)') dd
            read(Line(7:10),'(i10)') yy
            read(Line(12:13),'(i10)') hh
            read(Line(15:16),'(i10)') mi
            read(Line(18:19),'(i10)') ss
            read(Line(21:len_trim(Line)),*) wse
            write(21,200) yy, mo, dd, hh, mi, ss, wse
          enddo
999     close(11)
      enddo
200   Format(6i7,F14.4)      
      close(21)
      close(22)
      close(23)
      close(24)      

   end subroutine SLOSH2IMEDS   

!----------------------------------------------------------------------
!     29
!----------------------------------------------------------------------
   subroutine CSI2IMEDS

      implicit none

      integer                       :: nsta, i, j, yy, mo, dd, hh
      integer                       :: dir, mi=0, ss=0, X1, X2
      double precision              :: hs, tm, tp
      double precision, allocatable :: lat(:), lon(:), conv(:)
      character(120), allocatable   :: filename(:), datum(:), stnid(:)
      character(142)               :: tempC, templine, filestation
      character(120)                :: junkC

      write(6,*) "Enter the name of the CSI station list"
      read(5,*) filestation
      !filestation = 'CSI_Stations.txt'

      ! Read Station Data
      open(unit=20,file=trim(filestation),status='old',action='read')
      read(20,*) nsta
      read(20,*)
      allocate(lat(1:nsta), lon(1:nsta),filename(1:nsta), conv(1:nsta))
      allocate(stnid(1:nsta),datum(1:nsta))
      do i=1,nsta
        read(20,*) stnid(i), filename(i), lon(i), lat(i), datum(i), conv(i)
      enddo
      close(20)
        
      !Write Data files (hs, tm, tp, dir)
      open(unit=21,file="hs_CSI.IMEDS",status='replace',action='write')
      write(21,*) "% IMEDS generic format - Wave Height"
      write(21,*) "% ", "year ","month ","day ","hour ","min ", "sec ", "hs (m)"
      write(21,*) "CSI " , "UTC", " NAVD"

      open(unit=22,file="dir_CSI.IMEDS",status='replace',action='write')
      write(22,*) "% IMEDS generic format - Wave Direction"
      write(22,*) "% ", "year ","month ","day ","hour ","min ","sec ", "dir (degrees)"
      write(22,*) "CSI " , "UTC", " NAVD"
      
      open(unit=23,file="tm_CSI.IMEDS",status='replace',action='write')
      write(23,*) "% IMEDS generic format - Average Wave Period"
      write(23,*) "% ", "year ","month ","day ","hour ","min ","sec ", "tp (s)"
      write(23,*) "CSI " , "UTC", " NAVD"

      open(unit=24,file="tp_CSI.IMEDS",status='replace',action='write')
      write(24,*) "% IMEDS generic format - Peak Wave Period"
      write(24,*) "% ", "year ","month ","day ","hour ","min ","sec ", "tp (s)"
      write(24,*) "CSI " , "UTC", " NAVD"

      do i=1,nsta
        write(6,*) "Processing station Id: ", i, trim(stnid(i))      
          write(6,*) filename(i)
        open(unit=11,file=trim(filename(i)),status='old',action='read')
          read(11,*) 
          do j = 21,24
            write(j,*) trim(stnid(i)), lat(i), lon(i)
          enddo       
          do
990         read(11,'(a)',END=999) tempC
            if (tempC(1:1).eq.' ') goto 990
            X1 = index(tempC,"/")
            X2 = index(tempC(X1+1:len_trim(tempC)),"/")+X1
            read(tempC(1:X1-1),'(i10)') mo
            read(tempC(X1+1:X2-1),'(i10)') dd
            read(tempC(X2+1:X2+4),'(i10)') yy
            read(tempC(X2+5:len_trim(tempC)),'(a)') templine
            read(templine,*) hh, hs, tm, tp, dir
            if (hs.le.-1) hs = -999
            if (dir.le.-1) dir = -999
            if (tm.le.-1) tm = -999
            if (tp.le.-1) tp = -999
            if (hh.ge.100) hh = hh/100
            write(21,220) yy, mo, dd, hh, mi, ss, hs
            write(22,221) yy, mo, dd, hh, mi, ss, dir
            write(23,220) yy, mo, dd, hh, mi, ss, tm
            write(24,220) yy, mo, dd, hh, mi, ss, tp
          enddo        
999     close(11)
      enddo
220   Format(6i7,f14.4)
221   Format(6i7,i14)
      close(21)
      close(22)
      close(23)
      close(24)

   end subroutine CSI2IMEDS   
!----------------------------------------------------------------------
!     30
!----------------------------------------------------------------------
   subroutine NDBC2IMEDS   


      implicit none

      character(120),allocatable       :: filename(:),stnid(:)
      character(120)                   :: filestation
      integer                          :: i, j, k
      integer                          :: nsta,yy,mo,dd,hh,mi,ss=0
      double precision                 :: dir,hs,tp,tm, junkR
      double precision, allocatable    :: lat(:),lon(:)
      character(120)                   :: junkC

! Open Station list file and output files
      write(6,*) "Enter the name of the NDBC station list"
      read(5,*) filestation
      !filestation = 'NDBC_Stations.txt
      open(unit=10,file=trim(filestation),status='old',action='read')
        read(10,*) nsta
        allocate(stnid(1:nsta),lat(1:nsta),lon(1:nsta),filename(1:nsta))
        do i=1,nsta
          read(10,*) stnid(i), lon(i), lat(i), filename(i)
        enddo
      close(10)

! Open Each file for Output
      open(unit=21,file="dir_NDBC.IMEDS")
      open(unit=22,file="hs_NDBC.IMEDS")
      open(unit=23,file="tp_NDBC.IMEDS")
      open(unit=24,file="tm_NDBC.IMEDS")      

! Create Header For Each Output
      write(21,*) '% IMEDS generic format ver. 1.0 wind-direction'
      write(21,*) '% year month day hour min sec dir (degrees)'
      write(21,*) 'NDBC UTC NAVD'
      write(22,*) '% IMEDS generic format ver. 1.0 wave-height'
      write(22,*) '% year month day hour min sec hs (m)'
      write(22,*) 'NDBC UTC NAVD'
      write(23,*) '% IMEDS generic format ver. 1.0 peak wave-period'
      write(23,*) '% year month day hour min sec tp (s)'
      write(23,*) 'NDBC UTC NAVD'
      write(24,*) '% IMEDS generic format ver. 1.0 mean wave-period'
      write(24,*) '% year month day hour min sec tp (s)'
      write(24,*) 'NDBC UTC NAVD'     

! Process Stations Files
      do i=1,nsta

        write(21,*) "NDBC_"//trim(stnid(i)), lat(i), lon(i)
        write(22,*) "NDBC_"//trim(stnid(i)), lat(i), lon(i)
        write(23,*) "NDBC_"//trim(stnid(i)), lat(i), lon(i)
        write(24,*) "NDBC_"//trim(stnid(i)), lat(i), lon(i)

        open(11,file=trim(filename(i)),status='old',action='read')
          write(6,*) "Processing Station: ", stnid(i)
          do 
            read(11,*,END=999) junkC, yy, mo, dd, hh, mi, (junkR,k=1,14), hs, junkR, tp, tm, dir
            if (hs.eq.-999.0d0) hs = -999.0d0
            if (tp.eq.-999.0d0) tp = -999.0d0
            if (tm.eq.-999.0d0) tm = -999.0d0
            if (dir.eq.-999.0d0) then
               dir = -999.0d0
            else
               dir = 270.0d0 - dir
               if (dir.lt.0.0d0) dir = dir + 360.0d0
            endif   
            write(21,200) yy, mo, dd, hh, mi, ss, dir
            write(22,200) yy, mo, dd, hh, mi, ss, hs
            write(23,200) yy, mo, dd, hh, mi, ss, tp
            write(24,200) yy, mo, dd, hh, mi, ss, tm            
          enddo
999      close(11)
      enddo
200   Format(6i7,F14.4)
      close(21)
      close(22)
      close(23)
      close(24)

   end subroutine NDBC2IMEDS   
!----------------------------------------------------------------------
!     31
!----------------------------------------------------------------------
   subroutine IMEDScntrl(os)

      implicit none
      character(120)                :: filecntrl, RunID
      character(10), intent(inout)  :: os
      integer                       :: extension=5, cntrl=1, X1
      integer                       :: i, NumFiles, dtype
      integer                       :: dateu(5),j, NumTypes
      character(120)                :: filestationu,fort14u
      character(120)                :: fort63u,fort61u     
      character(120)                :: filestation
      character(120)                :: dfolder 

      write(6,*) "Enter the name of the input cntrl file"
      call choosefile(filecntrl,extension,os)
      write(6,*)
      open(unit=10,file=filecntrl,status='old',action='read')
        read(10,*) dateu(1), dateu(2), dateu(3), dateu(4), dateu(5)
        read(10,*) fort14u
        read(10,*) RunID
        read(10,*) dfolder
        read(10,*) NumTypes
        do j = 1,NumTypes
          read(10,*) NumFiles
          read(10,*) fort63u
          read(10,*) dtype
          do i = 1,NumFiles
            read(10,'(A)') filestation
            filestationu = trim(filestation)
             X1 = len_trim(filestationu)
             if (filestationu(X1-7:X1).eq.'Stations') then
               fort61u = filestationu(1:X1-9)//'.'//trim(RunID)//'.61'
             else
               fort61u = trim(filestationu)//'.'//trim(RunID)//'.61'
            endif
            call TSextract(os,cntrl,filestationu,fort14u,fort63u,RunID,dfolder)
            call IMEDS61(os, cntrl, dateu, filestationu, fort61u,dtype,dfolder)
          enddo
        enddo
      close(10)

   end subroutine IMEDScntrl
!----------------------------------------------------------------------
!     41
!----------------------------------------------------------------------
   subroutine ExportOOB(os)

      implicit none
      integer                       :: nnode, nelem, nope, neta
      double precision, allocatable :: xy(:,:)
      character(10), intent(inout)  :: os      
      character(120)                :: fort14, OBfile, ptfile
      integer                       :: iyyyy,imo,idd,ihh,imi,iss,iddd,fddd
      integer                       :: fyyyy,fmo,fdd,fhh,fmi,fss,ts, dt=30
      double precision, allocatable :: dlat(:), dlon(:), elon(:)
      integer                       :: i, j, n, nvdll, nbdv, icnt, k
!
      write(6,*) 
      write(6,*) '----------------------------------------------------------'
      write(6,*) 'This Program Exports Open Boundary Nodal Coordinates'
      write(6,*) 'to be used as input for the OSU Tidal Model Driver Program'
      write(6,*) '----------------------------------------------------------'
      write(6,*)       
      write(6,*) 'Enter name of grid file' 
      read(5,*) fort14
      write(6,*) 'Enter name of the output points file' 
      read(5,*) ptfile      
      write(6,*)
      write(6,*) 'Enter the Initial DateTime (Year Month Day Hour Min Second)'
      write(6,*) 'YYYY MM DD HH MM SS'
      read (5,*) iyyyy, imo, iddd, ihh, imi, iss
      write(6,*)      
      write(6,*) 'Enter the Final DateTime (Year Month Day Hour Min Second)'
      write(6,*) 'YYYY MM DD HH MM SS'
      read (5,*) fyyyy, fmo, fddd, fhh, fmi, fss
      write(6,*)
!
      fdd = fddd
      idd = iddd
      select case (fmo)
        case(02)
          fdd = fdd + 31
        case(03)
          fdd = fdd + 31 + 28                
        case(04)
          fdd = fdd + 31 + 28 + 31             
        case(05)
          fdd = fdd + 31 + 28 + 31 + 30                
        case(06)
          fdd = fdd + 31 + 28 + 31 + 30 + 31
        case(07)
          fdd = fdd + 31 + 28 + 31 + 30 + 31 + 30
        case(08)
          fdd = fdd + 31 + 28 + 31 + 30 + 31 + 30 + 31
        case(09)
          fdd = fdd + 31 + 28 + 31 + 30 + 31 + 30 + 31 + 31
        case(10)
          fdd = fdd + 31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30
        case(11)
          fdd = fdd + 31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31
        case(12)
          fdd = fdd + 31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30
      end select
      select case (imo)
        case(02)
          idd = idd + 31
        case(03)
          idd = idd + 31 + 28                
        case(04)
          idd = idd + 31 + 28 + 31             
        case(05)
          idd = idd + 31 + 28 + 31 + 30                
        case(06)
          idd = idd + 31 + 28 + 31 + 30 + 31
        case(07)
          idd = idd + 31 + 28 + 31 + 30 + 31 + 30
        case(08)
          idd = idd + 31 + 28 + 31 + 30 + 31 + 30 + 31
        case(09)
          idd = idd + 31 + 28 + 31 + 30 + 31 + 30 + 31 + 31
        case(10)
          idd = idd + 31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30
        case(11)
          idd = idd + 31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31
        case(12)
          idd = idd + 31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30
      end select      
      ts = (fyyyy-iyyyy)*(365*24*3600+6*3600)
      ts = ts + (fdd-idd)*24*3600
      ts = ts + (fhh-ihh)*3600
      ts = ts + (fmi-imi)*60
      ts = ts + (fss-iss)
!      
      open(14,file=fort14,status='old',action='read')
      open(60,file='tides.OSU',status='replace',action='write')
      open(50,file=ptfile,status='replace',action='write')
      allocate(dlat(1:1000),dlon(1:1000),elon(1:1000))
!
      k = 0
      read(14,*)
      read(14,*) nelem, nnode
      allocate( xy(2,nnode) )
      do i = 1, nnode
        read(14,*) n, xy(1,n), xy(2,n)
      enddo
      do i = 1, nelem
        read(14,*)
      enddo
      read(14,*) nope
      read(14,*) neta
      do i = 1, nope
        read(14,*) nvdll
        do j = 1, nvdll
          read(14,*) nbdv
          k = k + 1
          dlat(k)=xy(2,nbdv)
          dlon(k)=360d0+xy(1,nbdv)
          elon(k)=xy(1,nbdv)
        enddo
      enddo
      icnt = k
      write(60,'(a)')'   lat   lon            yy   mm  dd hh mi sec dt(min) TSLength'
      write(50,*) icnt
      do k = 1, icnt
        write(60,'(2f10.5,i7,2x,6(i2.2,x),2x,i15)') dlat(k), dlon(k), iyyyy, imo, iddd, ihh, imi, iss, dt, ts
        write(50,'(2f10.5)') elon(k), dlat(k)
      enddo

      close(14)
      close(50)
      close(60)

      end subroutine ExportOOB

!----------------------------------------------------------------------
!     42
!----------------------------------------------------------------------
   subroutine OSUextract

      implicit none
      character(120)                :: OSUfile, OSUout, junkC
      integer                       :: nc=13, i, j, ns
      integer, allocatable          :: cindex(:)
      double precision              :: junkD
      double precision, allocatable :: amp(:,:), ph(:,:)
      character(120), allocatable   :: cname(:)

      write(6,*) "Enter the name of the OSU file"
      read(5,*) OSUfile
      write(6,*) "Enter the number of stations"
      read(5,*) ns      
      write(6,*) "Enter the name of the extract file"
      read(5,*) OSUout

      allocate(amp(1:nc,1:ns), ph(1:nc,1:ns),cindex(1:nc),cname(1:nc))

! Read OSU file
      open(unit=20,file=trim(OSUfile),status='old',action='read')
      read(20,*)
      read(20,*)
      do j = 1,ns
        do i=1,nc
          read(20,*) junkD, junkD, junkC, junkC, amp(i,j), ph(i,j)
        enddo
      enddo
      close(20)
! m2, s2, n2, k2, k1, o1, p1, q1, mf, mm, m4, ms4, mn4
! 1 ,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11,  12,  13
      cindex = (/ 8,6,7,5,3,1,2,4,9,10,11,12,13 /)
      cname = (/'Q1 ','O1 ','P1 ','K1 ','N2 ','M2 ','S2 ','K2 ',&
      & 'Mf ','Mm ','M4 ','MS4','MN4'/)
! Write Extract File

      open(unit=30,file=trim(OSUout),status='replace',action='write')
      do i = 1,nc
        write(30,*) trim(cname(i))
        do j = 1,ns
          write(30,200) amp(cindex(i),j), ph(cindex(i),j), i, j
        enddo
      enddo
      close(30)
200   Format(F14.4, F14.2,2i10) 

   end subroutine OSUextract

!----------------------------------------------------------------------
!     43
!----------------------------------------------------------------------
   subroutine FES2004extract

      implicit none
      character(120)                :: FESfile, FESout, junkC
      integer                       :: nc=14, i, j, ns
      integer, allocatable          :: cindex(:)
      double precision              :: junkD
      double precision, allocatable :: amp(:,:), ph(:,:)
      character(120), allocatable   :: cname(:)

      write(6,*) "Enter the name of the FES2004 file"
      read(5,*) FESfile
      write(6,*) "Enter the number of stations"
      read(5,*) ns      
      write(6,*) "Enter the name of the extract file"
      read(5,*) FESout

      allocate(amp(1:nc,1:ns), ph(1:nc,1:ns),cindex(1:nc),cname(1:nc))

! Read OSU file
      open(unit=20,file=trim(FESfile),status='old',action='read')
      do i = 1,7
        read(20,*)
      enddo
      do j = 1,ns
        read(20,*)      
        read(20,*)              
        do i=1,nc
          read(20,*) junkC, amp(i,j), ph(i,j)
        enddo
      enddo
      close(20)
! m2, s2, k1, k2, n2, 2n2, o1, p1, q1, m4, mf, mm, mtm, msqm
! 1 ,  2,  3,  4,  5,   6,  7,  8,  9, 10, 11, 12,  13,   14
      cindex = (/ 9, 7, 8, 3, 5, 1, 2, 4, 6, 10, 11, 12, 13, 14 /)
      cname = (/'Q1  ','O1  ','P1  ','K1  ','N2  ','M2  ','S2  ','K2  ',&
      & '2N2 ','M4  ','Mf  ','Mm  ','Mtm ','Msqm'/)
! Write Extract File

      open(unit=30,file=trim(FESout),status='replace',action='write')
      do i = 1,nc
        write(30,*) trim(cname(i))
        do j = 1,ns
          write(30,200) amp(cindex(i),j), ph(cindex(i),j), i, j
        enddo
      enddo
      close(30)
200   Format(F14.6, F14.6,2i10) 


   end subroutine FES2004extract

!----------------------------------------------------------------------
!     EC2001extract
!----------------------------------------------------------------------
   subroutine EC2001extract

!      implicit none

      PARAMETER(MNP=400000)
      PARAMETER(MNHARF=10)
      DIMENSION X(MNP),Y(MNP),XOUT(MNP),YOUT(MNP)
      DIMENSION NM1(2*MNP),NM2(2*MNP),NM3(2*MNP)
      DIMENSION FREQ(MNHARF),NFACT(MNHARF),EQARG(MNHARF)
      DIMENSION NAME(MNHARF)
      DIMENSION EAMP(3,MNHARF),EPHA(3,MNHARF),ETAMP(MNHARF,MNP)
      DIMENSION ETPHA(MNHARF,MNP)
      DIMENSION UAMP(3,MNHARF),UPHA(3,MNHARF),UTAMP(MNHARF,MNP)
      DIMENSION UTPHA(MNHARF,MNP)
      DIMENSION VAMP(3,MNHARF),VPHA(3,MNHARF),VTAMP(MNHARF,MNP)
      DIMENSION VTPHA(MNHARF,MNP)
      DIMENSION NN(3),NO(3)
      CHARACTER*10 NAME
      integer                       :: cindex(1:10),cindex2(1:10)
      REAL NFACT
      integer                       :: i,j,n,nvdll,nbdv,cnt, R
      double precision, allocatable :: xyd(:,:), lat(:), lon(:)
      integer                       :: nnode, nelem, nope, neta
      character(120)                :: fort14     


      write(6,*) 
      write(6,*) '-----------------------------------------------------'
      write(6,*) 'This Program Extracts Harmonic Constiutents from the '
      write(6,*) 'EC2001 Tidal Database'
      write(6,*) '-----------------------------------------------------'
      write(6,*)        
      write(6,*) 'Enter the name of the grid file to interpolate'
      read(5,*) fort14
      write(6,*)

      open(16,file=fort14,status='old',action='read')
      open(17,file='tides.in',status='replace',action='write')
      read(16,*)
      read(16,*) nelem, nnode
      allocate( xyd(3,nnode) )
      do i = 1, nnode
        read(16,*) n, xyd(1,n), xyd(2,n), xyd(3,n)
      enddo
      do i = 1, nelem
        read(16,*)
      enddo
      read(16,*) nope
      read(16,*) neta
      cnt = 0
      do i = 1, nope
        read(16,*) nvdll
        do j = 1, nvdll
          read(16,*) nbdv
          write(17,502) xyd(1,nbdv),xyd(2,nbdv)
          cnt = cnt + 1
        enddo
      enddo
      close(17)
      close(16)

 501  format(i10)            
 502  format(2f18.10)

      allocate(lat(1:cnt),lon(1:cnt))

      open(17,file='tides.in',status='old',action='read')
      do i = 1,cnt
        read(17,*) lat(i),lon(i)
      enddo
      close(17)
      
      open(17,file='tides.in',status='replace',action='write')
      write(17,501) cnt
      do i = 1,cnt
        write(17,502) lat(i),lon(i)
      enddo
      close(17)
      write(6,*)
      write(6,*) 'ADCIRC Grid has been read'
      write(6,*) 'Number of Open Boundary Nodes is ', cnt
      write(6,*)

      OPEN(1,FILE='tides.out',status='replace',action='write')
      OPEN(12,FILE='tides.in',status='old',action='read')
      OPEN(14,FILE='ec2001.grd',status='old',action='read')

!...DEFINE PI
      PI=3.141592653589793

!...OUTPUT FILE HEADER
      WRITE(1,3900)
 3900 FORMAT(//,8X,'Constituent',13x,'Elevation',16x,'East Velocity', &
      &          12X,'North Velocity')
      WRITE(1,3901)
 3901 FORMAT(11x,'Name/',8x,3(6x,'Amplitude     Phase'))
      WRITE(1,3902)
 3902 FORMAT(5X,'Lon',10X,'Lat',12X,'(m)',8X,'(deg)',2(9x,'(m/s)',6X, &
      &                                                  '(deg)'),/)

!...READ INPUT LOCATIONS FOR HARMONIC ANALYSIS OUTPUT
      READ(12,*) NOUT
      DO I=1,NOUT
        READ(12,*) XOUT(I),YOUT(I)
 1000   FORMAT(//,' OUTPUT WILL BE GENERATED FOR THE POSITION: ',/, &
      &            5X,F11.6,' E LONGITUDE , ',F11.6,' N LATITUDE',//)
      END DO

!...INPUT GRID INFORMATION FROM UNIT 14
      READ(14,'(A24)') AGRID
      READ(14,*) NE,NP

!...NODAL COORDINATES
      DO I=1,NP
        READ(14,*) JKI,X(JKI),Y(JKI)
        END DO

!....CONNECTIVITY TABLE 
      DO I=1,NE
        READ(14,*) JKI,NHY,NM1(JKI),NM2(JKI),NM3(JKI)
      END DO

!....COMPUTE ELEMENT IN WHICH EACH LOCATION LIES.
      write(6,*) 'Interpolating Nodes ...'
      DO I=1,NOUT
        write(6,*) I
        AEMIN=1.0E+25
        KMIN=0
        X4=XOUT(I)
        Y4=YOUT(I)
        DO K=1,NE
          X1=X(NM1(K))
          X2=X(NM2(K))
          X3=X(NM3(K))
          Y1=Y(NM1(K))
          Y2=Y(NM2(K))
          Y3=Y(NM3(K))
          AREAS=(X1-X3)*(Y2-Y3)-(X3-X2)*(Y3-Y1)
          A1=(X4-X3)*(Y2-Y3)-(X3-X2)*(Y3-Y4)
          A2=(X1-X3)*(Y4-Y3)-(X3-X4)*(Y3-Y1)
          A3=(X1-X4)*(Y2-Y4)-(X4-X2)*(Y4-Y1)
          AA=ABS(A1)+ABS(A2)+ABS(A3)
          AE=ABS(AA-AREAS)/AREAS
          IF(AE.LT.AEMIN) THEN
            AEMIN=AE
            NNE=K
            N1=NM1(K)
            N2=NM2(K)
            N3=NM3(K)
            AREA=AREAS
            ENDIF
          END DO

!......PRINT WARNING IF NODE OUTSIDE THE DOMAIN
        IF(AEMIN.GT.1.0E-5) THEN                  !OUTSIDE AN ELEMENT
          WRITE(*,2000) AEMIN
 2000     FORMAT(///,' WARNING -  SPECIFIED LOCATION DOES NOT LIE', &
      &           ' WITHIN ANY ELEMENT IN THE DOMAIN.',/,' CHECK THE', &
      &           ' LONGITUDE AND LATITUDE FOR THIS LOCATION', &
      &           ' PROGRAM WILL ESTIMATE NEAREST ELEMENT',/,' THE', &
      &           ' PROXIMITY INDEX FOR THIS LOCATION EQUALS ',E15.6)
        ENDIF

!......COMPUTE INFORMATION REQUIRED TO INTERPOLATE AT OUTPUT LOCATION
        X1=X(N1)
        X2=X(N2)
        X3=X(N3)
        Y1=Y(N1)
        Y2=Y(N2)
        Y3=Y(N3)
        STA2=((X4-X1)*(Y3-Y1)-(Y4-Y1)*(X3-X1))/AREA
        STA3=(-(X4-X1)*(Y2-Y1)+(Y4-Y1)*(X2-X1))/AREA
        STA1=1.0-STA2-STA3

!......DETERMINE ASCENDING ORDER OF NODES
        IF((N1.LT.N2).AND.(N2.LT.N3)) THEN
          NN(1)=N1
          NN(2)=N2
          NN(3)=N3
          NO(1)=1
          NO(2)=2
          NO(3)=3
          ENDIF
        IF((N2.LT.N1).AND.(N1.LT.N3)) THEN
          NN(1)=N2
          NN(2)=N1
          NN(3)=N3
          NO(2)=1
          NO(1)=2
          NO(3)=3
          ENDIF
        IF((N1.LT.N3).AND.(N3.LT.N2)) THEN
          NN(1)=N1
          NN(2)=N3
          NN(3)=N2
          NO(1)=1
          NO(3)=2
          NO(2)=3
          ENDIF
        IF((N2.LT.N3).AND.(N3.LT.N1)) THEN
          NN(1)=N2
          NN(2)=N3
          NN(3)=N1
          NO(3)=1
          NO(1)=2
          NO(2)=3
          ENDIF
        IF((N3.LT.N1).AND.(N1.LT.N2)) THEN
          NN(1)=N3
          NN(2)=N1
          NN(3)=N2
          NO(2)=1
          NO(3)=2
          NO(1)=3
          ENDIF
        IF((N3.LT.N2).AND.(N2.LT.N1)) THEN
          NN(1)=N3
          NN(2)=N2
          NN(3)=N1
          NO(3)=1
          NO(2)=2
          NO(1)=3
          ENDIF

!......WRITE THE ELEMENT CONTAINING THE OUTPUT LOCATION
 2100   FORMAT(' SPECIFIED LOCATION WAS FOUND IN ELEMENT ',I8,/, &
      &      5X,' WHICH IS MADE UP OF NODES ',3I8,//)

!......READ HARMONIC CONSTITUENT INFORMATION
        OPEN(105,FILE='ec2001.tdb')
        READ(105,*) NHC
        DO J=1,NHC
          READ(105,3100) FREQ(J),NFACT(J),EQARG(J),NAME(J)
 3100     FORMAT(E21.10,F11.7,F13.8,2X,A10)
          END DO
        READ(105,*) NPP
 3200   FORMAT(//)
        DO K=1,3
          IBEG=1
          IF(K.GE.2) IBEG=NN(K-1)+1
          DO II=IBEG,NN(K)-1
            READ(105,3200)
            END DO
          READ(105,*) NODE,(EAMP(NO(K),J),EPHA(NO(K),J),J=1,NHC)
          READ(105,*)      (UAMP(NO(K),J),UPHA(NO(K),J),J=1,NHC)
          READ(105,*)      (VAMP(NO(K),J),VPHA(NO(K),J),J=1,NHC)
          IF(NODE.NE.NN(K)) THEN
          WRITE(*,*) 'ERROR FINDING NODE ',NN(K),' PROGRAM HAS FOUND', &
      &                ' NODE ',NODE,' IN ITS PLACE'
            WRITE(*,*)' PROGRAM WILL NOW BE TERMINATED'
            STOP
            ENDIF
          DO J=1,NHC
            EPHA(NO(K),J)=PI*EPHA(NO(K),J)/180.
            UPHA(NO(K),J)=PI*UPHA(NO(K),J)/180.
            VPHA(NO(K),J)=PI*VPHA(NO(K),J)/180.
            END DO
         END DO
       CLOSE(105)

!......COMPUTE HARMONIC CONSTITUENTS AT THE OUTPUT LOCATION
        DO J=1,NHC
          E1R=EAMP(1,J)*COS(EPHA(1,J))
          E1I=EAMP(1,J)*SIN(EPHA(1,J))
          E2R=EAMP(2,J)*COS(EPHA(2,J))
          E2I=EAMP(2,J)*SIN(EPHA(2,J))
          E3R=EAMP(3,J)*COS(EPHA(3,J))
          E3I=EAMP(3,J)*SIN(EPHA(3,J))
          U1R=UAMP(1,J)*COS(UPHA(1,J))
          U1I=UAMP(1,J)*SIN(UPHA(1,J))
          U2R=UAMP(2,J)*COS(UPHA(2,J))
          U2I=UAMP(2,J)*SIN(UPHA(2,J))
          U3R=UAMP(3,J)*COS(UPHA(3,J))
          U3I=UAMP(3,J)*SIN(UPHA(3,J))
          V1R=VAMP(1,J)*COS(VPHA(1,J))
          V1I=VAMP(1,J)*SIN(VPHA(1,J))
          V2R=VAMP(2,J)*COS(VPHA(2,J))
          V2I=VAMP(2,J)*SIN(VPHA(2,J))
          V3R=VAMP(3,J)*COS(VPHA(3,J))
          V3I=VAMP(3,J)*SIN(VPHA(3,J))
          ETR=E1R*STA1+E2R*STA2+E3R*STA3
          ETI=E1I*STA1+E2I*STA2+E3I*STA3
          UTR=U1R*STA1+U2R*STA2+U3R*STA3
          UTI=U1I*STA1+U2I*STA2+U3I*STA3
          VTR=V1R*STA1+V2R*STA2+V3R*STA3
          VTI=V1I*STA1+V2I*STA2+V3I*STA3
          ETAMP(J,I)=SQRT(ETR*ETR+ETI*ETI)
          UTAMP(J,I)=SQRT(UTR*UTR+UTI*UTI)
          VTAMP(J,I)=SQRT(VTR*VTR+VTI*VTI)
          IF(ETAMP(J,I).EQ.0.) THEN
             ETPHA(J,I)=0.
             ELSE
             ETPHA(J,I)=180.*ACOS(ETR/ETAMP(J,I))/PI
             IF(ETI.LT.0.) ETPHA(J,I)=360.-ETPHA(J,I)
             ENDIF
          IF(UTAMP(J,I).EQ.0.) THEN
             UTPHA(J,I)=0.
             ELSE
             UTPHA(J,I)=180.*ACOS(UTR/UTAMP(J,I))/PI
             IF(UTI.LT.0.) UTPHA(J,I)=360.-UTPHA(J,I)
             ENDIF
          IF(VTAMP(J,I).EQ.0.) THEN
             VTPHA(J,I)=0.
             ELSE
             VTPHA(J,I)=180.*ACOS(VTR/VTAMP(J,I))/PI
             IF(VTI.LT.0.) VTPHA(J,I)=360.-VTPHA(J,I)
             ENDIF
          END DO

        END DO

!......WRITE THE RESULTS
      OPEN(20,FILE=trim(fort14)//'.EC2001',status='replace',action='write')
      cindex = (/ 4, 3, 2, 7, 5, 6, 8, 9, 10, 1/)
      cindex2 = (/ 1, 2, 4, 5, 6, 7, 8, 9, 10, 11 /)      
!      cname = (/ 'STEADY', 'K1', 'O1','Q1','M2','S2','N2','K2','M4','M6'/)
      DO R=1,NHC
        J = cindex(R)
        WRITE(1,3999) NAME(J)
        WRITE(20,'(A)') NAME(J)
 3999   FORMAT(//,1X,A10,/)
        DO I=1,NOUT
          WRITE(20,503) ETAMP(j,I),ETPHA(J,I),cindex2(R),I
          WRITE(1,4000) XOUT(I),YOUT(I),ETAMP(J,I),ETPHA(J,I), &
      &                  UTAMP(J,I),UTPHA(J,I),VTAMP(J,I),VTPHA(J,I)
 4000     FORMAT(1X,2(F11.6,2X),3(E12.5,2X,F8.3,3X))
        END DO
      END DO
      CLOSE(20)

      WRITE(*,*) ' '                   
      WRITE(*,4100)
 4100 FORMAT('**** RESULTS HAVE BEEN STORED IN FILE: TIDES.OUT ****')
 503  format(2e18.10, 2i10)
      CLOSE(14)
      CLOSE(1)

   end subroutine EC2001extract

!----------------------------------------------------------------------
!     45. Fort.53 Decomposition
!----------------------------------------------------------------------
   subroutine Decompose53

      implicit none
! Part I
      character(120)                :: agrid, fort14
      integer                       :: ne, np
      double precision, allocatable :: xyd(:,:)
      integer,          allocatable :: nm(:,:)
      integer                       :: nope, neta, nvdl_max
      integer,          allocatable :: nvdll(:), nbdv(:,:)
      integer                       :: nbou, nvel, nvel_max
      integer,          allocatable :: nvell(:), ibtype(:),  nbvv(:,:), ibconn(:,:)
      double precision, allocatable :: bar(:,:,:)
      integer                       :: nodemax
      integer,          allocatable :: nsequencer(:)
      integer                       :: i, j, k
! Part II
      character(120)                :: fort53
      integer                       :: ncon
      character(120),allocatable    :: conname(:)
      double precision              :: junkD
      double precision              :: amp, pha


!-------------------------------Part I---------------------------------------      
      write(6,*)       
      write(6,*) 'Enter name of Grid File' 
      read(5,*) fort14
      write(6,*)
      write(6,*) 'Enter name of Fort.53 File' 
      read(5,*) fort53
      write(6,*)      
!
!Open(14), for dynamic allocation
      write(6,*) 'START! READING of Global data:  ', fort14(1:len_trim(fort14))
      open(14,file=fort14,status='old',action='read')
         call read14_alloc ( 14, ne, np, nope, nbou, nvdl_max, nvel_max, nodemax )
      close(14)
! Dymanic Memory Allocation for Global Grid
         allocate( xyd(3,np) )
         allocate( nm(ne,3) )
         allocate( nvdll(nope)  )
         allocate( nbdv(nope,nvdl_max) )
         allocate( nvell(nbou), ibtype(nbou)  )
         allocate( nbvv(nbou,nvel_max), ibconn(nbou,nvel_max), bar(3,nbou,nvel_max) )
         allocate( nsequencer(nodemax) )
!Re-open & read global grid data(14)
      open(14,file=fort14,status='old',action='read')
         call read14 ( 14, ne, np, nope, nbou, nvdl_max, nvel_max, nodemax, nsequencer,      &
                       agrid, xyd, nm, neta, nvdll, nbdv, nvel, nvell, ibtype, nbvv, ibconn, &
                       bar )
      close(14)
      write(6,*) 'GRID SUCCESSFULLY READ!'

      open(53,file=fort53,status='old',action='read')
        read(53,*) ncon
        allocate(conname(1:ncon))
        do i=1,ncon
          read(53,*) junkD, junkD, junkD, conname(i)
        enddo
        read(53,*)
        do i=1,ncon
          open(530+i,file=trim(fort53)//'.'//trim(conname(i))//'.amp.63',status='replace',action='write')
          open(630+i,file=trim(fort53)//'.'//trim(conname(i))//'.pha.63',status='replace',action='write')          
          write(530+i,*) agrid
          write(530+i,*) ne, np
          write(630+i,*) agrid
          write(630+i,*) ne, np          
        enddo
        do j=1,np
          read(53,*)
          do i=1,ncon
            read(53,*) amp, pha
            write(530+i,*) j, xyd(1,j),xyd(2,j), amp
            write(630+i,*) j, xyd(1,j),xyd(2,j), pha
          enddo
        enddo

        do i = 1,ncon
          do j = 1,ne
            write(530+i,*) j, 3, ( nm(j,k), k = 1, 3 )
            write(630+i,*) j, 3, ( nm(j,k), k = 1, 3 )
          enddo
          do j = 1,10
            write(530+i,*) 0
            write(630+i,*) 0            
          enddo
          close(530+i)
          close(630+i)
        enddo
      close(53)

      end subroutine decompose53
!----------------------------------------------------------------------
!     Add Wave Refraction to Fort.13
!----------------------------------------------------------------------
   subroutine Addwave13

      IMPLICIT NONE

      INTRINSIC           :: INDEX
      INTRINSIC           :: TRIM

      CHARACTER(LEN=100)  :: CoarseFile
      CHARACTER(LEN=100)  :: Fort14File
      CHARACTER(LEN=1)    :: JunkC
      CHARACTER(LEN=1000) :: Line
      CHARACTER(LEN=100)  :: New13File
      CHARACTER(LEN=100)  :: Orig13File
      CHARACTER(LEN=1)    :: UserSel

      INTEGER,ALLOCATABLE :: CoarseConn(:,:)
      INTEGER             :: CoarseNumElems
      INTEGER             :: CoarseNumVerts
      INTEGER             :: Elem
      INTEGER             :: IA
      INTEGER             :: IE
      INTEGER             :: IV
      INTEGER             :: JunkI
      INTEGER             :: NumAttrs
      INTEGER             :: NumElems
      INTEGER             :: NumNonDefault
      INTEGER             :: NumVerts
      INTEGER,ALLOCATABLE :: Refrac(:)

      REAL(8),ALLOCATABLE :: CoarseLat(:)
      REAL(8),ALLOCATABLE :: CoarseLon(:)
      REAL(8)             :: DefaultValue
      REAL(8),ALLOCATABLE :: Lat(:)
      REAL(8),ALLOCATABLE :: Lon(:)
      REAL(8),ALLOCATABLE :: N(:)
      REAL(8)             :: SubArea1
      REAL(8)             :: SubArea2
      REAL(8)             :: SubArea3
      REAL(8)             :: TotalArea
      REAL(8)             :: UserLatMax
      REAL(8)             :: UserLatMin
      REAL(8)             :: UserLonMax
      REAL(8)             :: UserLonMin
      REAL(8)             :: Value
      REAL(8)             :: X1
      REAL(8)             :: X2
      REAL(8)             :: X3
      REAL(8)             :: Y1
      REAL(8)             :: Y2
      REAL(8)             :: Y3

      WRITE(*,'(A)',ADVANCE='YES') ""
      WRITE(*,'(A)',ADVANCE='YES') "This program develops the nodal attribute for wave refraction in SWAN,"
      WRITE(*,'(A)',ADVANCE='YES') "  and then adds it to an existing nodal attributes file."
      WRITE(*,'(A)',ADVANCE='NO' ) "Enter name of the mesh (fort.14): "
      READ(*,'(A)') Fort14File

      OPEN(UNIT=14,FILE=TRIM(Fort14File),ACTION='READ')
      READ(14,'(A)') JunkC
      READ(14,*) NumElems,NumVerts
      ALLOCATE(Lon(1:NumVerts))
      ALLOCATE(Lat(1:NumVerts))
      DO IV=1,NumVerts
         READ(14,*) JunkI,Lon(IV),Lat(IV)
      ENDDO
      CLOSE(UNIT=14,STATUS='KEEP')
      WRITE(*,'(A)',ADVANCE='YES') "The mesh was read successfully."

      WRITE(*,'(A)',ADVANCE='YES') ""
      WRITE(*,'(A)',ADVANCE='YES') "How do you want to apply the wave refraction?"
      WRITE(*,'(A)',ADVANCE='YES') "  1. Enable within a specified lon/lat box."
      WRITE(*,'(A)',ADVANCE='YES') "  2. Enable within the coverage defined by a coarse mesh."
      WRITE(*,'(A)',ADVANCE='NO')  "Enter selection (1/2): "
      READ(*,'(A)') UserSel

      IF(UserSel.EQ."1")THEN

         WRITE(*,'(A)',ADVANCE="YES") "You have selected option 1."
         WRITE(*,'(A)',ADVANCE="YES") ""
         WRITE(*,'(A)',ADVANCE="YES") "How do you want to specify the lon/lat box?"
         WRITE(*,'(A)',ADVANCE="YES") "  A. Enable wave refraction in the northern"
         WRITE(*,'(A)',ADVANCE="YES") "     Gulf of Mexico (lon < -82, lat > 25)."
         WRITE(*,'(A)',ADVANCE="YES") "  B. Specify the lon/lat box inside which"
         WRITE(*,'(A)',ADVANCE="YES") "     wave refraction will be enabled."
         WRITE(*,'(A)',ADVANCE="NO" ) "Enter selection (A/B): "
         READ(*,'(A)') UserSel

         IF(UserSel=="A")THEN
            WRITE(*,'(A)',ADVANCE="YES") "You have selected option A."
            UserLonMax =   -82.0
            UserLonMin = -9999.0
            UserLatMax =  9999.0
            UserLatMin =    25.0
         ELSEIF(UserSel=="B")THEN
            WRITE(*,'(A)',ADVANCE="YES") "You have selected option B."
            WRITE(*,'(A)',ADVANCE="YES") ""
            WRITE(*,'(A)',ADVANCE="NO")  "Enter maximum longitude:"
            READ(*,*) UserLonMax
            WRITE(*,'(A)',ADVANCE="NO")  "Enter minimum longitude:"
            READ(*,*) UserLonMin
            WRITE(*,'(A)',ADVANCE="NO")  "Enter maximum latitude:"
            READ(*,*) UserLatMax
            WRITE(*,'(A)',ADVANCE="NO")  "Enter minimum latitude:"
            READ(*,*) UserLatMin
         ELSE
            WRITE(*,'(A)',ADVANCE="YES") "FATAL ERROR: That is not an acceptable selection."
            STOP
         ENDIF

         ALLOCATE(Refrac(1:NumVerts))
         Refrac = 0
         DO IV=1,NumVerts
            IF((UserLonMin.LE.Lon(IV)).AND.(Lon(IV).LE.UserLonMax).AND. &
               (UserLatMin.LE.Lat(IV)).AND.(Lat(IV).LE.UserLatMax))THEN
               Refrac(IV) = 1
            ENDIF
         ENDDO

      ELSEIF(UserSel.EQ."2")THEN

         WRITE(*,'(A)',ADVANCE='YES') ""
         WRITE(*,'(A)',ADVANCE='YES') "You want to enable the wave refraction within the coverage defined by a coarse mesh."
         WRITE(*,'(A)',ADVANCE='NO' ) "Enter name of the coarse mesh (GOMEX.grd): "
         READ(*,'(A)') CoarseFile

         OPEN(UNIT=14,FILE=TRIM(CoarseFile),ACTION='READ')
         READ(14,'(A)') JunkC
         READ(14,*) CoarseNumElems,CoarseNumVerts
         ALLOCATE(CoarseLon(1:CoarseNumVerts))
         ALLOCATE(CoarseLat(1:CoarseNumVerts))
         DO IV=1,CoarseNumVerts
            READ(14,*) JunkI,CoarseLon(IV),CoarseLat(IV)
         ENDDO
         ALLOCATE(CoarseConn(1:CoarseNumElems,1:3))
         DO IE=1,CoarseNumElems
            READ(14,*) JunkI,JunkI,CoarseConn(IE,1),CoarseConn(IE,2),CoarseConn(IE,3)
         ENDDO
         CLOSE(UNIT=14,STATUS='KEEP')
         WRITE(*,'(A)',ADVANCE='YES') "The coarse mesh was read successfully."

         ALLOCATE(Refrac(1:NumVerts))
         Refrac = 0
         DO IV=1,NumVerts
            Elem = 0
            Inner1: DO IE=1,CoarseNumElems
               X1 = Lon(IV)
               X2 = CoarseLon(CoarseConn(IE,2))
               X3 = CoarseLon(CoarseConn(IE,3))
               Y1 = Lat(IV)
               Y2 = CoarseLat(CoarseConn(IE,2))
               Y3 = CoarseLat(CoarseConn(IE,3))
               SubArea1 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))
               X1 = CoarseLon(CoarseConn(IE,1))
               X2 = Lon(IV)
               X3 = CoarseLon(CoarseConn(IE,3))
               Y1 = CoarseLat(CoarseConn(IE,1))
               Y2 = Lat(IV)
               Y3 = CoarseLat(CoarseConn(IE,3))
               SubArea2 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))
               X1 = CoarseLon(CoarseConn(IE,1))
               X2 = CoarseLon(CoarseConn(IE,2))
               X3 = Lon(IV)
               Y1 = CoarseLat(CoarseConn(IE,1))
               Y2 = CoarseLat(CoarseConn(IE,2))
               Y3 = Lat(IV)
               SubArea3 = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))
               X1 = CoarseLon(CoarseConn(IE,1))
               X2 = CoarseLon(CoarseConn(IE,2))
               X3 = CoarseLon(CoarseConn(IE,3))
               Y1 = CoarseLat(CoarseConn(IE,1))
               Y2 = CoarseLat(CoarseConn(IE,2))
               Y3 = CoarseLat(CoarseConn(IE,3))
               TotalArea = ABS((X2*Y3-X3*Y2)-(X1*Y3-X3*Y1)+(X1*Y2-X2*Y1))
               IF((SubArea1+SubArea2+SubArea3).LT.(1.001*TotalArea))THEN
                  Elem = IE
                  EXIT Inner1
               ENDIF
            ENDDO Inner1
            IF(Elem.NE.0)THEN
               Refrac(IV) = 1
            ENDIF
         ENDDO
         DEALLOCATE(CoarseLon)
         DEALLOCATE(CoarseLat)
         DEALLOCATE(CoarseConn)

      ELSE

         WRITE(*,'(A)',ADVANCE="YES") "FATAL ERROR: That is not an acceptable selection."
         STOP

      ENDIF

      WRITE(*,'(A)',ADVANCE='YES') "The nodal attribute for wave refraction was developed successfully."
      DEALLOCATE(Lon)
      DEALLOCATE(Lat)

      WRITE(*,'(A)',ADVANCE='YES') ""
      WRITE(*,'(A)',ADVANCE='NO' ) "Enter name of the original nodal attributes file (orig.13): "
      READ(*,'(A)') Orig13File
      WRITE(*,'(A)',ADVANCE='NO' ) "Enter name of the new nodal attributes file (new.13): "
      READ(*,'(A)') New13File

      OPEN(UNIT=13,FILE=TRIM(Orig13File),ACTION='READ')
      OPEN(UNIT=23,FILE=TRIM(New13File),ACTION='WRITE')

!... Process the header information.
      READ(13,'(A)') Line
      WRITE(23,'(A)') TRIM(Line)
      READ(13,'(A)') Line
      WRITE(23,'(A)') TRIM(Line)
      READ(13,'(A)') Line
      READ(Line,*) NumAttrs
      WRITE(23,'(I12)') NumAttrs + 1

!... Process the default information for the nodal attributes.
      DO IA=1,NumAttrs
         READ(13,'(A)') Line
         WRITE(23,'(A)') TRIM(Line)
         READ(13,'(A)') Line
         WRITE(23,'(A)') TRIM(Line)
         READ(13,'(A)') Line
         WRITE(23,'(A)') TRIM(Line)
         READ(13,'(A)') Line
         WRITE(23,'(A)') TRIM(Line)
      ENDDO
      WRITE(23,'(A)') "wave_refraction_in_swan"
      WRITE(23,'(A)') "unitless"
      WRITE(23,*)     1
      WRITE(23,*)     1.D0

!... Process the non-default vertices.
      DO IA=1,NumAttrs
         READ(13,'(A)') Line
         WRITE(23,'(A)') TRIM(Line)
         READ(13,'(A)') Line
         WRITE(23,'(A)') TRIM(Line)
         READ(Line,*) NumNonDefault
         DO IV=1,NumNonDefault
            READ(13,'(A)') Line
            WRITE(23,'(A)') TRIM(Line)
         ENDDO
      ENDDO

!... Add the wave refraction as an additional nodal attribute.
      NumNonDefault = 0
      DO IV=1,NumVerts
         IF(Refrac(IV).EQ.0)THEN
            NumNonDefault = NumNonDefault + 1
         ENDIF
      ENDDO
      WRITE(23,'(A)') "wave_refraction_in_swan"
      WRITE(23,*) NumNonDefault
      DO IV=1,NumVerts
         IF(Refrac(IV).EQ.0)THEN
            WRITE(23,*) IV,0.D0
         ENDIF
      ENDDO

      CLOSE(UNIT=13,STATUS='KEEP')
      CLOSE(UNIT=23,STATUS='KEEP')
      WRITE(*,'(A)',ADVANCE='YES') "The new nodal attributes file was written successfully."
      WRITE(*,'(A)',ADVANCE='YES') ""

   end subroutine
!----------------------------------------------------------------------
!     River Flux Extraction
!----------------------------------------------------------------------
   subroutine RiverExtract

      implicit none
! Part I
      character(120)                :: agrid, fort14
      integer                       :: ne, np
      double precision, allocatable :: xyd(:,:)
      integer,          allocatable :: nm(:,:)
      integer                       :: nope, neta, nvdl_max
      integer,          allocatable :: nvdll(:), nbdv(:,:)
      integer                       :: nbou, nvel, nvel_max
      integer,          allocatable :: nvell(:), ibtype(:),  nbvv(:,:), ibconn(:,:)
      double precision, allocatable :: bar(:,:,:)
      integer                       :: nodemax
      integer,          allocatable :: nsequencer(:)
      integer                       :: i, j, k
! Part II
      integer                       :: iriv, iptmax
      integer, allocatable          :: inode(:,:), icnt(:)
      character(120), allocatable   :: irivname(:)      
      double precision              :: dglat, dglon, junk, ddlat, ddlon
      double precision, allocatable :: idist(:,:), iflux(:,:), iflow(:)
! Part III
      character(120)                :: infile
      character(120), allocatable   :: orivname(:)
      integer                       :: flow, rivopt,optmax, oriv, flowopt
      integer, allocatable          :: ocnt(:)
      double precision              :: cflow, adjflux, orat, irat
      double precision, allocatable :: olat(:,:),olon(:,:),oflux(:,:)
      double precision, allocatable :: oflow(:),odist(:,:)
! Part IV
      character(120)                :: outfile


!-------------------------------Part I---------------------------------------      
      write(6,*) 
      write(6,*) '----------------------------------------------------------'
      write(6,*) 'River Boundary Extraction Tool'
      write(6,*) '----------------------------------------------------------'
      write(6,*)       
      write(6,*) 'Enter name of Grid File' 
      read(5,*) fort14
      write(6,*)
!
!Open(14), for dynamic allocation
      write(6,*) 'START! READING of Global data:  ', fort14(1:len_trim(fort14))
      open(14,file=fort14,status='old',action='read')
         call read14_alloc ( 14, ne, np, nope, nbou, nvdl_max, nvel_max, nodemax )
      close(14)
! Dymanic Memory Allocation for Global Grid
         allocate( xyd(3,np) )
         allocate( nm(ne,3) )
         allocate( nvdll(nope)  )
         allocate( nbdv(nope,nvdl_max) )
         allocate( nvell(nbou), ibtype(nbou)  )
         allocate( nbvv(nbou,nvel_max), ibconn(nbou,nvel_max), bar(3,nbou,nvel_max) )
         allocate( nsequencer(nodemax) )
!Re-open & read global grid data(14)
      open(14,file=fort14,status='old',action='read')
         call read14 ( 14, ne, np, nope, nbou, nvdl_max, nvel_max, nodemax, nsequencer,      &
                       agrid, xyd, nm, neta, nvdll, nbdv, nvel, nvell, ibtype, nbvv, ibconn, &
                       bar )
      close(14)
      write(6,*) 'GRID SUCCESSFULLY READ!'

!-------------------------------Part II---------------------------------------      
!Identify Number of Rivers and Maximum Number of Nodes per River
      iriv = 0
      iptmax = 0
      do k = 1,nbou
        if (ibtype(k).eq.52) then
          iriv = iriv + 1
          if (iriv==1) then
            iptmax = nvell(k)
          else
            if (nvell(k).gt.iptmax) then
              iptmax = nvell(k)
            end if
          end if
        end if
      enddo
      allocate(idist(1:iriv,1:iptmax),inode(1:iriv,1:iptmax),icnt(1:iriv),irivname(1:iriv))
      allocate(iflux(1:iriv,1:iptmax),iflow(1:iriv))
!Store Nodes 
      iriv = 0
      do k = 1,nbou
        if (ibtype(k).eq.52) then
          iriv = iriv + 1
          iriv = int(iriv)
          icnt(iriv) = nvell(k)
          do j = 1,nvell(k)
            inode(iriv,j)=nbvv(k,j)
          enddo
        end if
      enddo
!Store Node Distances
      dglat = 110859.57d0  ! meters
      dglon = 96075.195d0  ! meters
      junk = 0.0d0
      do i = 1,iriv
        idist(i,1) = 0.0d0
        write(6,*) i, iriv
        do j = 2,icnt(i)
          ddlon = (xyd(1,inode(i,j))-xyd(1,inode(i,j-1)))*dglon
          ddlat = (xyd(2,inode(i,j))-xyd(2,inode(i,j-1)))*dglat
          idist(i,j) = (ddlon**2.0d0+ddlat**2.0d0)**0.5d0 + idist(i,j-1)
        enddo
      enddo
!-------------------------------Part III---------------------------------------
      write(6,*)
      write(6,*) 'Useful Information (Avg Flows)'
      write(6,*) 
      write(6,*) '(3-Days Prior to LandFall to end of Wind Data)'
      write(6,*) 'Gustav  -  8552 cms (302,000 cfs)'
      write(6,*) 'Ike     - 12006 cms (424,000 cfs)'
      write(6,*) 'Katrina -  5069 cms (179,000 cfs)'
      write(6,*) 'Rita    -  5635 cms (199,000 cfs)'
      write(6,*) 
      write(6,*) '(Average of 7-Days Prior to LandFall)'
      write(6,*) 'Gustav  -  8807 cms (311,000 cfs)'
      write(6,*) 'Ike     - 12318 cms (435,000 cfs)'
      write(6,*) 'Katrina -  4842 cms (171,000 cfs)'
      write(6,*) 'Rita    -  5437 cms (192,000 cfs)'
      write(6,*)  
      write(6,*) 'Select Flow Option (1 or 2):'
      write(6,*) ' 1) Mississippi Flow (Atchafalaya = 30% of Total)'
      write(6,*) ' 2) Manual Entry'
      read(5,*) flowopt
      select case(flowopt)
      case(1)
        irivname(1) = 'Mississippi River'
        irivname(2) = 'Atchafalaya River'
        write(6,*) 'Enter Flowrate in Mississippi River (cms)'
        read(5,*) iflow(1)
        iflow(2) = iflow(1)/0.7d0*0.3d0
        write(6,*)
        write(6,505) 'River'      , 'Flowrate (cms)'  , 'Flowrate (cfs)'
        write(6,506) 'Mississippi', iflow(1)          , iflow(1)/0.3048d0**3.0d0
        write(6,506) 'Atchafalya' , iflow(2)          , iflow(2)/0.3048d0**3.0d0
        write(6,506) 'Total'      , iflow(1)+iflow(2) , (iflow(1)+iflow(2))/0.3048d0**3.0d0
  505 format(3a20)      
  506 format(a20,2f20.2)
      case(2)
        do i = 1,iriv
          write(6,*) 'Enter Flowrate (cms) for River:', i
  507 format(a30,1x,i3)          
          read(5,*) iflow(i)
          write(6,*)
          write(6,*) 'Enter River Identification:'
          read(5,*) irivname(i)
        enddo
      end select
      write(6,*)
      write(6,*) 'Select River Forcing Option (1 or 2):'
      write(6,*) '1) Uniform'
      write(6,*) '2) Input from File'
      write(6,*) '3) SL15'
      read(5,*) rivopt
      write(6,*)
!
      select case(rivopt)
      case(1) ! Uniform Flux
        do i = 1,iriv
          do j = 1,icnt(i)
            iflux(i,j)=iflow(i)/idist(i,icnt(i))
          enddo
        enddo
      case(2) ! Input-Specified Flux
        write(6,*) 'Enter name of River Flow Input File:'
        read(5,*) infile
        write(6,*)      
        open(10,file=infile,status='old',action='read')      
        read(10,*) 
        read(10,*) 
        read(10,*) oriv
        read(10,*) optmax
        allocate(ocnt(1:oriv),olat(1:oriv,1:optmax),olon(1:oriv,1:optmax),oflux(1:oriv,1:optmax),orivname(1:oriv))
        allocate(oflow(1:oriv),odist(1:oriv,1:optmax))
        write(6,*) 'Reading Source File'
        write(6,*) 
        write(6,505) 'River'      , 'Flowrate (cms)'  , 'Florate (cfs)'
        do i=1,oriv
          read(10,'(a)') orivname(i)
          read(10,*) ocnt(i)
          do j=1,ocnt(i)
            read(10,*) olon(i,j),olat(i,j),oflux(i,j)
          enddo
          odist(i,1) = 0.0d0
          do j = 2,ocnt(i)
            ddlon = (olon(i,j) - olon(i,j-1)) * dglon
            ddlat = (olat(i,j) - olat(i,j-1)) * dglat
            odist(i,j) = (ddlon**2.0d0+ddlat**2.0d0)**0.5d0 + odist(i,j-1)
          enddo
          oflow(i) = oflux(i,1)*odist(i,2)*0.5d0
          oflow(i) = oflow(i) + 0.5d0*(odist(i,ocnt(i))-odist(i,ocnt(i)-1))*oflux(i,ocnt(i))
          do j = 2,ocnt(i)-1
            oflow(i) = oflow(i) + oflux(i,j)*(odist(i,j+1)-odist(i,j-1))*0.5d0
          enddo
          write(6,506) orivname(i), oflow(i), oflow(i)/0.3048d0**3.0d0          
        enddo
        close(10)        
      case(3)
        oriv = 2
        optmax = 2
        allocate(ocnt(1:oriv),olat(1:oriv,1:optmax),olon(1:oriv,1:optmax),oflux(1:oriv,1:optmax),orivname(1:oriv))
        allocate(oflow(1:oriv),odist(1:oriv,1:optmax))        
        orivname(1) = 'Mississippi River'
        orivname(2) = 'Atchafalaya River'
        ocnt(1) = 13
        ocnt(2) = 13        
        olon(1,1) = -91.1941234000d0
        olon(1,2) = -91.1949560000d0
        olon(1,3) = -91.1957467000d0
        olon(1,4) = -91.1965052000d0
        olon(1,5) = -91.1972709000d0
        olon(1,6) = -91.1980693000d0
        olon(1,7) = -91.1988772000d0
        olon(1,8) = -91.1996913000d0
        olon(1,9) = -91.2005160000d0
        olon(1,10) = -91.2013430000d0
        olon(1,11) = -91.2021879000d0
        olon(1,12) = -91.2030450000d0
        olon(1,13) = -91.2038990000d0
        olon(2,1) = -91.7997400000d0
        olon(2,2) = -91.8004550000d0
        olon(2,3) = -91.8011700000d0
        olon(2,4) = -91.8018850000d0
        olon(2,5) = -91.8027340000d0
        olon(2,6) = -91.8037000000d0
        olon(2,7) = -91.8048000000d0
        olon(2,8) = -91.8059000000d0
        olon(2,9) = -91.8070000000d0
        olon(2,10) = -91.8078700000d0
        olon(2,11) = -91.8087400000d0
        olon(2,12) = -91.8096100000d0
        olon(2,13) = -91.8104800000d0
        olat(1,1) = 30.42215590000d0
        olat(1,2) = 30.42260730000d0
        olat(1,3) = 30.42311210000d0
        olat(1,4) = 30.42365300000d0
        olat(1,5) = 30.42418610000d0
        olat(1,6) = 30.42468180000d0
        olat(1,7) = 30.42516590000d0
        olat(1,8) = 30.42564220000d0
        olat(1,9) = 30.42610470000d0 
        olat(1,10) = 30.42656400000d0
        olat(1,11) = 30.42699830000d0
        olat(1,12) = 30.42741430000d0
        olat(1,13) = 30.42783500000d0
        olat(2,1) = 30.98025000000d0
        olat(2,2) = 30.98023500000d0
        olat(2,3) = 30.98022000000d0
        olat(2,4) = 30.98020500000d0
        olat(2,5) = 30.98018000000d0       
        olat(2,6) = 30.98011700000d0
        olat(2,7) = 30.98004500000d0      
        olat(2,8) = 30.97997300000d0        
        olat(2,9) = 30.97990000000d0        
        olat(2,10) = 30.97987000000d0       
        olat(2,11) = 30.97984000000d0       
        olat(2,12) = 30.97981000000d0
        olat(2,13) = 30.97978000000d0 
        oflux(1,1) = 2.360d0
        oflux(1,2) = 4.250d0
        oflux(1,3) = 5.860d0
        oflux(1,4) = 6.610d0
        oflux(1,5) = 6.890d0
        oflux(1,6) = 6.800d0
        oflux(1,7) = 5.760d0
        oflux(1,8) = 4.820d0
        oflux(1,9) = 2.830d0
        oflux(1,10) = 1.320d0
        oflux(1,11) = 0.660d0
        oflux(1,12) = 0.660d0
        oflux(1,13) = 0.570d0
        oflux(2,1) = 2.180d0
        oflux(2,2) = 2.180d0
        oflux(2,3) = 2.180d0
        oflux(2,4) = 2.180d0
        oflux(2,5) = 2.040d0
        oflux(2,6) = 2.000d0
        oflux(2,7) = 1.900d0
        oflux(2,8) = 1.810d0
        oflux(2,9) = 1.730d0 
        oflux(2,10) = 1.640d0
        oflux(2,11) = 1.580d0
        oflux(2,12) = 1.530d0
        oflux(2,13) = 1.460d0
        do i=1,oriv
          odist(i,1) = 0.0d0
          do j = 2,ocnt(i)
            ddlon = (olon(i,j) - olon(i,j-1)) * dglon
            ddlat = (olat(i,j) - olat(i,j-1)) * dglat
            odist(i,j) = (ddlon**2.0d0+ddlat**2.0d0)**0.5d0 + odist(i,j-1)
          enddo
          oflow(i) = oflux(i,1)*odist(i,2)*0.5d0
          oflow(i) = oflow(i) + 0.5d0*(odist(i,ocnt(i))-odist(i,ocnt(i)-1))*oflux(i,ocnt(i))
          do j = 2,ocnt(i)-1
            oflow(i) = oflow(i) + oflux(i,j)*(odist(i,j+1)-odist(i,j-1))*0.5d0
          enddo
        enddo        
      end select
      if (rivopt.ge.2) then
        write(6,*) 'Linear Interpolation of Fluxes gives an adjustment factor of:'        
        do i = 1 , iriv
          iflux(i,1) = oflux(i,1)
          iflux(i,icnt(i)) = oflux(i,ocnt(i))
          irat = odist(i,ocnt(i))/idist(i,icnt(i))
          do j = 2,icnt(i)-1
            do k = 1,ocnt(i)
              if( odist(i,k) >= idist(i,j)*irat ) goto 600
            enddo
  600     orat = (odist(i,k)-odist(i,k-1))/(oflux(i,k)-oflux(i,k-1))
          iflux(i,j) = oflux(i,k-1)+(idist(i,j)*irat-odist(i,k-1))/orat
          enddo
          !Sum Flow and Adjust
          cflow = iflux(i,1)*idist(i,2)*0.5d0
          cflow = cflow + 0.5d0*(idist(i,icnt(i))-idist(i,icnt(i)-1))*iflux(i,icnt(i))
          do j = 2,icnt(i)-1
            cflow = cflow + iflux(i,j)*(idist(i,j+1)-idist(i,j-1))*0.5d0
          enddo
          adjflux = iflow(i)/cflow
          write(6,603) orivname(i),adjflux
          do j = 1, icnt(i)
            iflux(i,j) = iflux(i,j)*adjflux
          enddo
        enddo
  603 format(a30,f10.8)
      endif
!-------------------------------Part IV---------------------------------------
      write(6,*)
      write(6,*) 'Enter name of River Flow Output File'
      read(5,*) outfile
      write(6,*)
      open(11,file=outfile,status='replace',action='write')
      do i = 1,iriv
        do j = 1,icnt(i)          
          write(11,510) iflux(i,j), junk, '! ' // irivname(i), inode(i,j)
        enddo
      enddo
  510 format(f18.10,5x,f4.2,10x,a20,i10)      
      close(11)
!
   end subroutine
  
!----------------------------------------------------------------------
!     TOOLS-1
!----------------------------------------------------------------------
   subroutine getmonthday(dy,yy,mo,dd)

      implicit none
      integer, intent(out)           :: mo, dd
      integer, intent(in)            :: yy, dy
      integer                       :: i
      integer, allocatable          :: cd(:)

      allocate(cd(1:13))
      cd = (/ 0,31,59,90,120,151,181,212,243,273,304,334,365 /)
      if( mod(yy,4) == 0 ) then
        cd = (/ 0,31,60,91,121,152,182,213,244,274,305,335,366 /)
      endif
      do i = 1,12
        if (dy.gt.cd(i)) then
          mo = i
          dd = dy-cd(i)
        endif
      enddo

   end subroutine
!----------------------------------------------------------------------
!     TOOLS-2
!----------------------------------------------------------------------
   subroutine TZconvert(yy,mo,dd,hh,TZC)

      implicit none
      integer, intent(inout)      :: yy, mo, dd, hh, TZC
      integer                     :: i 
      integer, allocatable          :: calen(:)

      allocate(calen(1:12))
      calen = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
      if( mod(yy,4) == 0 ) then
        calen(2) = 29
      endif
      hh = hh + TZC
      do while (hh.gt.23)
        hh = hh-24
        dd = dd + 1
        if (dd.gt.calen(mo)) then
          dd = 1
          mo = mo + 1
          if (mo.gt.12) then
            mo = 1
            yy = yy + 1
          endif
        endif
      enddo
      do while (hh.lt.0)
        hh = hh+24
        dd = dd - 1
        if (dd.eq.0) then
          mo = mo - 1
          dd = calen(mo)
        endif
      enddo      

   end subroutine TZconvert
!----------------------------------------------------------------------
!     TOOLS-3
!----------------------------------------------------------------------
   subroutine choosefile(filename,extnum,os)

      implicit none
      character(10)                     :: extension
      character(10), intent(inout)      :: os
      character(120)                    :: filename
      character(120)                    :: systemcommand
      character(120), allocatable       :: filelist(:)
      integer                           :: i, n, iopt
      integer, intent(in)               :: extnum

      select case(extnum)
        case(1)
          extension='IMEDS'
        case(2) 
          extension='Stations'
        case(3)
          extension='61'
        case(4)
          extension='txt'
        case(5)
          extension='cntrl'
      end select
      if (trim(os) == "unix") then
        CALL system('ls -1 -p *.'//TRIM(extension)//'> choosefilelist')
      else
        CALL system('dir /b *.'//TRIM(extension)//'> choosefilelist')
      endif

      open(15,file="choosefilelist")
        n = 0
        do
          read(15,'(A)',END=10)
          n = n + 1
        enddo
10    close(15)
      allocate(filelist(1:n))
      open(15,file="choosefilelist")
        do i = 1,n
          read(15,'(A)') filelist(i)
          write(6,*) i, ") ", trim(filelist(i)) 
        enddo
      close(15)
      write(6,*) n+1,") ", "User-Defined"
999   write(6,'(A)',ADVANCE="NO") "Enter Selection: "
      read(5,'(i10)',ERR=999) iopt
      if (iopt<n+1) then
        filename = trim(filelist(iopt))
      else
        write(6,'(A)',ADVANCE="NO") "Enter User-Defined Filename:"
        read(5,*) filename
      endif
      write(6,*) filename
   end subroutine

   
    
!-----+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
! READ GRID DATA FOR MEMORY ALLOCATION
   subroutine read14_alloc ( iunit, ne, np, nope, nbou, nvdl_max, nvel_max, nodemax )
!-----+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
      implicit none
      integer, intent(in) :: iunit
      integer, intent(out) :: ne, np, nope, nbou, nvdl_max, nvel_max
      integer, intent(out) :: nodemax ! for Sequencer of node
      integer :: i, j, k
!
      nvdl_max = 0
      nvel_max = 0
      nodemax = 0
      read(iunit,*)
      read(iunit,*) ne, np
        do k = 1, np
          read(iunit,*) i
          nodemax = max(i,nodemax)
        enddo
        do k = 1, ne
          read(iunit,*)
        enddo
        write(6,*) '  |'
        read(iunit,*) nope
        read(iunit,*) 
        do k = 1, nope
          read(iunit,*) i
          if( i >= nvdl_max ) nvdl_max = i
          do j = 1, i
            read(iunit,*)
          enddo
        enddo
        read(iunit,*) nbou
        read(iunit,*)
        do k = 1, nbou
          read(iunit,*) i
          if( i >= nvel_max ) nvel_max = i
          do j = 1, i
            read(iunit,*)
          enddo
        enddo
!
   end subroutine read14_alloc

!-----+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
! READ GRID DATA
   subroutine read14 ( iunit, ne, np, nope, nbou, nvdl_max, nvel_max,  nodemax, nsequencer,  &
                       agrid, xyd, nm, neta, nvdll, nbdv, nvel, nvell, ibtype, nbvv, ibconn, &
                       bar )
!-----+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
      implicit none
      integer, intent(in) :: iunit, nodemax
      integer, intent(inout) :: ne, np, nope, nbou, nvdl_max, nvel_max
      double precision, intent(out) :: xyd(3,np)
      integer, intent(out) :: nm(ne,3)
      character(120), intent(out) :: agrid
      integer, intent(out) :: neta, nvdll(nope), nbdv(nope,nvdl_max)
      integer, intent(out) :: nvel, nvell(nbou), nbvv(nbou,nvel_max), ibtype(nbou)
      integer, intent(out) :: ibconn(nbou,nvel_max)
      double precision, intent(out) :: bar(3,nbou,nvel_max)
!
      integer, intent(out) :: nsequencer(nodemax)
!
     integer :: i, j, k, jn, je, nhy
!
      nsequencer(:) = 0
      bar(:,:,:) = 0.0d0
      ibconn(:,:) = 0
!
        read(iunit,'(a)') agrid
         read(iunit,*) ne, np
         do k = 1, np
            read(iunit,*) jn, (xyd(j,k), j=1,3)
            nsequencer(jn) = k
         enddo
         write(6,*) '  + '
         do k = 1, ne
            read(iunit,*) je, nhy, ( nm(k,j), j = 1, 3 )
            do j = 1, 3
               if( nm(k,j) <= 0 ) write(6,*) k,j, nm(k,j)
               nm(k,j) = nsequencer(nm(k,j))
            enddo
         enddo
         read(iunit,*) nope
         read(iunit,*) neta
         do k = 1, nope
            read(iunit,*) nvdll(k)
            do j = 1, nvdll(k)
               read(iunit,*) nbdv(k,j)
               nbdv(k,j) = nsequencer(nbdv(k,j))
            enddo
         enddo
         read(iunit,*) nbou
         read(iunit,*) nvel
         do k = 1, nbou
            read(iunit,*) nvell(k), ibtype(k)
            select case(ibtype(k))
               case(0,1,2,10,11,12,20,21,22,30,52)
                  do j = 1, nvell(k)
                     read(iunit,*) nbvv(k,j)
                     nbvv(k,j) = nsequencer(nbvv(k,j))
                  enddo
               case(3, 13, 23)
                  do j = 1, nvell(k)
                     read(iunit,*) nbvv(k,j), (bar(i,k,j), i=1,2)
                     nbvv(k,j) = nsequencer(nbvv(k,j))
                  enddo
               case(4, 24)
                  do j = 1, nvell(k)
                     read(iunit,*) nbvv(k,j), ibconn(k,j), (bar(i,k,j), i=1,3)
                     nbvv(k,j) = nsequencer(nbvv(k,j))
                     ibconn(k,j) = nsequencer(ibconn(k,j))
                  enddo
               case default
                  write(6,*) ' IBTYPE is not allowed', ibtype(k)
                  write(6,*)
                  write(6,*) '**** Hit the Enter-Key to stop ****'
                  read(5,*)
                  stop
            end select
         enddo
!
!      deallocate( nsequencer )
!
   end subroutine read14
