c============= sample program for reading L1B radiance ===============
c
c Compile:
c  on IRIX
c  f77 -64 -bytereclen -C read_GLI_L1B_sample.f \
c      -o  read_GLI_L1B_sample \
c          -I/"your HDF directory"/include \
c          -L/"your HDF directory"/lib \
c          -lmfhdf -ldf -lz
c
c  on Solaris
c  f77 -64 -C  read_GLI_L1B_sample.f \
c      -o  read_GLI_L1B_sample \
c          -I/"your HDF directory"/include \
c          -L/"your HDF directory"/HDF4.1r4/lib \
c          -lmfhdf -ldf -ljpeg -lz -lnsl
c  
c==============================// by NASDA EORC H.Murakami, 2003/04/20
c
c--------------- chage this line for 250m data ----------
      parameter(nscan=420,npixel=1354,ndet=12)
c      parameter(nscan=145,npixel=5104,ndet=48)
c--------------------------------------------------------
c----- default input parameter size -----
      parameter(nspixel=npixel/ndet+2,nsdet=2)
      parameter(nsline=nscan*nsdet,nline=nscan*ndet)
c----- read valuables -------
      character     infile*120
      character     sds*20
      real*4        lon(nspixel*nscan*nsdet)
      real*4        lat(nspixel*nscan*nsdet)
      real*4        gsys(27),gsys_2km(2)
      integer*2     l2datai2(npixel*nline)
      integer*4     l2datai4(nscan*2)
      real*8        l2datar8(nspixel*nsline)
      integer*2     odet(nsdet),pxl(nspixel)
      integer*2     gch_v(19),gch_s(6),gch_m(7),gch_q(10)
      real*4        wt(2,2)
c----- global attributes ------
      character     Product_Name*100
      character*5   Data_Sub_type
      integer*4     Pixels_per_Scan_Line
      integer*4     Pixels_per_Scan_Line_for_SWIR_2km
      integer*4     Number_of_Scan_Lines
      integer*4     Lines_per_Scan
      integer*4     Lines_per_Scan_for_SWIR_2km
c----- general HDF arrays -----
      integer*4     sfstart,sffattr,ga_id
      integer*4     sfgainfo,sfrattr,sfn2index
      integer*4     sfselect,sfginfo,sfrdata,sfendacc,sfend
      integer       retn, count
c      integer*4     attr_index,sffinfo
      character     attr_name*40
      integer*4     sd_id,data_type,sds_index
      character     sds_name*20
      integer*4     sds_id,nattrs
      integer*4     rank,dims(6)
      integer*4     start(6),stride(6)
c      integer*4     edges(6)
      character*7   para_type(25)
c----------------------------------------------------------
c
      para_type( 4)='char8  '
      para_type( 5)='float32'
      para_type( 6)='float64'
      para_type(20)='int8   '
      para_type(21)='uint8  '
      para_type(22)='int16  '
      para_type(23)='uint16 '
      para_type(24)='int32  '
      para_type(25)='uint32 '
c
      data gch_v/1,2,3,4,6,8,9,11,13,14,15,16,17,18,19,20,21,22,23/
      data gch_s/1,2,3,4,5,6/
      data gch_m/1,2,3,4,5,6,7/
      data gch_q/1,2,3,4,0,0,0,0,5,6/
c
      sds='                    '
      infile( 1: 40)='                                        '
      infile(46: 80)='                                        '
      infile(91:120)='                                        '
c
 4000 format(a4,i4,1x,a20,a7,6i6)
c     ------------------------------------------------------
c
      write(6,'(a23)') '>>Write Input File Name'
      read(5,'(a)') infile
c
      i1=0
      i2=0
      do i=1,256
       if(infile(i:i+3).eq.'A2GL') i1=i
       if((i2.lt.1).and.(infile(i:i).eq.' ')) i2=i-1
      enddo
c
      sd_id=sfstart(infile,1)
c
      if(sd_id.lt.0) then
       write(6,*) 'No data=',infile(i1:i2)
       goto 777
      else
       write(6,*) 'HDF_open=',infile(i1:i2)
      endif
c
c     --------------- read global attributes ---------------
c
      ga_id=sffattr(sd_id,'Product Name')
      if(ga_id.gt.-1) then
       retn=sfgainfo(sd_id,ga_id,attr_name,data_type,count)
       retn=sfrattr(sd_id,ga_id,Product_Name)
       write(6,*) attr_name,' = ',Product_Name(1:count)
      endif
c
      ga_id=sffattr(sd_id,'Data Sub-type')
      if(ga_id.gt.-1) then
       retn=sfgainfo(sd_id,ga_id,attr_name,data_type,count)
       retn=sfrattr(sd_id,ga_id,Data_Sub_type)
      else
       Data_Sub_type='QUAT'
      endif
      write(6,*) attr_name,' = ',Data_Sub_type
c
      ga_id=sffattr(sd_id,'Pixels per Scan Line')
      if(ga_id.gt.-1) then
       retn=sfgainfo(sd_id,ga_id,attr_name,data_type,count)
       retn=sfrattr(sd_id,ga_id,Pixels_per_Scan_Line)
       write(6,*) attr_name,' = ',Pixels_per_Scan_Line
       npixels=Pixels_per_Scan_Line
      endif
c
      ga_id=sffattr(sd_id,'Number of Scan Lines')
      if(ga_id.gt.-1) then
       retn=sfgainfo(sd_id,ga_id,attr_name,data_type,count)
       retn=sfrattr(sd_id,ga_id,Number_of_Scan_Lines)
       write(6,*) attr_name,' = ',Number_of_Scan_Lines
       nscans=Number_of_Scan_Lines
      endif
c
      ga_id=sffattr(sd_id,'Lines per Scan')
      if(ga_id.gt.-1) then
       retn=sfgainfo(sd_id,ga_id,attr_name,data_type,count)
       retn=sfrattr(sd_id,ga_id,Lines_per_Scan)
       write(6,*) attr_name,' = ',Lines_per_Scan
       ndets=Lines_per_Scan
       nlines=nscans*ndets
      endif
c
      ga_id=sffattr(sd_id,'Pixels per Scan Line for SWIR 2km')
      if(ga_id.gt.-1) then
       retn=sfgainfo(sd_id,ga_id,attr_name,data_type,count)
       retn=sfrattr(sd_id,ga_id,Pixels_per_Scan_Line_for_SWIR_2km)
       write(6,*) attr_name,' = ',Pixels_per_Scan_Line_for_SWIR_2km
       npixels_2k=Pixels_per_Scan_Line_for_SWIR_2km
      else
       nlines_2k=0
      endif
c
      ga_id=sffattr(sd_id,'Lines per Scan for SWIR 2km')
      if(ga_id.gt.-1) then
       retn=sfgainfo(sd_id,ga_id,attr_name,data_type,count)
       retn=sfrattr(sd_id,ga_id,Lines_per_Scan_for_SWIR_2km)
       write(6,*) attr_name,' = ',Lines_per_Scan_for_SWIR_2km
       ndets_2k=Lines_per_Scan_for_SWIR_2km
       nlines_2k=nscans*ndets_2k
      else
       nlines_2k=0
      endif
c     ----------------------------------------------------------
c
c     ------------------- read SDS ----------------------
      do ir=1,6
       start(ir)=0
       stride(ir)=1
      enddo
c     ---------------------------------------------------------------
      sds_index=sfn2index(sd_id,'l1b_pos_samp')
      if(sds_index.gt.-1) then
       sds_id=sfselect(sd_id,sds_index)
       retn=sfginfo(sds_id,sds_name,rank,dims,data_type,nattrs)
       retn=sfrdata(sds_id,start,stride,dims,l2datai4)
       do n=1,dims(1)
        pxl(n)=l2datai4(n)
       enddo
c       write(6,*) (pxl(k),k=1,2)
       write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
       retn=sfendacc(sds_id)
      endif
c
c     -----------------------------------
      sds_index=sfn2index(sd_id,'l1b_pos_line')
      if(sds_index.gt.-1) then
       sds_id=sfselect(sd_id,sds_index)
       retn=sfginfo(sds_id,sds_name,rank,dims,data_type,nattrs)
       retn=sfrdata(sds_id,start,stride,dims,l2datai4)
       odet(1)=l2datai4(1)
       odet(2)=l2datai4(2)
c       write(6,*) (odet(k),k=1,2)
       write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
       retn=sfendacc(sds_id)
      endif
c
c     ------------------------------------------------------------
      sds_index=sfn2index(sd_id,'l1b_blk_lon')
      if(sds_index.gt.-1) then
       sds_id=sfselect(sd_id,sds_index)
       retn=sfginfo(sds_id,sds_name,rank,dims,data_type,nattrs)
       write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
       nspixels=dims(1)
       nslines=dims(2)
       retn=sfrdata(sds_id,start,stride,dims,l2datar8)
       do nm=1,dims(1)*dims(2)
        lon(nm)=float(l2datar8(nm))
       enddo
       retn=sfendacc(sds_id)
      endif
c
c     ------------------------------------------------------------
      sds_index=sfn2index(sd_id,'l1b_blk_lat')
      if(sds_index.gt.-1) then
       sds_id=sfselect(sd_id,sds_index)
       retn=sfginfo(sds_id,sds_name,rank,dims,data_type,nattrs)
       write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
       retn=sfrdata(sds_id,start,stride,dims,l2datar8)
       do nm=1,dims(1)*dims(2)
        lat(nm)=float(l2datar8(nm))
       enddo
       retn=sfendacc(sds_id)
      endif
c
c     ------------------------------------------------------------
      sds_index=sfn2index(sd_id,'gsys')
      if(sds_index.lt.0) sds_index=sfn2index(sd_id,'c1')
      if(sds_index.gt.-1) then
       sds_id=sfselect(sd_id,sds_index)
       retn=sfginfo(sds_id,sds_name,rank,dims,data_type,nattrs)
       retn=sfrdata(sds_id,start,stride,dims,gsys)
       retn=sfendacc(sds_id)
       nbl=dims(1)
      else
       nbl=0
      endif
c
c     ------------------------------------------------------------
      sds_index=sfn2index(sd_id,'gsys_2km')
      if(sds_index.gt.-1) then
       retn=sfginfo(sds_id,sds_name,rank,dims,data_type,nattrs)
       retn=sfrdata(sds_id,start,stride,dims,gsys_2km)
       retn=sfendacc(sds_id)
       gsys(nbl+1)=gsys_2km(1)
       gsys(nbl+2)=gsys_2km(2)
       nbl=nbl+2
      endif
c     ------------------------------------------------------------
c
c     ---------------- Specify SD data for reading -------------
      if(Data_Sub_type(1:1).eq.'V') then
       write(6,'(a44)') '>>Specify SD_name (e.g.,l1b_ch1_data;  1~19)'
      elseif(Data_Sub_type(1:1).eq.'S') then
       write(6,'(a44)') '>>Specify SD_name (e.g.,l1b_ch24_data;24~29)'
      elseif(Data_Sub_type(1:1).eq.'M') then
       write(6,'(a44)') '>>Specify SD_name (e.g.,l1b_ch30_data;30~36)'
      elseif(Data_Sub_type(1:1).eq.'Q') then
       write(6,'(a44)') '>>Specify SD_name (e.g.,l1b_ch20_data;20~29)'
      else
       write(6,'(a10,a4)') 'Invalid = ',Data_Sub_type(1:4)
       goto 999
      endif
      read(5,'(a)') sds
c
      if(sds(1:6).eq.'l1b_ch') then
       if(sds(8:8).eq.'_') then
        read(sds(7:7),'(i1)') nbc
       else
        read(sds(7:8),'(i2)') nbc
       endif
      else
       nbc=0
      endif
c
      if((Data_Sub_type(1:1).eq.'S').and.
     &   ((nbc.eq.28).or.(nbc.eq.29))) then
       sds(14:17)='_2km'
      endif
c
      if(sds(15:17).eq.'2km') then
       ndets=ndets_2k
       nlines=nlines_2k
       npixels=npixels_2k
      endif
c
c     ------------------------------------- Search band index
      ipwb=0
      if(Data_Sub_type(1:1).eq.'V') then
       nb=nbc
       nbg=gch_v(nbc)
       slope1=gsys(nbg)
       slope2=gsys(nbg+1)
       if((nb.eq.4).or.(nb.eq.5).or.(nb.eq.7).or.(nb.eq.8)) 
     &  ipwb=1
      elseif(Data_Sub_type(1:1).eq.'S') then
       nb=nbc-23
       if(nb.lt.1) then
        write(6,*) 'No data in ',Data_Sub_type
        goto 999
       endif
       nbg=gch_s(nbc-23)
       slope1=gsys(nbg)
       slope2=gsys(nbg)
      elseif(Data_Sub_type(1:1).eq.'M') then
       nb=nbc-29
       if(nb.lt.1) then
        write(6,*) 'No data in ',Data_Sub_type
        goto 999
       endif
       nbg=gch_m(nbc-29)
       slope1=gsys(nbg)
       slope2=gsys(nbg)
      elseif(Data_Sub_type(1:1).eq.'Q') then
       if(nbc.lt.20) then
        write(6,*) 'No data in ',Data_Sub_type
        goto 999
       endif
       nb=gch_q(nbc-19)
       nbg=gch_q(nbc-19)
       slope1=gsys(nbg)
       slope2=gsys(nbg)
      endif
c     ---------------------------------------------------
c
      sds_index=sfn2index(sd_id,sds)
      if(sds_index.lt.0) then
       write(6,*) 'No SDS=',sds
       goto 999
      endif
      sds_id=sfselect(sd_id,sds_index)
      retn=sfginfo(sds_id,sds_name,rank,dims,data_type,nattrs)
      retn=sfrdata(sds_id,start,stride,dims,l2datai2)
      retn=sfendacc(sds_id)
      write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
c
      write(6,'(a33)') '>>Specify pixel,line for reading'
      read(5,*) nc0,mc0
c
      if(sds(15:17).eq.'2km') then
       nc=(nc0-1)*2+1
       mc=(mc0-1)*2+1
      else
       nc=nc0
       mc=mc0
      endif
c
      if((mc.gt.0).and.(nc.gt.0).and.
     &   (mc.le.nlines).and.(nc.le.npixels)) then
c
c       -------------- calculate lat lon -------------
       rng=1.0*(nc-pxl(1))/(pxl(2)-pxl(1))
       ng=int(rng)
       rng2=rng-ng
       rng1=1.-rng2
       rmg=1.0*(mc-odet(1))/(odet(2)-odet(1))
       mg=int(rmg)
       rmg2=rmg-mg
       rmg1=1.-rmg2
       wt(1,1)=rng1*rmg1
       wt(2,1)=rng2*rmg1
       wt(1,2)=rng1*rmg2
       wt(2,2)=rng2*rmg2
       wts=0.
       av_lon=0.
       av_lat=0.
       do mm=1,2
        do nn=1,2
         wts=wts+wt(nn,mm)
         nmg=(mg+mm-1)*nspixels+ng+nn
         av_lon=av_lon+lon(nmg)*wt(nn,mm)
         av_lat=av_lat+lat(nmg)*wt(nn,mm)
        enddo
       enddo
       av_lon=av_lon/wts
       av_lat=av_lat/wts
c      ----------------------------------------------
c
       nm=(mc-1)*npixels+nc
       irad=l2datai2(nm)
       if(ipwb.gt.0) then
        if(irad.lt.0) then
         irad=mod(irad+65536,8192)
         if(irad.gt.4095) then
          slope=slope1
          rad=(irad-4096)*slope
         else
          slope=slope2
          rad=-999
         endif
        else
         if(irad.gt.16383) irad=irad-16384
         if(irad.gt.4095) then
          slope=slope1
          rad=(irad-4096)*slope
         else
          slope=slope2
          rad=irad*slope
         endif
        endif
       else
        if(irad.lt.0) then
         slope=slope1
         rad=-999
        else
         slope=slope1
         if(irad.gt.16383) irad=irad-16384
         rad=irad*slope
        endif
       endif
c
       write(6,'(a12,f8.3,a1,f8.3)') 'Lon, Lat  = ',av_lon,',',av_lat
       write(6,'(a12,i6,a1,i6)')     'pixel,line= ',nc0,',',mc0
       write(6,'(a12,i6,a1,f10.7)')
     &                        'DN, slope = ',l2datai2(nm),',',slope
       write(6,'(a12,f10.5,a15)')
     &                        'Radiance  = ',rad,' [W/m^2/str/um]'
c
      else
       write(6,'(a18,i5,a1,i5)') 'Out of range n,m =',nc,',',mc
       goto 999
      endif
c     -----------------------------------------------------------------

c     -------------------------- Close HDF file
  999  retn=sfend(sd_id)
c     ------------------------------------------
c
 777  continue
      stop
      end
c
