program gpull

c - Change 9/15/1997 by D. Smith 
c - Change forces the output to have a header that makes sense
c - with regard to the original data set, rather than blindly
c - putting in the values which the user has indicated.

*** extract a subgrid from a grid

      implicit double precision(a-h,o-z)
      character*300 fname
      integer hrec(1000000)
      real*4  zrec(1000000)
      integer*2 hrec2(1000000)
      equivalence(hrec(1),zrec(1),hrec2(1))

      common/static/zrec

      lin =1
      lout=2

      write(*,*) 'program gpull, 6/22/2000'

      write(*,2)
    2 format('Enter input file:           ',$)
      read(*,1) fname
    1 format(a)
      open(lin,file=fname,status='old',form='unformatted')

      write(*,3)
    3 format('Enter subgrid output file:  ',$)
      read(*,1) fname
      open(lout,file=fname,status='new',form='unformatted')

*** display header contents

      read(lin) glamn,glomn,dgla,dglo,nla,nlo,ikind
      glamx=glamn+(nla-1)*dgla
      glomx=glomn+(nlo-1)*dglo
      write(*,*) 'kind=',ikind
      write(*,*) 'LAT  min=',glamn
      write(*,*) '     del=',dgla,'    # lat=',nla
      write(*,*) '     max=',glamx
      write(*,*) 'LON  min=',glomn
      write(*,*) '     del=',dglo,'    # lon=',nlo
      write(*,*) '     max=',glomx
      write(*,*)

   15 write(*,5)
    5 format('Enter minimum latitude  of subgrid: ',$)
      read(*,*) glamn2
      if(glamn2.lt.glamn.or.glamn2.gt.glamx) go to 15

      iz = 0
      del = glamn2 - glamn
      iq = nint(del/dgla)
      xq = del/dgla
      dq = abs(float(iq) - xq)
      if(dq.gt.1d-8)then
        i = (del / dgla) + 1
        glamn2 = glamn + i*dgla
        iz = 1
c     else
c       glamn2 = glamn2
      endif

   16 write(*,6)
    6 format('Enter maximum latitude  of subgrid: ',$)
      read(*,*) glamx2
      if(glamx2.lt.glamn2.or.glamx2.gt.glamx) go to 16

      del = glamx2 - glamn
      iq = nint(del/dgla)
      xq = del/dgla
      dq = abs(float(iq) - xq)
      if(dq.gt.1d-8)then
        i = (del / dgla)
        glamx2 = glamn + i*dgla
        iz = 1
c     else
c       glamx2 = glamx2
      endif

      iy = 0
   17 write(*,7)
    7 format('Enter minimum longitude of subgrid: ',$)
      read(*,*) glomn2
      if(glomn2.lt.glomn.or.glomn2.gt.glomx) go to 17

      del = glomn2 - glomn
      iq = nint(del/dglo)
      xq = del/dglo
      dq = abs(float(iq) - xq)
      if(dq.gt.1d-8)then
        i = (del / dglo) + 1
        glomn2 = glomn + i*dglo
        iy = 1
c     else
c       glomn2 = glomn2
      endif

   18 write(*,8)
    8 format('Enter maximum longitude of subgrid: ',$)
      read(*,*) glomx2
      if(glomx2.lt.glomn2.or.glomx2.gt.glomx) go to 18

      del = glomx2 - glomn
      iq = nint(del/dglo)
      xq = del/dglo
      dq = abs(float(iq) - xq)
      if(dq.gt.1d-8)then
        i = (del / dglo)
        glomx2 = glomn + i*dglo
        iy = 1
c     else
c       glomx2 = glomx2
      endif

      if(iz.eq.1)write(6,50)
      if(iy.eq.1)write(6,51)
   50 format(1x,'*Warning: Input LAT limits do not align with grid'
     * ' spacing.',/,1x,' Output grid will align with grid and fit'
     * ,' inside input LAT limits')
   51 format(1x,'*Warning: Input LON limits do not align with grid'
     * ' spacing.',/,1x,' Output grid will align with grid and fit'
     * ,' inside input LON limits')

      nla2=idnint((glamx2-glamn2)/dgla)+1
      nlo2=idnint((glomx2-glomn2)/dglo)+1

      ila0=idnint((glamn2-glamn)/dgla)
      ila1=ila0+1
      ila2=ila0+nla2
      ilo0=idnint((glomn2-glomn)/dglo)
      ilo1=ilo0+1
      ilo2=ilo0+nlo2

      write(*,*) 'OUTPUT:'
      write(*,*) 'kind=',ikind
      write(*,*) 'LAT  min=',glamn2
      write(*,*) '     del=',dgla,'    # lat=',nla2
      write(*,*) '     max=',glamx2
      write(*,*) 'LON  min=',glomn2
      write(*,*) '     del=',dglo,'    # lon=',nlo2
      write(*,*) '     max=',glomx2
      write(*,*)


      write(lout) glamn2,glomn2,dgla,dglo,nla2,nlo2,ikind

*** skip southernmost rows

      if(ikind.eq.0) then
        do 10 irow=1,ila0
   10   read(lin) (hrec(i),i=1,nlo)
      elseif(ikind.eq.1)then
        do 20 irow=1,ila0
   20   read(lin) (zrec(i),i=1,nlo)
      elseif(ikind.eq.2)then
        do 110 irow=1,ila0
  110   read(lin) (hrec2(i), i=1,nlo)
      endif

*** read and write the subgrid rows

      if(ikind.eq.0) then
        do 30 irow=1,nla2
          read (lin) (hrec(i),i=1,nlo)
          write(lout)(hrec(i),i=ilo1,ilo2)
   30   continue
      elseif(ikind.eq.1)then
        do 40 irow=1,nla2
          read (lin) (zrec(i),i=1,nlo)
          write(lout)(zrec(i),i=ilo1,ilo2)
   40   continue
      elseif(ikind.eq.2)then
        do 130 irow=1,nla2
          read (lin) (hrec2(i),i=1,nlo)
          write(lout)(hrec2(i),i=ilo1,ilo2)
  130   continue
      endif

      stop
      end