;
;  This program plots a daily summary of the cloud retrieval 
;  2d variables. 
;
pro plot_average_daily_2d_cloud2,site,date,avg_int,ret_type

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

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

;  Date in the form yyyymmdd 
;date='20010511'

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

;  Retrieval version
;ret_type='mine'
;ret_type='mine'

;**************************************
;  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_grp/data/arm/'+site+'/'

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

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

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

;  Variables to get from the file
varname=['base_time','time_offset','height',$
         'ice_tau_solar',$
         'IWC_best_estimate',$
	 'iwc_source',$
         'liquid_tau_solar',$
	 'lwc_best_estimate',$
	 'lwc_source',$
         're_ice_best_estimate',$
         're_liquid_best_estimate']

if m_or_d eq 'month' then begin
  searchdate=strmid(date,0,6)
endif else searchdate=date

;  Get the average data
get_average_data_list,varname,data,$
  site,searchdate,avg_int,ret_type,m_or_d,fill_flag

;  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  

;  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
;*******************************************

;  Array of Julian_days
julian_day=double(day1+(base_time+time_offset)/sec_per_day)

;  Subset the data if this is a monthly file
if m_or_d eq 'month' then begin
  
  sjday=julday(strmid(date[0],4,2),strmid(date[0],6,2),$
       strmid(date[0],0,4),0,0,0)
  ejday=sjday+1
  caldat,sjday,emm,edd,eyyyy,ehour,eminute,esecond
  print,emm,edd,eyyyy,ehour,eminute,esecond

  result=where(julian_day ge sjday and julian_day le ejday,count)
  julian_day=julian_day[result]
  ice_tau_solar=ice_tau_solar[result,*]
  IWC_best_estimate=IWC_best_estimate[result,*]
  iwc_source=iwc_source[result,*]
  liquid_tau_solar=liquid_tau_solar[result,*]
  lwc_best_estimate=lwc_best_estimate[result,*]
  lwc_source=lwc_source[result,*]
  re_ice_best_estimate=re_ice_best_estimate[result,*]
  re_liquid_best_estimate=re_liquid_best_estimate[result,*] 

  numtimes=n_elements(julian_day)
endif  

;  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]

;*************************************************
;  Calculate Optical depth
;*************************************************

;  Get the dimensions of the arrays

s=size(ice_tau_solar,/dimensions)
n_times=s[0]
n_height=s[1]

;  Array to hold the total optical depth
total_ice_tau_solar=make_array(n_times,/float,value=0)
total_liquid_tau_solar=make_array(n_times,/float,value=0)

for i=0,n_times-1 do begin

  temp_ice=reform(ice_tau_solar[i,*])
  temp_liq=reform(liquid_tau_solar[i,*])

  result=where(temp_ice gt -8888,count_ice)
  if count_ice gt 0 then begin
    total_ice_tau_solar[i]=total(temp_ice[result])
  endif

  result=where(temp_liq gt -8888,count_liq)
  if count_liq gt 0 then begin
    total_liquid_tau_solar[i]=total(temp_liq[result])
  endif
endfor

ice_od=total_ice_tau_solar
liq_od=total_liquid_tau_solar
total_od=ice_od+liq_od

;*********************************************
;  Calculate total water
;  Add lwc and iwc
;*********************************************

missing=-8888
temp_iwc=IWC_best_estimate
result=where(temp_iwc le missing,count)
if count gt 0 then temp_iwc[result]=0
temp_lwc=lwc_best_estimate
result=where(temp_lwc le missing,count)
if count gt 0 then temp_lwc[result]=0
total_water=temp_iwc+temp_lwc
result=where(total_water eq 0,count)
total_water[result]=-9999
result=where(IWC_best_estimate eq missing,count)
if count gt 0 then total_water[result]=missing
wat_max=0.9
wat_min=-6.0

;********************************
;  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   

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

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

;********************************
;  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
!p.charsize=0.8

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

xl=0.05  &  xr=0.94
yb=0.05  &  yt=0.93
sx=0.08        
sy=0.00 
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 - ice_tau_solar
;**********************************

;  Plot location
pnum=0

;  Find the valid data
result=where(ice_tau_solar 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(ice_tau_solar[result])
  dmin=min(ice_tau_solar[result])

;  If there is one really large value, find a lower dmax
  gdata=ice_tau_solar[result]
  rlt=where(gdata lt 0.90*dmax,cnt)
  while float(cnt)/float(count) gt .95 do begin
    rlt=where(gdata lt 0.90*dmax,cnt)
    dmax=max(gdata[rlt])
  endwhile

  if dmin eq dmax then begin
    dmax=dmin+1.0
  endif
endelse

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

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

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

;  Set up the plot space
contour,ice_tau_solar,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,ice_tau_solar,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,'ice_tau_solar',$
  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 1 - liquid_tau_solar
;**********************************

;  Plot location
pnum=1

;  Find the valid data
result=where(liquid_tau_solar 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(liquid_tau_solar[result])
  dmin=min(liquid_tau_solar[result])

;  If there is one really large value, find a lower dmax
  gdata=liquid_tau_solar[result]
  rlt=where(gdata lt 0.90*dmax,cnt)
  while float(cnt)/float(count) gt .95 do begin
    rlt=where(gdata lt 0.90*dmax,cnt)
    dmax=max(gdata[rlt])
  endwhile

  if dmax eq dmin then begin
    dmax=1.0+dmin
  endif
endelse

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

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

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

;  Set up the plot space
contour,liquid_tau_solar,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,liquid_tau_solar,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='label_date',/nodata,$
  position=pos[pnum,*],/noerase

xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'liquid_tau_solar',$
  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 2 - Ice optical depth
;**********************************

;  Plot location
pnum=2

;  Restrict to good data
result=where(total_od gt 0,count_total)
if count_total gt 0 then begin
  total_od_total=total_od[result]
  julian_day_total=julian_day[result]

  result=where(ice_od gt 0,count_ice)
  if count_ice gt 0 then begin
    ice_od=ice_od[result]
    julian_day_ice=julian_day[result]
  endif

;  Put up the axes
  plot,julian_day_total,total_od_total,xstyle=1,ystyle=1,$
    color=0,xtickunits='time',xtickformat='(A1)',$
    position=pos[pnum,*],/noerase,yrange=[0.1,120.],/ylog,$
    xrange=[julian_day[0],julian_day[n_elements(julian_day)-1]]
  if count_ice gt 0 then oplot,julian_day_ice,ice_od,color=250,psym=4,symsize=0.4

;  Title
  xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'total optical depth',$
    color=0,/normal,charsize=tsz
  xyouts,pos[pnum,0]+0.21,pos[pnum,3]-dt,'ice optical depth',$
    color=250,/normal,charsize=tsz
endif

;**********************************
;  Plot 3 - Liquid optical depth
;**********************************

;  Plot location
pnum=3

;  Restrict to good data
result=where(total_od gt 0,count_total)
if count_total gt 0 then begin 
  total_od_total=total_od[result]
  julian_day_total=julian_day[result]

  result=where(liq_od gt 0,count_liq)
  if count_liq gt 0 then begin
    liq_od=liq_od[result]
    julian_day_liq=julian_day[result]
  endif

;  Put up the axes
  plot,julian_day_total,total_od_total,xstyle=1,ystyle=1,$
    color=0,xtickunits='time',xtickformat='(A1)',$
    position=pos[pnum,*],/noerase,yrange=[0.1,120.],/ylog,$
    xrange=[julian_day[0],julian_day[n_elements(julian_day)-1]]
  if count_liq gt 0 then oplot,julian_day_liq,liq_od,color=250,psym=4,symsize=0.4

;  Title
  xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'total optical depth',$
    color=0,/normal,charsize=tsz
  xyouts,pos[pnum,0]+0.21,pos[pnum,3]-dt,'liquid optical depth',$
    color=250,/normal,charsize=tsz
endif

;*********************************************
;  Plot 4 - re_ice_best_estimate
;*********************************************

do_plot4='yes'
if do_plot4 eq 'yes' then begin

;  Plot location
pnum=4

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

;  Convert from m to um
if count gt 0 then begin
  re_ice_best_estimate[result]=re_ice_best_estimate[result]*1e6  
endif

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

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

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

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

;  Set up the plot space
contour,re_ice_best_estimate,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,re_ice_best_estimate,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,'re_ice_best_estimate  um',$
  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,$
  ;title='counts',$
  ncolors=254,format='(F20.1)',$
  vertical=1,color=0,position=cbpos,charsize=cbsz;,$
  ;divisions=dmax
endif  ;end of do_plot4

;*********************************************
;  Plot 5 - re_liquid_best_estimate
;*********************************************

do_plot5='yes'
if do_plot5 eq 'yes' then begin

;  Plot location
pnum=5

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

;  Convert from m to um
if count gt 0 then begin
  re_liquid_best_estimate[result]=re_liquid_best_estimate[result]*1e6  
endif

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

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

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

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

;  Set up the plot space
contour,re_liquid_best_estimate,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,re_liquid_best_estimate,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',/nodata,$
  ;ytitle='Height (km)',$
  position=pos[pnum,*],/noerase

;  Title
xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'re_liquid_best_estimate  um',$
  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,$
  ;title='counts',$
  ncolors=254,format='(F20.1)',$
  vertical=1,color=0,position=cbpos,charsize=cbsz;,$
  ;divisions=dmax
endif

;**********************************************
;  Plot 6 - IWC_best_estimate
;**********************************************

;  Plot location
pnum=6

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

;  Take the log of the valid data
if count gt 0 then begin
  IWC_best_estimate[result]=alog10(IWC_best_estimate[result]) 
endif

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

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

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

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

;  Set up the plot space
contour,IWC_best_estimate,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,IWC_best_estimate,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',/nodata,$
  ytitle='Height (km)',position=pos[pnum,*],/noerase

;  Title
xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'IWC_best_estimate  log(g/m^3)',$
  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,$
  ;title='counts',$
  ncolors=254,format='(F6.2)',$
  vertical=1,color=0,position=cbpos,charsize=cbsz;,$
  ;divisions=dmax

;*********************************************
;  Plot 7 - lwc_best_estimate
;*********************************************

;  Plot location
pnum=7

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

;  Take the log of the valid data
if count gt 0 then begin
  lwc_best_estimate[result]=alog10(lwc_best_estimate[result]) 
endif

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

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

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

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

;  Set up the plot space
contour,lwc_best_estimate,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,lwc_best_estimate,julian_day,height,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='(A1)',/nodata,$
  ;ytitle='Height (km)',$
  position=pos[pnum,*],/noerase

;  Title
xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'lwc_best_estimate  log(g/m^3)',$
  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,$
  ;title='counts',$
  ncolors=254,format='(F6.2)',$
  vertical=1,color=0,position=cbpos,charsize=cbsz;,$
  ;divisions=dmax

;**********************************************
;  Plot 8 - iwc_source
;**********************************************

;  Plot location
pnum=8

;  Make an image array
image=iwc_source

;  Color the image array depending on the source 
dmax=8
dmin=1
source_colors=[0,indgen(8)*35]
for source=dmin,dmax do begin
  result=where(iwc_source eq source,count)
  if count gt 0 then image[result]=source_colors[source]
endfor

;  In the 5min data, no cloud has a source of 0.  put this to white.
if avg_int eq '5min' then begin
  result=where(iwc_source eq 0,count)
  if count gt 0 then image[result]=255
endif

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

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

;  Set up the plot space
contour,iwc_source,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,iwc_source,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,'iwc_source',$
  color=0,/normal,charsize=tsz

;**********************************************
;  Plot 9 - lwc_source
;**********************************************

do_lwc_source='no'
if do_lwc_source eq 'yes' then begin

;  Plot location
pnum=9

;  Make an image array
image=lwc_source

;  Color the image array depending on the source 
dmax=8
dmin=1
source_colors=[0,indgen(8)*35]
for source=dmin,dmax do begin
  result=where(lwc_source eq source,count)
  if count gt 0 then image[result]=source_colors[source]
endfor

;  In the 5min data, no cloud has a source of 0.  put this to white.
if avg_int eq '5min' then begin
  result=where(lwc_source eq 0,count)
  if count gt 0 then image[result]=255
endif

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

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

;  Set up the plot space
contour,lwc_source,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,lwc_source,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,'lwc_source',$
  color=0,/normal,charsize=tsz
endif

;**********************************************
;  Plot 9 - total water
;**********************************************

do_total_water='yes'
if do_total_water eq 'yes' then begin

;  Plot location
pnum=9

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

;  Take the log of the valid data
if count gt 0 then begin
  total_water[result]=alog10(total_water[result]) 
endif

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

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

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

;  Put all missing data to gray
result=where(total_water eq missing,count)
if count gt 0 then image[result]=254
;  Set up the plot space

contour,total_water,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,total_water,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,'total_water',$
  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,$
  ;title='counts',$
  ncolors=254,format='(F6.2)',$
  vertical=1,color=0,position=cbpos,charsize=cbsz;,$
  ;divisions=dmax


endif

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

;  Create a code for the ice sources
dcode=0.013
ccode=0.8
xyouts,xl,1.0-dcode,'1=mace et al 1998',/normal,color=source_colors[1],charsize=ccode 
xyouts,xl,1.0-2.0*dcode,'2=mace et al 2001',/normal,color=source_colors[2],charsize=ccode 
xyouts,xl,1.0-3.0*dcode,'3=ccm3 Parameterization (Kiehl et al 1998)',/normal,color=source_colors[3],charsize=ccode
xyouts,xl,1.0-4.0*dcode,'4=Liu and Illingworth, 2000',/normal,color=source_colors[4],charsize=ccode
;  Create a code for the liquid sources
dcode=0.013
ccode=0.8
xc=pos[1,0]
xyouts,xc,1.0-dcode,'5=Dong et al 1998 parameterization',/normal,color=source_colors[5],charsize=ccode
xyouts,xc,1.0-2.0*dcode,'6=Dong and Mace 2001 retrieval',/normal,color=source_colors[6],charsize=ccode
xyouts,xc,1.0-3.0*dcode,'7=frisch et al 1998 lwc retrieval',/normal,color=source_colors[7],charsize=ccode
xyouts,xc,1.0-4.0*dcode,'8=mwr physical parameterization (liq) '+$
  'or aircraft empirical (ice)',/normal,color=source_colors[8],charsize=ccode 

;  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='TWPC1'
endif else if site eq 'twpc2' then begin
  capsite='TWPC2'
endif else if site eq 'twpc3' then begin
  capsite='TWPC3'
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.43,0.97,align=0.5,charsize=1.5,/normal,color=0,$
title_str

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

;  Image directory
if avg_int eq '5min' and ret_type eq 'column' then begin
  imagedir='/data/mace20/arm_data/sgp/sgp.Integrate_Column_5min-Mace/images_daily/'
endif else if avg_int ne '5min' and ret_type eq 'column' then begin
  imagedir=arm_dir+site+'.Integrate_Column_'+avg_int+'-Mace/images_daily/'
endif else if avg_int ne '5min' 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' then begin
  imagedir=arm_dir+site+'.average.5min/images_daily/'+str_year+'/'
endif else if avg_int eq '5min' and ret_type ne 'mine' 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_cloud'
endif else begin
  imagename=site+'.average.'+avg_int+'.'+ret_type+'.'+date+'.2d_cloud'
endelse

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

end