c (c) copyright 2003 David G. Long, Brigham Young University
	program	egg
c
c	Example program to illustrate how to use BYU egg shape tables
c
c	Written: D. Long   23 June 2003
c
c       note: this code was written on a machine running HPUX 11 using
c       the HP f90 compiler
c       modifications may be required to get it to run on your machine
c
c       code should be linked to local HDF installation
c
	parameter (version = 1.0)
	character*30 prog_name
	parameter (prog_name = 'egg')
c
c L1B Variables
c	
        parameter (MAX_VAR_DIMS = 3)
	parameter (MAX_BUF_SIZE = 100*8*8*2)
c
	integer         sd_id,nsets,nattr
        integer         nrank,dims(MAX_VAR_DIMS)
	integer*2       buffer(MAX_BUF_SIZE)
c
c	Data product SDS items to be extracted are in l1b_frame.inc
c
ccccccccccccccccccccccccccccccccc l1b.inc cccccccccccccccccccccccccccccccccccccccc
c
c     Declarations for contents of the QuikSCAT Level 1B product      
c
c       (1) UTC Timetags are in 1 Vdata dataset
	
c       (2) Data in global attributes to be extracted:
c       
c           l1b_actual_frames   (#30)
c           sigma0_kpc_b(2)     (#34)
c           slice_kpc_b(12)     (#35)
c           sigma0_kpc_c(2)     (#36)
c           slice_kpc_c(12)     (#37)
c           receiver_gain_ratio (#38)
c       
c       (3) SDS Data Set info:
c	
c  Dataset/Name          Rank      Dimensions        Scale      Offset   Type Conv
c ---------------------------------------------------------------------------------
c       orbit_time         1   11250     0     0     1.00000    0.00000  uint4 int    
c    0  frame_inst_status  1   11250     0     0     1.00000    0.00000  uint4 int
c    1  frame_err_status   1   11250     0     0     1.00000    0.00000  uint4 int
c    2  frame_qual_flag    1   11250     0     0     1.00000    0.00000  uint2 int
c    3  num_pulses         1   11250     0     0     1.00000    0.00000  byte8 int
c    4  sc_lat             1   11250     0     0     1.00000    0.00000  real4 real
c    5  sc_lon             1   11250     0     0     1.00000    0.00000  real4 real
c    6  sc_alt             1   11250     0     0     1.00000    0.00000  real4 real
c    7  x_pos              1   11250     0     0     1.00000    0.00000  real4 real
c    8  y_pos              1   11250     0     0     1.00000    0.00000  real4 real
c    9  z_pos              1   11250     0     0     1.00000    0.00000  real4 real
c   10  x_vel              1   11250     0     0     1.00000    0.00000  real4 real
c   11  y_vel              1   11250     0     0     1.00000    0.00000  real4 real
c   12  z_vel              1   11250     0     0     1.00000    0.00000  real4 real
c   13  roll               1   11250     0     0     0.00100    0.00000  int*2 real
c   14  pitch              1   11250     0     0     0.00100    0.00000  int*2 real
c   15  yaw                1   11250     0     0     0.00100    0.00000  int*2 real
c   16  bandwidth_ratio    1   11250     0     0     0.00100    0.00000  int*2 real
c   17  x_cal_A            1   11250     0     0     0.01000    0.00000  int*2 real
c   18  x_cal_B            1   11250     0     0     0.01000    0.00000  int*2 real
c   19  cell_lat           2     100 11250     0     1.00000    0.00000  real4 real
c   20  cell_lon           2     100 11250     0     1.00000    0.00000  real4 real
c   21  sigma0_mode_flag   2     100 11250     0     1.00000    0.00000  uint2 int
c   22  sigma0_qual_flag   2     100 11250     0     1.00000    0.00000  uint2 int
c   23  slice_qual_flag    2     100 11250     0     1.00000    0.00000  uint4 int
c   24  sigma0             2     100 11250     0     0.01000    0.00000  int*2 real
c   25  pr_pt_ratio        2     100 11250     0     0.01000    0.00000  int*2 real
c   26  frequency_shift    2     100 11250     0     1.00000    0.00000  int*2 real
c   27  cell_area          2     100 11250     0                0.00000  int*2 real
c   28  cell_azimuth       2     100 11250     0     0.01000    0.00000  uint2 real
c   29  cell_incidence     2     100 11250     0     0.01000    0.00000  int*2 real
c   30  antenna_azimuth    2     100 11250     0     0.01000    0.00000  uint2 real
c   31  snr                2     100 11250     0     0.01000    0.00000  int*2 real
c   32  kpc_a              2     100 11250     0     0.00010    0.00000  int*2 real
c   33  slice_lat          3       8   100 11250     0.00100    0.00000  int*2 real
c   34  slice_lon          3       8   100 11250     0.00100    0.00000  int*2 real
c   35  slice_sigma0       3       8   100 11250     0.01000    0.00000  int*2 real
c   36  x_factor           3       8   100 11250     0.01000    0.00000  int*2 real
c   37  slice_azimuth      3       8   100 11250     0.01000    0.00000  uint2 real
c   38  slice_incidence    3       8   100 11250     0.01000    0.00000  int*2 real
c   39  slice_snr          3       8   100 11250     0.01000    0.00000  int*2 real
c   40  slice_kpc_a        3       8   100 11250     0.00010    0.00000  int*2 real
c ---------------------------------------------------------------------------------

	parameter (MAX_FRAMES = 2)
	parameter (MAX_PULSES = 100)
	parameter (MAX_SLICES = 8)

	character*21 TimeTags(11250) ! in vdata

	integer nframes                 ! in global attributes
	
        integer orbit_time(MAX_FRAMES)
	integer frame_inst_status(MAX_FRAMES) ! in SDS datasets
	integer frame_err_status(MAX_FRAMES)
	integer*2 frame_qual_flag(MAX_FRAMES)
	byte num_pulses(MAX_FRAMES)
c	real sc_lat(MAX_FRAMES),sc_lon(MAX_FRAMES),sc_alt(MAX_FRAMES)
c	real x_pos(MAX_FRAMES),y_pos(MAX_FRAMES),z_pos(MAX_FRAMES)
c	real x_vel(MAX_FRAMES),y_vel(MAX_FRAMES)
	real z_vel(MAX_FRAMES)
c	real roll(MAX_FRAMES),pitch(MAX_FRAMES),yaw(MAX_FRAMES)
c	real bandwidth_ratio(MAX_FRAMES)
c	real x_cal_A(MAX_FRAMES)
c	real x_cal_B(MAX_FRAMES)
	real cell_lat(MAX_PULSES,MAX_FRAMES)
	real cell_lon(MAX_PULSES,MAX_FRAMES)
	integer*2 sigma0_mode_flag(MAX_PULSES,MAX_FRAMES)
	integer*2 sigma0_qual_flag(MAX_PULSES,MAX_FRAMES)
c	integer slice_qual_flag(MAX_PULSES,MAX_FRAMES)
	real sigma0(MAX_PULSES,MAX_FRAMES)
cc	real pr_pt_ratio(MAX_PULSES,MAX_FRAMES)
c	real frequency_shift(MAX_PULSES,MAX_FRAMES)
c 	real cell_area(MAX_PULSES,MAX_FRAMES)
 	real cell_azimuth(MAX_PULSES,MAX_FRAMES)
 	real cell_incidence(MAX_PULSES,MAX_FRAMES)
 	real antenna_azimuth(MAX_PULSES,MAX_FRAMES)
c 	real snr(MAX_PULSES,MAX_FRAMES)
c 	real kpc_a(MAX_PULSES,MAX_FRAMES)
c 	real slice_lat(MAX_SLICES,MAX_PULSES,MAX_FRAMES)
c 	real slice_lon(MAX_SLICES,MAX_PULSES,MAX_FRAMES)
c 	real slice_sigma0(MAX_SLICES,MAX_PULSES,MAX_FRAMES)
c 	real x_factor(MAX_SLICES,MAX_PULSES,MAX_FRAMES)
c  	real slice_azimuth(MAX_SLICES,MAX_PULSES,MAX_FRAMES)
c  	real slice_incidence(MAX_SLICES,MAX_PULSES,MAX_FRAMES)
c  	real slice_snr(MAX_SLICES,MAX_PULSES,MAX_FRAMES)
c 	real slice_kpc_a(MAX_SLICES,MAX_PULSES,MAX_FRAMES)
c ---------------------------------------------------------------------------------

	character    tmpstr*24
	character    product*8
	character    rev_number*8
	character    begin_date*8,begin_time*12
	character    end_date*8,end_time*12
	character*80 values(20)
c
	real         sigma,thetai,azang
c
c       local variables
c
	integer dstart,dend,mstart,mend,year
	integer iregion, file_id
	integer iday,iyear,ihour,imin,ipolar
c
c local variables
c
        integer     ibeam,iasc,ib
	integer     irec,jrec,krec,firstrec,ios,n
c
c more local variables
c
	integer       argc,argn
       	character*180 fname,mname,line
	character*180 pathshape, innershapefile, outershapefile
	real          coef(36)
	real          response_array(300),clon_array(300),clat_array(300)
c
c flag masks
c
	integer*2     frame_qual_mask, sigma_mask, qual_mask
	integer*4     frame_mask, inst_mask, slice_mask
c
c *************************************
c
c	begin program
c
	irev = 0 ! rev counter
        irec = 0 ! input cells counter
        jrec = 0 ! output record counter
c
c       set mask flags
c
	call set_masks(frame_mask, frame_qual_mask, inst_mask, sigma_mask, qual_mask, slice_mask)
c
	write (*,*) 'SeaWinds Egg Demo Program'
	write (*,*) 'Program: ',prog_name,' Version ',version
	write (*,*)
c
c     get meta file name
c
	argc = iargc()		! number of command line arguments
	argn = 1		! command line argument counter
c
	if (argc.lt.argn) then
           write (*,1) 'Input L1B name?'
 1         format(1x,a,1x,$)
           read (*,'(a120)') fname
	else
	  call getarg(argn,fname)
	  argn = argn + 1
	endif
c
c       set instrument mode storage variable to force read of egg shape table
c
	instmode_store=-1
c
c       get environment variable pointing to director of tabular cell shape files
c
	call getenv('BYU_SHAPES_DIR',pathshape)
	print *,'BYU_SHAPES_DIR environment variable=',pathshape
	if (length(pathshape).lt.2) pathshape='./'
c
	krec=0              ! total frames read
c
c     open new input L1B file
c
	write (*,*) 'Opening input file "',fname(1:length(fname)),'"'
	sd_id=0
	call open_hdf(fname(1:length(fname)),sd_id,nsets,nattr)
	print *,'past open_hdf call',sd_id,nsets,nattr
	if (sd_id.lt.0) then
	   print *,'*** error opening hdf file ***'
	   stop
	endif
c
c     check file type
c
	call read_attrib_byname(sd_id,'ShortName',ntype,nval,product)
	if (nval.lt.0) then
	   print *,'*** Defective input file: not a QuikSCAT file'
	   print *,'*** skipping file'
	   call close_hdf(sd_id)
	   stop
	endif
	if (product .ne. 'QSCATL1B') then
	   print *,'*** The input file is not a QuikSCAT Level 1B file'
	   print *,'*** product: "',product,'"'
	   print *,'*** skipping file'
	   call close_hdf(sd_id)
	   stop
	endif
c
c       Get data from global attributes:
c
	call read_attrib_byname(sd_id,'rev_number',ntype,nval,rev_number)
	call read_attrib_byname(sd_id,'l1b_actual_frames',ntype,nval,values)
c
c       decode number of available frames
c
	read(values(1),'(i9)') nframes
c
c       read and decode orbit period
c
	call read_attrib_byname(sd_id,'rev_orbit_period',ntype,nval,values)
	read(values(1),'(f16.9)') rev_orbit_period
	if (rev_orbit_period.lt.100) rev_orbit_period=6061.0
c
	call read_attrib_byname(sd_id,'RangeBeginningDate',ntype,nval,begin_date)
	call read_attrib_byname(sd_id,'RangeBeginningTime',ntype,nval,begin_time)
c
	call read_attrib_byname(sd_id,'RangeEndingDate',ntype,nval,end_date)
	call read_attrib_byname(sd_id,'RangeEndingTime',ntype,nval,end_time)
c	
	print *
	print *,'Product type: ',product,'  Frames:',nframes
	print *,'Rev number:   ',rev_number,'  Orbit period:',rev_orbit_period
	tmpstr(1:21) = begin_date(1:8)//'T'//begin_time(1:12)
	print *,'Start time:   ',tmpstr(1:21)
	call convert_time(tmpstr(1:22),iyear,iday,ihour,imin)
	tmpstr(1:21) = end_date(1:8)//'T'//end_time(1:12)
	call convert_time(tmpstr(1:22),iyeare,idaye,ihoure,imine)
	print *,'End time:     ',tmpstr(1:21)
c
c       get UTC timetages
c
	call read_timetags(fname(1:length(fname)), nframes1, TimeTags)
	print *,'nframes from read_timetags:',nframes1
c
c       open output file
c
	open(unit=1,name='egg.DAT',status='unknown',form='formatted')
c
	firstrec=1
	iend=nframes
	nrec=0
c
c       for each frame in file
c
	print *,'Begin processing frames',firstrec,iend
	do 350 iframe=firstrec,iend
c
c       extract SDS data from input HDF file one frame at a time
c
	   call extract_int_sds(sd_id,'orbit_time',iframe,1,nrank,dims,orbit_time)
	   call extract_int_sds(sd_id,'frame_inst_status',iframe,1,nrank,dims,frame_inst_status)
	   call extract_int_sds(sd_id,'frame_err_status',iframe,1,nrank,dims,frame_err_status)
	   call extract_int_sds(sd_id,'frame_qual_flag',iframe,1,nrank,dims,frame_qual_flag)
	   call extract_int_sds(sd_id,'num_pulses',iframe,1,nrank,dims,num_pulses)
c	   call extract_sds2(sd_id,'sc_lat',iframe,1,nrank,dims,sc_lat,buffer)
c	   call extract_sds2(sd_id,'sc_lon',iframe,1,nrank,dims,sc_lon,buffer)
c	   call extract_sds2(sd_id,'sc_alt',iframe,1,nrank,dims,sc_alt,buffer)
c	   call extract_sds2(sd_id,'x_pos',iframe,1,nrank,dims,x_pos,buffer)
c	   call extract_sds2(sd_id,'y_pos',iframe,1,nrank,dims,y_pos,buffer)
c	   call extract_sds2(sd_id,'z_pos',iframe,1,nrank,dims,z_pos,buffer)
c	   call extract_sds2(sd_id,'x_vel',iframe,1,nrank,dims,x_vel,buffer)
c	   call extract_sds2(sd_id,'y_vel',iframe,1,nrank,dims,y_vel,buffer)
	   call extract_sds2(sd_id,'z_vel',iframe,1,nrank,dims,z_vel,buffer)
c	   call extract_sds2(sd_id,'roll',iframe,1,nrank,dims,roll,buffer)
c	   call extract_sds2(sd_id,'pitch',iframe,1,nrank,dims,pitch,buffer)
c	   call extract_sds2(sd_id,'yaw',iframe,1,nrank,dims,yaw,buffer)
c	   call extract_sds2(sd_id,'bandwidth_ratio',iframe,1,nrank,dims,bandwidth_ratio,buffer)
c	   call extract_sds2(sd_id,'x_cal_A',iframe,1,nrank,dims,x_cal_A,buffer)
c	   call extract_sds2(sd_id,'x_cal_B',iframe,1,nrank,dims,x_cal_B,buffer)
	   call extract_sds2(sd_id,'cell_lat',iframe,1,nrank,dims,cell_lat,buffer)
	   call extract_sds2(sd_id,'cell_lon',iframe,1,nrank,dims,cell_lon,buffer)
	   call extract_int_sds(sd_id,'sigma0_mode_flag',iframe,1,nrank,dims,sigma0_mode_flag)
	   call extract_int_sds(sd_id,'sigma0_qual_flag',iframe,1,nrank,dims,sigma0_qual_flag)
c	   call extract_int_sds(sd_id,'slice_qual_flag',iframe,1,nrank,dims,slice_qual_flag)
	   call extract_sds2(sd_id,'cell_sigma0',iframe,1,nrank,dims,sigma0,buffer)
c	   call extract_sds2(sd_id,'pr_pt_ratio',iframe,1,nrank,dims,pr_pt_ratio,buffer)
c	   call extract_sds2(sd_id,'frequency_shift',iframe,1,nrank,dims,frequency_shift,buffer)
c	   call extract_sds2(sd_id,'cell_area',iframe,1,nrank,dims,cell_area,buffer)
	   call extract_sds2(sd_id,'cell_azimuth',iframe,1,nrank,dims,cell_azimuth,buffer)
	   call extract_sds2(sd_id,'cell_incidence',iframe,1,nrank,dims,cell_incidence,buffer)
	   call extract_sds2(sd_id,'antenna_azimuth',iframe,1,nrank,dims,antenna_azimuth,buffer)
c
	   krec=krec+1		! count total frames read
	   nrec=nrec+1		! count frames read in file
c
	   if (mod(krec,100).eq.0) write (*,446) krec,irec,jrec
 446	   format('Frames ',I7,1x,'| Pulses ',I9,1x,'| Eggs ',I9)
c
c       convert time
c
	   call convert_time(TimeTags(iframe),iyear,iday,ihour,imin)
c
c       check instrument status, frame_err, and frame_qual
c
	   if (iand(frame_err_status(1),frame_mask).ne.0.or. ! eliminate all frame
     *		iand(frame_inst_status(1),inst_mask).ne.0.or. ! errors conservatively 
     *		iand(frame_qual_flag(1),frame_qual_mask).ne.0) goto 350	! skip frame
c
c       get instrument gate width mode from flag for start of frame
	   inst_mode=mod(frame_inst_status(1)/16,8)
	   if (inst_mode.ne.instmode_store) then
	      print *,'*** gate width mode change to:',inst_mode
	      instmode_store=inst_mode
	      write(innershapefile,'(A,''eggp_in'',I1,''.dat'')') pathshape(1:length(pathshape)),instmode_store+1
	      write(outershapefile,'(A,''eggp_out'',I1,''.dat'')') pathshape(1:length(pathshape)),instmode_store+1
c       
c	Load egg polynomial shape data
c
	      print *,'Reading egg shape data files for mode ',instmode_store
	      print *,innershapefile(1:length(innershapefile))
	      print *,outershapefile(1:length(outershapefile))
	      call read_eggshape_tables(innershapefile,outershapefile)
	      print *,'Shape files read'
	   endif
c
c       set asc/dsc flag for measurements
c
	   if (z_vel(1).ge.0.0) then
	      iasc=1
	   else
	      iasc=0
	   endif
c
c       for each pulse in frame
c
	   np=num_pulses(1)
	   do 3401 i=1,np
c
c     check sigma0_mode, slice quality
c
	      if (iand(sigma0_mode_flag(i,1),sigma_mask).ne.0) goto 3401
	      irec=irec+1
c
c       for demo, keep only measurements near equator
c
	if (abs(cell_lat(i,1)).lt.1.0) then
c
c       egg center
c
	      cen_lat=cell_lat(i,1)
	      cen_lon=cell_lon(i,1)
	      if (cen_lon.eq.0.0.and.cen_lat.eq.0.0) goto 3401 ! discard bad measurements
c
c       get beam (and polarization) of pulse
c
	      ibeam=mod(sigma0_mode_flag(i,1)/4,2) ! 0=inner, 1=outer
c
c       get info required to use tabularized cell shape
c
	      ant_az=antenna_azimuth(i,1)
c
c       note the orbit timer is a 32 Hz clock so period is 31.25 ms/tic
c
	      time_orb=orbit_time(1)*0.03125 ! convert to seconds
c
	      ib=ibeam+1	! 1=H (inner), 2=V (outer) beam
c
c       ibb odd=aft-looking, ibb even=forward-looking
	      ibb=2*ibeam
	      if (ant_az.gt.90.0.and.ant_az.lt.270.0) ibb=ibb+1
c
c       get coefficients of 2d polynomial which describes egg response
c
	      call get_egg_coef(time_orb,ant_az,ib,coef,morder,msize,rev_orbit_period)
c
c       compute normalized response function at each point of a grid in lat/lon
c       surrounding the center point of the egg.
c
c       set the grid spacing to be about 6 km
	      dellat=6.0/110.0
	      dellon=dellat/cosd(cen_lat);
c       and about 30 km in radius
	      isize=5
	      nsize=(2*isize+1)*(2*isize+1)
c
	      count=0
	      do iy=-isize,isize
		 dlat=dellat*iy
		 do ix=-isize,isize
		    dlon=dellon*ix
c
c     compute egg spatial response (only defines the main lobe of response)
c
		    sum=0.0
		    inc=0
		    do im=1,morder
		       do in=1,morder
			  if (im+in-1.le.morder) then
			     inc=inc+1
			     sum=sum+coef(inc)*dlat**(im-1)*dlon**(in-1)
			  endif
		       end do
		    end do
c		    sum=10.**(0.1*sum) ! convert to normal space
		    count=count+1
		    response_array(count)=sum
c
c       set lat/lon of grid center and store
c
		    clat=cen_lat+dellat*iy
		    clon=mod(cen_long+dellon*ix+360.0,360.0)
		    clon_array(count)=clon
		    clat_array(count)=clat
		 enddo
	      end do
c
c       write out to a file (standard unit 1)
c
	      write (1,1400) iframe,iasc,i-1,ant_az,ibeam,
	1	   cen_lon,cen_lat,
	2	   (response_array(j),clon_array(j),clat_array(j),j=1,nsize)
 1400	      format(i5,1x,i1,1x,i2,1x,f6.2,1x,i1,1x,f7.2,1x,f7.2,
	1	   <nsize>(1x,f7.2,1x,f7.2,1x,f7.2))
c       
	      jrec=jrec+1	! a count of total records written
c
c       get sign of sigma-0 value
c       
	      if (btest(sigma0_qual_flag(i,1),2).eq.1) then
		 isgn=-1
	      else
		 isgn=1
	      endif
	      sigma=sigma0(i,1)
	      thetai=cell_incidence(i,1)
	      azang=cell_azimuth(i,1)
c
        endif
c
 3401	   continue		! i loop: pulses
c
c     process next record
c
 350	continue		! iframe loop: frames
c
c     input file has been processed, close it
c
	call close_hdf(sd_id)
c
c       close output file
c
	close(1)

	write (*,*)
	write (*,*) 'Total input frames:  ',krec,'  Total input pulses: ',irec
	write (*,*) 'Output Eggs:',jrec
	write (*,*)
c
	write (*,*) prog_name,' completed successfully'
c
	stop
	end	
c
c ***************************************************************************
c
      subroutine convert_time(time_tag,iyear,iday,ihour,imin)
c
c       convert Seawinds ascii time tag into year, day, hour, minute
c
      character*(*) time_tag

      read(time_tag,10,err=20) iyear,iday,ihour,imin
 10   format(I4,1x,I3,1x,I2,1x,I2)
c      write (*,*) 'time_tag="',time_tag,'" ',iyear,iday,ihour,imin
      return
 20   iyear=0    ! encountered a header record or other error
      write (*,*) '*** time_tag err="',time_tag,'" ',iyear
      return
      end

c
c ***************************************************************************
c
c       set mask flags
c
	subroutine set_masks(frame_mask, frame_qual_mask, inst_mask, 
     $	    sigma_mask, qual_mask,slice_mask)
	integer*2     frame_qual_mask, sigma_mask, qual_mask
	integer*4     frame_mask, inst_mask, slice_mask
c
c       allow all bits
c
	frame_mask=-1  ! set all bits
	frame_qual_mask=-1  ! set all bits
	inst_mask=-1   ! set all bits
	sigma_mask=-1  ! set all bits
	qual_mask=-1   ! set all bits
	slice_mask=-1  ! set all bits
c
c       mask off the don't care bits
c
c frame_err_status
c	frame_mask=ibclr(frame_mask,12)    ! clear SES Data Overrun Flag #?#
c	frame_mask=ibclr(frame_mask,13)    ! clear SES Data Underrun Flag #?#
c	frame_mask=ibclr(frame_mask,14)    ! clear Payload Bus Interface (PBI) Flag #?#
	frame_mask=ibclr(frame_mask,23)    ! clear Attitude Data Flag #?#
	frame_mask=ibclr(frame_mask,24)    ! clear Ephemeris Data Flag #?#
c
c frame_qual_flag
	frame_qual_mask=ibclr(frame_qual_mask,0)    ! clear frame filler bit 1
	frame_qual_mask=ibclr(frame_qual_mask,1)    ! clear frame filler bit 2
	frame_qual_mask=ibclr(frame_qual_mask,2)    ! clear frame CRC flag bit 1
	frame_qual_mask=ibclr(frame_qual_mask,3)    ! clear frame CRC flag bit 2
c
c frame_inst_status
	inst_mask=ibclr(inst_mask,2)    ! clear first pulse count bit
	inst_mask=ibclr(inst_mask,4)    ! clear gate width bits
	inst_mask=ibclr(inst_mask,5)    ! clear gate width bits
	inst_mask=ibclr(inst_mask,6)    ! clear gate width bits
	inst_mask=ibclr(inst_mask,7)    ! clear hires/egg acquisition flag
	inst_mask=ibclr(inst_mask,8)    ! clear begin calibration bit
	inst_mask=ibclr(inst_mask,9)    ! clear SES select bit
	inst_mask=ibclr(inst_mask,10)   ! clear SAS select bit
	inst_mask=ibclr(inst_mask,11)   ! clear TWTA select bit
	inst_mask=ibclr(inst_mask,23)   ! clear TWTA monitor enable bit #?#
c
c sigma0_mode_flag
	sigma_mask=ibclr(sigma_mask,2)  ! clear beam flag bit
	sigma_mask=ibclr(sigma_mask,3)  ! clear location bit
	sigma_mask=ibclr(sigma_mask,6)  ! clear gate width bits
	sigma_mask=ibclr(sigma_mask,7)  ! clear gate width bits
	sigma_mask=ibclr(sigma_mask,8)  ! clear gate width bits
	sigma_mask=ibclr(sigma_mask,10) ! clear SES select bit
c	sigma_mask=ibclr(sigma_mask,11) ! clear spin rate bit
c
c sigma0_qual_flag
c	qual_mask=ibclr(qual_mask,0)  ! clear sigma0 useable bit #?#
	qual_mask=ibclr(qual_mask,1)  ! clear low snr bit #?#
	qual_mask=ibclr(qual_mask,2)  ! clear negative sigma-0 bit
	qual_mask=ibclr(qual_mask,3)  ! clear sigma0 out of range bit #?#
c
c slice_qual_flag 
	do i=1,8
	   j=(i-1)*4
	   slice_mask=ibclr(slice_mask,j+1) ! clear sign flag bit
	   slice_mask=ibclr(slice_mask,j+0) ! clear gain threshold bit #?#
	   slice_mask=ibclr(slice_mask,j+2) ! clear low SNR flag bit #?#
c	   slice_mask=ibclr(slice_mask,j+3) ! clear location bit #?#
	end do
c
	print 40,inst_mask,sigma_mask,slice_mask
 40	format(' Masks: ',z9,2x,z6,2x,z9)
	return
	end
c
c ******************************************************************
c
c       egg shape routines
c
c ******************************************************************
c
      subroutine read_eggshape_tables(f1,f2)
      character*(*) f1,f2
c
c     read QuikScat egg shape tables into common blocks
c
c       DGL 11/11/98
c
c       inputs
c     f1: file name for inner beam egg polynomial shape table file
c     f2: file name for outer beam egg polynomial shape table file
c
      dimension b(16)
      common /shape_table/norbs,nazs,morder(2),msize(2),orbtimes(36),
     $     coefs(16,36,36,2)
c
      norbs=36
      nazs=36
      ib=1
      open(unit=4,name=f1,status='old',form='formatted')
      do iorb=1,norbs
         do iaz=1,nazs
            read(4,*) rlong,rdlat,gmt,az_nom,clon,clat,morder(ib),msize(ib),
     $           (b(i),i=1,msize(ib))
            do i=1,msize(ib)
               coefs(i,iaz,iorb,ib)=b(i)
            end do
         end do
         orbtimes(iorb)=gmt
      end do
      close(4)
c
      ib=2
      open(unit=4,name=f2,status='old',form='formatted')
      do iorb=1,norbs
         do iaz=1,nazs
            read(4,*) rlong,rdlat,gmt,az_nom,clon,clat,morder(ib),msize(ib),
     $           (b(i),i=1,msize(ib))
            do i=1,msize(ib)
               coefs(i,iaz,iorb,ib)=b(i)
            end do
         end do
         orbtimes(iorb)=gmt
      end do
      close(4)
c
      return
      end
c
c
c
      subroutine get_egg_coef(time,az,ib,coef,morder1,msize1,orb_period)
      dimension coef(*)
c
c     given the orbit time, azimuth angle, beam and slice number,
c     compute QuikScat shape using stored table
c
c       DGL 11/11/98
c       revised DGL 7/12/99 to incorporate orbit period as passed parameter
c                           and to normalize orbit period by table period and
c                           to fix orbit time near end of orbit
c
c       inputs: 
c     time:    orbit time (in sec) from orbit crossing
c     az:      antenna azimuth angl
c     ib:      antenna beam number 1=inner, 2=outer
c     orbit_period: orbit period (in sec)
c
c       outputs:
c     coef:    2d polynomial coefficients
c     morder1: polynomial order
c     msize1:  number of coefficients
c
      common /shape_table/norbs,nazs,morder(2),msize(2),orbtimes(36),
     $     coefs(16,36,36,2)
c
      data table_orbit_period/6061.0/
c
c     linear interpolation function
c
      rinterp(xx,yy,r1,r2,r3,r4)=(r2-r1)*yy+r1+((r4-r3)*yy+r3
     $     -(r2-r1)*yy-r1)*xx
c
      t=table_orbit_period*(time/orb_period)
      if (t.lt.0.0) t=t+table_orbit_period
      if (t.gt.table_orbit_period) t=t-table_orbit_period
      if (t.gt.table_orbit_period) t=t-table_orbit_period
      do i=1,norbs
         j=i+1
         if (j.gt.norbs) then
	    j=j-norbs
	    if (t.ge.orbtimes(i).and.t.lt.orbtimes(j)+table_orbit_period) then
	       io1=i
	       io2=j
	       goto 10
	    endif
	 else
	    if (t.ge.orbtimes(i).and.t.lt.orbtimes(j)) then
	       io1=i
	       io2=j
	       goto 10
	    endif
	 endif
      end do
      print *,'*** orbit time error in get_egg_coef',time,t,orb_period,norbs
      stop
c
 10   continue
      delo=orbtimes(io2)-orbtimes(io1)
      if (io2.lt.io1) delo=delo+orb_period
      delo=(t-orbtimes(io1))/delo
c
      a=mod(az+720.0,360.0)/10.0
      ia1=ifix(a)+1
      ia2=ia1+1
      if (ia2.gt.nazs) ia2=ia2-nazs
      dela=a-ia1+1
c
      morder1=morder(ib)
      msize1=msize(ib)
c
c     bilinearly interpolate the various coefficients
c
      do ic=1,msize1
         coef(ic)=rinterp(dela,delo,coefs(ic,ia1,io1,ib),
     $        coefs(ic,ia1,io2,ib),coefs(ic,ia2,io1,ib),
     $        coefs(ic,ia2,io2,ib))
      end do
c
      return
      end
c
c  *************************************************************************
c
c  HDF reading code 
c
c **************************************************************************
c
      subroutine open_hdf(filename,sd_id,n_datasets,n_file_attr)
      
      ! open an HDF file to read and return the SD_ID, number of
      ! SD datasets, and number of file attributes to the calling
      ! program
      
      ! 05-13-98  rsd
      
      parameter (DFACC_RDONLY = 1)
      character filename*(*)
      integer sd_id
      integer n_datasets,n_file_attr
      integer sfstart
      
      sd_id = sfstart(filename,DFACC_RDONLY)
      call sffinfo(sd_id, n_datasets, n_file_attr)
      
      return
      end


      subroutine extract_sds(sd_id,in_name,irec,nrec,nrank,dims,real_array)
      
      ! extract the contents of an SDS from an HDF file
      ! 
      ! 5-12-98  rsd
      !
      parameter (MAX_VAR_DIMS = 3, MAX_NC_NAME=256)
      integer index,nrank,dims(MAX_VAR_DIMS)
c      parameter (MAX_BUF_SIZE = 500000)
      parameter (MAX_BUF_SIZE = 10000000)
      
      integer irec,nrec
      integer sds_id, sd_id, retn
      integer sfselect, sfendacc,sfginfo,sfgcal
      integer sfn2index
      character*(*) in_name
      integer start(MAX_VAR_DIMS),edge(MAX_VAR_DIMS),stride(MAX_VAR_DIMS)
      integer dim_sizes(MAX_VAR_DIMS)
      character name*(MAX_NC_NAME)
      integer rank,attributes,num_type
      double precision cal,cal_error,off,off_err

      integer*2 buffer(MAX_BUF_SIZE)
      real real_array(MAX_BUF_SIZE)
      
      index = sfn2index(sd_id, in_name)
      sds_id = sfselect(sd_id, index)   ! get SDS id
      retn = sfginfo(sds_id,name,rank,dim_sizes,num_type,attributes)
      retn = sfgcal(sds_id,cal,cal_error,off,off_err,num_type1)

      nrank = rank
      do i=1,rank
         edge(i)   = dim_sizes(i)
         dims(i)   = dim_sizes(i)
         start(i)  = 0
         stride(i) = 1
      enddo
      
      if (nrec .ne. 0) then
         edge(rank) = nrec
         dims(rank) = nrec
         start(rank) = irec-1
      endif
     
      iprod = 1
      do i=1,rank
         iprod = iprod*dims(i)
      enddo
      
      ntotal = iprod
      
      if (num_type .ne. 5) then
         
         retn = sfrdata(sds_id,start,stride,edge,buffer)
         do i=1,ntotal
            real_array(i) = cal*buffer(i) + off
            if (num_type .eq. 23 .and. real_array(i) .lt. 0.) 
     $           real_array(i) = real_array(i) + float(65536)*cal
         enddo
         
      else
         
         retn = sfrdata(sds_id,start,stride,edge,real_array)
         do i=1,ntotal
            real_array(i) = cal*real_array(i) + off
         enddo
         
      endif
      
      retn = sfendacc(sds_id)
      return
      end
      

      subroutine extract_int_sds(sd_id,in_name,irec,nrec,nrank,dims,buffer)
      
      ! extract the contents of an SDS from an HDF file
      ! 
      ! 5-12-98  rsd
      !
      parameter (MAX_VAR_DIMS = 3, MAX_NC_NAME=256)
      integer index,nrank,dims(MAX_VAR_DIMS)
c      parameter (MAX_BUF_SIZE = 500000)
      parameter (MAX_BUF_SIZE = 10000000)
      
      integer irec,nrec
      integer sds_id, sd_id, retn
      integer sfselect, sfendacc,sfginfo,sfgcal
      integer sfn2index
      character*(*) in_name
      integer start(MAX_VAR_DIMS),edge(MAX_VAR_DIMS),stride(MAX_VAR_DIMS)
      integer dim_sizes(MAX_VAR_DIMS)
      character name*(MAX_NC_NAME)
      integer rank,attributes,num_type
      double precision cal,cal_error,off,off_err

      integer*2 buffer(MAX_BUF_SIZE)

      index = sfn2index(sd_id, in_name)
      sds_id = sfselect(sd_id, index)   ! get SDS id
      retn = sfginfo(sds_id,name,rank,dim_sizes,num_type,attributes)
      retn = sfgcal(sds_id,cal,cal_error,off,off_err,num_type1)

      nrank = rank
      do i=1,rank
         edge(i)   = dim_sizes(i)
         dims(i)   = dim_sizes(i)
         start(i)  = 0
         stride(i) = 1
      enddo
      
      if (nrec .ne. 0) then
         edge(rank) = nrec
         dims(rank) = nrec
         start(rank) = irec-1
      endif
      
      retn = sfrdata(sds_id,start,stride,edge,buffer)
      iprod = 1
      do i=1,rank
         iprod = iprod*dims(i)
      enddo
      ntotal = iprod
      
      retn = sfendacc(sds_id)
      return

      end


      subroutine read_attrib_byname(sd_id,in_attr_name,
     $     num_type,n_values,fvalues)
      !
      ! read the name and value(s) of a global attribute (metadata)
      !  referenced by its name
      ! 5-14-98 rsd
      
      parameter (MAX_NC_NAME=256)
      integer sd_id
      integer attr_index
      character*(*) in_attr_name
      character attr_name*(MAX_NC_NAME)
      character attr_data*512
      character*(MAX_NC_NAME) values(20)
      character*(*) fvalues(*)
      character cr
      integer n_values
      integer n,oldn
      integer count
      integer retn,sfgainfo,sfrattr
      integer sffattr
 
      attr_index = sffattr(sd_id,in_attr_name)     
      retn = sfgainfo(sd_id,attr_index,attr_name,num_type,count)
      retn = sfrattr(sd_id,attr_index,attr_data)
      
      cr = char(10)
      ival = 0
      oldn = 1
 5    continue
      n = index(attr_data(oldn:(count-1)),cr)
      if (n .eq. 0) then
         ival = ival + 1
         values(ival) = attr_data(oldn:(count-1))
         goto 99
      else
         ival = ival + 1
         values(ival) = attr_data(oldn:(oldn+n-2))
      endif
      oldn = n + oldn
      goto 5
      
 99   continue
      n_values = ival - 2
      do i=1,n_values
          fvalues(i) = values(i+2)
      enddo
      return
      end


      subroutine extract_sds2(sd_id,in_name,irec,nrec,nrank,dims,real_array,buffer)
      
      ! extract the contents of an SDS from an HDF file
      ! 
      ! 5-12-98  rsd
      !
      parameter (MAX_VAR_DIMS = 3, MAX_NC_NAME=256)
      integer index,nrank,dims(MAX_VAR_DIMS)
      
      integer irec,nrec
      integer sds_id, sd_id, retn
      integer sfselect, sfendacc,sfginfo,sfgcal
      integer sfn2index
      character*(*) in_name
      integer start(MAX_VAR_DIMS),edge(MAX_VAR_DIMS),stride(MAX_VAR_DIMS)
      integer dim_sizes(MAX_VAR_DIMS)
      character name*(MAX_NC_NAME)
      integer rank,attributes,num_type
      double precision cal,cal_error,off,off_err

      integer*2 buffer(*)
      real real_array(*)
      
      index = sfn2index(sd_id, in_name)
      if (index.lt.0) print *,'*** sfn2index error extract_sds ***'

      sds_id = sfselect(sd_id, index)   ! get SDS id
      if (sds_id.lt.0) print *,'*** sfselect error extract_sds ***'

      retn = sfginfo(sds_id,name,rank,dim_sizes,num_type,attributes)
      if (retn.lt.0) print *,'*** sfginfo error extract_sds ***'

      retn = sfgcal(sds_id,cal,cal_error,off,off_err,num_type1)
      if (retn.lt.0) print *,'*** sfgcal error extract_sds ***'

      nrank = rank
      do i=1,rank
         edge(i)   = dim_sizes(i)
         dims(i)   = dim_sizes(i)
         start(i)  = 0
         stride(i) = 1
      enddo
      
      if (nrec .ne. 0) then
         edge(rank) = nrec
         dims(rank) = nrec
         start(rank) = irec-1
      endif
     
      iprod = 1
      do i=1,rank
         iprod = iprod*dims(i)
      enddo
      
      ntotal = iprod
      
      if (num_type .ne. 5) then
         
         retn = sfrdata(sds_id,start,stride,edge,buffer)
         if (retn.lt.0) print *,'*** sfrdata error extract_sds ***',in_name
         do i=1,ntotal
            real_array(i) = cal*buffer(i) + off
            if (num_type .eq. 23 .and. real_array(i) .lt. 0.) 
     $           real_array(i) = real_array(i) + float(65536)*cal
         enddo
         
      else
         
         retn = sfrdata(sds_id,start,stride,edge,real_array)
         if (retn.lt.0) print *,'*** sfrdata error extract_sds ***',in_name
         do i=1,ntotal
            real_array(i) = cal*real_array(i) + off
         enddo
         
      endif
      
      retn = sfendacc(sds_id)
      if (retn.lt.0) print *,'*** sfendacc error extract_sds ***'

      return
      end
      

      subroutine read_timetags(filename,nframes,timetags)
      !
      ! read SWS frame timetag info contained in HDF vdata
      !
      
      character*80 filename
      character vdata_name*30
      character*60 fields
      integer hopen,vsfatch,hclose,vsfgid,vsfinq
      integer vsfrd,vsfdtch
      integer retn,file_access_mode,n_dds
      integer file_id,vdata_id,vdata_ref
      integer n_records,interlace,vdata_size
      character*21 timetags(*)  ! changed from 22 for new l1b DGL 10/5/98
      integer DFACC_RDONLY, FULL_INTERLACE
      parameter (DFACC_RDONLY = 1,
     $           FULL_INTERLACE = 0)
      
      file_access_mode = 1
      n_dds = 0
      
      file_id = hopen(filename,DFACC_RDONLY,0)
      call vfstart(file_id)
      vdata_ref = -1
      vdata_ref = vsfgid(file_id,vdata_ref)
      vdata_id = vsfatch(file_id,vdata_ref,'r')
      
      retn = vsfinq(vdata_id,n_records,interlace,fields,
     $     vdata_size,vdata_name)
      if (vdata_name(1:10) .ne. 'frame_time') then  ! changed by DGL
         nframes = 0
         print *,'*** Problem reading frame_time vdata in read_timetags ',vdata_name
c         return
      endif
c      print *,'n_records,interlace,vdata_size = ',
c     $     n_records,interlace,vdata_size
      retn = vsfrd(vdata_id,timetags,n_records,interlace)
      nframes = n_records
      
c      do i=1,nframes
c         print *,timetags(i)
c      enddo
      
      retn =  vsfdtch(vdata_id)
      call vfend(file_id)
      retn = hclose(file_id)
      return
      end


      subroutine close_hdf(sd_id)
      
      !
      ! close an HDF SD interface previously opened by open_hdf()
      !
      ! 5-13-98  rsd
      
      integer sd_id
      integer sfend
      integer retn
      
      retn = sfend(sd_id)
      return
      end
      


C
C **************************************************************************
C
	INTEGER FUNCTION LENGTH(LINE)
C
C	COMPUTES THE LENGTH (NON-BLANK) OF A STRING
C
	CHARACTER*(*) LINE
	N=LEN(LINE)
	DO I=N,1,-1
		IF (LINE(I:I).NE.' '.AND.LINE(I:I).NE.CHAR(0)) THEN
			LENGTH=I
			GOTO 10
		ENDIF
	END DO
	LENGTH=1
10	RETURN
	END
