c=====================================================================
c   Read GLI L1B/L2 scene data and mapping to binary data
c
c   Compile&Run:
c     run_glismap_v## $GLIFileName $SdsName 
c                 ($lon_min $lon_max $lat_min $lat_max $resolution)
c      * $GLIFileName: GLI-L1B/Ocean-L2 file name with path
c      * $SdsName: GLI-L1B/Ocean-L2 SDS name; chla, sst..
c      * $lon_min...: output lat/lon area specification (option)
c      * if no area specification; mapping full-area
c      * $resolution: output lat/lon grid interval by degree (option)
c      * if no area specification, dd0 is used
c
c   Examples of "$SdsName":
c             chla, sst, cdom, ss, nLw_443,...
c             l1b_all, l1b_ch4_data, l1b_ch28_data_2km,
c             nLw, l1b_rgb  (<- 1byte x 3ch data for RGB browse)
c
c   Output data format
c      pixel & line size: you can find from output file name;
c                        "A2GL1**_pixel_line"
c      one-ch data: header+ 2-byte binary x pixel x line
c      RGB data   : header+ 1-byte binary x pixel x line
c      header: its size corresponds to the one line (2/1-byte x pixel)
c              weitten by ascii, see "output header" in this code.
c              including following information
c              "pixel,line,lon_min,lat_max,resolution,slope,offset"
c      parameter value [each unit]: DN*slope
c
c   When you want process 1-km data, change two lines below,
c                              "----- for 1km data -----".
c   When 250m data, change the two lines and comment out 1-km ones
c                              "----- for 250m data -----".
c   * if your computer has not enough memory, please tune these 
c     dimension size or "unlimit stacksize" of your environment.
c
c   NOTES:
c     Please refer GLI Ocean or CAL HP about GLI data: 
c       http://suzaku.eorc.jaxa.jp/GLI/ocean/index.html
c       http://suzaku.eorc.jaxa.jp/GLI/cal/index.html
c     The HDF library must be installed before this process. 
c     The HDF library can be obtained from below.
c       http://hdf.ncsa.uiuc.edu
c
c   Copyright:
c     2003, JAXA Earth Observation Research and Application Center
c       You can use this only for Research or application study, 
c       not for commercial application.
c
c     last modified by JAXA/EORC H.Murakami, 2003/12/15
c
c=====================================================================
c     ----- for 1km data -----
c      parameter(nscan=420,npixel=1354,ndet=12)   ! 1-km NRT
      parameter(nscan=150,npixel=1354,ndet=12)   ! 1-km STD    *
      parameter(nl=5271,ml=4800,dd0=0.01)   ! Full resolution  *
c     ----- for 250m data -----
c      parameter(nscan=145,npixel=5104,ndet=48)
c      parameter(nl=9904,ml=7000,dd0=0.0025)   ! Full resolution  *
c      parameter(nl=5000,ml=4000,dd0=0.0025)   ! Full resolution  *
c      parameter(nl=6000,ml=4200,dd0=0.005)    ! middle resolution
c      parameter(nl=9608,ml=2400,dd0=0.01)   ! 1-km resolution
c      parameter(nl=5000,ml=4200,dd0=0.04)     ! for browse
c----- default input parameter size -----
      parameter(nspixel=npixel/ndet+2,nsdet=2)
      parameter(nsline=nscan*nsdet,nline=nscan*ndet)
c
c----- global attributes ------
      character     Product_Name*120
      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----------------------------------------
      real*4        lon0(nspixel*nscan*nsdet)
      real*4        lat0(nspixel*nscan*nsdet)
      real*4        lon1(nspixel*nscan*nsdet)
      real*4        lat1(nspixel*nscan*nsdet)
      real*4        gsys(40),gsys_2km(2)
      integer*2     vnirb(19),glivn
      integer*4     map(nl,ml)
      real*4        slope,slope2,offset
      integer*2     iflag(npixel,nline)
      integer*2     i2_out(nl,ml)
      character     i1_out(nl,ml,3)*1
      integer*2     odet(nsdet),pxl(npixel)    ! npixel for mud
c
c----- SDS Array -----
      integer*2     l2datai2(npixel*nline)
c      integer*4     l2datai4(npixel*nline)
      integer*4     l2datai4(1354*4000)   ! only used for 1km
      real*4        l2datar4(nspixel*nsline)
      real*8        l2datar8(nspixel*nsline)
c----- general HDF arrays -----
      character     infile1*120,infile*50
      integer*4     sfsdtstr,sfgdtstr
      integer*4     hopen,vfstart,vfatch,vfdtch,vfend,hclose
      integer*4     vfsnam,vfscls,sfid2ref,vfadtr,sds_ref
      integer*4     sfstart,sffinfo,sffattr,ga_id
      integer*4     sfgainfo,sfrattr,sfn2index
      integer*4     sfselect,sfginfo,sfrdata,sfendacc,sfend
      integer*4     sfcreate,sfwdata
      integer*4     sfsnatt,sfscatt
      integer       retn, count
      character     attr_name*40
      integer*4     sd_id,data_type,sds_index
      character     sds_name*20
      integer*4     sds_id,data_type,nattrs
      integer*4     rank,dims(6)
      integer*4     start(6),stride(6),edges(6)
      character*7   para_type(25)
      integer*4     sffattr,sfrattr,sfgainfo
      integer*4     attr_index
      integer*4     sd_id2
c----------------------------------------------------------
      character*1   bafac1(50)
      integer*2     bafai2(50)
      integer*4     bafai4(50)
      real*4        bafar4(50)
      real*8        bafar8(50)
c
      character     outfile*55,head*125
      real*4        lat_max,lat_min,lon_max,lon_min
      character     sds0*17,sds1*17
      character     sds_rgb(3)*13
      real*4        rng_rgb(3)
      character     a0*1,a9*1,a8*1,a7*1,a6*1
      real*4        gains(2,3)
      real*4        d_max(3),d_min(3),d_avg(3),d_cld(3)
      character     pal(256,3)*1
      real*4        f_max(3),f_min(3),f_avg(3)
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 vnirb/1,2,3,4,6,8,9,11,13,14,15,16,17,18,19,20,21,22,23/
c
      data sds0/'                 '/
      a0=char(0)
c      a7=char(150)                         ! for gray background
      a7=char(1)                           ! for black background
      a6=char(0)                           ! for black background
      a8=char(254)
      a9=char(255)
      infile1( 1: 40)='                                        '
      infile1(46: 80)='                                        '
      infile1(91:120)='                                        '
      Product_Name( 1: 40)='                                        '
      Product_Name(41: 80)='                                        '
      Product_Name(81:100)='                                        '
      rloc4=-999.
      rloc1=-999.
      resol=-999.
c     ------------------ RGB Scaling parameters ---------------
      data d_max/130.,140.,155./
      data d_min/ 24., 45., 80./
c      data d_avg/ 36., 50.,70./    ! for ocean dark target
      data d_avg/ 47., 65.,91./    ! for ocean middle
c      data d_avg/ 60., 80.,110./   ! for ocean bright
c      data d_avg/ 90.,120.,165./   ! for land  <--generally recommneded
      data d_cld/432.,488.,536./   ! for land  <--generally recommneded
c      data d_cld/540.,610.,670./   ! for extremely bright target
c
c     ===================================
c      data f_max/0.,0.,0./
c      data f_min/0.,0.,0./
c      data f_avg/0.,0.,0./
      data f_max/432.   ,488.   ,536.   /
      data f_min/ 14.950, 28.530, 45.095/
      data f_avg/ 90.   ,120.   ,165.   /
c     ===================================
c
c     =======
      idsp=1
c     =======
c
 4000 format(a4,i4,1x,a20,a7,6i6)
c
c     =================================== Input HDF file name
      open(10,file='glismap.name',status='old')
      read(10,'(a)') infile1
      read(10,'(a)') sds0
      read(10,*,end=111,err=111) rloc1
      read(10,*,end=111,err=111) rloc2
      read(10,*,end=111,err=111) rloc3
      read(10,*,end=111,err=111) rloc4
      read(10,*,end=111,err=111) resol
      close(10)
  111 continue
c     =======================================================
c
      if(resol.gt.0) then
       dd=resol
      elseif((rloc1.gt.0).and.(rloc4.eq.-999.)) then
       dd=rloc1
      else
       dd=dd0
      endif
c
      if(rloc4.eq.-999.) then
       rloc1=-999.
       rloc2=-999.
       rloc3=-999.
       rloc4=-999.
      endif
c
      sd_id=sfstart(infile1,1)
c
      i1=0
      i2=0
      do i=1,120
       if((infile1(i:i+3).eq.'A2GL').or.
     &    (infile1(i:i+2).eq.'MOD')) i1=i
       if((i2.lt.1).and.(infile1(i:i).eq.' ')) i2=i-1
      enddo
      infile=infile1(i1:i2)
c
      if(sd_id.lt.0) then
       write(6,*) 'No data=',infile1(1:i2)
       goto 777
      else
       write(6,*) 'HDF_open=',infile1(1: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)
       i3=0
       do i=1,120-30
        if(Product_Name(i:i+3).eq.'A2GL') i3=i
       enddo
       if(idsp.gt.0) write(6,*) attr_name,' = ',Product_Name(i3:i3+32)
      else
       i3=0
       if(idsp.gt.0) write(6,*) attr_name,' = ',Product_Name
      endif
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)
       if(idsp.gt.0) write(6,*) attr_name,' = ',Pixels_per_Scan_Line
       npixels=Pixels_per_Scan_Line
      else
       npixels=0
      endif
c
c     --------------------- dimension check npixels
      if(npixel.lt.npixels) then
       write(6,*) 'Err, npixel < npixels=',npixels
       goto 777
      endif
c     -----------------------------------------------------
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)
       if(idsp.gt.0)
     &  write(6,*) attr_name,' = ',Pixels_per_Scan_Line_for_SWIR_2km
       npixels_2k=Pixels_per_Scan_Line_for_SWIR_2km
      else
       npixels_2k=0
      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)
       if(idsp.gt.0) write(6,*) attr_name,' = ',Number_of_Scan_Lines
       nscans=Number_of_Scan_Lines
      else
       nscans=0
      endif
c
c     -------------------- dimension check nscans
      if(nscan.lt.nscans) then
       write(6,*) 'Err, nscan < nscans ',nscans
       goto 777
      endif
c     -----------------------------------------------------
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)
       if(idsp.gt.0) write(6,*) attr_name,' = ',Lines_per_Scan
       ndets=Lines_per_Scan
       nlines=nscans*ndets
      else
       nlines=0
      endif
c
      if(Product_Name(32:33).eq.'mu') then
       ndets=1
       nscans=nlines
      endif
c
c     --------------------  dimension check ndets
      if(nsdet.lt.nsdets) then
       write(6,*) 'Err, ndet < ndets ',ndets
       goto 777
      endif
c     --------------------  dimension check nlines
      if(nline.lt.nlines) then
       write(6,*) 'Err, nline < nlines ',nlines
       goto 777
      endif
c     -----------------------------------------------------
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)
       if(idsp.gt.0) 
     &  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
      if(sds0(15:17).eq.'2km') then
       ndets=ndets_2k
       nlines=nlines_2k
       npixels=npixels_2k
      endif
c
      ga_id=sffattr(sd_id,'MODIS_VNIR_CHANELS')
      if(ga_id.gt.-1) then
       modcnv=1
       write(6,*) 'MODIS CNV data'
      else
       modcnv=0
      endif
c
c     -------------------- read SDS data --------------------
      do ir=1,6
       start(ir)=0
       stride(ir)=1
      enddo
c
c     ----------------- flag ------------------
      sds_index=sfn2index(sd_id,'l2_flags')
      if(sds_index.lt.0) then
       sds_index=sfn2index(sd_id,'flag')
      endif
c
      if(sds0(1:6).eq.'RcRefl') sds_index=-1
c
      if(sds_index.gt.-1) then
       sds_id=sfselect(sd_id,sds_index)
       retn=sfginfo(sds_id,sds_name,rank,dims,data_type,nattrs)
       npixels=dims(1)
       nlines=dims(2)
       if(data_type.eq.23) then
        retn=sfrdata(sds_id,start,stride,dims,l2datai2)
        do m=1,nlines
        do n=1,npixels
         nm=(m-1)*npixels+n
         iflag(n,m)=mod(l2datai2(nm),2)
        enddo
        enddo
       elseif(data_type.eq.25) then
        retn=sfrdata(sds_id,start,stride,dims,l2datai4)
        do m=1,nlines
        do n=1,npixels
         nm=(m-1)*npixels+n
         iflag(n,m)=mod(l2datai4(nm)/32,2)
        enddo
        enddo
       endif
       write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
       retn=sfendacc(sds_id)
      else
       do m=1,nlines
        do n=1,npixels
         nm=(m-1)*npixels+n
         iflag(n,m)=0
        enddo
       enddo
      endif
c     --------------------------------------------------------------
c
c     ---------------------------------------------------------------
      pxl(1)=-9999
      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,dims(1))
       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)
       if(l2datai4(1).gt.0) then
        odet(1)=l2datai4(1)
        odet(2)=l2datai4(2)
       else
        odet(1)=l2datai4(1)+1
        odet(2)=l2datai4(2)+1
       endif
       write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
       retn=sfendacc(sds_id)
      endif
c
c     ------------------------------------------------------------
      do nm=1,nspixel*nscan*nsdet
       lon0(nm)=-999
      enddo
      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)
        lon0(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)
        lat0(nm)=float(l2datar8(nm))
       enddo
       retn=sfendacc(sds_id)
      endif
c
c     ------------------------------------------------------------
      if_scanlat=0
      sds_index=sfn2index(sd_id,'lat')
      if(sds_index.lt.0) then
       sds_index=sfn2index(sd_id,'Latitude')
      endif
      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,l2datar4)
       retn=sfendacc(sds_id)
       write(6,*) 'Read scan-point geolocations' 
       if(idsp.gt.0)
     & write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
       if(rank.eq.3) then
        nspixels= dims(1)
        nsdets=   dims(3)
        nscans=   dims(2)
        nslines=  dims(2)*dims(3)
        if_scanlat=1
       elseif(rank.eq.2) then
        nspixels= dims(1)
        nslines=  dims(2)
        nsdets=   2
        nscans=   dims(2)/nsdets
        if_scanlat=1
       elseif(rank.gt.3) then
        nbands=dims(4)
        if_scanlat=2
       endif
       do nm=1,nspixels*nslines
        lat0(nm)=l2datar4(nm)
       enddo
       if(pxl(1).eq.-9999) then
        do np=1,nspixels
         pxl(np)=(np-1)*ndet
        enddo
       endif
       if(nsline.lt.nslines) then
        write(6,*) 'Err, nsline < nslines=',nslines
        goto 777
       endif
      endif
c
      sds_index=sfn2index(sd_id,'lon')
      if(sds_index.lt.0) then
       sds_index=sfn2index(sd_id,'Longitude')
      endif
      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,l2datar4)
       retn=sfendacc(sds_id)
       do nm=1,nspixels*nslines
        lon0(nm)=l2datar4(nm)
       enddo
       if(idsp.gt.0)
     & write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
      endif
c      write(6,*) nscans,nsdets,nspixels,nslines
c
c     -------------- Longitude/Latitude correction -------------
      if(if_scanlat.eq.0) then
       nsdets=2
       nslines=nscans*nsdets
       slon=lon0(nscans*nspixels/2+1)
       do nsc=1,nscans
        do nd=1,nsdets
         nsd=(nsc-1)*nsdets+nd
         do np=1,nspixels
          nsdp1=(nsc+nd-nsdets)*nspixels+np
          nsdp=(nsd-1)*nspixels+np
          if(lon0(nsdp1).gt.-541.) then
           if(slon.gt.-541.) then
            dif=lon0(nsdp1)-slon
            if(dif.gt.180.) then
             lon1(nsdp)=lon0(nsdp1)-360.
            elseif(dif.lt.-180.) then
             lon1(nsdp)=lon0(nsdp1)+360.
            else
             lon1(nsdp)=lon0(nsdp1)
            endif
           else
            lon1(nsdp)=lon0(nsdp1)
           endif
           lat1(nsdp)=lat0(nsdp1)
          else
           lon1(nsdp)=-999.
           lat1(nsdp)=-999.
          endif
         enddo
        enddo
       enddo
      else
c       write(6,*) nlines,nspixels
       slon=lon0(nslines*nspixels/2+1)
       do nsc=1,nscans
        do nd=1,nsdets
         nsd=(nsc-1)*nsdets+nd
         do np=1,nspixels
          nsdp=(nsd-1)*nspixels+np
          if(lon0(nsdp).gt.-541.) then
           if(slon.gt.-541.) then
            dif=lon0(nsdp)-slon
            if(dif.gt.180.) then
             lon1(nsdp)=lon0(nsdp)-360.
            elseif(dif.lt.-180.) then
             lon1(nsdp)=lon0(nsdp)+360.
            else
             lon1(nsdp)=lon0(nsdp)
            endif
           else
            lon1(nsdp)=lon0(nsdp)
           endif
           lat1(nsdp)=lat0(nsdp)
          else
           lon1(nsdp)=-999.
           lat1(nsdp)=-999.
          endif
         enddo
        enddo
       enddo
      endif
c     ----------------------------------------------------------
c
c     ------------------------------ min,max of the area
      if(rloc4.eq.-999.) then
       lat_max= -90.
       lat_min=  90.
       lon_max=-541.
       lon_min= 541.
c       write(6,*) nslines,nspixels,dd
       do m=1,nslines
        do n=1,nspixels
         nm=(m-1)*nspixels+n
         if((lat1(nm).gt.-89.82).and.(lat1(nm).lt.89.82)) then
          if(lat1(nm).gt.lat_max) then
           lat_max=lat1(nm)
          elseif(lat1(nm).lt.lat_min) then
           lat_min=lat1(nm)
          endif
         endif
         if(lon1(nm).gt.lon_max) then
          lon_max=lon1(nm)
         elseif((lon1(nm).gt.-541.).and.(lon1(nm).lt.lon_min)) then
          lon_min=lon1(nm)
         endif
        enddo
       enddo
       lat_min=nint(lat_min/dd)*dd
       lat_max=nint(lat_max/dd)*dd
       lon_min=nint(lon_min/dd)*dd
       lon_max=nint(lon_max/dd)*dd
      else
       lon_min=rloc1
       lon_max=rloc2
       lat_min=rloc3
       lat_max=rloc4
      endif
      mll=int((lat_max-lat_min)/dd+1)
      nll=int((lon_max-lon_min)/dd+1)
      if(mll.gt.ml) then
       write(6,'(a8,i6,a3,i6)') 'Warning ',mll,' > ',ml
       lat_max=lat_max-((mll-ml)/2)*dd
       lat_min=lat_max-(ml-1)*dd
       mll=ml
      endif
      if(nll.gt.nl) then
       write(6,'(a8,i6,a3,i6)') 'Warning ',nll,' > ',nl
       lon_min=lon_min+((nll-nl)/2)*dd
       lon_max=lon_min+(nl-1)*dd
       nll=nl
      endif
      write(6,'(a7,4f8.2)') 
     &             ' Area =',lon_min,lon_max,lat_min,lat_max
      write(6,'(a11,2i7)') ' [nll,mll]=',nll,mll
c
c     ===================================================
      CALL LINEARMAP
     &(npixel,nspixel,nsline,nspixels,nslines,npixels,nlines,
     & nscans,ndets,lat_max,lat_min,lon_max,lon_min,dd,
     & nsdet,nsdets,nl,ml,nll,mll,pxl,odet,lat1,lon1,map)
c     ===================================================
c
c     ------------------------------------------------------------
      if(sds0(1:3).eq.'sst') then
        slope =0.01
        offset=263.15
      elseif(sds0(1:3).eq.'chl') then
        slope =0.0015
        offset=0.
      elseif(sds0(1:4).eq.'cdom') then
        slope =0.0010
        offset=0.
      elseif(sds0(1:4).eq.'k490') then
        slope =0.0010
        offset=0.
      elseif(sds0(1:2).eq.'ss') then
        slope =0.0100
        offset=0.
      elseif(sds0(1:6).eq.'RcRefl') then
        slope =0.00002
        offset=0.
      elseif(sds0(1:3).eq.'tau') then
        slope =0.0001
        offset=0.
      elseif(sds0(1:3).eq.'par') then
        slope =0.002
        offset=0.
      elseif(sds0(1:3).eq.'alp') then
        slope =0.0001
        offset=-3.0
      elseif(infile(21:22).eq.'NL') then
        slope =0.001
        offset=0.
      else
        slope =1.
        offset=0.
      endif
c     ------------------------------------------------------
c
c     ------------------------------------ mapping data
      iotype=0
c     --------------- for RGB in NL products ---------------
      if(((sds0(1:4).ne.'nLw_').and.(sds0(1:3).eq.'nLw')).or.
     &   (sds0(1:6).eq.'RcRefl')) then
       iotype=1
       nbl=1
       if(sds0(1:6).eq.'RcRefl') then
        sds_rgb(1)='RcRefl_678'
        sds_rgb(2)='RcRefl_545'
        sds_rgb(3)='RcRefl_460'
        rng_rgb(1)=0.5
        rng_rgb(2)=0.5
        rng_rgb(3)=0.5
       elseif(sds0(1:3).eq.'nLw') then
c        sds_rgb(1)='nLw_666'      !    for MODIS CNV
        sds_rgb(1)='nLw_625'      !    for GLI
        sds_rgb(2)='nLw_545'
        sds_rgb(3)='nLw_460'
        rng_rgb(1)=2.2
        rng_rgb(2)=3.0
        rng_rgb(3)=3.3
       endif
c
       do k=1,3
        sds_index=sfn2index(sd_id,sds_rgb(k))
        if(sds_index.lt.0) then
         write(6,*) 'no SDS data=',sds0
         goto 777
        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)
c
        do m1=1,mll
         do n1=1,nll
          nm=map(n1,m1)
          m=(nm-1)/npixels+1
          n=nm-(m-1)*npixels
          if((m.gt.0)  .and.(n.gt.0).and.
     &       (m.le.nlines).and.(n.le.npixels)) then
           if(iflag(n,m).lt.1) then
            if(l2datai2(nm).ge.0) then
             if(l2datai2(nm).gt.offset/slope) then
              rdat=((slope*l2datai2(nm)+offset)/rng_rgb(k))**0.75*250
             else
              rdat=0.
             endif
            else
             rdat=((slope*(65536+l2datai2(nm))+offset)
     &                 /rng_rgb(k))**0.75*250
            endif
            if(rdat.gt.254.) then
             rdat=254.
            elseif(rdat.lt.0.) then
             rdat=0.
            endif
           else
            rdat=100.
           endif
          else
c           rdat=128.
           rdat=1.                  !   *** set Black ***
          endif
          i1_out(n1,m1,k)=char(int(rdat))
         enddo
        enddo
        write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &                (dims(ir),ir=1,rank)
       enddo
c
c     --------------------- for RGB in L1B --------------------
      elseif(sds0(5:7).eq.'rgb') then
       nbl=1
c
       if(sds0(1:3).ne.'l1b') then
        write(6,'(a30)') 'Input VNIR file (path/A2GL1..)'
        read(5,'(a)') infile1
        sd_id2=sfstart(infile1,1)
        if(sd_id2.gt.-1) then
         iotype=3
c
         ga_id=sffattr(sd_id2,'MODIS_VNIR_CHANELS')
         if(ga_id.gt.-1) then
          modcnv=1
          write(6,*) 'MODIS CNV data'
         else
          modcnv=0
         endif
c
        else
         write(6,*) 'No data= ',infile1
         sd_id2=sd_id
         iotype=2
        endif
       else
        sd_id2=sd_id
        iotype=2
       endif
c
       sds_index=sfn2index(sd_id2,'gsys')
       sds_id=sfselect(sd_id2,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)
c
       if(Product_Name(i3:i3+4).eq.'A2GL2') then
        write(6,*) 'making 250m RGB data...'
        sds_rgb(1)='l1b_ch22_data'
        gains(1,1)=gsys(3)
        gains(2,1)=gsys(3)
        sds_rgb(2)='l1b_ch21_data'
        gains(1,2)=gsys(2)
        gains(2,2)=gsys(2)
        sds_rgb(3)='l1b_ch20_data'
        gains(1,3)=gsys(1)
        gains(2,3)=gsys(1)
c
       else
c
        sds_rgb(1)='l1b_ch13_data'
c        sds_rgb(1)='l1b_ch15_data'
c        sds_rgb(1)='l1b_ch17_data'
        gains(1,1)=gsys(vnirb(13))
        gains(2,1)=gsys(vnirb(13))
c
        sds_rgb(2)='l1b_ch8_data'
        gains(1,2)=gsys(vnirb(8))
        gains(2,2)=gsys(vnirb(8)+1)
c
        sds_rgb(3)='l1b_ch5_data'
        gains(1,3)=gsys(vnirb(5))
        gains(2,3)=gsys(vnirb(5)+1)
c
c      ------ for modis cnv -----
        if(modcnv.gt.0) then
         sds_rgb(1)='l1b_ch10_data'
         gains(1,1)=gsys(vnirb(10))
         gains(2,1)=gsys(vnirb(10))
         sds_rgb(2)='l1b_ch9_data'
         gains(1,2)=gsys(vnirb(9))
         gains(2,2)=gsys(vnirb(9))
        endif
c      --------------------------
c
       endif
c
       i_high=0
       do k=3,1,-1
        sds_index=sfn2index(sd_id2,sds_rgb(k))
        if(sds_index.lt.0) then
         write(6,*) 'no SDS data=',sds_rgb(k)
         goto 777
        endif
        sds_id=sfselect(sd_id2,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
        nrmax=0
        if(sds_rgb(k).eq.'l1b_ch22_data') then
         rmax=150./gains(1,k)
         do nm=1,npixels*nlines
          if(l2datai2(nm).gt.rmax) then
           l2datai2(nm)=-16383
           nrmax=nrmax+1
          endif
         enddo
         write(6,'(a11,f5.1,a2)') 
     &             'Saturation=',100.*nrmax/npixels/nlines,' %'
        endif
c
        dmin=100000.
        dmax=0.
        davg=0.
        dnum=0.
        do m1=1,mll
         do n1=1,nll
          nm=map(n1,m1)
          m=(nm-1)/npixels+1
          n=nm-(m-1)*npixels
          if((m.gt.0)  .and.(n.gt.0).and.
     &       (m.le.nlines).and.(n.le.npixels)) then
           if(l2datai2(nm).gt.0) then
            if(l2datai2(nm).gt.16383) l2datai2(nm)=l2datai2(nm)-16384
            if(l2datai2(nm).gt.8191)  l2datai2(nm)=l2datai2(nm)-8192
            if(l2datai2(nm).gt.4095) then
             rdat=(l2datai2(nm)-4096)*gains(1,k)
            else
             rdat=l2datai2(nm)*gains(2,k)
            endif
            if(rdat.lt.dmin) dmin=rdat
            if(rdat.gt.dmax) dmax=rdat
            davg=davg+rdat
            dnum=dnum+1
            rdat=rdat/gains(1,k)
            if(rdat.gt.32767) rdat=32767
            i2_out(n1,m1)=nint(rdat)
           elseif(l2datai2(nm).lt.0) then
            rdat=mod(l2datai2(nm)+65536,8192)
            if(rdat.gt.4095) then
             i2_out(n1,m1)=nint(rdat-4096)
            else
             i2_out(n1,m1)=32767
            endif
           else
            i2_out(n1,m1)=0
           endif
          else
           i2_out(n1,m1)=-1
          endif
         enddo
        enddo
        if(dmin.lt.d_min(k)) d_min(k)=dmin
        if(dmax.gt.d_max(k)) d_max(k)=dmax
        if(d_max(k).gt.d_cld(k)) then
         d_max(k)=d_cld(k)
         i_high=1
        elseif(i_high.gt.0) then
         d_max(k)=d_cld(k)
        endif
        if(nrmax.gt.10) d_max(k)=d_cld(k)        ! for 250m
c
        if(k.eq.1) d_max(1)=d_cld(1)*(d_max(2)/d_cld(2)) ! for saturation R
c
        if(dnum.gt.0) davg=davg/dnum
        if((davg.gt.d_min(k)).and.(davg.lt.d_max(k))) then
c         if(davg.gt.d_avg(k)) d_avg(k)=davg
c         if(davg.lt.d_avg(k)) d_avg(k)=davg
c         if((davg.lt.d_avg(k)*1.5).and.(davg.gt.d_avg(k)*0.5))
c     &       d_avg(k)=davg
c         d_avg(k)=davg
        else
         d_avg(k)=(d_min(k)+d_max(k))/2.
        endif
c
c     -------------------------------- fix scale
        if(f_max(k).gt.1.) then
         d_max(k)=f_max(k)
         d_min(k)=f_min(k)
         d_avg(k)=f_avg(k)
        endif
c     --------------------------------
c
        write(6,'(6f10.3)') dmin,davg,dmax,d_min(k),d_avg(k),d_max(k)
c
c        rh1=0.50
        rh1=0.66       ! for land
        do m1=1,mll
         do n1=1,nll
          if(i2_out(n1,m1).gt.-1) then
           rdat=i2_out(n1,m1)*gains(1,k)
           if(rdat.lt.d_avg(k)) then
            rdat=(rdat-d_min(k))/(d_avg(k)-d_min(k))*rh1
           else
            rdat=(rdat-d_avg(k))/(d_max(k)-d_avg(k))*(1.-rh1)+rh1
           endif
           rdat=rdat*254.
           if(rdat.gt.254.) then
            rdat=254.
           elseif(rdat.lt.0.) then
            rdat=0.
           endif
           i1_out(n1,m1,k)=char(int(rdat))
          else
           i1_out(n1,m1,k)=a7
          endif
         enddo
        enddo
c
       enddo
c
c     --------------------- for L1B band data ---------------------
      elseif(sds0(1:3).eq.'l1b') then
       iotype=4
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
       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     ------------------- for CHLA and SST --------------------
      else
       nbl=1
       sds_index=sfn2index(sd_id,sds0)
       if(sds_index.lt.0) then
        write(6,*) 'no SDS data=',sds0
        goto 777
       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)
       do m1=1,mll
        do n1=1,nll
         nm=map(n1,m1)
         m=(nm-1)/npixels+1
         n=nm-(m-1)*npixels
         if((m.gt.0)  .and.(n.gt.0).and.
     &      (m.le.nlines).and.(n.le.npixels)) then
          if(iflag(n,m).lt.1) then
           nm=(m-1)*npixels+n
           i2_out(n1,m1)=l2datai2(nm)
          else
           i2_out(n1,m1)=-1
          endif
         else
          i2_out(n1,m1)=-2
         endif
        enddo
       enddo
       write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
c     --------------------------------------------------------------
      endif
c
c     ------------------- overlay RGB --------------------
      if(iotype.eq.3) then
       nbl=1
c
       i4=0
       do i=1,4
        if(sds0(i+1:i+3).eq.'rgb') i4=i
       enddo
       if(sds0(i4:i4).eq.'_') i4=i4-1
c
c       open(30,file='GRAD.pal',READONLY,
       open(30,file='GRAD_black.pal',READONLY,
     & status='old',form='unformatted',access='direct',recl=256*3)
       read(30,rec=1) ((pal(i,k),i=1,256),k=1,3)
       close(30)
c
       sds_index=sfn2index(sd_id,sds0(1:i4))
       if(sds_index.lt.0) then
        write(6,*) 'no SDS data=',sds0(1:i4)
        goto 777
       endif
c
       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)
c
       if(sds0(1:i4).eq.'sst') then
        sstmax=-2.
        sstmin=60.
        do m1=1,mll
         do n1=1,nll
          nm=map(n1,m1)
          m=(nm-1)/npixels+1
          n=nm-(m-1)*npixels
          if((m.gt.0)  .and.(n.gt.0).and.
     &       (m.le.nlines).and.(n.le.npixels)) then
           nm=(m-1)*npixels+n
           if((iflag(n,m).lt.1).and.(l2datai2(nm).gt.0)) then
            sst=l2datai2(nm)*slope+offset-273.15
            if(sst.gt.sstmax) sstmax=sst
            if(sst.lt.sstmin) sstmin=sst
           endif
          endif
         enddo
        enddo
        if(sstmax.gt.36.25) sstmax=36.25
        if(sstmin.lt.-2.00) sstmin=-2.00
        rng=sstmax-sstmin
        sstmax=sstmax+rng*0.05
        sstmin=sstmin-rng*0.05
        if(rng.lt.3.6) then
         rng=4.
         sstmax=(sstmin+sstmax+rng)/2.
         sstmin=(sstmin+sstmax-rng)/2.
        endif
        sstmin=int(sstmin)
        rng=sstmax-sstmin
        dnc=rng/255.
        write(6,'(3f8.4)') sstmin,sstmax,dnc
       endif
c
       do m1=1,mll
        do n1=1,nll
         nm=map(n1,m1)
         m=(nm-1)/npixels+1
         n=nm-(m-1)*npixels
         if((m.gt.0)  .and.(n.gt.0).and.
     &      (m.le.nlines).and.(n.le.npixels)) then
          nm=(m-1)*npixels+n
          if((iflag(n,m).lt.1).and.(l2datai2(nm).gt.0)) then
           if(sds0(1:i4).eq.'sst') then
            sst=l2datai2(nm)*slope+offset-273.15
            idata=nint((sst-sst_min)/dnc)
           elseif(sds0(1:i4).eq.'chla') then
            idata=nint((log10(l2datai2(nm)*slope+offset)+2.)/0.015)
           else
            idata=nint(l2datai2(nm)/256.)
           endif
           if(idata.gt.255) idata=255
           if(idata.lt.  0) idata=0
           do k=1,3
            i1_out(n1,m1,k)=pal(idata+1,k)
           enddo
          endif
         endif
        enddo
       enddo
       write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &               (dims(ir),ir=1,rank)
c
       if(sds0(1:i4).eq.'sst') then
        slope=dnc
        offset=sstmin
       elseif(sds0(1:i4).eq.'chla') then
        slope=0.015
        offset=-2.
       endif
      endif
c     ---------------------------------------------------------------
c
      nbg=0
      nbgain=0
      glivn=0
      do nb=1,nbl
c     -------------------------------------L1B channels
       if((nbl.gt.1).or.(sds0(1:6).eq.'l1b_ch')) then
        glivn=1
        nbgain=nbgain+1
c 
c     ------------------------------------- Search band
c
  888   nbg=nbg+1
        if((nbg.gt.36).or.(nbgain.gt.nbl)) then
c         write(6,*) nbg,nbgain
         goto 999
        endif
       
        if(nbg.lt.10) then
         write(sds1,'(a6,i1.1,a5)') 'l1b_ch',nbg,'_data'
        elseif((npixels.lt.1400).and.
     &         ((nbg.eq.28).or.(nbg.eq.29)) ) then
         write(sds1,'(a6,i2.2,a9)') 'l1b_ch',nbg,'_data_2km'
        else
         write(sds1,'(a6,i2.2,a5)') 'l1b_ch',nbg,'_data'
        endif
        sds_index=sfn2index(sd_id,sds1)
        if(sds_index.lt.0) then
c         write(6,*) sds_index
         goto 888
        endif
c     ---------------------------------------------------
c
        slope=gsys(nbgain)
c
        if((nbg.eq.4).or.(nbg.eq.5).or.
     &     (nbg.eq.7).or.(nbg.eq.8)) then
         nbgain=nbgain+1
         slope2=gsys(nbgain)
        else
         slope2=slope
        endif
c        write(6,*) nb,nbg,slope,slope2
c
        if((sds0(1:6).eq.'l1b_ch').and.
     &     (sds1(1:8).ne.sds0(1:8))) then
c         write(6,*) sds0(1:8),sds1(1:8)
         goto 666
        endif
c
        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)
c
        do m1=1,mll
         do n1=1,nll
          nm=map(n1,m1)
          m=(nm-1)/npixels+1
          n=nm-(m-1)*npixels
          if((m.gt.0).and.(n.gt.0).and.
     &       (m.le.nlines).and.(n.le.npixels)) then
c
           if(sds1(15:17).eq.'2km') then
            m=nint(m*0.5)
            n=nint(n*0.5)
            nm=(m-1)*(npixels_2k)+n
           endif
c
           if(l2datai2(nm).lt.0) then
c           if(l2datai2(nm).le.0) then
            rdat=mod(l2datai2(nm)+65536,8192)
            if((rdat.gt.4095).and.(l2datai2(nm).lt.-16384)) then
             i2_out(n1,m1)=nint(rdat-4096)
c             write(6,*) l2datai2(nm)
            else
             i2_out(n1,m1)=-2
c             i2_out(n1,m1)=32767
            endif
           else
            if(l2datai2(nm).gt.16383) l2datai2(nm)=l2datai2(nm)-16384
            if(glivn.lt.1) then
             i2_out(n1,m1)=l2datai2(nm)
            else
             if(l2datai2(nm).gt.4095) then
              i2_out(n1,m1)=l2datai2(nm)-4096
             else
              idata=nint(l2datai2(nm)*slope2/slope)
              if(idata.gt.32767) then
               i2_out(n1,m1)=idata-65536
              else
               i2_out(n1,m1)=idata
              endif
             endif
            endif
           endif
c
          else
           i2_out(n1,m1)=-1
          endif
         enddo
        enddo
        write(6,4000) 'SDS=',sds_index,sds_name,para_type(data_type),
     &                (dims(ir),ir=1,rank)
       else
        sds1=sds0
       endif
c     -----------------------------------------------------------------
c
c     -------------------------- Output ------------------------------
       outfile( 1: 25)='                         '
       outfile(26: 55)='                              '
       head( 1: 40)='                                        '
       head(41: 80)='                                        '
       head(81:125)='                                             '
c
       ll=0
       do l=1,17
        if(sds1(l:l).ne.' ') ll=l
       enddo
c
       if(i3.lt.1) then
        write(outfile,'(a24,a1,i5.5,a1,i5.5,a1)')
     &         infile(1:24),'_',nll,'_',mll,'_'
       else
        write(outfile,'(a24,a1,i5.5,a1,i5.5,a1)')
     &         Product_Name(i3:i3+23),'_',nll,'_',mll,'_'
       endif
       outfile(38:37+ll)=sds1(1:ll)
c
c     ---------------------------------------- output header
       if((sds1(1:4).eq.'nLw_').or.
     &    (sds1(1:3).eq.'par').or.(sds1(1:3).eq.'tau').or.
     &    (sds1(1:3).eq.'alp')) then
        write(head,'(2i6,3f8.2,2f9.4,a1,a8,a1,a55)') 
     &     nll,mll,lon_min,lat_max,dd,slope,offset,',',
     &     sds1(1:8),',',outfile(1:37+15)
       elseif(infile(21:22).eq.'NL') then
        write(head,'(2i6,3f8.2,2f9.4,a1,a8,a1,a55)') 
     &     nll,mll,lon_min,lat_max,dd,1.,offset,',',
     &     sds1(1:8),',',outfile(1:37+ll)
       elseif((sds1(1:3).eq.'l1b').and.(sds1(5:7).ne.'rgb')) then
        if(dd.gt.0.0099) then
         write(head,'(2i6,3f8.2,2e12.5,a1,a8,a1,a55)') 
     &      nll,mll,lon_min,lat_max,dd,slope,offset,',',
     &      sds1(1:8),',',outfile(1:37+ll)
        else
         write(head,'(2i6,2f8.2,f8.4,2e12.5,a1,a8,a1,a55)') 
     &      nll,mll,lon_min,lat_max,dd,slope,offset,',',
     &      sds1(1:8),',',outfile(1:37+ll)
        endif
       else
        if(dd.gt.0.0099) then
         write(head,'(2i6,3f8.2,2f9.4,a1,a8,a1,a55)') 
     &      nll,mll,lon_min,lat_max,dd,slope,offset,',',
     &      sds1(1:8),',',outfile(1:37+ll)
        else
         write(head,'(2i6,2f8.2,f8.4,2f9.4,a1,a8,a1,a55)') 
     &      nll,mll,lon_min,lat_max,dd,slope,offset,',',
     &      sds1(1:8),',',outfile(1:37+ll)
        endif
       endif
c
       if(nb.eq.1) write(6,*) 'head info=',head(1:62)
       write(6,*) 'OutPut=',outfile(1:37+ll)
c     ----------------------------------------
c
       if((sds1(1:3).eq.'nLw').and.(sds1(4:4).ne.'_')) then
        do m=1,mll
         do n=1,nll
c
          ifw=0
c        -- nLw-R <0; R->0 --
          if((i1_out(n,m,1).eq.a0).and.
     &       (ichar(i1_out(n,m,2)).gt.63)) ifw=1
          do k=2,3
           if(i1_out(n,m,k).eq.a0) ifw=1
          enddo
          if(ifw.gt.0) then
           do k=1,3
            i1_out(n,m,k)=a9
           enddo
          endif
c        ---------------------
c
c        --------------------- gray to black **
          ifb=0
          do k=1,3
           if(i1_out(n,m,k).eq.a7) ifb=ifb+1
          enddo
          if(ifb.gt.2) then
           do k=1,3
            i1_out(n,m,k)=a6
           enddo
          endif
c        -------------------------------------
         enddo
        enddo
c
       elseif((sds1(1:6).eq.'RcRefl').or.
     &        (sds1(5:7).eq.'rgb')) then
c
c     --------------------------------------------- for 250m
        if(Product_Name(i3:i3+4).eq.'A2GL2') then
         do m=1,mll
          do n=1,nll
           if(i1_out(n,m,1).eq.a8) then
            i1_out(n,m,1)=i1_out(n,m,2)
           endif
          enddo
         enddo
        endif
c     --------------------------------------------- 
c
        do m=1,mll
         do n=1,nll
c
c     -------------------------------------------- for error
          ifw=0
          do k=1,3
           if(i1_out(n,m,k).eq.a0) ifw=ifw+1
          enddo
          if(ifw.gt.2) then
           do k=1,3
            i1_out(n,m,k)=a9
           enddo
          endif
c
c        ------------------------------------ gray to black **
          ifb=0
          do k=1,3
           if(i1_out(n,m,k).eq.a7) ifb=ifb+1
          enddo
          if(ifb.gt.2) then
           do k=1,3
            i1_out(n,m,k)=a6
           enddo
          endif
c        --------------------------------------------
c
c     -------------------------------------------- for saturation
          ifw=0
          do k=1,3
           if(i1_out(n,m,k).eq.a8) ifw=1
          enddo
          if(ifw.gt.0) then
           do k=1,3
            i1_out(n,m,k)=a8
           enddo
          endif
         enddo
        enddo
       endif
c
       if(((sds0(1:4).ne.'nLw_').and.(sds0(1:3).eq.'nLw')).or.
     &    (sds0(1:6).eq.'RcRefl').or.(sds1(5:7).eq.'rgb')) then
        open(20,file=outfile(1:37+ll),status='unknown'
     &  ,access='direct',recl=nll,form='unformatted')
        write(20,rec=1) head(1:115)
        mc=1
        do k=1,3
         do m=1,mll
          mc=mc+1
          write(20,rec=mc) (i1_out(n,m,k),n=1,nll)
         enddo
        enddo
        close(20)
       else
        open(20,file=outfile(1:37+ll),status='unknown'
     &  ,access='direct',recl=nll*2,form='unformatted')
        write(20,rec=1) head(1:115)
        do m=1,mll
         write(20,rec=m+1) (i2_out(n,m),n=1,nll)
        enddo
        close(20)
       endif
c
c     ---------------------------------------------------------------
  666 continue
      enddo
c
c     -------------------------- Close HDF file
  999  retn=sfend(sd_id)
c     ------------------------------------------
c
 777  continue
      stop
      end
c
c
c==========================================================
      SUBROUTINE LINEARMAP
     & (npixel,nspixel,nsline,nspixels,nslines,npixels,nlines,
     &  nscans,ndets,lat_max,lat_min,lon_max,lon_min,dd,
     &  nsdet,nsdets,nl,ml,nll,mll,pxl,det,lat0,lon0,map)
c==========================================================
      real*4       lat_max,lat_min,lon_max,lon_min
      real*4       lat0(nspixel*nsline)
      real*4       lon0(nspixel*nsline)
      integer*2    pxl(npixel),det(nsdet)    ! npixel for mud
      integer*4    map(nl,ml)
c
      real*4       rxlonr(2),rxlatr(2)
      real*8       rxlon(0:3),rxlat(0:3),rxpix(0:3),rxlin(0:3)
      real*8       aa1,aa2,bb1,bb2,aa0,aaa1,aaa2,bbb1,bbb2
      real*8       rlati,rlong
      integer*4    ixlonr(2),ixlatr(2)
      real*4       rlon0(0:1,0:1),rlat0(0:1,0:1)
      real*4       rlon(0:1,0:1),rlat(0:1,0:1)
      real*4       slon
c
      do m=1,mll
      do n=1,nll
       map(n,m)=-1
      enddo
      enddo
      n_error=0
c
c      ddx=120*dd
      ddx=160*dd
      ddy=80*dd
      do nsc=1,nscans
c      do nsc=74,74
       m=(nsc-1)*nsdets+1
       do n=1,nspixels-1
c       do n=90,90
        nm0=(m-1)*nspixels+n
c
        ifg=0
        do md=0,1
        do nd=0,1
         nm=(m+md-1)*nspixels+n+nd
         if((lon0(nm).lt.lon_min-ddx).or.
     &      (lon0(nm).gt.lon_max+ddx).or.
     &      (lat0(nm).lt.lat_min-ddy).or.
     &      (lat0(nm).gt.lat_max+ddy)) then
          ifg=1
         endif
        enddo
        enddo
c
        if(ifg.lt.1) then
c
c     -------------------- define subareas -------------------
c
         slon=lon0((m-1)*nspixels+n)
         do md=0,1
          do nd=0,1
           nm=(m+md-1)*nspixels+n+nd
           if(lon0(nm)-slon.gt.180.) then
            rlon0(nd,md)=lon0(nm)-360.
           elseif(lon0(nm)-slon.lt.-180.) then
            rlon0(nd,md)=lon0(nm)+360.
           else
            rlon0(nd,md)=lon0(nm)
           endif
           rlat0(nd,md)=lat0(nm)
          enddo
         enddo
c
         rt1=1.*(   1.-det(1))/(det(nsdets)-det(1))
         rt2=1.*(ndets-det(1))/(det(nsdets)-det(1))
c
         rlon(0,0)=(1.-rt1)*rlon0(0,0)+rt1*rlon0(0,1)
         rlon(1,0)=(1.-rt1)*rlon0(1,0)+rt1*rlon0(1,1)
         rlon(0,1)=(1.-rt2)*rlon0(0,0)+rt2*rlon0(0,1)
         rlon(1,1)=(1.-rt2)*rlon0(1,0)+rt2*rlon0(1,1)
c
         rlat(0,0)=(1.-rt1)*rlat0(0,0)+rt1*rlat0(0,1)
         rlat(1,0)=(1.-rt1)*rlat0(1,0)+rt1*rlat0(1,1)
         rlat(0,1)=(1.-rt2)*rlat0(0,0)+rt2*rlat0(0,1)
         rlat(1,1)=(1.-rt2)*rlat0(1,0)+rt2*rlat0(1,1)
c
         rxlonr(1)= 999.
         rxlonr(2)=-999.
         rxlatr(1)= 999.
         rxlatr(2)=-999.
         do md=0,1
          do nd=0,1
           if(rlon(nd,md).lt.rxlonr(1)) rxlonr(1)=rlon(nd,md)
           if(rlon(nd,md).gt.rxlonr(2)) rxlonr(2)=rlon(nd,md)
           if(rlat(nd,md).lt.rxlatr(1)) rxlatr(1)=rlat(nd,md)
           if(rlat(nd,md).gt.rxlatr(2)) rxlatr(2)=rlat(nd,md)
          enddo
         enddo
c
         ixlonr(1)=nint((rxlonr(1)-lon_min)/dd+0.5)
         ixlonr(2)=nint((rxlonr(2)-lon_min)/dd+1.5)
         ixlatr(1)=nint((lat_max-rxlatr(2))/dd+0.5)
         ixlatr(2)=nint((lat_max-rxlatr(1))/dd+1.5)
c         write(6,*) ixlatr(1),ixlatr(2),ixlonr(1),ixlonr(2)
c
c     ------------------ calculate coefficients ---------------
         do md=0,1
          do nd=0,1
           rxlon(md*2+nd)=1.D0*(rlon0(nd,md)-rlon0(0,0))
           rxlat(md*2+nd)=1.D0*(rlat0(nd,md)-rlat0(0,0))
           rxpix(md*2+nd)=1.D0*(pxl(n+nd)-pxl(n))
           rxlin(md*2+nd)=1.D0*(det(1+md)-det(1))
          enddo
         enddo
c
         aaa1=0.D0
         aaa2=0.D0
         bbb1=0.D0
         bbb2=0.D0
         numa=0
         do icn=1,3
          k1=icn
          k2=icn+1
          if(k2.gt.3) k2=1
          aa0=rxlon(k1)*rxlat(k2)-rxlon(k2)*rxlat(k1)
          if(aa0.ne.0.D0) then
           aa1=(rxlat(k2)*rxpix(k1)-rxlat(k1)*rxpix(k2))/aa0
           aa2=(rxlon(k1)*rxpix(k2)-rxlon(k2)*rxpix(k1))/aa0
           bb1=(rxlat(k2)*rxlin(k1)-rxlat(k1)*rxlin(k2))/aa0
           bb2=(rxlon(k1)*rxlin(k2)-rxlon(k2)*rxlin(k1))/aa0
           numa=numa+1
           aaa1=aaa1+aa1
           aaa2=aaa2+aa2
           bbb1=bbb1+bb1
           bbb2=bbb2+bb2
          endif
c          write(6,'(i2,4f9.3)') icn,aa1,aa2,bb1,bb2
         enddo
         if(numa.gt.0) then
          aaa1=aaa1/numa
          aaa2=aaa2/numa
          bbb1=bbb1/numa
          bbb2=bbb2/numa
         else
c          write(6,*) 'error (n,m)=', n,m
          n_error=n_error+1
         endif
c         write(6,*) aaa1,aaa2,bbb1,bbb2,numa
c     ----------------------------------------------------------
c
         do mm=ixlatr(1),ixlatr(2)
         do nn=ixlonr(1),ixlonr(2)
          rlati=1.D0*(lat_max-(mm-1.)*dd)-rlat0(0,0)
          rlong=1.D0*(lon_min+(nn-1.)*dd)-rlon0(0,0)
          rnd=float(aaa1*rlong+aaa2*rlati)
          rmd=float(bbb1*rlong+bbb2*rlati)
c          rnc=pxl(n)+rnd
          rnc=pxl(n)+rnd+1
          rmc=(nsc-1.)*ndets+det(1)+rmd
          nc=nint(rnc)
          mc=nint(rmc)
          dnscn1=(rmd+0.5)/ndets
          dnpix1=(rnd+0.5)/(pxl(n+1)-pxl(n))
c          write(6,*) dnscn1,dnpix1
c
          if((mm.gt.0).and.(nn.gt.0).and.
     &       (mm.le.mll).and.(nn.le.nll).and.
     &       (mc.gt.0).and.(nc.gt.0).and.
     &       (mc.le.nlines).and.(nc.le.npixels)) then
c
           if((dnpix1.ge.0.0).and.(dnpix1.le.1.0).and.
     &        (dnscn1.ge.0.0).and.(dnscn1.le.1.0)) then
             nscn1=nsc
           elseif(map(nn,mm).lt.1) then
            if((dnpix1.gt.-0.05).and.(dnpix1.lt.1.05).and.
     &         (dnscn1.gt.-0.08).and.(dnscn1.lt.1.08)) then
             nscn1=nsc
            else
             nscn1=0
            endif
           else
            nscn1=0
           endif
c
           if(nsc.eq.nscn1) then
            nm=(mc-1)*npixels+nc
            map(nn,mm)=nm
           endif
c
          endif
         enddo
         enddo
        endif
       enddo
      enddo
      if(n_error.gt.0) write(6,*) 'Mapping error =',n_error
c
      return
      end
c==========================================================
