;
;  This program plots the merged sounding for any site
;
pro plot_merged_sounding_nsa,input_flag
;
;  Common to hold data
;
common merged_sounding_data, base_time,time_offset,height,temp, $
	press,mxrat,sflag,numtimes,numheights
;
;  find date and site depending on input flag
;
input_flag=strcompress(input_flag,/remove_all)
if input_flag eq 'one' then begin
  site=' '
  read,site,prompt='Enter ARM site (nsa,sgp,twpc1,twpc2): '
  site=strcompress(site,/remove_all)
  curdate=' '
  read,curdate,prompt='Enter date (yyyymmdd): '
  curdate=strcompress(curdate,/remove_all)
endif else if input_flag eq 'eos' then begin
  fname=''
  openr,1,'/data/mace4/arm_data/eos/eos_netcdf_file.dat'
  readf,1,fname
  close,1
  fname=strcompress(fname,/remove_all)
  check=strmid(fname,0,29)
  if check eq '/data/mace2/arm_data/nsa/nsaC' then begin
    site='nsa'
    curdate=strmid(fname,69,15)
    openw,12,'/data/mace4/arm_data/eos/eos_page_date.dat'
    printf,12,strmid(fname,69,8)
    close,12
  endif else if check eq '/data/mace3/arm_data/sgp/sgpC' then begin
    site='sgp'
    curdate=strmid(fname,69,15)
    openw,12,'/data/mace4/arm_data/eos/eos_page_date.dat'
    printf,12,strmid(fname,69,8)
    close,12
  endif else if check eq '/data/mace/arm_data/twp/Manus' then begin
    site='twpc1'
    curdate=strmid(fname,77,15)
    openw,12,'/data/mace4/arm_data/eos/eos_page_date.dat'
    printf,12,strmid(fname,77,8)
    close,12
  endif else if check eq '/data/mace/arm_data/twp/Nauru' then begin
    site='twpc2'
    curdate=strmid(fname,77,15)
    openw,12,'/data/mace4/arm_data/eos/eos_page_date.dat'
    printf,12,strmid(fname,77,8)
    close,12
  endif
endif  
;
;  Format the date to process
;
print,'date  ',curdate,'site',site
curdate=strcompress(string(curdate),/remove_all)
date_string_fm_adtg, curdate, date_string
;
;  Get the data
;
get_merged_sounding_nsa, curdate,site,datastream,input_flag
;
;  Calculate fractional hour and Julian day
;
ijday_ihrfrac_fm_numsec,base_time,time_offset,hrfrac,jday,yy,mm,dd,hh,mi,ss
;
;  Calculate IDL Julian day for plotting
;
julian_day=make_array(n_elements(jday),/double,value=-9999)
julian_day=julday(mm,dd,yy,hh,mi,ss)
;
;  Calculate relative humidity
;
svp=6.112*exp((17.67*(temp-273.15))/((temp-273.15)+243.5))
smxrat=1000.*(0.622*(svp/((press/100.)-svp)))
rh=100.*(mxrat/smxrat)
result=where( rh gt 100.,count ) 
if count gt 0 then rh[result]=100.
;
;  Calculate potential temperature
;
theta=temp*((press/100000.)^(-287.04/1004.))
;
;  Convert height to km
;
height=height/1000.
;
;  Transpose the data for plotting
;
theta=transpose(theta)
temp=transpose(temp)
rh=transpose(rh)
;
;  Test to make sure there is data
;
tmax=max(temp)
tmin=min(temp)
if tmax eq tmin then goto,no_data_to_plot
;
;  Set up the plot
;
!p.background=!d.n_colors-1
loadct,5 ;39
!p.charsize=2.2
window, 0, xsize=600, ysize=650, TITLE='merged sounding',/pixmap
!p.multi = [0,1,3]
!y.omargin=[4,6]
!y.margin=[1,4]
!x.margin=[8,12]  ;L,R
!x.ticklen=-0.03  ;change the xticks to out of the plot
!y.ticklen=-0.02  ;change the yticks to out of the plot
if input_flag eq 'eos' then $
  dummy=label_date(date_format=['%N/%D/%Z!C%H:%I:%S']) $
else $
  dummy=label_date(date_format=['%H:%I'])
;
;  Plot 1 - Relative Humidity
;
dmax=100.
dmin=0.
;
;  Set up the plot coordinates
;
contour,rh,julian_day,height,/nodata,xstyle=1,ystyle=1,color=0,$
  xtickformat='label_date',xtickunits='time',title='Relative Humidity'
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]
;
;  Scale and create the image
;
image=bytscl(rh,max=dmax,min=dmin)
tv,congrid(image,sx,sy),px[0]+1,py[0]
;
;  Set up the axes
;
axis,julian_day[0],height[0],xaxis=0,xtickformat='label_date',$
  xtickunits='time',xstyle=1,color=0,$
  xrange=[julian_day[0],julian_day[n_elements(julian_day)-1]]
axis,julian_day[0],height[n_elements(height)-1],xaxis=1,$
  xtickformat='(A1)',xtickunits='time',color=0,xstyle=1,$
  xrange=[julian_day[0],julian_day[n_elements(julian_day)-1]]   
axis,julian_day[0],height[0],yaxis=0,ytitle='Height (km)',color=0,$
  ystyle=1,yrange=[height[0],height[n_elements(height)-1]]
axis,julian_day[n_elements(julian_day)-1],height[0],yaxis=1,$
  ytickformat='(A1)',color=0,ystyle=1,$
yrange=[height[0],height[n_elements(height)-1]]
;
;  Add a color bar
; 
colorbar,image,dmax,dmin,'RH (%)',px,py,sx,sy
;
;  Plot 2 - Temperature
;
dmax=330.
dmin=190.
;
;  make no data white color
;
result=where(temp eq -9999,count)
if count gt 0 then temp[result]=1e6
;
;  Set up the plot coordinates
;
contour,temp,julian_day,height,/nodata,xstyle=1,ystyle=1,color=0,$
  xtickformat='label_date',xtickunits='time',title='Temperature'
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]
;
;  Scale and create the image
;
image=bytscl(temp,max=dmax,min=dmin)
tv,congrid(image,sx,sy),px[0]+1,py[0]
;
;  Set up the axes
;
axis,julian_day[0],height[0],xaxis=0,xtickformat='label_date',$
  xtickunits='time',color=0,xstyle=1,$
  xrange=[julian_day[0],julian_day[n_elements(julian_day)-1]]
axis,julian_day[0],height[n_elements(height)-1],xaxis=1,$
  xtickformat='(A1)',xtickunits='time',color=0,xstyle=1,$
  xrange=[julian_day[0],julian_day[n_elements(julian_day)-1]]
axis,julian_day[0],height[0],yaxis=0,ytitle='Height (km)',color=0,$
  yrange=[height[0],height[n_elements(height)-1]],ystyle=1
axis,julian_day[n_elements(julian_day)-1],height[0],yaxis=1,$
  ytickformat='(A1)',color=0,$
  yrange=[height[0],height[n_elements(height)-1]],ystyle=1
;
;  Add a color bar
; 
colorbar,image,dmax,dmin,'Temperature (K)',px,py,sx,sy
;
;  Plot 3 - Potential Temperature
;
dmax=390.
dmin=250.
;
;  Set up the plot coordinates
;
contour,theta,julian_day,height,xstyle=1,ystyle=1,/nodata,color=0,$
  xtickformat='label_date',xtickunits='time',title='Potential Temperature'
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]
;
;  Scale and create the image
;
image=bytscl(theta,max=dmax,min=dmin)
tv,congrid(image,sx,sy),px[0]+1,py[0]
;
;  Set up the axes
;
axis,julian_day[0],height[0],xaxis=0,xtickformat='label_date',$
  xtickunits='time',color=0,xtitle='Time (UTC)',xstyle=1,$
  xrange=[julian_day[0],julian_day[n_elements(julian_day)-1]]
axis,julian_day[0],height[n_elements(height)-1],xaxis=1,$
  xtickformat='(A1)',xtickunits='time',color=0 
  xrange=[julian_day[0],julian_day[n_elements(julian_day)-1]]
axis,julian_day[0],height[0],yaxis=0,ytitle='Height (km)',color=0,$
  ystyle=1,yrange=[height[0],height[n_elements(height)-1]]
axis,julian_day[n_elements(julian_day)-1],height[0],yaxis=1,$
  ytickformat='(A1)',color=0,$
  ystyle=1,yrange=[height[0],height[n_elements(height)-1]]
;
;  Add a color bar
; 
colorbar,image,dmax,dmin,'Theta (K)',px,py,sx,sy


;
;  Add title to plot page
;
if site eq 'twpc1' then capsite='TWP C1'
if site eq 'twpc2' then capsite='TWP C2'
if site eq 'nsa' then capsite='NSA'
if site eq 'sgp' then capsite='SGP'
a_string=capsite+' Merged Sounding, '+date_string
xyouts, 0.5, 0.94, alignment=0.5,a_string, /normal, color=0, charsize=2.2
xyouts, 0.5, 0.89, alignment=0.5,datastream, /normal, color=0, charsize=1.9
;
;  Create the gif file
;
if input_flag eq 'eos' then begin

  gif_file='/web/mace/homepages/research/eos/images/'+$
    site+'/'+datastream+'.'+curdate
  write_gif,gif_file+'.gif',tvrd()
  unix_command='rm '+gif_file+'.png'
  spawn,unix_command

  openw,13,'/data/mace4/arm_data/eos/eos_image_name.dat'
  printf,13,datastream+'.'+curdate+'.gif'
  close,13

endif else $
  write_gif,'/data/mace4/arm_data/temp_image/'+$
  datastream+'.'+curdate+'.gif',tvrd()

;********************
;  Plot the Merged Sounding Status Flag
;********************
;
;  Set up plot
;
window,0,xsize=500,ysize=700,title='Merged Sounding Status Flag',/pixmap
!p.multi=[0,1,1]
!x.ticklen=0.02
!y.ticklen=0.02
if site eq 'twpc1' then capsite='TWP C1'
if site eq 'twpc2' then capsite='TWP C2'
if site eq 'nsa' then capsite='NSA'
if site eq 'sgp' then capsite='SGP'
if input_flag eq 'eos' then $
  dummy=label_date(date_format=['%N/%D/%Z!C%H:%I:%S']) $
else $
  dummy=label_date(date_format=['%H:%I'])
;
;  Determine how many sounding flag categories depending on 
;  site
;
if site eq 'sgp' then begin
  smax=17
  smin=1
endif else if site eq 'nsa' then begin
  smax=9
  smin=1
endif else if site eq 'twpc1' or site eq 'twpc2' then begin
  smax=7
  smin=1
endif 
;
;  Plot the sounding flag data
;
plot, julian_day, sflag, $
  background=!d.n_colors-1,color=0,charsize=1.3, $
  ;xtitle="Time (UTC)", $
  xrange=[julian_day[0],julian_day[n_elements(julian_day)-1]],$
  xstyle=1,psym=6,symsize=0.6,xtickunits='time',xtickformat='label_date',$
  yrange=[smin,smax],ystyle=1,yticks=smax-1,$
  position=[0.1,0.5,0.9,0.9],/normal

xyouts,0.5,0.95,alignment=0.5,capsite+'  '+date_string,/normal,color=0,charsize=2

a_string='Merged Sounding Status Flag'
xyouts,0.5,0.92,alignment=0.5,a_string,/normal,color=0,charsize=2

a_string='Time (UTC)'
xyouts,0.5,0.42,alignment=0.5,a_string,/normal,color=0,charsize=1.3

dy=0.42/float(smax+1)

if site eq 'twpc1' or site eq 'twpc2' or site eq 'nsa' then begin
  a_string=' 1. Sonde, MWR vapor'
endif else if site eq 'sgp' then begin
  a_string=' 1. CART Sonde, MWR vapor, AERI profile'
endif
;xyouts,0.3,0.28,alignment=0.0,a_string,/normal,color=0,charsize=1.25
xyouts,0.2,1*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2

if site eq 'twpc1' or site eq 'twpc2' or site eq 'nsa' then begin
  a_string=' 2. Sonde'
endif else if site eq 'sgp' then begin
  a_string=' 2. CART Sonde, MWR vapor'
endif
;xyouts,0.3,0.25,alignment=0.0,a_string,/normal,color=0,charsize=1.25
xyouts,0.2,2*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2

if site eq 'twpc1' or site eq 'twpc2' then begin
  a_string=' 3. ECMWF, MWR vapor'
endif else if site eq 'nsa' then begin
  a_string=' 3. Barrow NWS Sonde, MWR vapor'
endif else if site eq 'sgp' then begin
  a_string=' 3. CART Sonde, AERI profile'
endif
;xyouts,0.3,0.22,alignment=0.0,a_string,/normal,color=0,charsize=1.25
xyouts,0.2,3*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2

if site eq 'twpc1' or site eq 'twpc2' then begin
  a_string=' 4. ECMWF'
endif else if site eq 'nsa' then begin
  a_string=' 4. Barrow NWS Sonde'
endif else if site eq 'sgp' then begin
  a_string=' 4. CART Sonde'
endif
;xyouts,0.3,0.19,alignment=0.0,a_string,/normal,color=0,charsize=1.25
xyouts,0.2,4*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2

if site eq 'twpc1' or site eq 'twpc2' then begin
  a_string=' 5. Standard Atmosphere, MWR vapor'
endif else if site eq 'nsa' then begin
  a_string=' 5. NCEP MOLTS Sounding, MWR vapor'
endif else if site eq 'sgp' then begin
  a_string=' 5. OUN Sonde, MWR vapor, AERI profile'
endif
;xyouts,0.3,0.16,alignment=0.0,a_string,/normal,color=0,charsize=1.25
xyouts,0.2,5*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2

if site eq 'twpc1' or site eq 'twpc2' then begin
  a_string=' 6. Standard Atmosphere'
endif else if site eq 'nsa' then begin
  a_string=' 6. NCEP MOLTS Sounding'
endif else if site eq 'sgp' then begin
  a_string=' 6. OUN Sonde, MWR vapor'
endif
;xyouts,0.3,0.13,alignment=0.0,a_string,/normal,color=0,charsize=1.25
xyouts,0.2,6*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2

if site eq 'twpc1' or site eq 'twpc2' then begin
  a_string=' 9. No Sounding'
endif else if site eq 'nsa' then begin
  a_string=' 7. ECMWF Sounding, MWR vapor'
endif else if site eq 'sgp' then begin
  a_string=' 7. OUN Sonde, AERI profile'
endif
;xyouts,0.3,0.10,alignment=0.0,a_string,/normal,color=0,charsize=1.25
xyouts,0.2,7*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2

if site eq 'nsa' then begin
  a_string=' 8. ECMWF Sounding'
endif else if site eq 'sgp' then begin
  a_string=' 8. MOLTS, MWR vapor, AERI profile'
endif
if site eq 'nsa' or site eq 'sgp' then $
;xyouts,0.3,0.07,alignment=0.0,a_string,/normal,color=0,charsize=1.25
xyouts,0.2,8*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2

if site eq 'nsa' then begin
  a_string=' 9. No Sounding'
endif else if site eq 'sgp' then begin
  a_string=' 9. MOLTS, MWR Vapor'
endif
if site eq 'nsa' or site eq 'sgp' then $
;xyouts,0.3,0.04,alignment=0.0,a_string,/normal,color=0,charsize=1.25
xyouts,0.2,9*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2

if site eq 'sgp' then begin
  a_string='10. MOLTS, AERI Profile'
  xyouts,0.2,10*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2
endif

if site eq 'sgp' then begin
  a_string='11. ECMWF, MWR vapor, AERI profile'
  xyouts,0.2,11*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.2
endif

if site eq 'sgp' then begin
  a_string='12. ECMWF, MWR vapor'
  xyouts,0.2,12*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.25
endif

if site eq 'sgp' then begin
  a_string='13. ECMWF, AERI profile'
  xyouts,0.2,13*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.25
endif

if site eq 'sgp' then begin
  a_string='14. RUC, MWR vapor, AERI profile'
  xyouts,0.2,14*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.25
endif

if site eq 'sgp' then begin
  a_string='15. RUC, MWR vapor'
  xyouts,0.2,15*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.25
endif

if site eq 'sgp' then begin
  a_string='17. No Sounding'
  xyouts,0.2,16*dy,alignment=0.0,a_string,/normal,color=0,charsize=1.25
endif

;
;  Create gif file
;
if input_flag eq 'eos' then begin
  gif_file='/web/mace/homepages/research/eos/images/'+$
    site+'/'+site+'.'+curdate+'.sounding_flag'
  write_gif,gif_file+'.gif',tvrd()
  unix_command='rm '+gif_file+'.png'
  spawn,unix_command  

  openw,13,'/data/mace4/arm_data/eos/eos_image_name.dat'
  printf,13,site+'.'+curdate+'.sounding_flag'+'.gif'
  close,13

endif else $
write_gif,'/data/mace4/arm_data/temp_image/'+$
  site+'.'+curdate+'.sounding_flag.gif',tvrd()

no_data_to_plot:print,'no data to plot'
end


;****************************************************
;  colorbar.pro creates a colorbar for the DBZ images
;****************************************************

pro colorbar, image,dmax,dmin,cb_title,px,py,sx,sy

min_image=0;min(image)
max_image=255;max(image)

cb_arr=fltarr(2,max_image-min_image+1)
for j=0,max_image-min_image do $
  cb_arr(*,j)=dmin+((dmax-dmin)/(max_image-min_image))*j

cb_img=bindgen(2,max_image-min_image+1)

for j=0,max_image-min_image do cb_img(*,j)=min_image+j

tv, congrid(cb_img,(sx/35),sy),px[1]+sx/8,py[0]
axis,px[1]+sx/8,py[0],yaxis=0,color=0,$
  yrange=[cb_arr[0,0],cb_arr[0,max_image-min_image]],$
  ystyle=1,/device,charsize=2.0,ytitle=cb_title,ticklen=-0.02

return 
end