c
c "getbin.f"	is run as part of the "step-2.csh" MODIS processing in
c		/Users/a2/data/modis/test" directory
c
c The embedded IDL/SeaDas command "out" creates 
c		4-byte integer table of signals
c		4-byte real table for latitude longitude
c
c Feb.-11, 2010 amue
c
	parameter (ny=7840, nx=5416, nt=1)
	parameter (nn=84922880, nnav=339691520)
	integer*4 a1(nx,ny),a2(nx,ny)
	real*4 xlat(nx,ny),xlon(nx,ny)
	real xmin,xmax,ymin,ymax
	character*14 input
	integer iyear,iyday,ihr,imin,iresolution
	double precision dx,dy,ds,ds2,dx2,dy2,deg2rad,arg
	double precision dx3,dy3,ds3
c
	ymin = 80.
	ymax = 83.
	xmin = -72.
	xmax = -52.
c
c output ascii data within the frame specified by [xmin,ymin,xmax,ymax]
c
	open (unit=16,file='loc')
c
c input control file created in step-2.csh via "hdfdump -h filename.L1B_QKM"
c
	open (unit=10, file='log.dat')
	read(10,*)input,refl1,refl2,rad1,rad2,ngrid
	write(6,110)input,refl1,refl2,rad1,rad2
110	format(a14,1x,e13.7,1x,e13.7,1x,f9.7,1x,f9.7)
	close (10)
c
c Reading FLAT binary files of a single L1B MODIS scene
c
c band-1:
	   open(unit=10,file='bin1',
     1		 form='unformatted',access='direct', recl=nx*4)
c
c i = [1,nx=5416] is the along-scan number (scan number) for 250-m data
c j = [1,ny=7840] is the along-track number (line number) for 250-m data
c
	   do 10 j=1,ny
	      read(10,rec=j)(a1(i,j),i=1,nx)
10	   continue
	   close (10)
c Location:
	   open(unit=10,file='nav',
     1		 form='unformatted',access='direct',recl=nx*4)
c
c Latitude:
c
	   do 12 j=1,ny
	      read(10,rec=j)(xlat(i,j),i=1,nx)
12	   continue
c
c Longitude:
c
	   do 11 j=ny+1,2*ny
	      read(10,rec=j)(xlon(i,j-ny),i=1,nx)
11	   continue
	   close (10)
c
	   write(6,*)input,': Done reading grid into memory'
c
c get resolution ds as a function of scan-number "i" at "j=1,ny/2,ny"
c
	iresolution = 0
	if (iresolution .eq. 1) then
	  deg2rad = atan2(1.,1.)/45.
	  do 33 i=2,nx-1
	     j=ny/2
	     arg = deg2rad*xlat(i,j)
	     dx = (xlon(i-1,j)-xlon(i+1,j))*dcos(arg) 
	     dy = (xlat(i-1,j)-xlat(i+1,j))
	     ds = 1852.*dsqrt(dx*dx+dy*dy)*0.5*60.
	     j=1
	     arg = deg2rad*xlat(i,j)
	     dx2 = (xlon(i-1,j)-xlon(i+1,j))*dcos(arg) 
	     dy2 = (xlat(i-1,j)-xlat(i+1,j))
	     ds2 = 1852.*dsqrt(dx2*dx2+dy2*dy2)*0.5*60.
	     j = ny
	     arg = deg2rad*xlat(i,j)
	     dx3 = (xlon(i-1,j)-xlon(i+1,j))*dcos(arg) 
	     dy3 = (xlat(i-1,j)-xlat(i+1,j))
	     ds3 = 1852.*dsqrt(dx3*dx3+dy3*dy3)*0.5*60.
	     write(9,109)i-nx/2,j,xlon(i,j),xlat(i,j),ds,ds2,ds3
109	     format(2i5,2f9.3,3f10.3)
33	  continue
	end if
c
c select sub-setting frame
c
	k = 0
c	do 31 i=1,nx
c
c Limited to better than 500-m resolution (determined from above)
c
	do 31 i=nx/2-2000,nx/2+2000
	   do 32 j=1,ny
		if (xlat(i,j).gt.ymin) then
		  if (xlat(i,j).lt.ymax) then
		    if (xlon(i,j).gt.xmin) then
		      if (xlon(i,j).lt.xmax) then
c
c Scale dimensionless integers to dimensionless reflectance in [0,1]
c	(use refl1 and refl2 to get reflectance) 
c	(use rad1 and rad2 to get radiances in W/m^2/u/sr)
c
		xband1 = a1(i,j)*refl1
		write(16,101)xlat(i,j),xlon(i,j),xband1
		k = k+1
		      end if
		    end if
		   end if
		end if
32	   continue
31	continue
	write(99,103)input,k,xlon(nx/2,ny/2),xlat(nx/2,ny/2)
101	format(f8.3,f9.3,f7.4)
102	format(a18,1x,i4,i3,2i2)
103	format(a18,i12,2f9.3)
	stop
	end