program venus_stereo2xyz
c
c  program to extract a sub array from an integer*2 Mercator grid
c
      implicit real*8 (a-h,o-z)
      common /grid/nlt,nln,dlt,dln,iproj
      common /bounds/rlt0,rltf,rln0,rlnf,rlnm
      common/info/lcrt(2),lin(10),nin,lout(10),nout
      integer*2 iin(21600) 
c
c  open the input file
c
      call inioc(lu1)
c
c open the output file
c
      call outio(lu2)
c
c setup the mapping from pixel to lon, lat
c
      call lambertd(rlt,rln,i,j,0)
c 
c read the file
c
      ntot=0
      do 200 i=1,nlt
      call rddiscb(lu1,iin,2*nln,istat)
      if(istat.lt.0) then
        write(lcrt(2),'(a)')' Problem reading input file '
        stop
      endif
      do 150 j=1,nln
c
c  skip empty cells
c
      idat=iin(j)
c     if(idat.eq.-32768) go to 150
      if(idat.le.0.) go to 150
      call lambertd(rltt,rlnt,i,j,-1)
      topo=(idat-11877)/1000.
c
c  write out latitude, longitude and gravity anomaly (mgal) or depth
c
      ntot=ntot+1
      write(lu2,907)rlnt,rltt,topo
  907 format(2f10.4,f8.3)
  150 continue
 200  continue
      call closio(-1)
      write(*,906)ntot
 906  format(' # of records in output file ',i8)
      stop
      end
c
c
      subroutine lambertd(rlt,rln,i,j,icall)
c
c  routine to compute the indices i,j associated with the
c  lambert projection of rlt,rln.
c
c  input
c   rlt   -   latitude (deg)
c   rln   -   longitude (deg)
c   icall -   0-set up grid parameters
c             1-calculate index
c
c output
c   i     -   row of matrix for rlt
c   j     -   column of matrix for rln
c
      implicit real*8 (a-h,o-z)
      common /grid/nlt,nln,dlt,dln,iproj
      common /bounds/rlt0,rltf,rln0,rlnf,rlnm
      common /info/lcrt(2),lin(10),nin,lout(10),nout
      data rad /.0174533/
      save rad
c
c  if icall equals 0 then get location parameters
c
      if(icall.eq.0) then
      write(lcrt(2),900)
  900 format(' # of lat, # of lon (both even): ',$)
      read(lcrt(1),*)nlt,nln
      write(lcrt(2),901)
  901 format(' minimum lat, latitude spacing (deg): ',$)
      read(lcrt(1),*)rlt0,dlt
      rltf=rlt0+dlt*nlt
      write(lcrt(2),902)
  902 format(' central longitude for grid (deg): ',$)
      read(lcrt(1),*)rlnm
      drnf=nln/2*dlt/cos(rltf*rad)
      rn1=rlnm-drnf
      rn2=rlnm+drnf
      drn0=nln/2*dlt/cos(rlt0*rad)
      rn3=rlnm-drn0
      rn4=rlnm+drn0
c
c  check to see if the left corners of the box are negative
c  and add 360. if it's true.
c
      if(rn1.lt.0..or.rn3.lt.0.) then
      rn1=rn1+360.
      rn2=rn2+360.
      rn3=rn3+360.
      rn4=rn4+360.
      rlnm=rlnm+360.
      endif
c
c  print corners of area
c
      write(lcrt(2),903)
 903  format('  corners of area  ')
      write(lcrt(2),904)rltf,rn1,rltf,rn2
      write(lcrt(2),904)rlt0,rn3,rlt0,rn4
  904 format(2f7.2,6x,2f7.2)
      write(lcrt(2),905)
  905 format(' continue?  1-yes  0-no: ',$)
      read(lcrt(1),*)iyes
      if(iyes.eq.0)stop
      return
      endif
c
c compute the indices of the point
c
      if(icall.eq.1) then
        rln1=rln
        i=(rltf-rlt)/dlt+1
        if(i.lt.1.or.i.gt.nlt) i=-1
  20    jd=int((rln1-rlnm)*cos(rlt*rad)/dlt)
        if(rln1-rlnm.lt.0.)jd=jd-1
        j2=nln/2+jd+1
        j=j2
c
c  check to see if the point lies to the left of the box
c
        if(j2.lt.1) then
        rln1=rln1+360.
        if(rln1.le.rlnf)go to 20
        if(j.lt.1.or.j.gt.nln) j=-1
        endif
c
c  compute latitude and longitude
c
      else
        if(i.lt.1.or.i.gt.nlt) then
        rlt=-999.
        return
        endif
        if(j.lt.1.or.j.gt.nln) then
        rln=-999.
        return
        endif
        rlt=rltf+dlt*(.5-i)
        term=dlt*(j-nln/2-.5)/cos(rad*rlt)
        rln=rlnm+term
      endif
      return
      end
c
      block data luio
c
c  common blocks for logical unit # transfers.
c
      common/info/lcrt(2),lin(10),nin,lout(10),nout
c
c   lcrt(1)  -  lu for terminal read.
c   lcrt(2)  -  lu for terminal write.
c   lin(i)   -  i'th input file lu (10-19 only).
c   nin      -  # of input files opened.
c   lout(i)  -  i'th output file lu (20-29 only).
c   nout     -  # of output files opened.
c
c***** calls no other routines
c
      data lcrt/5,6/
      data lin,nin/10*5,0/
      data lout,nout/10*6,0/
      end