;================================================================
;g11_smoke_main_new4_firemv_case_HI.pro
;Retrieve biomass burning aerosols from GOES-11
;change 1: move fires out of the GASP image boundary in AK
;change 2: GASP data added scattering angle, new standard deviation,
;                   and surface reflectance, azimuth bug free
;change 3: Added hmshysplit fire duration
;change 4: new GASP output
;change5: in HI, loose aodstd=0.5, bulk_aod searching cirlcle limits to 10.
;created: 2009--03-06
;author: Jian Zeng (ERT)
;
; Last Modified:2010-05-20            by J. Yang
;               For Operational in SSD
;---------------------------------------------------------------------------
pro hist_az,aod,az,direc,bin,bulk_aod

    amin=0
    amax=360
    bin=10
    nn=fix((amax-amin)/bin)
    xx=indgen(nn+1)*bin+amin+bin/2.
    hist=histogram(az,min=amin,max=amax,bin=bin)
    result=max(hist,loc)
    ind0=where(hist eq result,cont)

    if cont eq 1 then begin
      direc=xx(loc)
      ind=where(az ge direc-bin/2. and az le direc+bin/2.)
      if max(aod(ind)) gt 0.2 then bulk_aod=max(aod(ind)) else bulk_aod=max(aod)
    endif else begin
      bulk_aod=max(aod,loc)
      direc=az(loc)
    endelse

end
;---------------------------------------------------------------------------------
pro find_bulk_aod,maxaod,maxaz,maxratio,maxdist,bulk_aod,direc

ind=where(maxratio gt 50.,cp)
bulk_aod=-999.0
direc=-999.0

if cp gt 0L then begin
  aod=maxaod(ind)
  az=maxaz(ind)
  if n_elements(ind) gt 10 then begin
    hist_az,aod,az,direc,bin,bulk_aod
  endif else begin
    bulk_aod=max(aod,loc)
    direc=az(loc)
  endelse

endif else begin  ;add else on Jan 15, 2007
  ;bulk_aod=max(maxaod(0:10),loc)
  masz=size(maxaod,/dim)
  if masz(0) gt 5 then begin
     if masz(0) gt 10 then cc=10 else cc=masz(0)-1
     bulk_aod=max(maxaod(0:cc),loc)
     direc=maxaz(loc)
  endif
endelse

end
;--------------------------------------------------------------------------
pro find_sqr,i0,j0,outi,outj,lev,img,nsz,aod,dstep

  sz=size(img,/dim)
  outi=lonarr(sz(0)*sz(1))
  outj=lonarr(sz(0)*sz(1))
  num=0
  
  for ii=1,sz(0)-2 do begin
  for jj=1,sz(1)-2 do begin
  
    if img(ii,jj) lt -9998.0 then begin
      s1=(ii-1) > 0
      s2=(ii+1) < (sz(0)-1)
      e1=jj-1 > 0
      e2=jj+1 < (sz(1)-1)
      temp=img(s1:s2,e1:e2)
      ind=where(temp gt -9998.0, num0)
      if num0 gt 0 then begin
        aa=mean(temp(ind))
      endif else aa=img(ii,jj)
    endif else aa=img(ii,jj)

    if aa ge lev then begin
       
        outi(num)=(ii+i0-dstep) > 0 < nsz(0)
        outj(num)=(jj+j0-dstep) > 0 < nsz(1)
        num=num+1
    endif
  endfor
  endfor

  if num gt 0 then begin
    outi=outi(0:num-1)
    outj=outj(0:num-1)
  endif else begin
    outi=-1
    outj=-1
  endelse
  
end
;-------------------------------------------------------------
pro find_plume,lev,ilon0,ilat0,img,ind_out

  sz0=size(img,/dimension)

   ;dilute method
   doind=[ilon0+ilat0*sz0(0)]
   oldind=[ilon0+ilat0*sz0(0)]

   con=1
   circle=0
   
   while (con) do begin
     nnew=0
     newind=[-1]
     circle=circle+1
     for ido=0,n_elements(doind)-1 do begin
       ind0=doind(ido)

       if circle eq 1 then dstep=5 else dstep=2
       
       result=array_indices(img,ind0)
       i0=result(0)
       j0=result(1)
       s1=(i0-dstep) > 0
       s2=(i0+dstep) < (sz0(0)-1)
       e1=j0-dstep > 0
       e2=j0+dstep < (sz0(1)-1)
       temp=img(s1:s2,e1:e2)
       
       find_sqr,i0,j0,outi,outj,lev,temp,sz0,img,dstep

       if max(outi) gt 0L then begin
         outind=outi+outj*sz0(0)

         if min(outind) lt 0L then goto,nonew
         result=where(outind eq ind0,ncompl=nnew,compl=new)
         if nnew gt 0 then begin
           for inew=0,nnew-1 do begin
             ind=new(inew)
             result=where(oldind eq outind(ind),ndl)
             if ndl eq 0 then begin
               if (newind(0) eq -1) then $
                 newind=[outind(ind)] $
               else begin
                 result=where(oldind eq outind(ind),ndl2)
                 if ndl2 eq 0 then $
                   newind=[newind,outind(ind)]
               endelse
             endif
           endfor
         endif
         nonew:
       endif

     endfor
    
     if newind(0) ne -1 then begin
        result=sort(newind)
        doind=newind(uniq(newind,result))
        oldind=[oldind,doind]
     endif else con=0


   endwhile

  ind_out=oldind

end
;--------------------------------------------------------------
function doublefind,old_ind,new_ind

  all=[old_ind,new_ind]
  result=all(uniq(all,sort(all)))

  return, result

end
;--------------------------------------------------------------------------------
pro convert_hourly,jtime,dir_temp,dir_grib,pref,nlon,nlat,yr,month,nday,$
                   dlon,dlat,clon,clat,lon,lat,$
                   glon=glon,glat=glat,num0_smokeaod=num0_smokeaod,$
                   sum0_smokeaod=sum0_smokeaod,$
                   num0_cloud=num0_cloud,num0_nonsmoke=num0_nonsmoke,$
                   num0_nofire=num0_nofire,num0_nodata=num0_nodata,grib=grib,$
                   data_suffix=data_suffix,hysplit_data=hysplit_data
  
  shr=0 ;11
  
  for kk=shr,jtime-1 do begin
    if kk le 3 or kk ge 11 then begin
    outfl_f77_NAM=dir_temp+pref+'.'+yr+month+nday+string(kk,format='(i2.2)')+$
    data_suffix+'.NAM3.dat'

    remap_GASPsmoke_NAM,nlat,nlon,outdata=comb_NAM,$
                        num_smokeaod=reform(num0_smokeaod(*,*,kk)),$
                        sum_smokeaod=reform(sum0_smokeaod(*,*,kk)),$
                        num_cloud=reform(num0_cloud(*,*,kk)),$
                        num_nonsmoke=reform(num0_nonsmoke(*,*,kk)),$
                        num_nofire=reform(num0_nofire(*,*,kk)),$
                        num_nodata=reform(num0_nodata(*,*,kk))
    combconc_NAM=conv_smoke_combaod_combconc(comb_NAM)
    combconc=reform(combconc_NAM(*,*,0))
    
    if grib eq 1 then begin
        outfl_combaod_NAM=dir_temp+pref+'.'+yr+month+nday+$
        string(kk,format='(i2.2)')+data_suffix+'.NAM3.dat'

        write_binary_all,outfl_combaod_NAM,nlon,nlat,lon,lat,$
                         comb_NAM,combconc
        print,'created----',outfl_combaod_NAM
        outfl_grib_NAM=dir_grib+pref+'.'+yr+month+nday+$
            string(kk,format='(i2.2)')+data_suffix+'.NAM3.grd'
       
        ndate=yr+month+nday+string(kk,format='(i2.2)')+'00'

        openw,1,'infile0'
        printf,1,outfl_combaod_NAM
        printf,1,outfl_grib_NAM
        printf,1,strmid(ndate,2,10)
        close,1
        spawn,'../GRIB/grib_main_NAM3'
        ;spawn,'/home/jzeng/GRIB/GRIB/grib_main_NAM3'
        spawn,'rm -f infile0'
        print,'created '+outfl_grib_NAM

        ;file_move,outfl_grib_NAM,dir_grib,/overwrite
        ;print,'file_move '+outfl_grib_NAM+' to '+dir_grib
        
    endif ;if grib eq 1

    
    if hysplit_data eq 1 then begin
       write_binary_f77,outfl_f77_NAM,combconc
    endif

    endif ; if kk le 3  
  endfor

  
end
 
;----------------------------------------------------------------------------------------------
pro find_pattern_boundary,paod,out_ind,nbulkaod=nbulkaod,bulkaod=bulkaod
    
    out_ind=[-1]
    if max(paod) gt 0. then begin
      if keyword_set(nbulkaod) then begin
        if nbulkaod gt 0 then begin
          maod=mean(bulkaod(0:nbulkaod-1))
          if maod gt 0.5 then pnt=0.2 else pnt=0.4
          thresh=maod*pnt
        endif else thresh=max(paod)*0.05  ;0.01 - 05/05/2009
        out_ind=where(paod ge thresh)
      endif ;else thresh=0  ;max(paod)*0.05
     
    endif
end

;----------------------------------------------------------------------------------
pro disp_aircon,ch1,lat,lon,outfile=outfile,data=data,$
                title=tit,fires_lat=fires_lat,fires_lon=fires_lon,limit=limit
    if max(ch1 gt 0.) then begin
    img=hist_equal(ch1) / 2. + 131 < 255
    img_range=[0,1.2]
    ind0=where(data gt img_range(0),cp)
    if cp gt 0 then img(ind0) = bytscl((data(ind0)) < img_range(1))/2. < 127 >1

    set_plot,'Z'
    loadct,0,/silent
    tvlct,r0,g0,b0,/get
    loadct,39,/silent
    tvlct,r1,g1,b1,/get

    r=bytarr(256)+0B
    g=bytarr(256)+0B
    b=bytarr(256)+0B
    bot=50
    top=253
    ncolors=top-bot+1
    r(1:127) = congrid(r1(bot:top),127) 
    g(1:127) = congrid(g1(bot:top),127)
    b(1:127) = congrid(b1(bot:top),127)
    ncolors=127
    r(0) = 0
    g(0) = 0
    b(0) = 0

    r(128) = 255
    g(128) = 130
    b(128) = 0
    ;plume color=129
    r(129) = 255
    g(129) = 100
    b(129) = 255
     r(130) = 10
    g(130) = 10
    b(130) = 10

; the upper half(129:255) is the gray scale from colortable=0

    r(130:255) = congrid(r0,126)
    g(130:255) = congrid(g0,126)
    b(130:255) = congrid(b0,126)


set_plot,'ps'
DEVICE,/inches,xoffset=0.5,yoffset=1,xsize=7.5,ysize=7,/COLOR,$
       filename=outfile+'.ps',bits_per_pixel=8
TVLCT, 0, 0, 0, !p.background
!p.charthick=3.0
!p.charsize=1.0
!p.thick=3.0
!p.font=2

a=findgen(16)*(!PI*2/16.)
usersym,cos(a),sin(a),/fill

tvlct,r,g,b
imagemap_unscl,img,lat,lon,limit=limit,missing=0,mapcolor=0,$
         position=[0.05,0.18,0.95,0.93],/noborder
map_continents,/usa,/cont,/countries,/noerase,/HIRE,color=133;129
colorbar,position=[0.1,0.08,0.9,0.11],minrange=img_range(0),$
         maxrange=img_range(1),format='(f5.1)',div=6,$
         col=130,bottom = 1,ncolors=ncolors,chars = 1.1,charthick=2.0,$
         title='Aerosol Optical Thickness'
xyouts,0.916,0.05,'>',/norm,charsize=1.2,col=133

if n_elements(fires_lon) ne 0 then begin
    oplot,fires_lon,fires_lat,psym=5,symsize=0.8,color=133
    
endif

xyouts,0.05,0.95,/norm,tit,color=133,charsize=1.2

device,/close
set_plot,'x'
ps2,outfile+'.ps','jpg',/overwrite,/printn

endif
end
;==================================================================
;main program
;==================================================================
@read_GASP_sca_loose.pro
@pattern_main_new4.pro
@convert_julianday.pro
@get_hysplitgrid.pro
@calc_hysplit_latlon.pro
@write_binary_f77.pro
@read_binary.pro
@write_binary_all.pro
@write_binary.pro
@remap_GASPsmoke_NAM.pro
@conv_smoke_combaod_combconc.pro
@generate_combsmoke.pro
@read_HMS_fire_west_2.pro
@merge_HMS_fires_west_mv_AK_HI_2.pro
@imagemap_unscl.pro
@colorbar.pro
@ps2.pro
@make_png.pro
@plot_check.pro
@disp_smoke_nocld_opr.pro
@ind2xy.pro

lon_infile='lon.dat'
lat_infile='lat.dat'

dir_aod='/data/GASP/GASP-WEST/US/aod/'
dir_fire='/data/new_SMOKE/arrH/'
;dir_out='../output/'
dir_out='/data/new_SMOKE_WEST/SMOKE_HAWAII/output/'
dir_temp=dir_out+'NAM/'
dir_grib=dir_out+'grib/'

yr0=[strmid(systime(0),strlen(systime(0))-4,4)]
date0=[fix(systime(/jul)-julday(1,0,yr0)-0.5)]
;date0=[fix(systime(/jul)-julday(1,0,yr0)-309.5)]
;date0=[160]
if date0(0) le 0 then begin
   yr0(0)=string(fix(yr0(0))-1,format='(i4.4)')
   if (fix(yr0(0)) mod 4 eq 0 and fix(yr0(0)) mod 100 ne 0) or (fix(yr0(0)) mod 400 eq 0) then date0(0)=[366] else date0(0)=[365]
endif

ncep=1
if ncep eq 1 then begin
  grib=1   ;0--do not generate GRIB files    1--generate GRIB files
  skip_NAM3=0  ;1---skip the NAM3 grid convert, 0--NAM3 convert and GRIB is grib=1
endif else begin
  grib=0   ;0--do not generate GRIB files    1--generate GRIB files
  skip_NAM3=1  ;1---skip the NAM3 grid convert, 0--NAM3 convert and GRIB is grib=1
endelse

process=1   ; 1--new process   0--read from smoke files already generated
hysplit_data=1
;check=1
check=0    ; 1--some detailed images will be generated and it will take longer to run, no need in operational
;plot_img=1  ;plot AOD and fire image
plot_img=0  
limit=[18,-162,23,-154 ]
smoke_plot=1 ;1-plot smoke image

nx_g=2500L
ny_g=912L

stime=0
etime=2400
scale=0.74   ;biomass-burning AOD ratio

;pref='G11.HI'
pref='GWHI'
instr='i14'

screen=0
maxang='170'
minang='70'
if screen eq 1 then cha='.screen'+minang+'-'+maxang else cha=''

aodstd=0.5

;data_suffix='.2smoke'+cha+'.mv.new4.nu.aodstd'+string(aodstd,format='(f3.1)')
;data_suffix='.combaod2'
;data_suffix='.2smoke.combaod.hmshysplitcomb2.new'
data_suffix='.2smoke.combaod.hmshysplitcomb2'

;get grid information
get_hysplitgrid,nlat,nlon,dlat,dlon,clat,clon,loc='HI'
calc_hysplit_latlon,nlon,nlat,clat,clon,dlat,dlon,newlon,newlat
read_GASP_latlon,lat_infile,lon_infile,nx_g,ny_g,glat,glon

for idate=0,n_elements(date0)-1 do begin
  yr=string(yr0(idate),format='(i4.4)')
  jday=string(date0(idate),format='(i3.3)')
  convert_julianday,yr,jday,mon,day
  month=string(mon,format='(i2.2)')
  nday=string(day,format='(i2.2)')

  dir=dir_aod

  if jday eq '001' then begin
    yr_last=string(fix(yr0(idate))-1,format='(i4.4)')
    month_last='12'
    nday_last='31'
    jday_last=string(julday(month_last,nday_last,long(yr_last))-julday(1,1,long(yr_last))+1,format='(i3.3)')
    dir_last=dir_out+yr_last+month_last+nday_last+'/'
    file_last=file_search(dir_aod+yr_last+jday_last+'*'+instr+'*US.all.aod.gz',count=nfll)
  endif else begin
    jday_last=string(jday-1,format='(i3.3)')
    convert_julianday,yr,jday_last,mon,day
    yr_last=yr
    month_last=string(mon,format='(i2.2)')
    nday_last=string(day,format='(i2.2)')
    file_last=file_search(dir_aod+yr+jday_last+'*'+instr+'_US.all.aod.gz',count=nfll)

    dir_last=dir_out+yr_last+month_last+nday_last+'/'

  endelse

  fire_infile=file_search(dir_fire+'hmshysplit'+yr+month+nday+'.txt',count=nfire)
  fire_infile_last=file_search(dir_fire+'hmshysplit'+yr_last+month_last+nday_last+'.txt',count=nfire_last)

  if nfire_last ne 0 then begin
    merge_HMS_fires_west_mv_AK_HI_2,fire_infile_last,lat_out_last_hi,lon_out_last_hi,ilat_out_last_hi,ilon_out_last_hi,/yester,dur_out=dur_out_hi,region='HI'
  endif
  

  cd,dir,current=cdir

  file0=file_search(yr+jday+'*'+instr+'*all.aod.gz',count=nfl)
  cd,cdir

  if nfl eq 0 then $
    print,'no GASP data on day '+jday+' of the year '+yr $
   else begin
    jtime=24
    sum0_smokeaod=fltarr(nlon,nlat,jtime)+0.
    num0_smokeaod=intarr(nlon,nlat,jtime)+0L
    num0_cloud=intarr(nlon,nlat,jtime)+0L
    num0_nonsmoke=intarr(nlon,nlat,jtime)+0L
    num0_nofire=intarr(nlon,nlat,jtime)+0L
    num0_nodata=intarr(nlon,nlat,jtime)+0L
    
    for ifl=0,nfl-1 do begin
    ;for ifl=0,1 do begin
      goes_infile=file0(ifl)
      pos=strpos(goes_infile,'_'+instr+'_US.all.aod.gz')
      time=strmid(goes_infile,pos-6,4)
      hr=strmid(time,0,2)
      ihr=fix(hr)
      date=yr+month+nday+time
      
      if time ge stime and time le etime then begin

        cd,dir_out,curr=cdir
        result=file_search(yr+month+nday,/test_direc)
        if result eq '' then spawn,"mkdir "+yr+month+nday
        outdir=dir_out+yr+month+nday+'/'
        cd,cdir

        outfl_combaod=outdir+pref+'.'+date+data_suffix+'.dat'
        
        aod=read_GASP_sca_loose(dir+goes_infile,nx_g=nx_g,ny_g=ny_g,bcld=bcld,bch1=bch1,bsca=sca_ang,caodstd=aodstd)

        if screen eq 1 then begin
           if maxang ne '' then $
             ind=where((sca_ang le minang and sca_ang gt -998.0) or sca_ang ge maxang,cp) $
           else ind=where((sca_ang le minang and sca_ang gt -998.0),cp)

           if cp gt 0 then aod(ind)=-9999.0
        endif

        gind=where(glat gt -90. and glon gt -180.)
        lat1=min(glat(gind))
        lat2=max(glat(gind))
        lon1=min(glon(gind))
        lon2=max(glon(gind))
        fires_lat=fltarr(1)-9999.0
        fires_lon=fltarr(1)-9999.0
        if gind(0) eq -1 then begin
          print, 'GOES geoloc missing*****'

        endif else begin

         if process eq 1 then begin
            
          if nfire_last eq 0 then begin
             print,'yesterday--no HMSHYSPLIT file found on '+yr_last+month_last+nday_last
        
          endif else begin
             if  ilat_out_last_hi(0) eq -1 then begin
               print,'yesterday--no fire found in HI on '+yr_last+month_last+nday_last
        
             endif else begin
               a3=where(ilat_out_last_hi eq -3,nnfire)
               fires_lat=fltarr(nnfire)-9999.0
               fires_lon=fltarr(nnfire)-9999.0
               num_fires=0

               a3=where(dur_out_hi eq 2400,nn3)
               if nn3 gt 0 then begin
                  fires_lat[num_fires:num_fires+nn3-1]=lat_out_last_hi(a3)
                  fires_lon[num_fires:num_fires+nn3-1]=lon_out_last_hi(a3)
                  num_fires=num_fires+nn3 
               endif

               a3=where(dur_out_hi lt 2400 and dur_out_hi gt 0,nn3)
               if nn3 gt 0 then begin
                
                 ;;include all HI fires found yesterday
                  ;fires_lat[num_fires:num_fires+nn3-1]=lat_out_last_hi(a3)
                  ;fires_lon[num_fires:num_fires+nn3-1]=lon_out_last_hi(a3)
                  ;num_fires=num_fires+nn3
 
                 if nfire eq 0L then begin
                   print,'today---no HMS file found on '+yr+month+nday             
                 endif else begin

                   merge_HMS_fires_west_mv_AK_HI_2,fire_infile,lat_out_hi,$
                         lon_out_hi,ilat_out_hi,ilon_out_hi,$
                         time0=time,region='HI'
                            
                   if nn3 gt 0 then begin
                       ind2=where(ilat_out_hi eq -3,cp2)
                       if cp2 gt 0 then begin
                         for ihi=0,nn3-1 do begin
                           dist=abs(lat_out_hi(ind2)-lat_out_last_hi(a3(ihi)))+abs(lon_out_hi(ind2)-lon_out_last_hi(a3(ihi)))
                           mind=min(dist,locc)
                           if mind lt 0.1 then begin
                              fires_lat(num_fires)=lat_out_hi(ind2(locc))
                              fires_lon(num_fires)=lon_out_hi(ind2(locc))
                              num_fires=num_fires+1
                           endif  ;if min(dist) l
                         endfor  ;for ihi=0,nn3-1
                       endif  ;if cp2 gt 0 
                   endif  ;if nn3 gt 0
                   
                 endelse  ;if nfire eq 0L
               endif  ;if nn3 gt 0

               if num_fires gt 0 then begin
                  fires_lat=fires_lat(0:num_fires-1) 
                  fires_lon=fires_lon(0:num_fires-1) 
               endif
             endelse ;if ilat_out_last(0) eq -1
          endelse  ;if nfire_last eq 0L
  
          paod=fltarr(nx_g,ny_g)-9999.0

          if fix(time) gt 1000 or nfll eq 0 then allfile=dir+file0 else allfile=[file_last,dir+file0]

          pattern_main_new4,allfile,instr,dir_out,yr,jday,time,yr_last,month_last,nday_last,dir_last,$
                            nx_g,ny_g,paod,glat,glon,data_suffix,target_xs=4,target_ys=4,screen,maxang=maxang,minang=minang

          if plot_img eq 1 then begin
              bartit='Aerosol Optical Depth'
              tit='GOES-11 Observation ('+yr+' '+month+' '+nday+$
                 '  '+time+'Z)'
              outfile='../output/img/AOD/'+pref+data_suffix+'.'+date
              ind=where(glat gt limit(0)-1 and glat lt limit(2)+1 and glon gt limit(1)-1 and glon lt limit(3)+1)
              if ind(0) eq -1 then print, 'GOES geoloc missing*****' else begin

                   yy=ind/nx_g
                   xx=ind-ind/nx_g*nx_g
                   iy1=min(yy)
                   iy2=max(yy)
                   ix1=min(xx)
                   ix2=max(xx)

                   aod_=aod(ix1:ix2,iy1:iy2)
                   glat_=glat(ix1:ix2,iy1:iy2)
                   glon_=glon(ix1:ix2,iy1:iy2)
                   bch1_=bch1(ix1:ix2,iy1:iy2)
                   bcld_=bcld(ix1:ix2,iy1:iy2)

                   xmag=4
                   ymag=4
                   sz0=size(aod_,/dim)
                   nx=sz0(0)
                   ny=sz0(1)
                   img_=congrid(aod_,xmag*nx,ymag*ny)
                   ch1_=congrid(bch1_,xmag*nx,ymag*ny, /interp)
                   cld_=congrid(bcld_,xmag*nx,ymag*ny)
                   lat_=congrid(glat_,xmag*nx,ymag*ny, /interp)
                   lon_=congrid(glon_,xmag*nx,ymag*ny, /interp)
                   ind00=where(glat_ gt -90 and glon_ gt -180.,ncomp=nnull)
                   ind=where(lat_ ge min(glat_(ind00)) and lat_ le max(glat_(ind00)) and lon_ ge min(glon_(ind00)) and lon_ le max(glon_(ind00)))

                   print,max(fires_lat)
                   ;stop
                   if max(fires_lat) gt -998.0 then disp_aircon,ch1_(ind),lat_(ind),lon_(ind),outfile=outfile,data=img_(ind),title=tit,limit=limit,$
                      fires_lat=fires_lat,fires_lon=fires_lon $
                   else disp_aircon,ch1_(ind),lat_(ind),lon_(ind),outfile=outfile,data=img_(ind),title=tit,limit=limit

                 
                 endelse
          endif

          if max(fires_lat) lt -998.0 then begin
        
            print,'no fires detected.......'
                
            if max(paod) lt 0. then $
              combaod=generate_combsmoke(aod,bcld,nfires=0,scale=scale,ismoke=1) $
            else begin
              find_pattern_boundary,paod,out_ind

              if out_ind(0) ne -1 then begin
                smokeaod=fltarr(nx_g,ny_g)-9999.0
                smokeaod(out_ind)=paod(out_ind)
                combaod=generate_combsmoke(aod,bcld,nfires=-1,smokeaod=smokeaod,$
                                         scale=scale,ismoke=1)
              endif
            endelse

          endif else begin

            nfires=n_elements(fires_lat)
            plume_ind=[-1]
            smokeaod=[-9999.0]
        
            bulkaod=fltarr(n_elements(fires_lat))-99.0
            nbulkaod=0
        
            for ifire=0,n_elements(fires_lat)-1 do begin
              lat0=fires_lat(ifire)
              lon0=fires_lon(ifire)

              print,ifire+1,'  of',nfires,lat0,lon0

              if lat0 gt lat1 and lat0 lt lat2 and lon0 gt lon1 and lon0 lt lon2 then begin
                az=fltarr(nx_g,ny_g)-9999.0  ;azimuth angle from the north
                distance=fltarr(nx_g,ny_g)-9999.0
                ddlat =glat - lat0
                ddlon = glon - lon0
                diff = abs(ddlat)+abs(ddlon)
                del = min(diff,minloc)
                ind2xy,minloc,glat,x0,y0
                if (del gt 1.0 ) then begin
                 mini=-1
                 minj=-1
               endif else begin
                 mini=x0
                 minj=y0
               endelse

               if mini ne -1 and minj ne -1 then begin
                maxcir=10
                step=1
                maxaod=fltarr(maxcir)-99.0
                maxdist=fltarr(maxcir)-99.0
                maxratio=fltarr(maxcir)+0.
                maxaz=fltarr(maxcir)-99.

                if plume_ind(0) eq -1 and n_elements(plume_ind) eq 1 then plume_ind=[mini+minj*nx_g]     ;Jian Zeng, July 12, 2007
                maxaod(0)=aod(mini,minj)
        
                if maxaod(0) gt 0. then maxratio(0)=100.  ;change from -9998.0 to 0. Jan15  2007
                maxdist(0)=0.
                i1=mini
                i2=mini
                j1=minj
                j2=minj
        
                naod=1L
                for iaod=1,maxcir-1,step do begin
                  if i1 eq 0 and i2 eq nx_g-1 and j1 eq 0 and $
                    j2 eq ny_g-1 then goto,skipping $
                  else begin
                    i1=(mini-iaod) > 0
                    i2=(mini+iaod) < (nx_g-1)
                    j1=(minj-iaod) > 0
                    j2=(minj+iaod) < (ny_g-1)
                    if i2-i1 ge 2 then begin
                    dump_aod=[reform(aod(i1,j1:j2)),reform(aod(i2,j1:j2)),$
                              reform(aod(i1+1:i2-1,j1)),reform(aod(i1+1:i2-1,j2))]
                    dump_lat=[reform(glat(i1,j1:j2)),reform(glat(i2,j1:j2)),$
                              reform(glat(i1+1:i2-1,j1)),reform(glat(i1+1:i2-1,j2))]
                    dump_lon=[reform(glon(i1,j1:j2)),reform(glon(i2,j1:j2)),$
                              reform(glon(i1+1:i2-1,j1)),reform(glon(i1+1:i2-1,j2))]

                    nn=n_elements(dump_aod)

                    if nn gt 0L then begin
                      distance=fltarr(nn)-99.0
                      az=fltarr(nn)-99.0
                      for idump=0,nn-1 do begin
                        id=idump
                        if dump_lat(id) gt -90. and dump_lon(id) gt -180 then begin
                          dist0=map_2points(dump_lon(id),dump_lat(id),$
                                lon0,lat0,/meters)
                          aa=map_2points(dump_lon(id),dump_lat(id),lon0,lat0)
                          distance(idump)=dist0/1000.
                          az(idump)=aa(1)
                          if aa(1) lt 0 then az(idump)=360+aa(1)
                        endif

                      endfor

                      if max(dump_aod) gt 0. then begin
                        maxv=max(dump_aod,location)
                        maxaod(naod)=maxv
                        result=where(abs(dump_aod - maxv) lt 1.e-5,cp)
                        if cp eq 1L then loca=array_indices(dump_aod,location) $
                        else begin
                          if iaod-step le 0 then begin
                            loca=array_indices(dump_aod,location)
                          endif else begin
                            min_az=999.0
                            for iaz=0,cp-1 do begin
                              dist_az=abs(az(result(iaz)) - maxaz(iaod-step))
                              if dist_az lt min_az then begin
                                min_az=dist_az
                                loca=result(iaz)
                              endif
                            endfor
                          endelse
                        endelse

                        maxdist(naod)=distance(loca)
                        maxaz(naod)=az(loca)
                        ind0=where(dump_aod gt 0.04,cp)
                        maxratio(naod)=float(cp)/float(nn)*100.

                      endif else begin ;if max(dump_aod) gt 0.
                        maxaod(naod)=-99.
                        maxdist(naod)=max(distance)
                        maxaz(naod)=-99.0
                        maxratio(naod)=0.
                        print,'******',iaod,maxaod(naod),maxdist(naod),$
                            maxaz(naod),maxratio(naod)

                      endelse
                      
                      naod=naod+1
        
                    endif ;if nn gt 0L
        
                    endif ;added on Oct. 30 2007
                  endelse  ;if i1 eq 0 and i2 eq nx_g-1 ... else
                endfor ;for iaod=1,maxcir-1,step

                skipping:

                find_bulk_aod,maxaod(0:naod-1),maxaz(0:naod-1),$
                             maxratio(0:naod-1),maxdist(0:naod-1),$
                              bulk_aod,wind_direc
                print,'bulk_aod, wind_direction:',bulk_aod,wind_direc

                if bulk_aod lt 0.05 then begin
                  print,'This fire is possibly under cloud or a false fire'
                  i1=(mini-100) > 0
                  i2=(mini+100) < (nx_g-1)
                  j1=(minj-100) > 0
                  j2=(minj+100) < (ny_g-1)
                  paod_part=paod(i1:i2,j1:j2)
        
                  bulk=max(paod_part)
                  if bulk gt 0 then begin
                    find_pattern_boundary,paod,out_ind,nbulkaod=1,bulkaod=bulk
                    if out_ind(0) ne -1 then begin
                      result=doublefind(plume_ind,out_ind)
                      plume_ind=result
                    endif
                  endif
                  ;stop
                endif else begin
        
                  bulkaod(nbulkaod)=bulk_aod
                  nbulkaod=nbulkaod+1

                  nlev=1
                  if bulk_aod gt 0.5 then percent=[0.2] else percent=[0.4]
        
                  for ilev=0,nlev-1 do begin
                    
                    lev=bulk_aod*percent(ilev)
                    find_plume,lev,mini,minj,aod,ind_out
                    help,'in---',ind_out,plume_ind

                    result=doublefind(plume_ind,ind_out)
                    plume_ind=result
                    help,'out---',plume_ind
                    print,aod(mini,minj)

                    ;combine with paod
                    if paod(mini,minj) gt 0. then $
                      find_plume,lev,mini,minj,paod,ind_out $
                    else begin
                     if max(paod) gt 0. then begin
                      lim=min([10,mini,minj])  
                      il=mini
                      jl=minj
                      max_aod=0.
                      for ilim=mini-lim,mini+lim do begin
                      for jlim=minj-lim,minj+lim do begin
                        iilim=ilim > 0 < (nx_g-1)
                        jjlim=jlim > 0 < (ny_g-1)
                        if paod(iilim,jjlim) gt max_aod then begin
                          il=iilim
                          jl=jjlim
                          max_aod=paod(iilim,jjlim)
                        endif
                      endfor
                      endfor

                      if paod(il,jl) gt 0. then $
                       find_plume,lev,il,jl,paod,ind_out
                     endif
                    endelse
        
                    result=doublefind(plume_ind,ind_out)
                    plume_ind=result
        
                    if check eq 1 then begin
                       plot_check,yr,month,nday,time,mini,minj,aod,$
                               glat,glon,nx_g,ny_g,bulk_aod,lat0,lon0,$
                               plume_ind,maxaod,maxaz,maxratio,maxdist,maxcir
                    endif

                  endfor ;for ilev=0,nlev-1
        
                endelse ;if bulk_aod lt 0.05
                endif else print,'fire out of boundary....' ;if mini ne -1
              endif ;if lat0 gt 20 and lat0 lt 50 and
            endfor ;for ifire=0,n_elements(fires_lat)-1
        
            ;third check
            find_pattern_boundary,paod,out_ind,nbulkaod=nbulkaod,bulkaod=bulkaod
            if out_ind(0) ne -1 then begin
               result=doublefind(plume_ind,out_ind)
               plume_ind=result
            endif
        
            if plume_ind(0) ne -1 then begin
              aod_plume=fltarr(nx_g,ny_g)-9999.0

              aod_plume(plume_ind)=aod(plume_ind)

              ;generate smoke AOD combined with cloud screening
              combaod=generate_combsmoke(aod,bcld,nfires=nfires,$
                                         smokeaod=aod_plume,$
                                         scale=scale,ismoke=1)
            endif else begin
              print,'no smoke detected.......'

              combaod=generate_combsmoke(aod,bcld,nfires=0,scale=scale,ismoke=1,$
                                         smokeaod=smokeaod)
            endelse
        
            
          endelse ;if min(fires_lat) lt -998.0
          
          write_binary,outfl_combaod,combaod

         endif else begin ;if process eq 1
        
          combaod=fltarr(nx_g,ny_g)-9999.0
          print,outfl_combaod
          waiton=1
          while waiton do begin
            okfile=file_search(outfl_combaod,count=nok)
            if  nok eq 1 then begin
              waiton=0
              read_binary,outfl_combaod,combaod
            endif else begin
              print,'Data not ready yet, wait for 10 minutes....'
              wait,600
        
            endelse
          endwhile
        
        endelse ;if process eq 1

        if smoke_plot eq 1 then begin
           disp_smoke_nocld_opr,dir_out,glat,glon,combaod,yr,month,nday,time,limit=limit,data_suffix=data_suffix
         endif
         if skip_NAM3 eq 0 then begin
        
          sum_smokeaod=fltarr(nlon,nlat)+0.
          num_smokeaod=intarr(nlon,nlat)+0L
          num_cloud=intarr(nlon,nlat)+0L
          num_nonsmoke=intarr(nlon,nlat)+0L
          num_nofire=intarr(nlon,nlat)+0L
          num_nodata=intarr(nlon,nlat)+0L
        
          for i=0,nx_g-1 do begin
          for j=0,ny_g-1 do begin
            if glon(i,j) ge -180. and glat(i,j) ge -90. then begin
              ii=fix((glon(i,j)-clon+0.5*dlon)/dlon)
              jj=fix((glat(i,j)-clat+0.5*dlat)/dlat)
              if ii ge 0 and jj ge 0 and ii lt nlon and jj lt nlat then begin
                if combaod(i,j) lt -9998.0 then begin
                  num_nodata(ii,jj)=num_nodata(ii,jj)+1L
                  num0_nodata(ii,jj,ihr)=num0_nodata(ii,jj,ihr)+1L
                endif else begin
                  if combaod(i,j) eq -1. then begin
                    num_cloud(ii,jj)=num_cloud(ii,jj)+1L
                    num0_cloud(ii,jj,ihr)=num0_cloud(ii,jj,ihr)+1L
                  endif else begin
                   if combaod(i,j) eq -2. then begin
                     num_nofire(ii,jj)=num_nofire(ii,jj)+1L
                     num0_nofire(ii,jj,ihr)=num0_nofire(ii,jj,ihr)+1L
                   endif else begin
                     if combaod(i,j) eq -3. then begin
                       num_nonsmoke(ii,jj)=num_nonsmoke(ii,jj)+1L
                       num0_nonsmoke(ii,jj,ihr)=num0_nonsmoke(ii,jj,ihr)+1L
                     endif else begin
                       if combaod(i,j) ge -0.5 and combaod(i,j) le 2.05 then begin
                         sum_smokeaod(ii,jj)=sum_smokeaod(ii,jj)+$
                                                 combaod(i,j)
                         num_smokeaod(ii,jj)=num_smokeaod(ii,jj)+1L
                         sum0_smokeaod(ii,jj,ihr)=sum0_smokeaod(ii,jj,ihr)+$
                                                 combaod(i,j)
                         num0_smokeaod(ii,jj,ihr)=num0_smokeaod(ii,jj,ihr)+1L
                       endif
                     endelse
                   endelse
                  endelse
                endelse
              endif
            endif
          endfor
          endfor

         endif ; if skip_NAM3 eq 0
        

        endelse ;if ind(0) eq -1 then print, 'GOES geoloc missing*****'
        
      endif ;if time ge stime

    endfor ;for ifl=0,nfl-1
;print,grib,skip_NAM3
;stop

    if skip_NAM3 eq 0 then begin
      for iiii=0,23 do print,iiii,max(num0_cloud(*,*,iiii)),max(num0_nodata(*,*,iiii)),max(num0_smokeaod(*,*,iiii))
;stop
      convert_hourly,jtime,dir_temp,dir_grib,pref,nlon,nlat,yr,month,nday,dlon,$
                     dlat,clon,clat,newlon,newlat,glon=glon,glat=glat,$
                     num0_smokeaod=num0_smokeaod,data_suffix=data_suffix,$
                     sum0_smokeaod=sum0_smokeaod,hysplit_data=hysplit_data,$
                     num0_cloud=num0_cloud,num0_nonsmoke=num0_nonsmoke,$
                     num0_nofire=num0_nofire,num0_nodata=num0_nodata,grib=grib
       
    endif
   endelse ;if nfl eq 0
   
endfor ;for idate=0,n_elements(date0)-1


end