openr, 3, '/project/p1162/jcohen/snowdepth/ims_lsmask_nogland.txt'
openr, 5, '/project/p1162/jcohen/snowdepth/cmc_analysis_ps_lat_long.txt' 
openr, 6, '/project/p1162/jcohen/snowdepth/cmc_analysis_2009.txt' 

x=706
y=706
n=365

dateid = '2009020100'

sd = fltarr(x,y)
tmp = fltarr(x)
lon=fltarr(x,y)

lat=fltarr(x,y)
mask=intarr(x,y)
fltr=intarr(1024,1024)
ftmp=intarr(1024)
hdrtext=""

; read land/sea mask
for i=0,1024-1 do begin
  readf, 3, ftmp, format = '(1024i1.1)'
  fltr(i,*) = ftmp
endfor
close, 3
print, "finished reading land/sea mask"

; find the x by y array below the specified dateid
tmpdate = strarr(1)
readf, 6, tmpdate
tmpdate = strcompress(tmpdate, /REMOVE_ALL)
while(dateid ne tmpdate) do begin 
  skip = strarr(1)
  for i=0,x-1 do begin
    readf, 6, skip
  endfor
  tmpdate = strarr(1)
  readf, 6, tmpdate
  tmpdate = strcompress(tmpdate, /REMOVE_ALL)
endwhile
for i=0,x-1 do begin
  readf, 6, tmp, format = '(706f6.1)'
  sd(i,*) = tmp
endfor

close, 6
print, "finished reading snow depth"

; read lats/lons of polar stereographic grid
for k=0,8 do begin
  readf, 5, hdrtext, format='(a20)'
endfor

for i=0,x-1 do begin
  for j=0,y-1 do begin
    readf, 5, ii, jj, tmplat, tmplon, format = '(2(9x,i3),2(3x,G13.6))'
    lat(FIX(ii)-1, FIX(jj)-1)=tmplat
    lon(FIX(ii)-1, FIX(jj)-1)=tmplon
  endfor
endfor
close, 5

; debugging odd values over Lake Superior 
for i=190,390 do begin
  for j=75,275 do begin
;   print, "i j lat lon sd = ",i, " ", j, " ", lat(i,j), " ", lon(i,j), " ", sd(i,j)
  endfor
endfor

; use only the portion of the land/sea mask that corresponds to the snow depth grid
for i=0,x-1 do begin
  for j=0,y-1 do begin
     mask(i,j) = fltr(i+159,j+159)
   endfor
endfor

; set snow depths over water (using land/sea mask) to zero
for i=0,x-1 do begin
  for j=0,y-1 do begin
     if mask(i,j) eq 0 then sd(i,j)=0.
   endfor
endfor

; scale sdow depth values for plotting
                 sd_colors = 0.25*(sd(*,*))
    
                 idx = where(sd_colors ge 1.0 and sd_colors lt 2.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 0
                 idx = where(sd_colors ge 2.0 and sd_colors lt 3.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 1
                 idx = where(sd_colors ge 3.0 and sd_colors lt 4.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 2
                 idx = where(sd_colors ge 4.0 and sd_colors lt 5.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 3
                 idx = where(sd_colors ge 5.0 and sd_colors lt 10.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 4
                 idx = where(sd_colors ge 10.0 and sd_colors lt 20.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 5
                 idx = where(sd_colors ge 20.0 and sd_colors lt 30.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 6
                 idx = where(sd_colors ge 30.0 and sd_colors lt 40.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 7
                 idx = where(sd_colors ge 40.0 and sd_colors lt 60.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 8
                 idx = where(sd_colors ge 60.0 and sd_colors lt 80.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 9
                 idx = where(sd_colors ge 80.0 and sd_colors lt 100.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 10
                 idx = where(sd_colors ge 100.0 and sd_colors lt 200.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 11
                 idx = where(sd_colors ge 200.0 and sd_colors lt 400.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 12
                 idx = where(sd_colors ge 400.0 and sd_colors lt 600.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 13
                 idx = where(sd_colors ge 600.0 and sd_colors lt 800.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 14
                 idx = where(sd_colors ge 800.0 and sd_colors lt 1000.0)
                 if (size(idx, /N_DIM) gt 0) then sd_colors(idx) = 15
                 sd_colors = LONG(sd_colors)
                 offset=64
                 sd_colors=sd_colors+64 
                 print, "max(sd_colors) =", max(sd_colors)
                 print, "min(sd_colors) =", min(sd_colors)

; debugging odd values over Lake Superior 
for i=190,390 do begin
  for j=75,275 do begin
;   print, "i j lat lon sd_colors = ",i, " ", j, " ", lat(i,j), " ", lon(i,j), " ", sd_colors(i,j)
  endfor
endfor

print, "Snow depth min max = ", min(sd(*,*)), " ", max(sd(*,*))

lats = [30, 40, 50,60]
lons = [-125,-120,-110,-100,-90,-80,-70]


c1 = [1,2,3,4,5,10,20,30,40,60,80,100,200,400,600,800,1000]
c2 = [1,2,3,4,5,10,50,100,500,1000]
SET_PLOT, 'PS'
DEVICE, /COLOR,                 $
        FILENAME='./tmpfile.ps',        $
        /HELVETICA,             $
        /BOLD,                  $
        /INCHES,                $
        XSIZE=8.5,               $
        YSIZE=11,              $
        SCALE_FACTOR=1.0,       $
        XOFFSET=0,              $
        YOFFSET=0

!P.FONT   = 0
!P.REGION = [ 0.9, 0.9, 0.9, 0.9 ]
!X.MARGIN = [ 10, 10 ]
!Y.MARGIN = [ 30, 10 ]

!P.CHARSIZE = 1.25

CTload, 40

map_set, 90, -80, -20, /stereographic, /continent, /USA, /ISOTROPIC,$
limit=[32., -122., 55, -97., 32., -72., 27., -97.], $
title='!17 '+ 'Snow Depth '+ dateid, $
/noerase, mlinethick=3
map_grid, LATS=lats, LONS=lons

; map extents for N. Hemisphere plot
;limit=[19.6282, -170., 19.475, 100., 19.475, 10., 19.6282, -80.], $

index = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
FOR j=0,15 DO CTLoad, index[j], NCOLORS=16, BOTTOM=16*j, /BREWER

contour, sd(*,*), lon, lat, levels=c1,/OVERPLOT,/FILL,$
      c_colors=[0+offset, 1+offset, 2+offset, 3+offset,4+offset,5+offset,6+offset,7+offset,8+offset,9+offset,10+offset,11+offset,12+offset,13+offset,14+offset,15+offset]
contour, sd(*,*), lon, lat, levels=c2,/OVERPLOT,c_colors=15+offset

CTload, 40
plot,[0.0,1.0],[0.0,0.0],                     $
   position=[0.1,0.05,.9,.1],               $
   ystyle=4,                                $
   xstyle=1, xrange=[0,16],                 $
   xtickv=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16], $
   xticks=16,              $
   xtickname=['1','2','3','4','5','10','20','30','40','60','80','100','200','400','600','800','1000'], $
   /noerase

index = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
FOR j=0,15 DO CTLoad, index[j], NCOLORS=16, BOTTOM=16*j, /BREWER

polyfill,[0,0,1,1],[0,1,1,0],color=0+offset
polyfill,[1,1,2,2],[0,1,1,0],color=1+offset
polyfill,[2,2,3,3],[0,1,1,0],color=2+offset
polyfill,[3,3,4,4],[0,1,1,0],color=3+offset
polyfill,[4,4,5,5],[0,1,1,0],color=4+offset
polyfill,[5,5,6,6],[0,1,1,0],color=5+offset
polyfill,[6,6,7,7],[0,1,1,0],color=6+offset
polyfill,[7,7,8,8],[0,1,1,0],color=7+offset
polyfill,[8,8,9,9],[0,1,1,0],color=8+offset
polyfill,[9,9,10,10],[0,1,1,0],color=9+offset
polyfill,[10,10,11,11],[0,1,1,0],color=10+offset
polyfill,[11,11,12,12],[0,1,1,0],color=11+offset
polyfill,[12,12,13,13],[0,1,1,0],color=12+offset
polyfill,[13,13,14,14],[0,1,1,0],color=13+offset
polyfill,[14,14,15,15],[0,1,1,0],color=14+offset
polyfill,[15,15,16,16],[0,1,1,0],color=15+offset

CTload, 40
map_set, 90, -80, -20, /stereographic, /continent, /USA, /ISOTROPIC,$
limit=[32., -122., 55, -97., 32., -72., 27., -97.], $
title=' ', $
/noerase, mlinethick=3
map_grid, LATS=lats, LONS=lons, /NOERASE

device,/close_file
end