;
;  This program plots a daily summary of the arm 2d variables. 
;
pro plot_average_daily_2d_arm,site,date,avg_int,ret_type

;*************************************
;*************************************
;  INPUTS
;*************************************
;*************************************

;  Site is the following 'sgp','nsa','twpc1','twpc2','twpc3'
site='iphex'

;  Date in the form yyyymmdd 
;date='20000929' ;nsa
date='20140518' ;iphex

;  Averaging interval  '5min' '1hr'
avg_int='5min'

;  Retrieval version
ret_type='mine'  ;'basep','mine','temp'

;**************************************
;  Commons
;**************************************

;  Color common
common colors,r_orig,g_orig,b_orig,r_curr,g_curr,b_curr

;*****************************************
;  Constants
;*****************************************

;  Indicate which files you want to use, the monthly or daily files.
m_or_d='day'  ;day,month

;  Missing data value
missing=-8888

;  Nocloud value
nocloud=-9999

;  Mplflag value
mplflag=-9000

;  Seconds in a day
sec_per_day=double(24L*60L*60L)

;  Get IDL julian day of day 1
day1=julday(1,1,1970,0,0,0)

;  Arm data directory
arm_dir='/uufs/chpc.utah.edu/common/home/mace-group3/arm/'+site+'/'

;*****************************************
;  Read in the monthly averaged data
;*****************************************

;  Make sure date is a string
date=strcompress(date,/remove_all)

;  Don't need to fill daily data
fill_flag='no'

;  List of variables to read in
varname=['base_time','time_offset','height',$
  'ReflectivityBestEstimate',$
  'ReflectivityPossiblePoints',$
  'ReflectivityCloudyPoints',$
  'temp',$
  'MeanDopplerVelocity',$
  'pressure',$
  'SpectralWidth',$
  'mxrat',$
  'column_cfrac',$
  'column_cfrac_mpl']     

;  Get the variables in the list from the file
if site eq 'iphex' then begin
  get_average_data_list_iphex,varname,data,$
    site,date,avg_int,ret_type,m_or_d,fill_flag
endif else begin
  get_average_data_list,varname,data,$
    site,date,avg_int,ret_type,m_or_d,fill_flag
endelse

;  Pull out the individual data arrays
for v=0,n_elements(varname)-1 do begin
  xstr=varname[v]+'=data.'+varname[v]
  result=execute(xstr)
endfor  

if site eq 'iphex' then column_cfrac=column_cfrac*100.0

;  Get the size of the arrays
numtimes=n_elements(time_offset)
numheights=n_elements(height)
  
;  Convert height from meters to km
height=height/1000.0  

;***********************************
;  Convert time to julian day for plotting
;***********************************

;  Calculate julian_day
julian_day=double(day1+(base_time+time_offset)/sec_per_day)

;  Array of standard dates
caldat,julian_day,mm,dd,yyyy,hour,minute,second

;  Starting date and time
syear=yyyy[0]
smonth=mm[0]
sday=dd[0]
shour=hour[0]
smin=minute[0]
ssec=second[0]

;  Ending date and time
eyear=yyyy[numtimes-1]
emonth=mm[numtimes-1]
eday=dd[numtimes-1]
ehour=hour[numtimes-1]
emin=minute[numtimes-1]
esec=second[numtimes-1]

;********************************
;  Set up the plot device
;********************************

;  Set up the device
set_plot,'z'

;  xsize of the plot window
xwin=800  

;  ysize of the plot window
ywin=800  

;  Size the window
device,set_resolution=[xwin,ywin]

;********************************
;  Set up the colors
;  0-253 rainbow
;  254   gray
;  255   white
;********************************

;  Load the color table
loadct,39   

;  Top is the last color to scale
top_color=252

;  make hot pink for mplflag data in color 253
r_curr[253]=255
g_curr[253]=20
b_curr[253]=147

;  Make gray for missing data in color 253
r_curr[254]=200
g_curr[254]=200
b_curr[254]=200

;  Make the background white
!p.background=!d.n_colors-1  

;********************************
;  Set up the plot
;********************************

;  Date label format
dummy=label_date(date_format=['%H:%I'])

;  Xticks point out from plot
!x.ticklen=-0.02

;  Charsizes
tsz=1.2
cbsz=0.8
!y.charsize=0.8

;********************************
;  Position vectors
;********************************

xl=0.05  &  xr=0.94
yb=0.05  &  yt=0.94
sx=0.08        ;x buffer space
sy=0.00 ;0.05  ;y buffer space
numplots_x=2  ;number of plots in x direction
numplots_y=5  ;number of plots in y direction

position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos

;  String positioning
dr=0.09  ;string offset from the right xcoor
dl=0.025  ;string offset from the left xcoor
dt=0.025  ;string offset from the top ycoor

;**********************************
;  Plot 0 - reflectivity
;**********************************

;  Plot location
pnum=0

;  Find the valid data
result=where(ReflectivityBestEstimate gt missing,count)
if count gt 0 and site ne 'iphex' then begin
  ReflectivityBestEstimate[result]=ReflectivityBestEstimate[result]/100.0
endif

;  Filter with clutter 
;  Assign all bad clutter reflectivity data -9999
if (site eq 'sgp' or site eq 'nsa' or site eq 'twpc3') and $
   n_elements(qc_ReflectivityClutterFlagMode) ne 0 then begin
  result=where(qc_ReflectivityClutterFlagMode eq 0, count)
  if count gt 0 then ReflectivityBestEstimate[result]=nocloud
endif

;  Turn missing data to gray
result=where(ReflectivityPossiblePoints eq 0,count)
if count gt 0 then ReflectivityBestEstimate[result]=missing

;  Find the max and min data values
dmax=20.0
dmin=-60.0

;  Bytscl the data
image=bytscl(ReflectivityBestEstimate,top=top_color,max=dmax,min=dmin)

;  Put all no cloud values to white
result=where(ReflectivityBestEstimate eq nocloud,count)
if count gt 0 then image[result]=255

;  Put all missing data to gray
result=where(ReflectivityBestEstimate eq missing,count)
if count gt 0 then image[result]=254

;  Put all mplflag data to hot pink
result=where(ReflectivityBestEstimate eq mplflag,count)
if count gt 0 then image[result]=253

;  Set up the plot space
contour,ReflectivityBestEstimate,julian_day,height,/nodata,$
  xstyle=4,ystyle=4,position=pos[pnum,*]

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image,sx,sy),px[0],py[0]

;  Put up the axes
contour,ReflectivityBestEstimate,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='label_date',/nodata,$
  ytitle='Height (km)',position=pos[pnum,*],/noerase

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'reflectivity  DbZ',$
  color=0,/normal,charsize=tsz
xyouts,pos[pnum,0]+0.2,pos[pnum,3]-dt,'mplmask',$
  color=253,/normal,charsize=tsz

;  Color bar
cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=254,format='(F6.2)',divisions=4,$
  vertical=1,color=0,position=cbpos,charsize=cbsz

;**********************************
;  Plot 1 - temperature
;**********************************

;  Plot location
pnum=1

;  Find the valid data
result=where(temp gt missing,count)

;  Find the max and min data values
dmax=300.0
dmin=200.0

;  Bytscl the data
image=bytscl(temp,top=top_color,max=dmax,min=dmin)

;  Put all no cloud values to white
result=where(temp eq nocloud,count)
if count gt 0 then image[result]=255

;  Put all missing data to gray
result=where(temp eq missing,count)
if count gt 0 then image[result]=254

;  Set up the plot space
contour,temp,julian_day,height,/nodata,/noerase,$
  xstyle=4,ystyle=4,position=pos[pnum,*]

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image,sx,sy),px[0],py[0]

;  Put up the axes
contour,temp,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='label_date',$
  position=pos[pnum,*],/noerase,/nodata

;  Put up the axes
contour,temp,julian_day,height,xstyle=4,ystyle=4,$
  color=255,xtickunits='time',xtickformat='label_date',$
  position=pos[pnum,*],/noerase,$
  levels=[200,210,220,240,260,280,300,320,340,360,380],$
  c_labels=[1,1,1,1,1,1,1,1,1,1,1],c_charsize=1.0

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'temp  Kelvin',$
  color=255,/normal,charsize=tsz

;  Color bar
cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=254,format='(F6.2)',divisions=4,$
  vertical=1,color=0,position=cbpos,charsize=cbsz

;**********************************
;  Plot 2 - Velocity
;**********************************

;  Plot location
pnum=2

;  Find the valid data
result=where(MeanDopplerVelocity gt missing,count)
if count gt 0 and site ne 'iphex' then begin
  MeanDopplerVelocity[result]=MeanDopplerVelocity[result]/1000.0
endif
;if count gt 0 and site eq 'iphex' then begin
;  MeanDopplerVelocity[result]=MeanDopplerVelocity[result]*1000.0
;endif
;  Turn missing data to gray
result=where(ReflectivityPossiblePoints eq 0,count)
if count gt 0 then MeanDopplerVelocity[result]=missing

;  Find the max and min data values
result=where(MeanDopplerVelocity gt missing,count)
if count eq 0 then begin
  dmax=1.0
  dmin=0.0
endif else begin
  dmax=max(MeanDopplerVelocity[result])
  dmin=min(MeanDopplerVelocity[result])
  if dmin eq dmax then begin
    dmax=dmin+1.0
  endif
endelse

;  Bytscl the data
image=bytscl(MeanDopplerVelocity,top=top_color,max=dmax,min=dmin)

;  Put all no cloud values to white
result=where(MeanDopplerVelocity eq nocloud,count)
if count gt 0 then image[result]=255

;  Put all missing data to gray
result=where(MeanDopplerVelocity eq missing,count)
if count gt 0 then image[result]=254

;  Set up the plot space
contour,MeanDopplerVelocity,julian_day,height,/nodata,$
  xstyle=4,ystyle=4,position=pos[pnum,*],/noerase

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image,sx,sy),px[0],py[0]

;  Put up the axes
contour,MeanDopplerVelocity,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',/nodata,$
  ytitle='Height (km)',position=pos[pnum,*],/noerase

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'doppler velocity  m/s',$
  color=0,/normal,charsize=tsz

;  Color bar
cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=254,format='(F8.2)',$
  vertical=1,color=0,position=cbpos,charsize=cbsz

;**********************************
;  Plot 3 - pressure
;**********************************

;  Plot location
pnum=3

;  Find the valid data
result=where(pressure gt -8888,count)

;  Convert from pascals to millibars
if count gt 0 then begin
  pressure[result]=pressure[result]*0.01
endif

;  Find the max and min data values
dmax=1000.0
dmin=100.0

;  Bytscl the data
image=bytscl(pressure,top=253,max=dmax,min=dmin)

;  Put all no cloud values to white
result=where(pressure eq nocloud,count)
if count gt 0 then image[result]=255

;  Put all missing data to gray
result=where(pressure eq missing,count)
if count gt 0 then image[result]=254

;  Set up the plot space
contour,pressure,julian_day,height,/nodata,/noerase,$
  xstyle=4,ystyle=4,position=pos[pnum,*]

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image,sx,sy),px[0],py[0]

;  Put up the axes
contour,pressure,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',$
  position=pos[pnum,*],/noerase,/nodata;,$
  ;levels=[100,200,300,400,500,600,700,800,900,1000],$
  ;c_labels=[1,1,1,1,1,1,1,1,1,1],c_charsize=1.0

contour,pressure,julian_day,height,xstyle=4,ystyle=4,$
  color=255,xtickunits='time',xtickformat='(A1)',$
  position=pos[pnum,*],/noerase,$
  levels=[100,115,150,200,300,400,500,700,900,1000],$
  c_labels=[1,1,1,1,1,1,1,1,1,1],c_charsize=1.0

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'pressure  mb',$
  color=255,/normal,charsize=tsz

;  Color bar
cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=254,format='(F7.2)',divisions=4,$
  vertical=1,color=0,position=cbpos,charsize=cbsz

;*********************************************
;  Plot 4 - spectrum
;*********************************************

;  Plot location
pnum=4

;  Find the valid data
result=where(SpectralWidth gt missing,count)
if count gt 0 and site ne 'iphex' then begin
  SpectralWidth[result]=SpectralWidth[result]/10.0
endif

;  Turn missing data to gray
result=where(ReflectivityPossiblePoints eq 0,count)
if count gt 0 then SpectralWidth[result]=missing

;  Find the max and min data values
result=where(SpectralWidth gt missing,count)
if count eq 0 then begin
  dmax=1.0
  dmin=0.0
endif else begin
  dmax=max(SpectralWidth[result])
  dmin=min(SpectralWidth[result])
  if dmin eq dmax then begin
    dmax=dmin+1.0
  endif
endelse

;  Bytscl the data
image=bytscl(SpectralWidth,top=253,max=dmax,min=dmin)

;  Put all no cloud values to white
result=where(SpectralWidth eq nocloud,count)
if count gt 0 then image[result]=255

;  Put all missing data to gray
result=where(SpectralWidth eq missing,count)
if count gt 0 then image[result]=254

;  Set up the plot space
contour,SpectralWidth,julian_day,height,/nodata,$
  xstyle=4,ystyle=4,position=pos[pnum,*],/noerase

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image,sx,sy),px[0],py[0]

;  Put up the axes
contour,SpectralWidth,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',/nodata,$
  ytitle='Height (km)',position=pos[pnum,*],/noerase

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'spectrum width  m/s',$
  color=0,/normal,charsize=tsz

;  Color bar
cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=254,format='(F6.2)',$
  vertical=1,color=0,position=cbpos,charsize=cbsz
  
;**********************************
;  Plot 5 - mixing ratio
;**********************************

;  Plot location
pnum=5

;  Find the valid data
result=where(mxrat gt -8888,count)

;  Find the max and min data values
if count gt 0 then begin
  dmax=max(mxrat[result])
  dmin=min(mxrat[result])
endif

;  Bytscl the data
image=bytscl(mxrat,top=253,max=dmax,min=dmin)

;  Put all no cloud values to white
result=where(mxrat eq nocloud,count)
if count gt 0 then image[result]=255

;  Put all missing data to gray
result=where(mxrat eq missing,count)
if count gt 0 then image[result]=254

;  Set up the plot space
contour,mxrat,julian_day,height,/nodata,/noerase,$
  xstyle=4,ystyle=4,position=pos[pnum,*]

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image,sx,sy),px[0],py[0]

;  Put up the axes
contour,mxrat,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',$
  position=pos[pnum,*],/noerase,/nodata

contour,mxrat,julian_day,height,xstyle=4,ystyle=4,$
  color=255,xtickunits='time',xtickformat='(A1)',$
  position=pos[pnum,*],/noerase,$
  ;levels=[100,115,150,200,300,400,500,700,900,1000],$
  c_labels=[1,1,1,1,1,1,1,1,1,1],c_charsize=1.0

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'mixing ratio g/kg',$
  color=255,/normal,charsize=tsz

;  Color bar
cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=254,format='(F7.2)',divisions=4,$
  vertical=1,color=0,position=cbpos,charsize=cbsz  

;*********************************************
;  Plot 7 - reflectp
;*********************************************

;  Plot location
pnum=7

;  Find the valid data
result=where(ReflectivityPossiblePoints gt -8888,count)

;  Find the max and min data values
if count eq 0 then begin
  dmax=1.0
  dmin=0.0
endif else begin
  dmax=max(ReflectivityPossiblePoints[result])
  dmin=min(ReflectivityPossiblePoints[result])
  if dmin eq dmax then begin
    dmax=dmin+1.0
  endif
endelse

;  Bytscl the data
image=bytscl(ReflectivityPossiblePoints,top=top_color,max=dmax,min=dmin)

;  Put all no cloud values to white
result=where(ReflectivityPossiblePoints eq nocloud,count)
if count gt 0 then image[result]=255

;  Put all missing data to gray
result=where(ReflectivityPossiblePoints eq missing,count)
if count gt 0 then image[result]=254

;  Set up the plot space
contour,ReflectivityPossiblePoints,julian_day,height,/nodata,/noerase,$
  xstyle=4,ystyle=4,position=pos[pnum,*]

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image,sx,sy),px[0],py[0]

;  Put up the axes
contour,ReflectivityPossiblePoints,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',$
  position=pos[pnum,*],/noerase,/nodata

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'# possible points',$
  color=255,/normal,charsize=tsz

;  Color bar
cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=254,format='(F7.2)',divisions=4,$
  vertical=1,color=0,position=cbpos,charsize=cbsz

;*********************************************
;  Plot 9 - reflectc
;*********************************************

;  Plot location
pnum=9

;if site eq 'sgp' then ReflectivityCloudyPoints=$
;  qc_ReflectivityClutterFlagCloudyPoints

;  Find the valid data
result=where(ReflectivityCloudyPoints gt missing,count)

;  Find the max and min data values
if count eq 0 then begin
  dmax=1.0
  dmin=0.0
endif else begin
  dmax=max(ReflectivityCloudyPoints[result])
  dmin=min(ReflectivityCloudyPoints[result])
  if dmin eq dmax then begin
    dmax=dmin+1.0
  endif
endelse

;  Bytscl the data
image=bytscl(ReflectivityCloudyPoints,top=top_color,max=dmax,min=dmin)

;  Put all no cloud values to white
result=where(ReflectivityCloudyPoints eq nocloud,count)
if count gt 0 then image[result]=255

;  Put all missing data to gray
result=where(ReflectivityCloudyPoints eq missing,count)
if count gt 0 then image[result]=254

;  Set up the plot space
contour,ReflectivityCloudyPoints,julian_day,height,/nodata,/noerase,$
  xstyle=4,ystyle=4,position=pos[pnum,*]

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image,sx,sy),px[0],py[0]

;  Put up the axes
contour,ReflectivityCloudyPoints,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',$
  position=pos[pnum,*],/noerase,/nodata

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'# cloudy points',$
  color=255,/normal,charsize=tsz

;  Color bar
cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=254,format='(F7.2)',divisions=4,$
  vertical=1,color=0,position=cbpos,charsize=cbsz

;**********************************************
;  Plot 8 - column cloud fraction
;**********************************************

;  Plot location
pnum=8


;  Plot the data
plot,julian_day,column_cfrac,/noerase,$
  xstyle=1,ystyle=1,position=pos[pnum,*],color=0,psym=-4,$
  yrange=[0,120],xtickunits='time',xtickformat='(A1)',$
  ytitle='%'
oplot,julian_day,column_cfrac_mpl,color=250

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,$
  'column cloud fraction mmr',$
  color=0,/normal,charsize=tsz

xyouts,pos[pnum,0]+0.33,pos[pnum,3]-dt,$
  'mpl',$
  color=250,/normal,charsize=tsz


;**********************************************
;  Plot 6 - cloud fraction
;**********************************************

;  Plot location
pnum=6

;
;  Calculate cloud fraction
;

;  Percentage of points that are cloudy
  data=float(ReflectivityCloudyPoints)/float(ReflectivityPossiblePoints)*100.

;  Where possible points is 0, there is missing data.
;  Where possible points is -8888, there is missing data
;  Give these points the missing data value,-8888,gray
  result=where(ReflectivityPossiblePoints le 0.0 or ReflectivityPossiblePoints eq missing,count)
  print,'ReflectivityPossiblePoints is 0 for ',count,' points'
  if count ne 0 then data[result]=missing

;
;  End of cloud fraction calculation
;

;  Find the valid data
result=where(data gt -8888,count)

;  Find the max and min data values
if count eq 0 then begin
  dmax=1.0
  dmin=0.0
endif else begin
  dmax=max(data[result])
  dmin=min(data[result])
  if dmin eq dmax then begin
    dmax=dmin+1.0
  endif
endelse

;  Bytscl the data
image=bytscl(data,top=top_color,max=dmax,min=dmin)

;  Put all no cloud values to white
result=where(data eq -9999,count)
if count gt 0 then image[result]=255

;  Put all missing data to gray
result=where(data eq -8888,count)
if count gt 0 then image[result]=254

;  Set up the plot space
contour,data,julian_day,height,/nodata,/noerase,$
  xstyle=4,ystyle=4,position=pos[pnum,*]

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image,sx,sy),px[0],py[0]

;  Put up the axes
contour,data,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',$
  position=pos[pnum,*],/noerase,/nodata,ytitle='Height (km)'

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'cloud fraction  %',$
  color=255,/normal,charsize=tsz

;  Color bar
cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=top_color,format='(F7.2)',divisions=4,$
  vertical=1,color=0,position=cbpos,charsize=cbsz

;************************************
;  Plot title and keys
;************************************

;  Add a title
if smonth eq 1 then smon='Jan'
if smonth eq 2 then smon='Feb'
if smonth eq 3 then smon='Mar'
if smonth eq 4 then smon='Apr'
if smonth eq 5 then smon='May'
if smonth eq 6 then smon='Jun'
if smonth eq 7 then smon='Jul'
if smonth eq 8 then smon='Aug'
if smonth eq 9 then smon='Sep'
if smonth eq 10 then smon='Oct'
if smonth eq 11 then smon='Nov'
if smonth eq 12 then smon='Dec'
if site eq 'sgp' then begin
  capsite='SGP'
endif else if site eq 'nsa' then begin
  capsite='NSA'
endif else if site eq 'twpc1' then begin
  capsite='TWP  Manus'
endif else if site eq 'twpc2' then begin
  capsite='TWP  Nauru'
endif else if site eq 'twpc3' then begin
  capsite='TWP  Darwin'
endif else if site eq 'iphex' then begin
  capsite='IPHEX'
endif

str_year=string(syear,format='(I4)')
str_month=string(smonth,format='(I02)')
str_day=string(sday,format='(I02)')

title_str=capsite+' '+avg_int+' '+ret_type+$
  '!C'+str_day+' '+smon+' '+str_year
xyouts,0.5,0.97,align=0.5,charsize=1.7,/normal,color=0,$
title_str

;**********************************
;  Create gif file
;**********************************

if site eq 'iphex' then begin
  imagedir=''
endif else if avg_int eq '1hr' or avg_int eq '6hr' then begin
  imagedir=arm_dir+site+'.average.'+avg_int+'.'+ret_type+'/images_daily/'+str_year+'/'
endif else if avg_int eq '5min' and (ret_type eq 'mine' or ret_type eq 'temp') then begin
  imagedir=arm_dir+site+'.average.5min/images_daily/'+str_year+'/'
endif else if avg_int eq '5min' and (ret_type ne 'mine' and ret_type ne 'temp') then begin
  imagedir=arm_dir+site+'.average.5min.'+ret_type+'/images_daily/'+str_year+'/'
endif else begin
  imagedir=''
endelse

;if ret_type eq 'mine' then begin
  imagename=site+'.average.'+avg_int+'.'+date+'.2d_arm'
;endif else begin
;  imagename=site+'.average.'+avg_int+'.'+ret_type+'.'+date+'.2d_arm'
;endelse

print,imagedir+imagename
write_gif,imagedir+imagename+'.gif',tvrd( )

end