PROGRAM interp_ned
c  
c   Author: Steve Roecker, RPI  February, 2004
c
c   Spline inerpolation routines are from Numerical recipes.
c   Byte swapping routines are adapted from some hints I got from a few sites on the web
c   All other routines written by Steve Roecker
c
c   Read in a DEM from the USGS NED (generated by the seamless data distribution
c   Web site) and interpolate using bicubic splines to find values at requested locations
c
c   When you request data from the USGS, make sure to ask for GRIDFLOAT format.  That
c   data most likely will be in "Little Endian" form, and if it is this program will convert
c   it for you.
c
c   The USGS site will provide you with a number of files, two of which are used
c   here: the *.HDR file and the *.FLT file.  You will be asked for the basename of these
c   files and it will open the appropriate files for you.
c
c   Note:  this program may not compile with SUN's f77 compiler (it doesn't like the tricky
c   bits in the byte swapping routines) but WILL compile with g77.
c
c   g77 interp_ned.f -o interp_ned
c
c   should do it.
c
c   Some generic comments:
c
c	Input is a binary file with a sequence if real*4 floats.  Parameters
c	describing this file are read in through a separate ascii file.  The
c	variables that should be specified are:
c
c	nx	number of rows
c	ny	number of columns
c	yo	y origin (e.g., latitude of point 1)
c	xo	x origin (e.g., longitude of point 1)
c	dy	y increment between nodes
c	dx	x incremant between nodes
c
c	m is maximum number of columns
c	n is maximum number of rows

	integer m, n
	parameter (m=2000,n=2000,mn=m*n)
	integer  i,j
	real*4  f,xx1,xx2,x1(m),x2(n),y(mn),y2(mn)
	real*4  y2tmp(n),ytmp(n),yytmp(n)

	real*8 dx, dy, dx2, dy2, yo, xo
	real*8 xllcorner, yllcorner, cellsize
	real*8 xmin, xmax, ymin, ymax

	character*80  filename
	character*256 aline
	character*80 mapfile
	character*80 headerf
	character*132 byteord
	character*132 pval
	character*256 basename

	logical litend

	write(*,*) ' Enter base name for *.HDR, *.FLT files: '
	read(*,'(a)')  basename
	call dblank(basename, lenb)
	headerf = basename(1:lenb)//".HDR"
	mapfile = basename(1:lenb)//".FLT"

c----read in parms
	open(4,file=headerf)
	call getvars("ncols", pval, len, ierr)
	if (ierr.eq.0) read(pval(1:len),*) ncols
	call getvars("nrows", pval, len, ierr)
	if (ierr.eq.0) read(pval(1:len),*) nrows
	call getvars("xllcorner", pval, len, ierr)
	if (ierr.eq.0) read(pval(1:len),*) xllcorner
	call getvars("yllcorner", pval, len, ierr)
	if (ierr.eq.0) read(pval(1:len),*) yllcorner
	call getvars("NODATA_value", pval, len, ierr)
	if (ierr.eq.0) read(pval(1:len),*) nullval
	call getvars("byteorder", pval, leno, ierr)
	if (ierr.eq.0) byteord = pval(1:leno)
	leno = leno - 1
	call getvars("cellsize", pval, len, ierr)
	if (ierr.eq.0) read(pval(1:len),*) cellsize

	write(*,*) ' ncols = ', ncols
	write(*,*) ' nrows = ', nrows
	write(*,*) ' xllcorner = ', xllcorner
	write(*,*) ' yllcorner = ', yllcorner
	write(*,*) ' cellsize = ', cellsize
	write(*,*) ' null value = ', nullval
	write(*,'(2a)') ' byte order = ', byteord(1:leno)

c---decide if we need to swap bytes
	call endian(litend)
	iswap = 0
	if (litend) then
	  write(*,*) ' This machine is little-endian '
	  if (byteord.eq.'LSBFIRST') then
	    write(*,*) ' Data is litte-endian also '
	  else
	    write(*,*) ' Data is big-endian, swap required  '
	    iswap = 1
	  endif
	else
	  write(*,*) ' This machine is big-endian '
	  if (byteord(1:leno).eq."LSBFIRST") then
	    write(*,*) ' Data is litte-endian, swap required  '
	    iswap = 1
	  else
	    write(*,*) ' Data is big-endian also '
	  endif
	endif

c       xdim = 1.d0/(ncols-1)
c       ydim = 1.d0/(nrows-1)
	xdim = cellsize
	ydim = cellsize

c  lsize is the length in words (32 bits, so 2 16 bit integers per word)
	lsize = 4
	nbytes = nrows*ncols
	if (nbytes.gt.mn) then
	  write(*,*) ' Error: nbytes is ', nbytes
	  write(*,*) ' and it cannot be greater than ', mn
	  write(*,*) ' You must redimension this program to continue .. '
	  stop
	endif
	lrecl = lsize*nbytes
	ierr = 0

	xmin = xllcorner
	ymin = yllcorner
	xmax = xllcorner + xdim*(ncols-1)
	ymax = yllcorner + ydim*(ncols-1)

	nx = ncols
	ny = nrows
	xo = xllcorner
	yo = yllcorner + ydim*(nrows-1)
	dx = xdim
	dy = ydim

	if (nx.gt.m) then
	  write(*,*) ' Error; nx = ', nx
	  stop
	end if
	if (ny.gt.n) then
	  write(*,*) ' Error; ny = ', ny
	  stop
	end if

c----set up the nominal grid (x1, x2)
c  Note that we add onto yo here because of the requirements of
c	the spline subroutine (increasing x).  We will adjust for
c	this later in the program
	do 10 i = 1, ny
	  x2(i) = yo + dabs(dy)*(i-1)
10	continue
	do 11 i = 1, nx
	  x1(i) = xo + dabs(dx)*(i-1)
11	continue

c	write(9,125) (x2(i),i=1,ny)
c	write(9,125) (x1(i),i=1,nx)
c125	format(10f10.5)

	write(*,*) ' nx, ny = ', nx, ny
	write(*,*) ' xo, yo = ', xo, yo
	write(*,*) ' dx, dy = ', dx, dy
	write(*,*) yo, xo, dx2, dy2

c----now read in elevations
	nelev = nx*ny
	write(*,*) ' nelev = ', nelev

c---lsize = 4 for 4 byte real
	lsize = 4
	lrecl = lsize*nelev

	open (unit=7, file= mapfile, access='direct',
     +      recl=lrecl)
	read(7,rec=1) (y(i),i=1,nelev)
	close(7)

        if (iswap.eq.1) then
          do i = 1, nbytes
            call byteswapr4(y(i))
          enddo
        endif

	write(*,*) ' Generating splines in splie2 ...'

	call splie2(x1,x2,y,nx,ny,y2)

	write(*,*) ' Out of splie2'
	nelev = 0

	write(*,*) ' Enter a y (or latitude) and  x (or longitude) : '
	read(*,*) xlat, xlon
	if (xlat.gt.ymax.or.xlat.lt.ymin) then
	  write(*,*) ' WARNING: Your Latitude is outside the bounds of your model!'
	endif
	if (xlon.gt.xmax.or.xlon.lt.xmin) then
	  write(*,*) ' WARNING: Your Longitude is outside the bounds of your model!'
	endif
c---Here's where we make up for the requirement of increasing latitudes above
	xx2 = 2*yo - xlat
	xx1 = xlon
	do 112 j=1,nx
	  do 111 k=1,ny
	    nel = j + (k-1)*nx
	    ytmp(k) = y(nel)
	    y2tmp(k) = y2(nel)
111	  continue
	  call splint(x2,ytmp,y2tmp,ny,xx2,yytmp(j))
112	continue
     	call spline(x1,yytmp,nx,1.e30,1.e30,y2tmp)
	call splint(x1,yytmp,y2tmp,nx,xx1,f)

	write(*,*) ' Interpolated elevation is: ', f, ' meters or ', f/.3048,' feet'

c----As a check, output nearest neighbor elevations
	ilat = (yo-xlat)/ydim 
	ilon = (xlon-xo)/xdim + 1
	nwel = ncols*ilat + ilon
	neel = nwel + 1
	swel = nwel + ncols
	seel = swel + 1
	write(*,*) ' Northwest elevation = ', y(nwel), ' meters or ', y(nwel)/.3048,' feet'
	write(*,*) ' Northeast elevation = ', y(neel), ' meters or ', y(neel)/.3048,' feet'
	write(*,*) ' Southwest elevation = ', y(swel), ' meters or ', y(swel)/.3048,' feet'
	write(*,*) ' Southeast elevation = ', y(seel), ' meters or ', y(seel)/.3048,' feet'

	stop
	end

      SUBROUTINE splie2(x1a,x2a,ya,m,n,y2a)
      INTEGER m,n,NN
      REAL x1a(m),x2a(n),y2a(m,n),ya(m,n)
      PARAMETER (NN=2000)
CU    USES spline
      INTEGER j,k
      REAL y2tmp(NN),ytmp(NN)

c---dimension check
	if (n.gt.NN) then
	  write(*,*) ' Error in splie2: n, NN = ', n, NN
	  stop
	end if
c---test
c	write(*,*) ' IN SPLIE2: '
c	write(*,*) ' m,n = ', m, n
c	write(*,*)'  ya(1,1), ya(1,2) = ', ya(1,1), ya(1,2)
c	write(*,*)'  ya(2,1), ya(2,2) = ', ya(2,1), ya(2,2)
c	write(*,*)'  ya(3,1), ya(3,2) = ', ya(3,1), ya(3,2)
      do 13 j=1,m
	do 11 k=1,n
	  ytmp(k)=ya(j,k)
11      continue
	call spline(x2a,ytmp,n,1.e30,1.e30,y2tmp)
	do 12 k=1,n
	  y2a(j,k)=y2tmp(k)
12      continue
13    continue
      return
      END


      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      REAL yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=2000)
      INTEGER i,k
      REAL p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then
	y2(1)=0.
	u(1)=0.
      else
	y2(1)=-0.5
	u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
	sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
	p=sig*y2(i-1)+2.
	y2(i)=(sig-1.)/p
	u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
	qn=0.
	un=0.
      else
	qn=0.5
	un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
	y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END

      SUBROUTINE splin2(x1a,x2a,ya,y2a,m,n,x1,x2,y)
      INTEGER m,n,NN
      REAL x1,x2,y,x1a(m),x2a(n),y2a(m,n),ya(m,n)
      PARAMETER (NN=2000)
CU    USES spline,splint
      INTEGER j,k
      REAL y2tmp(NN),ytmp(NN),yytmp(NN)
      do 12 j=1,m
	do 11 k=1,n
	  ytmp(k)=ya(j,k)
	  y2tmp(k)=y2a(j,k)
11      continue
	call splint(x2a,ytmp,y2tmp,n,x2,yytmp(j))
12    continue
      call spline(x1a,yytmp,m,1.e30,1.e30,y2tmp)
      call splint(x1a,yytmp,y2tmp,m,x1,y)
      return
      END

      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      REAL x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
	k=(khi+klo)/2
	if(xa(k).gt.x)then
	  khi=k
	else
	  klo=k
	endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) pause 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.
      return
      END

        subroutine dblank(line, len)
        character*256 line

        do i = 1, 256
          j = 257 - i
          if (line(j:j).ne.' ') go to 1
        enddo
        len = 0
        return

1       do i = 1, j
          if (line(i:i).eq.' ') then
            do k = i, j+1
              line(k:k) = line(k+1:k+1)
            enddo
          endif
        enddo

        do i = 1, 256
          j = 257 - i
          if (line(j:j).ne.' ') go to 2
        enddo

2       len = j

        do i = 1, j
          if (line(i:i).eq.' ') go to 1
        enddo

        return
        end


	subroutine endian(litend)

c     checks if this is a little endian machine
c     returns litend=.true. if it is, litend=.false. if not

	integer*1 j(2)
	integer*2 i
	equivalence (i,j)
	logical litend

	i = 1
	if (j(1).eq.1) then
	   litend = .true.
	else
	   litend = .false.
	end if
	return

	end

	subroutine byteswapi2(k)

c     does a byteswap on integer2 number

	integer*1 ii(2), jj(2)
	integer*2 i, j, k
	equivalence (i,ii)
	equivalence (j,jj)

	i = k

	jj(1) = ii(2)
	jj(2) = ii(1)

	k = j

	return

	end


	subroutine byteswapr4(r)

c     does a byteswap on real*4 number

	integer*1 ii(4), jj(4)
	real*4 r, s, t
	equivalence (s,ii)
	equivalence (t,jj)

	s = r

	jj(1) = ii(4)
	jj(2) = ii(3)
	jj(3) = ii(2)
	jj(4) = ii(1)

	r = t

	return

	end


c----Subroutine to parse the parameter file to get values
c
c	Blank lines are ignored
c	Comment lines begin with a "#"
c	Each parameter specification line contains space/tab delimited fields that
c	 specify  the name of the parameter 
c	and the value specified.  Any text appearing after the last field is ignored 
c	and can be used to insert line-specific comments.
c
c	Syntax:
c
c	getvars(vname, parval, len, error)
c
c	vname		132 character max string with name of variable
c	parval		132 character max string with value of variable
c	len		length of parval
c	error		=0	if variable found; 1 variable not found
c

	subroutine getvars(vname, parval, len, error)

	character*132 vname, varname, parval
	integer error, len

	character*132 aline

	lunspc = 4

	ierr = 0
	error = 1

	rewind (lunspc)

	nline = 0
1	call getline(lunspc,aline,ierr)
	if (ierr.eq.1) go to 60
	nline = nline + 1
	if (ierr.ne.0) go to 1

c------recover the variable name
	ib = 1
	call getfield(lunspc,aline,ib,ie,varname,lenv,ierr)
	if (ierr.ne.0) go to 15
c-----see if this parameter name is valid
	if (varname(1:lenv).ne.vname(1:lenv)) go to 1

c-----recover the variable setting
	ib = ie
	call getfield(lunspc,aline,ib,ie,parval,nvl,ierr)
	if (ierr.ne.0) go to 15

c	write(*,*) ' j = ', j
c	write(*,*) ' parval:  ', parval(1:nvl)

	len = nvl
	error = 0
	return

15	write(*,*) ' getvars Warning:  getfield error at line ',nline
	go to 1

50	write(*,*) ' getvars Read error:  line = ', nline
	go to 1


60	error = 1
c	write(*,*) ' getvars warning:  EOF reached in pararmeter file'

	return
	end


c----Subroutine to retrieve a field from a line of text.
c	Fields are separated by tabs or spaces.
c
c	If the field is a "\" then the next line of text.
c	is read in from unit "lun", but blank lines and comment
c	lines beginning with a "#" are ignored.
c
c	If a field begins with double quotes (") then all 
c	text is read until another double quote is encountered.
c
c	Note that continuation of quoted strings is not allowed.
c
c
c	Input:
c	lun	the unit number
c	aline	the line to interpret
c	ib	start at the ib'th character of aline
c
c	Output:
c	field	the field
c	len	the length of outs
c	ierr	= 0
c	ie	= the position of aline at the end of the read
c	
	subroutine getfield(lun,aline,ib,ie,field,len,ierr)

	character*132 aline, field
	character*1 tab

	tab = char(9)

1	nch = 0
	do 10 i = ib, 132
c----ignore leading tabs or spaces
	  if (nch.eq.0.and.(aline(i:i).eq.' '.or.aline(i:i).eq.tab)) go to 10
	  if (aline(i:i).eq.' '.or.aline(i:i).eq.tab) go to 12
	  nch = nch + 1
c----test for a quoted string
	  if (nch.eq.1.and.aline(i:i).eq.'"') go to 14
c----test for a continuation
	  if (nch.eq.1.and.aline(i:i).eq.'\\') go to 22
	  field(nch:nch) = aline(i:i)
	  len = nch
10	continue
	write(*,*) ' Error:  EOL reached in getfield! '
	ierr = 1
	return

12	ierr = 0
	ie = i
	return

c----quoted string section
14	i1 = i+1
	nch = 0
	do 20 i = i1, 132
	  if (aline(i:i).eq.'"') go to 12
	  nch = nch + 1
	  field(nch:nch) = aline(i:i)
	  len = nch
20	continue
	write(*,*) ' Error:  End of quoted string not found in getfield! '
	ierr = 1
	return

c----continuation section
22	call getline(lun,aline,ierr)
	if (ierr.eq.1) go to 60
	if (ierr.ne.0) go to 50
	ib = 1
	go to 1

50	write(*,*) ' Read error in getfield! '
	ierr = 1
	return

60	write(*,*) ' Error:  EOF encountered in getfield! '
	ierr = 1

	return
	end


c----Subroutine to retrieve the next line of text that is
c	neither blank nor a comment (lines beginning with a "#")
c
c	Input:
c	lun	the unit number
c	aline	the line to interpret
c
c	Output:
c	aline	the line
c
c	ierr = 1 if EOF has been reached
c	
	subroutine getline(lun,aline,ierr)

	character*132 aline
	character*1 tab

	tab = char(9)

1	read(lun,'(a132)', err=50, end =60) aline
c------skip over comments
	if (aline(1:1).eq.'#') go to 1
	do 10 i = 1, 132
	  if (aline(i:i).ne.' '.and.aline(i:i).ne.tab) go to 2
10	continue
	go to 1

2	ierr = 0
	return

50	write(*,*) ' Read error in getline! '
	ierr = 2
	return

60	ierr = 1

	return
	end