;title:   acse_plotday.pro
;purpose: create a suite of data plots for operational ACSE measurements
;input:   datearg - string argument for the date to process "YYYYMMDD"
;output:  set of plots
;-----------------------------------------------------------------------

;pro acse_plotday,datearg
datearg='20140816'  ;The date for processing.  Note: To run as a procedure, comment out this line and uncomment the previous line.


;parameters
;-----------
archivedir='E:/archive'
wbanddir=archivedir+'/wband/mom'
motiondir=archivedir+'/wband/logs'
mwrdir=archivedir+'/mwr'
ceilodir=archivedir+'/ceilometer'
rwpdir=archivedir+'/windprofiler/QAWinds'
radiodir=archivedir+'/radiosonde'
shipdir=archivedir+'/shipdata'
hatprodir=archivedir+'/hatpro'
snakedir=archivedir+'/seasnake'
irtdir=archivedir+'/irt'
dopliddir=archivedir+'/dopplerlidar'
plotdir='C:/data/acse/plots'
colorfile='C:/data/acse/programs/colorscopy.tbl'
wfactor=50 ;factor for decreasing W-band data rate
cfactor=20 ;factor for decreasing W-band crossbow data rate
kfactor=10 ;factor for decreasing W-band kongsberg data rate
d_hfactor=4  ;factor for decreasing Doppler lidar height resolutoin
d_tfactor=10  ;factor for decreasing Doppler lidar time scale
;-----------
device,decomposed=0.
read_data=1
plot_data=1

;read data
;----------
;----------
if read_data then begin
sarg=datearg+'00'
earg=datearg+'24'
;WBAND
acse_readwband,sarg,earg,w_time,w_height,w_dbz,w_dopp,w_spwd,w_stn,directory=wbanddir,factor=wfactor
;WBAND Motion
acse_readcross,sarg,earg,cr_time,cr_pitch,cr_roll,cr_zacc,directory=motiondir,factor=cfactor
acse_readkongs,sarg,earg,ko_time,ko_pitch,ko_roll,ko_vel,directory=motiondir,factor=kfactor
;MWR
acse_readmwr,sarg,earg,m_time,m_lwp,m_pwv,tamb=m_tamb,rh=m_rh,pres=m_pres,tir=m_tir,cbase=m_cbase,height=m_height,tprof=m_tprof,qprof=m_qprof,directory=mwrdir
m_tamb=m_tamb-273.15
m_tir=m_tir-273.15
acse_readmwrbrt,sarg,earg,mb_time,directory=mwrdir,t22p234=mb_t22p234,t28p000=mb_t28p000,t52p280=mb_t52p280,t56p660=mb_t56p660
;Ceilo
acse_readceilo,sarg,earg,ctime,detstat,cbase1,cbase2,cbase3,vertvis,highsig,directory=ceilodir
;RWP
acse_readwind,sarg,earg,1,p_timehb,p_heighthb,p_speedhb,p_directionhb,directory=rwpdir,algorithm=1
acse_readwind,sarg,earg,0,p_timelb,p_heightlb,p_speedlb,p_directionlb,directory=rwpdir,algorithm=1
;Radiosonde
acse_readrds,sarg,earg,r_time,r_height,r_pres,r_temp,r_dewpt,r_rh,r_speed,r_direction,directory=radiodir,htlim=10.
numsonde=n_elements(r_time)
;Ship data
acse_readship,sarg,earg,o_time,o_lat,o_lon,o_head,o_cog,o_sog,o_wdir,o_wspd,o_tsea,o_ssea,o_tair,o_rhair,o_pres,directory=shipdir
;HATPRO
acse_readhatpro,sarg,earg,h_timei,h_lwp,h_pwv,timeh=h_timeh,height=h_height,tprof=h_tprof,qprof=h_qprof,directory=hatprodir
acse_readhatprobrt,sarg,earg,hb_time,directory=hatprodir,t22p24=hb_t22p24,t27p84=hb_t27p84,t52p28=hb_t52p28,t56p66=hb_t56p66
;Sea Snake
acse_readsnake,sarg,earg,s_time,s_temp,directory=snakedir
;irt
acse_readirt,sarg,earg,i_time1,i_irt1,i_time2,i_irt2,irtr1=i_irtr1,irtr2=i_irtr2,directory=irtdir
;Doppler lidar
acse_readdlidar,sarg,earg,d_time,d_height,d_back,d_dopp,d_intens,directory=dopliddir,jtime=d_jtime,hfactor=d_hfactor,tfactor=d_tfactor


endif ;read_data

;plot data
;----------
;----------
cd,plotdir
datestr=strmid(sarg,0,8)+' '+strmid(sarg,8,2)+':00 - '+strmid(earg,0,8)+' '+strmid(earg,8,2)+':00'
if plot_data then begin
loadct,44,file=colorfile
blk=strarr(10)+' '
th=1
cs=1.2
xra=[float(strmid(sarg,8,2)),float(strmid(earg,8,2))]
wyra=[0,6]
;Wband
;--------
cbase=interpol(cbase1,ctime,w_time)
wband_momentsplot,w_time,w_height,w_dbz,w_dopp,w_spwd,w_stn,datestr,XRA=xra,YRA=wyra,SITE='ACSE',ylog=ylog,cbase=cbase
ifile='acse.wband.'+sarg+'.'+earg+'.png'
tvlct,r,g,b,/get
write_png,ifile,tvrd(true=1),r,g,b

;Doppler lidar
;---------------
dopplerlidar_plot,d_time,d_height,d_back,d_dopp,d_intens,datestr,XRA=xra,SITE='ACSE'
ifile='acse.dopplerlidar.'+sarg+'.'+earg+'.png'
tvlct,r,g,b,/get
write_png,ifile,tvrd(true=1),r,g,b

loadct,44,file=colorfile
;Radar metaplot
;--------------
window,0,xsize=500,ysize=600
pos=[.12,.08,.87,.95]
pos1=[pos[0],pos[1]+(pos[3]-pos[1])*3./5,pos[2],pos[1]+(pos[3]-pos[1])*5./5]
pos2=[pos[0],pos[1]+(pos[3]-pos[1])*2./5,pos[2],pos[1]+(pos[3]-pos[1])*3./5]
pos3=[pos[0],pos[1]+(pos[3]-pos[1])*1./5,pos[2],pos[1]+(pos[3]-pos[1])*2./5]
pos4=[pos[0],pos[1]+(pos[3]-pos[1])*0./5,pos[2],pos[1]+(pos[3]-pos[1])*1./5]
;Vd
;Filter out data where S-to-N is lower than a threshold
;------------------------------------------------------------
pdopp=w_dopp & pdbz=w_dbz
threshold = -14.
iwh=where(w_stn le threshold)
if iwh[0] ne -1 then begin
   pdopp[iwh]=-9
   pdbz[iwh]=-99
endif
levels=[-10,findgen(25)/(24/5.4)-2.7]
nlev=n_elements(levels)
numc=256
colors=[0, fix(indgen(nlev-2)*(numc*.80/(nlev-3)))+fix(numc*.196), numc-1]
labels=strtrim(string(levels,format='(f4.1)'),2)
contour,pdopp,w_time,w_height,levels=levels,c_colors=colors,/cell_fill,position=pos1,/norm,font=-1,$
  /xstyle,/ystyle,yrange=wyra,xrange=xra,xtickname=blk,ytitle='Height AGL (km)',ylog=ylog,$
  title='ACSE, W-Band Meta, for '+datestr
colorbar_shupe,[pos1[2]+0.02,pos1[1]+0.02,pos1[2]+0.05,pos1[3]-0.05],levels,colors,$
  nlabels=6,title='[m/s]',labels=labels,charthick=th,color=1,charsize=1,laboffset=1
xyouts,pos1[0]+0.02,pos1[3]-0.04,/norm,'!6Doppler Velocity!5',font=-1,charsize=1.7,color=1,charthick=th
;vertical
plot,cr_time,cr_zacc,position=pos2,/noerase,xrange=xra,/xstyle,xtickname=blk,color=1,/nodata,ytitle='Vert (m/s)',$
    yrange=[-1,2],/ystyle
oplot,cr_time,cr_zacc,color=50
oplot,ko_time,ko_vel,color=254
xyouts,5,1.3,/data,'Zaccel (m/s!u2!n)',charsize=1.1,color=50
xyouts,5,0.4,/data,'Zvel (m/s)',charsize=1.1,color=254
;pitch
plot,cr_time,cr_pitch,position=pos3,/noerase,xrange=xra,/xstyle,xtickname=blk,color=1,/nodata,ytitle='Pitch (deg)',$
    yrange=[-1,1],/ystyle
oplot,cr_time,cr_pitch,color=50
oplot,ko_time,ko_pitch,color=254
;roll
plot,cr_time,cr_roll,position=pos4,/noerase,xrange=xra,/xstyle,xtitle='Time [UTC]',color=1,/nodata,ytitle='Roll (deg)',$
    yrange=[-2,2],/ystyle
oplot,cr_time,cr_roll,color=50
oplot,ko_time,ko_roll,color=254
xyouts,pos3[2]+0.03,pos3[1]-0.1,/norm,orient=90,'Antenna (Kongsberg)',color=254,charsize=1.1
xyouts,pos3[2]+0.03,pos3[3],/norm,orient=90,'Ship (Crossbow)',color=50,charsize=1.1
;write a png file
ifile='acse.wbandmeta.'+sarg+'.'+earg+'.png'
tvlct,r,g,b,/get
write_png,ifile,tvrd(true=1),r,g,b

;Combo plot #1
;--------------
window,0,xsize=600,ysize=600
pos=[.1,.09,.98,.95]
rs=pos[0]+(pos[2]-pos[0])*0.65
pos1=[pos[0],pos[1]+(pos[3]-pos[1])*2./4,rs,pos[3]]
pos3=[pos[0],pos[1]+(pos[3]-pos[1])*1./4,rs,pos[1]+(pos[3]-pos[1])*2./4]
pos4=[pos[0],pos[1],rs,pos[1]+(pos[3]-pos[1])*1./4]
pos5=[rs+.015,pos[1]+(pos[3]-pos[1])*2./4,pos[2],pos[3]]
;Plot reflectivity
levels=[-60.,findgen(29)*2.5-50.]
numc=256
colors=[0,fix(indgen(28)*(numc*.78/27))+fix(numc*.196),numc-1]
labels=strtrim(string(fix(levels)),2)
contour,pdbz,w_time,w_height,position=pos1,/norm,color=1,/xstyle,/ystyle,xrange=xra,yrange=wyra,$
  levels=levels,c_colors=colors,/fill,xtickname=blk,$
  title='            ACSE Combo Plot '+datestr,ytitle='Ht [km, AGL]',ylog=ylog
oplot,ctime,cbase1,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase2,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase3,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,vertvis,color=1,psym=2,symsize=0.1,min_value=0.01
colorbar_shupe,[pos3[2]+0.03,pos3[1]-0.02,pos3[2]+0.06,pos3[3]-0.1],levels,colors,$
  nlabels=6,title='dBZ',labels=labels,color=1,charsize=cs,laboffset=1
xyouts,pos1[0]+0.02,pos1[3]-0.03,/norm,'!6Reflectivity!5',charsize=1.8,color=1,charthick=2
;plot LWP
lyra=[0,600]
plot,m_time,smooth(m_lwp,1),position=pos3,/noerase,/norm,color=1,xtickname=blk,xrange=xra,/xstyle,$
  ytitle='LWP [g/m!u2!n]',yrange=lyra,/nodata
oplot,m_time,smooth(m_lwp,1),color=50
oplot,h_timei,smooth(h_lwp,1),color=254
xyouts,pos3[0]+0.02,pos3[3]-0.03,/norm,'!6LWP!5',charsize=1.8,color=1,charthick=2
xyouts,pos3[0]+0.12,pos3[3]-0.03,/norm,'MWR',charsize=1.6,charthick=2,color=50
xyouts,pos3[0]+0.22,pos3[3]-0.03,/norm,'HATPRO',charsize=1.6,charthick=2,color=254
;plot winds
len=(pos4[2]-pos4[0])/(xra[1]-xra[0])*2./3
plot,m_time,m_lwp,position=pos4,/noerase,/norm,/nodata,color=1,xtickformat='timeticks',$
  xtitle='Time [hours, UTC]',ytitle='Ht [km, AGL]',xrange=xra,/xstyle,yrange=[0,6],/ystyle
for i=0,n_elements(p_timehb)-1 do begin
  for j=0,n_elements(p_heighthb)-1 do begin
    if (p_speedhb[i,j] lt 50 and p_speedhb[i,j] gt 0.1) then windhatch,p_timehb[i],p_heighthb[j],p_speedhb[i,j],p_directionhb[i,j],length=len,clip=pos4,thick=2,color=(fix(7*p_speedhb[i,j])+50)<254
  endfor
endfor
xyouts,pos4[0]+0.02,pos4[3]-0.03,/norm,'!6Winds!5',charsize=1.8,color=1,charthick=2
labs=['25','20','15','10','5','0']
offset=0.06
for i=0,n_elements(labs)-1 do xyouts,pos4[2]+0.07,pos4[3]-0.02*(i+1)-offset,/norm,labs[i],color=fix(labs[i])*7+50,align=0.5
xyouts,pos4[2]+0.07,pos4[3]-0.02*(n_elements(labs)+1)-offset,/norm,color=1,'m/s',align=0.5
;plot soundings
sxra=[-50,10]
;sxra=[-35,5]
col=[60,100,200,250,150]
plot,r_temp[0,*],r_height[0,*],position=pos5,/noerase,/norm,color=1,ytickname=blk,yrange=wyra,/ystyle,$
  xrange=sxra,/xstyle,xtitle='T [C]',xminor=1,/nodata,xticks=4
for i=0,(numsonde<5)-1 do begin
  mi=round((r_time[i] mod 1)*60)
  if mi eq 60 then begin
    mi=0
    hrplus=1
  endif else hrplus=0
  if mi lt 10 then mi='0'+strtrim(mi,2) else mi=strtrim(mi,2)
  hr=fix(r_time[i])+hrplus
  if hr lt 10 then hr='0'+strtrim(hr,2) else hr=strtrim(hr,2)
  stst=hr+':'+mi
  oplot,r_temp[i,*],r_height[i,*],thick=2,color=col[i]
  oplot,r_dewpt[i,*],r_height[i,*],thick=2,color=col[i],linestyle=2
  xyouts,pos5[2]-0.07,pos5[3]-0.02-i*0.02,/norm,stst,color=col[i]
endfor
oplot,r_temp[0,*]*0,r_height[0,*],linestyle=1,color=1,thick=1
;add location info
latm=mean(o_lat)
lonm=mean(o_lon)
posstr='Ship Location (mean): '+string(latm,format="(f5.2)")+' N, '+string(lonm,format="(f6.2)")+' E'
xyouts,pos[2]-0.01,0.02,/norm,posstr,align=1,color=1
;write a png file
ifile='acse.combo.'+sarg+'.'+earg+'.png'
tvlct,r,g,b,/get
write_png,ifile,tvrd(true=1),r,g,b


;Thermo plot - MWR
;--------------------
window,0,xsize=500,ysize=600
pos=[.12,.08,.87,.95]
pos1=[pos[0],pos[1]+(pos[3]-pos[1])*2./3,pos[2],pos[1]+(pos[3]-pos[1])*3./3]
pos2=[pos[0],pos[1]+(pos[3]-pos[1])*1./3,pos[2],pos[1]+(pos[3]-pos[1])*2./3]
pos3=[pos[0],pos[1]+(pos[3]-pos[1])*0./3,pos[2],pos[1]+(pos[3]-pos[1])*1./3]
yra=[0,4]

;Plot reflectivity
levels=[-60.,findgen(29)*2.5-50.]
numc=256
colors=[0,fix(indgen(28)*(numc*.78/27))+fix(numc*.196),numc-1]
labels=strtrim(string(fix(levels)),2)
contour,pdbz,w_time,w_height,position=pos1,/norm,color=1,/xstyle,/ystyle,xrange=xra,yrange=yra,$
  levels=levels,c_colors=colors,/fill,xtickname=blk,$
  title='ACSE MWR-sounding, '+datestr,ytitle='Ht [km, AGL]'
oplot,ctime,cbase1,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase2,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase3,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,vertvis,color=1,psym=2,symsize=0.1,min_value=0.01
colorbar_shupe,[pos1[2]+0.03,pos1[1]+0.02,pos1[2]+0.06,pos1[3]-0.05],levels,colors,$
  nlabels=6,title='dBZ',labels=labels,color=1,charsize=cs,laboffset=1
xyouts,pos1[0]+0.02,pos1[3]-0.03,/norm,'!6Reflectivity!5',charsize=1.8,color=1,charthick=2
;Temperature
levels=findgen(38)*1.5-37.5
numc=256
colors=[0,fix(indgen(36)*(numc*.78/35))+fix(numc*.196),numc-1]
labels=strtrim(string(fix(levels)),2)
contour,m_tprof-273.16,m_time,m_height,position=pos2,/noerase,/norm,color=1,/xstyle,/ystyle,xrange=xra,yrange=yra,$
  levels=levels,c_colors=colors,/fill,xtickname=blk,$
  ytitle='Ht [km, AGL]'
oplot,ctime,cbase1,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase2,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase3,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,vertvis,color=1,psym=2,symsize=0.1,min_value=0.01
colorbar_shupe,[pos2[2]+0.03,pos2[1]+0.02,pos2[2]+0.06,pos2[3]-0.05],levels,colors,$
  nlabels=6,title='C',labels=labels,color=1,charsize=cs,laboffset=1
for i=0,n_elements(r_time)-1 do begin
  polyfill,(r_time[i]+[-0.5,-0.5,0.25,0.25])>0.01,[0,1,1,0]*yra[1]*0.99,color=0
  tcolors=interpol(colors,levels,r_temp[i,*])
  wh=where(r_height[i,*] lt yra[1],nm)
  if nm gt 0 then begin
    for j=0,nm-1 do begin
      polyfill,(r_time[i]+[-0.44,-0.44,0.19,0.19])>0.01,[r_height[i,wh[j]],r_height[i,wh[j]+1],r_height[i,wh[j]+1],r_height[i,wh[j]]]<yra[1]*0.99,color=tcolors[j]
    endfor
  endif
endfor
xyouts,pos2[0]+0.02,pos2[3]-0.03,/norm,'!6Temperature!5',charsize=1.8,color=1,charthick=2
;calculate radiosonde vapor density
estr=6.11  ;hPa
Ttr=273.15 ;K
Llv=2.54e6 ;J/kg
Rv=461.0   ;J/K/kg
eps=0.622
es=estr*exp(Llv/Rv*((1.0/Ttr)-(1.0/(r_temp+Ttr))))
ws=eps*es/(r_pres-es)
r_rho=r_pres/(r_temp+Ttr)/287.*100.  ;kg/m3
r_qv=r_rh/100.*ws/(1.0+r_rh/100.*ws)*1000.  ;g/kg
r_qq=r_qv*r_rho  ;g/m3
;plot vapor density
levels=[findgen(33)*3/10.]
numc=256
colors=[0,fix(indgen(31)*(numc*.78/30))+fix(numc*.196),numc-1]
labels=string(levels,format='(f3.1)')
contour,m_qprof,m_time,m_height,position=pos3,/noerase,/norm,color=1,/xstyle,/ystyle,xrange=xra,yrange=yra,$
  levels=levels,c_colors=colors,/fill,xtickformat='timeticks',xtitle='Time [hours, UTC]',$
  ytitle='Ht [km, AGL]'
oplot,ctime,cbase1,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase2,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase3,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,vertvis,color=1,psym=2,symsize=0.1,min_value=0.01
colorbar_shupe,[pos3[2]+0.03,pos3[1]+0.02,pos3[2]+0.06,pos3[3]-0.05],levels,colors,$
  nlabels=6,title='g/m!u3!n',labels=labels,color=1,charsize=cs,laboffset=1
for i=0,n_elements(r_time)-1 do begin
  polyfill,(r_time[i]+[-0.5,-0.5,0.25,0.25])>0.01,[0,1,1,0]*yra[1]*0.99,color=0
  wh=where(r_qq[i,*] gt max(levels),nm)
  if nm gt 0 then r_qq[i,wh]=max(levels)
  tcolors=interpol(colors,levels,r_qq[i,*])
  wh=where(r_height[i,*] lt yra[1],nm)
  if nm gt 0 then begin
    for j=0,nm-1 do begin
      polyfill,(r_time[i]+[-0.44,-0.44,0.19,0.19])>0.01,[r_height[i,wh[j]],r_height[i,wh[j]+1],r_height[i,wh[j]+1],r_height[i,wh[j]]]<yra[1]*0.99,color=tcolors[j]
    endfor
  endif
endfor
xyouts,pos3[0]+0.02,pos3[3]-0.03,/norm,'!6Vapor Density!5',charsize=1.8,color=1,charthick=2
;write a png file
ifile='acse.thermo_mwr.'+sarg+'.'+earg+'.png'
tvlct,r,g,b,/get
write_png,ifile,tvrd(true=1),r,g,b

;Thermo plot - HATPRO
;---------------------
window,0,xsize=500,ysize=600
pos=[.12,.08,.87,.95]
pos1=[pos[0],pos[1]+(pos[3]-pos[1])*2./3,pos[2],pos[1]+(pos[3]-pos[1])*3./3]
pos2=[pos[0],pos[1]+(pos[3]-pos[1])*1./3,pos[2],pos[1]+(pos[3]-pos[1])*2./3]
pos3=[pos[0],pos[1]+(pos[3]-pos[1])*0./3,pos[2],pos[1]+(pos[3]-pos[1])*1./3]
yra=[0,4]
;Plot reflectivity
levels=[-60.,findgen(29)*2.5-50.]
numc=256
colors=[0,fix(indgen(28)*(numc*.78/27))+fix(numc*.196),numc-1]
labels=strtrim(string(fix(levels)),2)
contour,pdbz,w_time,w_height,position=pos1,/norm,color=1,/xstyle,/ystyle,xrange=xra,yrange=yra,$
  levels=levels,c_colors=colors,/fill,xtickname=blk,$
  title='ACSE HATPRO-sounding, '+datestr,ytitle='Ht [km, AGL]'
oplot,ctime,cbase1,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase2,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase3,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,vertvis,color=1,psym=2,symsize=0.1,min_value=0.01
colorbar_shupe,[pos1[2]+0.03,pos1[1]+0.02,pos1[2]+0.06,pos1[3]-0.05],levels,colors,$
  nlabels=6,title='dBZ',labels=labels,color=1,charsize=cs,laboffset=1
xyouts,pos1[0]+0.02,pos1[3]-0.03,/norm,'!6Reflectivity!5',charsize=1.8,color=1,charthick=2
;Temperature
levels=findgen(38)*1.5-37.5
numc=256
colors=[0,fix(indgen(36)*(numc*.78/35))+fix(numc*.196),numc-1]
labels=strtrim(string(fix(levels)),2)
contour,h_tprof-273.16,h_timeh,h_height,position=pos2,/noerase,/norm,color=1,/xstyle,/ystyle,xrange=xra,yrange=yra,$
  levels=levels,c_colors=colors,/fill,xtickname=blk,$
  ytitle='Ht [km, AGL]'
oplot,ctime,cbase1,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase2,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase3,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,vertvis,color=1,psym=2,symsize=0.1,min_value=0.01
colorbar_shupe,[pos2[2]+0.03,pos2[1]+0.02,pos2[2]+0.06,pos2[3]-0.05],levels,colors,$
  nlabels=6,title='C',labels=labels,color=1,charsize=cs,laboffset=1
for i=0,n_elements(r_time)-1 do begin
  polyfill,(r_time[i]+[-0.5,-0.5,0.25,0.25])>0.01,[0,1,1,0]*yra[1]*0.99,color=0
  tcolors=interpol(colors,levels,r_temp[i,*])
  wh=where(r_height[i,*] lt yra[1],nm)
  if nm gt 0 then begin
    for j=0,nm-1 do begin
      polyfill,(r_time[i]+[-0.44,-0.44,0.19,0.19])>0.01,[r_height[i,wh[j]],r_height[i,wh[j]+1],r_height[i,wh[j]+1],r_height[i,wh[j]]]<yra[1]*0.99,color=tcolors[j]
    endfor
  endif
endfor
xyouts,pos2[0]+0.02,pos2[3]-0.03,/norm,'!6Temperature!5',charsize=1.8,color=1,charthick=2
;plot vapor density
levels=[findgen(33)*3/10.]
numc=256
colors=[0,fix(indgen(31)*(numc*.78/30))+fix(numc*.196),numc-1]
labels=string(levels,format='(f3.1)')
contour,h_qprof,h_timeh,h_height,position=pos3,/noerase,/norm,color=1,/xstyle,/ystyle,xrange=xra,yrange=yra,$
  levels=levels,c_colors=colors,/fill,xtickformat='timeticks',xtitle='Time [hours, UTC]',$
  ytitle='Ht [km, AGL]'
oplot,ctime,cbase1,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase2,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,cbase3,color=1,psym=2,symsize=0.1,min_value=0.01
oplot,ctime,vertvis,color=1,psym=2,symsize=0.1,min_value=0.01
colorbar_shupe,[pos3[2]+0.03,pos3[1]+0.02,pos3[2]+0.06,pos3[3]-0.05],levels,colors,$
  nlabels=6,title='g/m!u3!n',labels=labels,color=1,charsize=cs,laboffset=1
for i=0,n_elements(r_time)-1 do begin
  polyfill,(r_time[i]+[-0.5,-0.5,0.25,0.25])>0.01,[0,1,1,0]*yra[1]*0.99,color=0
  wh=where(r_qq[i,*] gt max(levels),nm)
  if nm gt 0 then r_qq[i,wh]=max(levels)
  tcolors=interpol(colors,levels,r_qq[i,*])
  wh=where(r_height[i,*] lt yra[1],nm)
  if nm gt 0 then begin
    for j=0,nm-1 do begin
      polyfill,(r_time[i]+[-0.44,-0.44,0.19,0.19])>0.01,[r_height[i,wh[j]],r_height[i,wh[j]+1],r_height[i,wh[j]+1],r_height[i,wh[j]]]<yra[1]*0.99,color=tcolors[j]
    endfor
  endif
endfor
xyouts,pos3[0]+0.02,pos3[3]-0.03,/norm,'!6Vapor Density!5',charsize=1.8,color=1,charthick=2
;write a png file
ifile='acse.thermo_hat.'+sarg+'.'+earg+'.png'
tvlct,r,g,b,/get
write_png,ifile,tvrd(true=1),r,g,b


;Make a wind comparison plot
;------------------------------
window,0,xsize=800,ysize=650
pos=[.08,.08,.98,.95]
posi=fltarr(4,5,2)
posi[*,0,0]=[pos[0]+(pos[2]-pos[0])*0./5,pos[1]+(pos[3]-pos[1])*1.1/2,pos[0]+(pos[2]-pos[0])*1./5,pos[1]+(pos[3]-pos[1])*2./2]
posi[*,0,1]=[pos[0]+(pos[2]-pos[0])*0./5,pos[1]+(pos[3]-pos[1])*0./2,pos[0]+(pos[2]-pos[0])*1./5,pos[1]+(pos[3]-pos[1])*0.9/2]
posi[*,1,0]=[pos[0]+(pos[2]-pos[0])*1./5,pos[1]+(pos[3]-pos[1])*1.1/2,pos[0]+(pos[2]-pos[0])*2./5,pos[1]+(pos[3]-pos[1])*2./2]
posi[*,1,1]=[pos[0]+(pos[2]-pos[0])*1./5,pos[1]+(pos[3]-pos[1])*0./2,pos[0]+(pos[2]-pos[0])*2./5,pos[1]+(pos[3]-pos[1])*0.9/2]
posi[*,2,0]=[pos[0]+(pos[2]-pos[0])*2./5,pos[1]+(pos[3]-pos[1])*1.1/2,pos[0]+(pos[2]-pos[0])*3./5,pos[1]+(pos[3]-pos[1])*2./2]
posi[*,2,1]=[pos[0]+(pos[2]-pos[0])*2./5,pos[1]+(pos[3]-pos[1])*0./2,pos[0]+(pos[2]-pos[0])*3./5,pos[1]+(pos[3]-pos[1])*0.9/2]
posi[*,3,0]=[pos[0]+(pos[2]-pos[0])*3./5,pos[1]+(pos[3]-pos[1])*1.1/2,pos[0]+(pos[2]-pos[0])*4./5,pos[1]+(pos[3]-pos[1])*2./2]
posi[*,3,1]=[pos[0]+(pos[2]-pos[0])*3./5,pos[1]+(pos[3]-pos[1])*0./2,pos[0]+(pos[2]-pos[0])*4./5,pos[1]+(pos[3]-pos[1])*0.9/2]
posi[*,4,0]=[pos[0]+(pos[2]-pos[0])*4./5,pos[1]+(pos[3]-pos[1])*1.1/2,pos[0]+(pos[2]-pos[0])*5./5,pos[1]+(pos[3]-pos[1])*2./2]
posi[*,4,1]=[pos[0]+(pos[2]-pos[0])*4./5,pos[1]+(pos[3]-pos[1])*0./2,pos[0]+(pos[2]-pos[0])*5./5,pos[1]+(pos[3]-pos[1])*0.9/2]

yra=[0,6]
for i=0,(numsonde-1)<4 do begin
  if i eq 0 then begin
    ytic=''
    ytit='Height [km]'
    noer=0 
  endif else begin
    ytic=strarr(10)+' '
    ytit=''
    noer=1
  endelse
  ;get closest RWP and ship time
  wi=closest(p_timehb,r_time[i],0)
  if r_time[i]-p_timehb[wi] lt 1 then wi=(wi-1)>0
  li=closest(p_timelb,r_time[i],0)
  if r_time[i]-p_timelb[li] lt 1 then li=(li-1)>0
  si=closest(o_time,r_time[i]-0.5,0)
  dell=abs(r_time[i]-p_timelb[li])
  delh=abs(r_time[i]-p_timehb[wi])
  
  mxs=(max([max(p_speedlb[li,*]),max(p_speedhb[wi,*]),max(r_speed[i,*])])+1)<25
  xras=[0,mxs]
  xrad=[0,360]

  plot,r_speed[i,*],r_height[i,*],/nodata,color=1,position=posi[*,i,0],noerase=noer,ytickname=ytic,ytitle=ytit,$
      xrange=xras,/xstyle,yrange=yra,/ystyle,thick=th,xthick=th,ythick=th,charthick=th,charsize=1.1
  oplot,r_speed[i,*],r_height[i,*],color=50
  if dell le 2 then oplot,p_speedlb[li,*],p_heightlb,color=210,psym=2,symsize=0.3
  if delh le 2 then oplot,p_speedhb[wi,*],p_heighthb,color=254,psym=2,symsize=0.3
;  oplot,p_speedla[li,*],p_heightla,color=100,psym=2,symsize=0.3
;  oplot,p_speedha[wi,*],p_heightha,color=150,psym=2,symsize=0.3 
  plots,o_wspd[si],0,color=100,psym=2,thick=th
  
  plot,r_direction[i,*],r_height[i,*],/nodata,color=1,position=posi[*,i,1],/noerase,ytickname=ytic,ytitle=ytit,$
      xrange=xrad,/xstyle,yrange=yra,/ystyle,xminor=2,thick=th,xthick=th,ythick=th,charthick=th,charsize=1.1
  oplot,r_direction[i,*],r_height[i,*],color=50
  if dell le 2 then oplot,p_directionlb[li,*],p_heightlb,color=210,psym=2,symsize=0.3
  if delh le 2 then oplot,p_directionhb[wi,*],p_heighthb,color=254,psym=2,symsize=0.3
;  oplot,p_directionla[li,*],p_heightla,color=100,psym=2,symsize=0.3
;  oplot,p_directionha[wi,*],p_heightha,color=150,psym=2,symsize=0.3
  plots,o_wdir[si],0,color=100,psym=2,thick=th
  xyouts,posi[0,i,1]+0.01,posi[3,i,1]-0.03,/norm,'Sonde: '+string(r_time[i],format="(i2)"),color=50,charthick=th,charsize=1.2
  xyouts,posi[0,i,1]+0.01,posi[3,i,1]-0.06,/norm,'RWP_H: '+string(p_timehb[wi],format="(i2)"),color=254,charthick=th,charsize=1.2
  xyouts,posi[0,i,1]+0.01,posi[3,i,1]-0.09,/norm,'RWP_L: '+string(p_timelb[li],format="(i2)"),color=210,charthick=th,charsize=1.2
  xyouts,posi[0,i,1]+0.01,posi[3,i,1]-0.12,/norm,'Ship: '+string(o_time[si],format="(f4.1)"),color=100,charthick=th,charsize=1.2
  
  xyouts,pos[0]+(pos[2]-pos[0])/2.,pos[3]+0.01,/norm,align=0.5,color=1,'ACSE, Wind comp, '+datestr,charsize=1.3,charthick=th
  xyouts,pos[0]+(pos[2]-pos[0])/2.,posi[1,0,0]-0.05,/norm,'Wind Speed [m/s]',color=1,align=0.5,charsize=1.1,charthick=th
  xyouts,pos[0]+(pos[2]-pos[0])/2.,pos[1]-0.05,/norm,'Wind Direction [deg]',color=1,align=0.5,charsize=1.1,charthick=th

endfor
;write a png file
ifile='acse.windcomp.'+sarg+'.'+earg+'.png'
tvlct,r,g,b,/get
write_png,ifile,tvrd(true=1),r,g,b

;Basic Met/state compare
;------------------------
window,0,xsize=500,ysize=600
pos=[.15,.08,.95,.95]
pos1=[pos[0],pos[1]+(pos[3]-pos[1])*3./4,pos[2],pos[1]+(pos[3]-pos[1])*4./4]
pos2=[pos[0],pos[1]+(pos[3]-pos[1])*2./4,pos[2],pos[1]+(pos[3]-pos[1])*3./4]
pos3=[pos[0],pos[1]+(pos[3]-pos[1])*0./4,pos[2],pos[1]+(pos[3]-pos[1])*2./4]
th=2
cs=1.2
;pres
mxint=(ceil(max([max(m_pres),max(o_pres),max(r_pres[*,0])]))+1)<1040
mnint=(floor(min([min(m_pres),min(o_pres),min(r_pres[*,0])]))-1)>970
plot,m_time,m_pres,position=pos1,/nodata,xrange=xra,/xstyle,$
  title='!6ACSE State, '+datestr, yrange=[mnint,mxint],/ystyle,xtickname=blk,ytitle='Pres [mb]',$
  thick=th,xthick=th,ythick=th,charthick=th,charsize=cs
oplot,m_time,smooth(m_pres,21),color=254
oplot,o_time,smooth(o_pres,21),color=48
plots,r_time,r_pres[*,0],psym=2,symsize=1.5,color=60,thick=2
xyouts,pos1[0]+0.02,pos1[3]-0.03,/norm,'Ship',color=48,charsize=cs
xyouts,pos1[0]+0.11,pos1[3]-0.03,/norm,'Sonde',color=60,charsize=cs
xyouts,pos1[0]+0.22,pos1[3]-0.03,/norm,'MWR',color=254,charsize=cs
;rh
mnint=(floor(min([min(m_rh),min(o_rhair),min(r_rh[*,0])]))-2)>20
plot,m_time,m_rh,position=pos2,/noerase,/nodata,xrange=xra,/xstyle,$
  yrange=[mnint,104],/ystyle,xtickname=blk,ytitle='RH [%]',$
  thick=th,xthick=th,ythick=th,charthick=th,charsize=cs
oplot,m_time,m_rh,color=254
oplot,o_time,smooth(o_rhair,31),color=48
plots,r_time,r_rh[*,0],psym=2,symsize=1.5,color=60,thick=2
xyouts,pos2[0]+0.02,pos2[3]-0.03,/norm,'Ship',color=48,charsize=cs
xyouts,pos2[0]+0.11,pos2[3]-0.03,/norm,'Sonde',color=60,charsize=cs
xyouts,pos2[0]+0.22,pos2[3]-0.03,/norm,'MWR',color=254,charsize=cs
;temp
mxint=(ceil(max([max(o_tair),max(o_tsea),max(m_tamb),max(r_temp[*,0]),max(s_temp),max(i_irt1),max(i_irt2)]))+1)<20>5
mnint=(floor(min([min(o_tair),min(o_tsea),min(m_tamb),min(r_temp[*,0]),min(s_temp),min(i_irt1),min(i_irt2)]))-1)>(-30)<0
plot,m_time,m_tamb,position=pos3,/noerase,/nodata,xrange=xra,/xstyle,$
  yrange=[mnint,mxint],/ystyle,xtickformat='timeticks',ytitle='T [C]',xtitle='Time [UTC]',$
  thick=th,xthick=th,ythick=th,charthick=th,charsize=cs
oplot,i_time1,smooth(i_irt1,71),color=200
oplot,i_time2,smooth(i_irt2,71),color=220
oplot,m_time,smooth(m_tamb,11),color=254
oplot,o_time,smooth(o_tair,21),color=48
oplot,s_time,s_temp,color=120,thick=2
oplot,o_time,o_tsea,color=85,thick=2
plots,r_time,r_temp[*,0],psym=2,symsize=1.5,color=60,thick=2
xyouts,pos3[0]+0.02,pos3[3]-0.03,/norm,'Ship',color=48,charsize=cs
xyouts,pos3[0]+0.11,pos3[3]-0.03,/norm,'Sonde',color=60,charsize=cs
xyouts,pos3[0]+0.22,pos3[3]-0.03,/norm,'MWR',color=254,charsize=cs
xyouts,pos3[0]+0.31,pos3[3]-0.03,/norm,'IRT1',color=200,charsize=cs
xyouts,pos3[0]+0.38,pos3[3]-0.03,/norm,'IRT2',color=220,charsize=cs
xyouts,pos3[0]+0.47,pos3[3]-0.03,/norm,'Snake',color=120,charsize=cs
xyouts,pos3[0]+0.58,pos3[3]-0.03,/norm,'Sea-8m',color=85,charsize=cs
;write a png file
ifile='acse.statecomp.'+sarg+'.'+earg+'.png'
tvlct,r,g,b,/get
write_png,ifile,tvrd(true=1),r,g,b

;Make microwave comparison plot
window,0,xsize=500,ysize=600
pos=[.12,.08,.95,.95]
pos1=[pos[0],pos[1]+(pos[3]-pos[1])*4./5,pos[2],pos[1]+(pos[3]-pos[1])*5./5]
pos2=[pos[0],pos[1]+(pos[3]-pos[1])*3./5,pos[2],pos[1]+(pos[3]-pos[1])*4./5]
pos3=[pos[0],pos[1]+(pos[3]-pos[1])*2./5,pos[2],pos[1]+(pos[3]-pos[1])*3./5]
pos4=[pos[0],pos[1]+(pos[3]-pos[1])*1./5,pos[2],pos[1]+(pos[3]-pos[1])*2./5]
pos5=[pos[0],pos[1]+(pos[3]-pos[1])*0./5,pos[2],pos[1]+(pos[3]-pos[1])*1./5]

xra=[0,24]
mnpt=min([n_elements(mb_time),n_elements(hb_time)/2.])
sm=31<(mnpt-1)


ymx=max([max(mb_t22p234),max(mb_t28p000),max(hb_t22p24),max(hb_t27p84)])<200
ymn=min([min(mb_t22p234),min(mb_t28p000),min(hb_t22p24),min(hb_t27p84)])
yra=[(ymn-2)>0,ymx+2]
plot,mb_time,mb_t22p234,position=pos1,/nodata,color=1,xrange=xra,/xstyle,xtickname=blk,$
   thick=th,xthick=th,ythick=th,charthick=th,ytitle='T!dB!n [K]',yminor=1,$
   title='!6ACSE Micro Compare, '+datestr,yrange=yra,/ystyle
oplot,mb_time,smooth(mb_t22p234,sm),color=50
oplot,mb_time,smooth(mb_t28p000,sm),color=50,linestyle=1,thick=2
oplot,hb_time,smooth(hb_t22p24,2*sm),color=250
oplot,hb_time,smooth(hb_t27p84,2*sm),color=250,linestyle=1,thick=2
xyouts,pos1[0]+0.3,pos1[3]-0.02,/norm,'Solid~22.24GHz, Dash~28GHz'

ymx=max([max(mb_t52p280),max(mb_t56p660-110.),max(hb_t52p28),max(hb_t56p66-110.)])<200
ymn=min([min(mb_t52p280),min(mb_t56p660-110.),min(hb_t52p28),min(hb_t56p66-110.)])
yra=[(ymn-2)>0,ymx+2]
plot,mb_time,mb_t52p280,position=pos2,/noerase,/nodata,color=1,xrange=xra,/xstyle,xtickname=blk,$
  thick=th,xthick=th,ythick=th,charthick=th,ytitle='T!dB!n [K]',yrange=yra,/ystyle
oplot,mb_time,smooth(mb_t52p280,sm),color=50
oplot,mb_time,smooth(mb_t56p660-110.,sm),color=50,linestyle=1,thick=2
oplot,hb_time,smooth(hb_t52p28,2*sm),color=250
oplot,hb_time,smooth(hb_t56p66-110.,2*sm),color=250,linestyle=1,thick=2
xyouts,pos2[0]+0.3,pos2[3]-0.02,/norm,'Solid=52.28GHz, Dash=56.66GHz -110K'

ymx=max([max(m_pwv),max(h_pwv/10000.)])
yra=[0,ymx+0.1]
plot,m_time,m_pwv,position=pos3,/noerase,/nodata,color=1,xrange=xra,/xstyle,xtickname=blk,$
  thick=th,xthick=th,ythick=th,charthick=th,ytitle='PWV [cm]',yrange=yra,/ystyle,yminor=1
oplot,m_time,smooth(m_pwv,sm),color=50
oplot,h_timei,smooth(h_pwv/10000.,2*sm),color=250
xyouts,pos3[0]+0.2,pos3[1]+0.02,/norm,'MWR',color=50,charsize=2
xyouts,pos3[0]+0.4,pos3[1]+0.02,/norm,'HATPRO',color=250,charsize=2

ymx=max([max(m_lwp),max(h_lwp)])<990
yra=[0,ymx]
plot,m_time,m_lwp,position=pos4,/noerase,/nodata,color=1,xrange=xra,/xstyle,xtickname=blk,$
  thick=th,xthick=th,ythick=th,charthick=th,ytitle='LWP',yminor=1,yrange=yra,/ystyle
oplot,m_time,smooth(m_lwp,sm),color=50
oplot,h_timei,smooth(h_lwp,2*sm),color=250

ymx=max([max(m_tir)])
ymn=min([min(m_tir)])
yra=[ymn-10,ymx+10]
plot,m_time,m_tir,position=pos5,/noerase,/nodata,color=1,xrange=xra,/xstyle,xtickformat='timeticks',$
  thick=th,xthick=th,ythick=th,charthick=th,ytitle='T!dIR!n [K]',xtitle='Time [hours,UTC]',yrange=yra,/ystyle
oplot,m_time,smooth(m_tir,sm/2.),color=50

;write a png file
ifile='acse.microcomp.'+sarg+'.'+earg+'.png'
tvlct,r,g,b,/get
write_png,ifile,tvrd(true=1),r,g,b

endif ;plot_data

end