;============================================================================ ;g13_smoke_main5.pro ;written by Jian Zeng (ERT) ;version 5 ;created:2006-12-22 ; ;20% for bulk_aod > 0.5 and 40% for bulk_AOD lt 0.5 version 2 ;variable threshold ; version 3 ;neglect gap-filled and back to 20-40% threhsold, $ ;change find_winddirection and bulk AOD version4 2006-7-7 ; ;2006-09-21: ; Change the maximum cirle from 50 in step of 1 to 100 in step of 2 ; Add 'std' to denote different screenings ;2006-09-28: ; change the NAM grid size to 801*534, lower_left=(-175,0) ; covert F77 data from 30-minutes files to hourly data files ;2006-12-05: ; change aodstd into variables. ; First only use January empirical formula ; Second 12 formulas for 12 months ; Also with yearly mean formula ;2006-12-22: ; combine source scheme---PACKAGE/generate_goessmoke_opr5_case.pro ; and pattern scheme---../Pattern/pattern_main.pro ;2007-01-12: ; use HMS fire product and predict the smoke movement for the second day ; and compare to HYSPLIT forecast. ;2007-02-15: ; check how long the previous day's fires last today. Compare to HYSPLIT ; 24hr forecast ;2007-03-12: ; version 5: keep the fires of two consecutive days at the same location. ;----------------------------------------------------------------------------- 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 end ;---------------------------------------- pro find_sqr,i0,j0,outi,outj,lev,img nxny_file='/data/GASP/US/ref/USnxny' openr,1,nxny_file readf,1,outnx,outny nx_g = outnx ny_g = outny close,1 ;print,'in find_sqr',nx_g,ny_g 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-2) > 0 < 2128 ; outj(num)=(jj+j0-2) > 0 < 880 outi(num)=(ii+i0-2) > 0 < nx_g outj(num)=(jj+j0-2) > 0 < ny_g 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) result=array_indices(img,ind0) i0=result(0) j0=result(1) s1=(i0-2) > 0 s2=(i0+2) < (sz0(0)-1) e1=j0-2 > 0 e2=j0+2 < (sz0(1)-1) temp=img(s1:s2,e1:e2) find_sqr,i0,j0,outi,outj,lev,temp 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,pref,nlon,nlat,yr,month,nday,lon,lat,$ 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=11 dir_ftp='/net/www/aftp/pub/smcd/spb/shobha/GOES_AOD/SMOKE5/' for kk=shr,jtime-1 do 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 outfl_grib_NAM=dir_temp+pref+'.'+yr+month+nday+$ string(kk,format='(i2.2)')+data_suffix+'.NAM3.grd' date=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(date,2,10) close,1 spawn,'/home/oper/new_SMOKE/PACKAGE5/grib_main_NAM3' spawn,'rm -f infile0' ; spawn,'rm -f '+outfl_combaod_NAM print,'created '+outfl_grib_NAM ;print,'move to ---',dir_ftp2 ; file_move,outfl_grib_NAM,dir_ftp,/overwrite endif ;if grib eq 1 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 endif else thresh=max(paod)*0.05 out_ind=where(paod ge thresh) endif end ;================================================================== ;main program ;================================================================== @read_GASP.pro @read_GASP_latlon.pro @pattern_main.pro @convert_julianday.pro @get_hysplitgrid3.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 @merge_HMS_fires3.pro @disp_smoke_nocld.pro @imagemap_unscl.pro @colorbar.pro @make_png.pro @ps2.pro @read_HMS_fire3.pro ;lat_infile='lat.dat' ;lon_infile='lon.dat' lat_infile='/data/GASP/US/ref/lat.dat' lon_infile='/data/GASP/US/ref/lon.dat' ;dir_fire='/net/orbit111l/disk2/pub/zeng/HMS/DATA/FIRE/HMS_HYSPLIT/hmshysplit/' ;dir_aod='/net/orbit111l/disk3/pub/cxu/GOES/SAVE_AOD/' ;dir_out='OUTPUT/' ;dir_ftp2='/net/www/aftp/pub/smcd/spb/shobha/GOES_AOD/' ;dir_temp='NAM/' ;dir_src='/home/oper/SMOKE5' dir_src='/home/oper/new_SMOKE/PACKAGE5' ;dir_fire='/home/oper/SMOKE5/arrH/' dir_fire='/data/new_SMOKE/arrH/' dir_data='/data/new_SMOKE' ;dir_fire=dir_data+'/SMOKE_input/arrH/' ;dir_aod=dir_data+'/SMOKE_input/aod/' ;dir_aod='/data/GASP/US/aod/' ;dir_aod='/home/oper/SMOKE5/AOD/' dir_aod='/data/new_SMOKE/AOD/' dir_out=dir_data+'/Thresh20_40/COMB_AOD_CLOUD_FIRE/2SCHEMES/' dir_temp=dir_data+'/Thresh20_40/COMB_AOD_CLOUD_FIRE/NAM/' dir_ftp2=dir_data+'/GOES_AOD/' dir_Image=dir_data+'/Image/' ;dir_out='/data/new_SMOKE/Thresh20_40/COMB_AOD_CLOUD_FIRE/2SCHEMES/' ;dir_ftp2='/data/new_SMOKE/GOES_AOD/' ;dir_temp='/data/new_SMOKE/Thresh20_40/COMB_AOD_CLOUD_FIRE/NAM/' 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)-1.5)] 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 ;date0=[fix(systime(/jul)-julday(1,0,yr0)-1.5)] ;yr0=['2009'] ;date0=[163] ;date0=[365] grib=1 ;0--do not generate GRIB files 1--generate GRIB files all=1 skip_NAM3=0 ;1---skip the NAM3 grid convert, 0--NAM3 convert and GRIB is grib=1 process=1 ; 1--new process 0--read from smoke files already generated hysplit_data=0 ;0---do not generate HYSPLIT binary data, 1--generate HYSPLIT binary data img_plot=1 ;1--plot smoke image 0-no smoke image plotting ;nx_g=2000L ;ny_g=850L nx_g=2128L ny_g=880L ;stime=1415 ;stime=1045 stime=1015 etime=2245 ;etime=1315 pref='G13' scale=0.74 ;biomass-burning AOD ratio std='' ;'.newstd'---for new aodstd ;get grid information get_hysplitgrid3,nlat,nlon,dlat,dlon,clat,clon calc_hysplit_latlon,nlon,nlat,clat,clon,dlat,dlon,xl,xr,yl,yu,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)') ; if jday eq '000' then begin ; yr=string(fix(yr0(idate))-1,format='(i4.4)') ; month='12' ; nday='31' ;; date0=[fix(systime(/jul)-julday(1,0,yr)-19.5)] ; date0=[fix(systime(/jul)-julday(1,0,yr)-0.5)] ; jday=string(date0(idate),format='(i3.3)') ; print,'test3', yr, month, nday, date0, jday ; endif if jday eq '001' then begin yr_last=string(fix(yr0(idate))-1,format='(i4.4)') month_last='12' nday_last='31' endif else begin yr_last=yr convert_julianday,yr_last,jday-1,mon_last,day_last month_last=string(mon_last,format='(i2.2)') nday_last=string(day_last,format='(i2.2)') endelse ; dir=dir_aod+yr+'_'+jday+'/' dir=dir_aod ; fire_infile=file_search(dir_fire+yr+'/hmshysplit'+yr+month+nday+'.txt',count=nfire) ; fire_infile_last=file_search(dir_fire+yr_last+'/hmshysplit'+yr_last+month_last+nday_last+'.txt',count=nfire_last) 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 $ merge_HMS_fires3,fire_infile_last,lat_out_last,lon_out_last,ilat_out_last,ilon_out_last,/yester print,'fire file '+fire_infile print,'last fire file'+fire_infile_last print,'yr mon dd'+dir_fire+yr+'/hmshysplit'+yr+month+nday+'.txt' cd,dir,current=cdir infile0=file_search(dir_aod+yr+jday+'*all.aod.gz',count=nfl) print,'infile all '+infile0 cd,cdir if nfl eq 0 then $ print,'no GASP data on day '+jday+' of the year '+yr $ else begin jtime=23 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 goes_infile=infile0(ifl) pos=strpos(goes_infile,'_i18_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 print,goes_infile outfl_allaod=dir_temp+pref+'.'+date+'.all.aod_conc.NAM3'+std+'.dat' 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 data_suffix='.2smoke.combaod'+std+'.hmshysplitcomb2' outfl_combaod=outdir+pref+'.'+date+data_suffix+'.dat' print,'goes_infile'+dir+goes_infile ; aod=read_GASP(dir+goes_infile,nx_g=nx_g,ny_g=ny_g,bcld=bcld,bch1=bch1) aod=read_GASP(goes_infile,nx_g=nx_g,ny_g=ny_g,bcld=bcld,bch1=bch1) 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 grib eq 1 and all eq 1 then begin aod_new=generate_combsmoke(aod,bcld) remap_GASPsmoke_NAM,nlat,nlon,dlat=dlat,dlon=dlon,clat=clat,clon=clon,$ glat=glat,glon=glon,$ indata=aod_new,outdata=aod_NAM3 allconc=conv_smoke_combaod_combconc(aod_NAM3) write_binary_all,outfl_allaod,nlon,nlat,newlon,newlat,$ aod_NAM3,reform(allconc(*,*,0)) print,'outfl_allaod'+outfl_allaod outfl_grib_NAM3=dir_temp+pref+'.'+date+'.all.aod_conc.NAM3.grd' print,'outfl_grib_NAM3'+dir_temp+pref+'.'+date+'.all.aod_conc.NAM3.grd' openw,1,'infile0' printf,1,outfl_allaod printf,1,outfl_grib_NAM3 printf,1,strmid(date,2,10) close,1 spawn, dir_src+'/grib_main_NAM3' ;spawn,'/home/oper/new_SMOKE5/PACKAGE5/grib_main_NAM3' spawn,'rm -f infile0' spawn,'rm -f '+outfl_allaod print,'created '+outfl_grib_NAM3 ;file_move,outfl_grib_NAM3,dir_ftp2+'ALL/',/overwrite endif ;if grib eq 1 if process eq 1 then begin if nfire_last eq 0 then begin print,'yesterday--no HMS file found on '+yr_last+month_last+nday_last endif else begin if ilat_out_last(0) eq -1 then begin print,'yesterday--not fire found on '+yr_last+month_last+nday_last endif else begin if nfire eq 0L then begin print,'today---no HMS file found on '+yr+month+nday endif else begin nnfire=n_elements(lat_out_last) fires_lat=fltarr(nnfire)-9999.0 fires_lon=fltarr(nnfire)-9999.0 merge_HMS_fires3,fire_infile,lat_out,lon_out,ilat_out,ilon_out,$ nst=nst,ndur=ndur,type=type,xgrid=xgrid,ygrid=ygrid,$ fires_lat=fires_lat_orig,fires_lon=fires_lon_orig,$ time0=time nnn=0L for innn=0,n_elements(ilat_out)-1 do begin indnnn=where(ilat_out_last eq ilat_out(innn) and $ ilon_out_last eq ilon_out(innn),cp) if cp gt 0 then begin fires_lat(nnn)=lat_out(innn) fires_lon(nnn)=lon_out(innn) nnn=nnn+1 endif endfor if nnn gt 0 then begin fires_lat=fires_lat(0:nnn-1) ; lat_out fires_lon=fires_lon(0:nnn-1) ; lon_out endif endelse ;if nfire eq 0L endelse ;if ilat_out_last(0) eq -1 endelse ;if nfire_last eq 0L paod=fltarr(nx_g,ny_g)-9999.0 print,'pattern_main'+dir pattern_main,dir,outdir,yr,jday,time,nx_g,ny_g,paod,glat,glon,data_suffix 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 mindist=99999.0 mini=-1 minj=-1 for i=0,nx_g-1 do begin for j=0,ny_g-1 do begin if glat(i,j) gt -90. and glon(i,j) gt -180. then begin dist0=map_2points(glon(i,j),glat(i,j),lon0,lat0,/meters) if dist0 lt mindist then begin mindist=dist0 mini=i minj=j endif endif endfor endfor if mini ne -1 and minj ne -1 then begin maxcir=100 step=2 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 -9998. then maxratio(0)=100. 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) 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(lon0,lat0,dump_lon(id),dump_lat(id)) 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 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 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 result=doublefind(plume_ind,ind_out) plume_ind=result ;combine with paod if paod(mini,minj) gt 0. then $ find_plume,lev,mini,minj,paod,ind_out $ else begin lim=min([100,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 endelse result=doublefind(plume_ind,ind_out) plume_ind=result endfor ;for ilev=0,nlev-1 endelse ;if bulk_aod lt 0. endif else print,'fire out of boundary....' ;if mini ne -1 endif ;if lat0 gt 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 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 img_plot eq 1 then begin ;file_mv=1, move smoke image to ftp webserver 0--no image moving ; disp_smoke_nocld,glat,glon,combaod,yr,month,nday,time,file_mv=0 disp_smoke_nocld,dir_Image,glat,glon,combaod,yr,month,nday,time,file_mv=0 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 if skip_NAM3 eq 0 then begin convert_hourly,jtime,dir_temp,pref,nlon,nlat,yr,month,nday,newlon,newlat,$ 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