c
c
c     This program extracts (topography/bathymetry) heights for a given region.
c     
c     Input: 
c
c     262.0 <= minlon < maxlon < 402.0 (equivalent to 98W - 42E)
c       0.0 <= minlat < maxlat <  75.0
c
c     Output:
c     One line per gridpoint:
c     longitude, latitude, height  (negative over water)
c
c
      program extract
      implicit real*8 (a-h, o-z)

      character*120 filename
      dimension igrid(0:3360),hd(3338,1562)

      x0 = 262.0d0
      y0 =   0.0d0
      c24 = 24.0d0
      c24div = 1.0d0/c24
      c12 = 12.0d0
      c12div = 1.0d0/c12
      delta  = 0.1d0*c24div
      c2div  = 0.50d0
      c4div  = 0.25d0

   50 write (6, 6001)
 6001 format ('Input minlon, maxlon, minlat, maxlat')
      write (6, 6002)
 6002 format ('262.0 <= minlon < maxlon < 402.0') 
      write (6, 6003)
 6003 format ('  0.0 <= minlat < maxlat <  75.0')

      read (5, *) aminlon, amaxlon, aminlat, amaxlat

      if (aminlon .lt. 262.0d0)  then
        write (6, 6111)
 6111   format ('minlon is less than 262.0; try again')
        go to 50
      else if (amaxlon .le. aminlon)  then
        write (6, 6112)
 6112   format ('maxlon is less than minlon; try again')
        go to 50
      else if (amaxlon .ge. 402.0d0)  then
        write (6, 6113)
 6113   format ('maxlon is larger than 402.0; try again')
        go to 50
      else if (aminlat .lt.   0.0d0)  then
        write (6, 6114)
 6114   format ('minlat is less than 0.0; try again')
        go to 50
      else if (amaxlat .le. aminlat)  then
        write (6, 6115)
 6115   format ('maxlat is less than minlat; try again')
        go to 50
      else if (amaxlat .ge.  75.0d0)  then
        write (6, 6116)
 6116   format ('maxlat is larger than 75.0; try again')
        go to 50
      end if

      imin = (aminlon-x0)*c24
      imax = (amaxlon-x0)*c24+1
      jmin = (aminlat-y0)*c24
      jmax = (amaxlat-y0)*c24+1

      write(*,*) ' imin,imax,jmin,jmax = ',imin,imax,jmin,jmax

      write (6, 6061)
 6061 format('Input filename of bathymetry master file')
      read (5, 5061) filename
 5061 format (a120)
      open (1, file=filename, status='old', form='formatted',
     $        access='sequential')
      rewind (1)

c      write (6, 6062)
 6062 format('Input filename of output file')
c      read (5, 5061) filename
c      open (2, file=filename, status='new', form='formatted')

      do j = 0, 1799
        if (j .gt. jmax)  then
          go to 999
        else if (j .ge. jmin  .and.  j .le. jmax)  then
          read (1, *) (igrid(i), i=0, 3359)
          do i = imin, imax
            hd(i-imin+1,j+1) = float(igrid(i))
c            print *,hd(i-imin+1,j+1),igrid(i)
c            write (2, 8101) dfloat(i)*c24div+x0, 
c     $           dfloat(j)*c24div+y0, igrid(i)
c 8101 format(2f12.5, i10)
          end do
        else
          read (1, *) jjj
        end if
      end do

      

 999  continue
      print *,'conversion ok'

      write(20,800)hd
800   format(10f10.0)
      close(1)
      close(2)

      stop

      end