c (c) copyright 2003 David G. Long, Brigham Young University program slice c c Example program to illustrate how to use BYU slice shape tables c c Written: D. Long 23 June 2003 c modified to compile with lf95 1 feb 2007, fix slice lon scaling 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 = 'slice') 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 = 1) parameter (MAX_PULSES = 100) parameter (MAX_SLICES = 8) character*21 TimeTags(11250) ! in vdata ! note change from 22 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) integer slice_qual_flag(MAX_PULSES,MAX_FRAMES) c 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) c real cell_azimuth(MAX_PULSES,MAX_FRAMES) c 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) real slice_lat(MAX_SLICES,MAX_PULSES,MAX_FRAMES) real slice_lon(MAX_SLICES,MAX_PULSES,MAX_FRAMES) real slice_sigma0(MAX_SLICES,MAX_PULSES,MAX_FRAMES) c real x_factor(MAX_SLICES,MAX_PULSES,MAX_FRAMES) real slice_azimuth(MAX_SLICES,MAX_PULSES,MAX_FRAMES) real slice_incidence(MAX_SLICES,MAX_PULSES,MAX_FRAMES) c real slice_snr(MAX_SLICES,MAX_PULSES,MAX_FRAMES) 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 real a(8),b(8),nlat(8),nlon(8) integer ibeam,iasc,ib real cx,cy 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 response_threshold 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 Slice 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 slice 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,file='slice.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) call extract_int_sds(sd_id,'slice_qual_flag',iframe,1,nrank,dims,slice_qual_flag) c call extract_sds2(sd_id,'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) c call extract_sds2(sd_id,'cell_azimuth',iframe,1,nrank,dims,cell_azimuth,buffer) c 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 call extract_sds2(sd_id,'snr',iframe,1,nrank,dims,snr,buffer) c call extract_sds2(sd_id,'kpc_a',iframe,1,nrank,dims,kpc_a,buffer) call extract_sds2(sd_id,'slice_lat',iframe,1,nrank,dims,slice_lat,buffer) call extract_sds2(sd_id,'slice_lon',iframe,1,nrank,dims,slice_lon,buffer) call extract_sds2(sd_id,'slice_sigma0',iframe,1,nrank,dims,slice_sigma0,buffer) c call extract_sds2(sd_id,'x_factor',iframe,1,nrank,dims,x_factor,buffer) call extract_sds2(sd_id,'slice_azimuth',iframe,1,nrank,dims,slice_azimuth,buffer) call extract_sds2(sd_id,'slice_incidence',iframe,1,nrank,dims,slice_incidence,buffer) c call extract_sds2(sd_id,'slice_snr',iframe,1,nrank,dims,slice_snr,buffer) call extract_sds2(sd_id,'slice_kpc_a',iframe,1,nrank,dims,slice_kpc_a,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,'| Slices ',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,''shape_in'',I1,''.tab'')') * pathshape(1:length(pathshape)),instmode_store+1 write(outershapefile,'(A,''shape_out'',I1,''.tab'')') * pathshape(1:length(pathshape)),instmode_store+1 c c Load slice shape data c print *,'Reading slice shape data files for mode ',instmode_store print *,innershapefile(1:length(innershapefile)) print *,outershapefile(1:length(outershapefile)) call read_shape_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.or. * iand(slice_qual_flag(i,1),slice_mask).ne.0) goto 3401 c c for demo, keep only measurements near equator c if (abs(cell_lat(i,1)).lt.1.0) then c c get beam (and polarization) of pulse c ibeam=mod(sigma0_mode_flag(i,1)/4,2) ! 0=inner, 1=outer c c for each slice c irec=irec+1 ! a count of pulses examined do 340 ii=1,8 ! use only inner 8 slices c c compute slice center c cen_long=mod(slice_lon(ii,i,1)/cosd(cell_lat(i,1))+cell_lon(i,1),360.0) cen_lat=slice_lat(ii,i,1)+cell_lat(i,1) 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 lat,lon of slice corners using table c call get_shape(time_orb,ant_az,ib,ii, * cen_lat,cen_long,nc,nlat,nlon,rev_orbit_period) ncor=nc c c write out to a file (standard unit 1) c write (1,1400) iframe,iasc,i-1,ii,ant_az,ibeam, * cell_lon(i,1),cell_lat(i,1), * (nlon(j),nlat(j),j=1,8) 1400 format(i5,1x,i1,1x,i2,1x,i1,1x,f6.2,1x,i1,1x,f7.2,1x,f7.2,16(1x,F7.2)) c jrec=jrec+1 ! a count of total records written c c get sign of sigma-0 value c c if (btest(slice_qual_flag(i,1),(i-1)*4+1).eq.1) then if (btest(slice_qual_flag(i,1),(i-1)*4+1)) then isgn=-1 else isgn=1 endif sigma=slice_sigma0(ii,i,1) thetai=slice_incidence(ii,i,1) azang=slice_azimuth(ii,i,1) c 340 continue ! ii loop: slices 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 Slices:',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 slice shape routines c c ****************************************************************** c subroutine read_shape_tables(f1,f2) character*(*) f1,f2 c c read QuikScat shape tables into common blocks c c DGL 9/26/98 c c inputs c f1: file name for inner beam shape table file c f2: file name for outer beam shape table file c common /shape_table/norbs,nazs,nss,orbtimes(36),az_angs(8,36,36,2), $ corners(16,8,36,36,2) c nss=8 ib=1 open(unit=39,file=f1,status='old',form='formatted') read(39,*) norbs,nazs read(39,*) orbtimes do iorb=1,norbs do iaz=1,nazs read(39,*) (az_angs(i,iaz,iorb,ib),i=1,8) end do end do do iorb=1,norbs do iaz=1,nazs read(39,*) ((corners(i,islice,iaz,iorb,ib),i=1,16),islice=1,8) end do end do close(39) c ib=2 open(unit=39,file=f2,status='old',form='formatted') read(39,*) norbs,nazs read(39,*) orbtimes do iorb=1,norbs do iaz=1,nazs read(39,*) (az_angs(i,iaz,iorb,ib),i=1,8) end do end do do iorb=1,norbs do iaz=1,nazs read(39,*) ((corners(i,islice,iaz,iorb,ib),i=1,16),islice=1,8) end do end do close(39) c return end c c c subroutine get_shape(time,az,ib,islice,cen_lat,cen_lon, $ nc,cor_lat,cor_lon,orb_period) dimension cor_lat(*),cor_lon(*) c c given the orbit time, azimuth angle, beam and slice number, c compute QuikScat shape using stored table c c DGL 9/26/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 islice: slice number (1-8 for inner 8 slices) c cen_lat: latitude of center of slice c cen_lon: longitude of center of slice c orbit_period: orbit period (in sec) c c outputs: c nc: number of corners (8) c cor_lat: latitude of corners c cor_lon: longitude of corners c common /shape_table/norbs,nazs,nss,orbtimes(36),az_angs(8,36,36,2), $ corners(16,8,36,36,2) c double precision x(8),y(8) double precision aearth,FLAT,temp,r,rsq data aearth,FLAT/6378.1363,3.3528131778969144e-3/ c data table_orbit_period/6061.0/ c rinterp(xx,yy,r1,r2,r3,r4)=(r2-r1)*yy+r1+((r4-r3)*yy+r3-(r2-r1)*yy-r1)*xx c nc=8 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_shape',time,t,orb_period,norbs print *,orbtimes call exit(-1) io1=1 io2=2 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)*0.1 ia1=ifix(a)+1 ia2=ia1+1 if (ia2.gt.nazs) ia2=ia2-nazs dela=a-ia1+1 c c bilinearly interpolate the azimuth angle c cx=rinterp(dela,delo,cosd(az_angs(islice,ia1,io1,ib)), $ cosd(az_angs(islice,ia1,io2,ib)), $ cosd(az_angs(islice,ia2,io1,ib)), $ cosd(az_angs(islice,ia2,io2,ib))) sx=rinterp(dela,delo,sind(az_angs(islice,ia1,io1,ib)), $ sind(az_angs(islice,ia1,io2,ib)), $ sind(az_angs(islice,ia2,io1,ib)), $ sind(az_angs(islice,ia2,io2,ib))) c r=sqrt(cx*cx+sx*sx) cx=cx/r cy=cy/r c c for each slice corner, interpolate and rotate c do i=1,8 x(i)=rinterp(dela,delo,corners(i,islice,ia1,io1,ib), $ corners(i,islice,ia1,io2,ib), $ corners(i,islice,ia2,io1,ib), $ corners(i,islice,ia2,io2,ib)) y(i)=rinterp(dela,delo,corners(i+8,islice,ia1,io1,ib), $ corners(i+8,islice,ia1,io2,ib), $ corners(i+8,islice,ia2,io1,ib), $ corners(i+8,islice,ia2,io2,ib)) temp= cx*x(i)+sx*y(i) y(i)=-sx*x(i)+cx*y(i) x(i)= temp end do c c convert x,y corners to lat,lon c r=(1.0-FLAT*sind(cen_lat)**2)*aearth rsq=r*cosd(cen_lat) do i=1,8 temp=x(i)/rsq if (temp.gt.1.0) temp=1.0 if (temp.lt.-1.0) temp=-1.0 temp=asind(temp) cor_lat(i)=asind((y(i)-(1.0-cosd(temp))* * sind(cen_lat)*rsq)/r)+cen_lat cor_lon(i)=temp+cen_lon 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