PRO hw4
; Read data.  File names
ufn = '/net/b/klotter/uwnd.mon.ltm.nc'

LatMin = -45. ; negative for South
LatMax = 15.
LonMin = -90. ; negative for West
LonMax = -30.

Pres = 700.  ; milibar level used
Month = 5 ; June

map_limit = [LatMin, LonMin, LatMax, LonMax]; for maps

xx_level = INDGEN(21)*1 -10.
CU_level = xx_level
CU_level = [-9999., CU_level, 9999]
CU_color =  ['ffffff'x, '000088'x, '0000aa'x, '0000cc'x, '0000ff'x, $
  '0055ff'x, '0088ff'X, '00aaff'x, '00ccff'x, '00ffff'x, '11ffdd'x, '44ffaa'x, $
  '66ee44'x, '99ff88'x, 'ddff22'x, 'ffff00'x, 'ffcc00'x, 'ffaa00'X, $
  'ff8800'x, 'ff5500'x, 'ff2200'X, '555555'x, '111111'x ]
CU_color = REVERSE(CU_color)
CU_thick = FLTARR(N_ELEMENTS(CU_level)) +.001
y = WHERE(xx_level MOD 5 EQ 0, amt) + 1
CU_line = CU_level[y]
CU_style = INTARR(amt)+1
i = WHERE(CU_line EQ 0)
CU_style[i] = 2

; Uwind data read
lun_uwnd = NCDF_OPEN(ufn)
; note that it is assumed that the files start at Jan 1948

; Read data.  Assume that the lon lat time and level are same for both U/V
NCDF_VARGET, lun_uwnd, "lon", lon
NCDF_VARGET, lun_uwnd, "lat", lat
NCDF_VARGET, lun_uwnd, "time", time
NCDF_VARGET, lun_uwnd, "level", level

NCDF_VARGET, lun_uwnd, "uwnd", uwind
; assuming that SF=1, add_offset=0
NCDF_CLOSE, lun_uwnd

; selection parameters used for plot.  limit to just those of interest here
; LONGITUDES.  Complicated due to wrap-around for Prime Meridian
lon_ind = [ WHERE(lon GE LonMin+360 AND lon LE LonMax+360), $
	WHERE(lon GE LonMin AND lon LE LonMax) ] ; index of longitudes used
check = WHERE(lon_ind GE 0, amt)
lon_ind = lon_ind[check]
lon_used = lon[lon_ind]
lon_ind2 = WHERE(lon_used GT 180, amt) ; wrap around longitudes to plot over PM
IF amt GT 0 THEN lon_used[lon_ind2] -= 360.  ; value of lon's used
n_lon = N_ELEMENTS(lon_used)

lat_ind = WHERE(lat GE LatMin AND lat LE LatMax, n_lat)
lat_used = lat[lat_ind]

prs_ind = WHERE(Pres EQ level, n_prs)

time_ind = Month

uwind_used = FLTARR(n_lon, n_lat, n_prs, 1)
uwind_used[*,*,*,*] = uwind[lon_ind[*], lat_ind[*], prs_ind[*], time_ind]
uwind_used = REFORM(uwind_used)

; test print
PRINT, 'ORIGINAL SIZE ', SIZE(uwind)
PRINT, 'NEW SIZE ', SIZE(uwind_used)

DEVICE, DECOMPOSED=1
WINDOW, 2, XSIZE=600, YSIZE=600
!P.background='ffffff'X ; white background, default black print
!P.color = 0
ERASE, 'ffffff'X

MAP_SET, 0, 0, 0, /CYLINDRICAL,  $
		/ISOTROPIC, LIMIT=map_limit, POSITION=[.1, .6, .4, .9], $
		/NOERASE, /NOBORDER, TITLE='June Zonal Wind Plot 700 mb', CHARSIZE=1.
CONTOUR, uwind_used, lon_used, lat_used, $
	/FILL, LEVELS=cu_level, $
	C_COLORS=cu_color, /OVERPLOT, COLOR=0 
CONTOUR, uwind_used, lon_used, lat_used,  $
	LEVELS=cu_line, $
	 /OVERPLOT, COLOR=0, C_LINESTYLE=cu_style, C_LABELS=cu_style
MAP_CONTINENTS, /COASTS
AXIS, XAXIS=0, XRANGE=[LonMin, LonMax],XSTYLE=1, $
	XTITLE='Longitude in Degrees',CHARSIZE=1.2, COLOR=0
AXIS, YAXIS=0, YRANGE=[LatMin, LatMax],YSTYLE=1, $
	YTITLE='Latitude in Degrees',CHARSIZE=1.2, COLOR=0

OPLOT, [LonMin, LonMax], [.999,.999]*LatMax, /NOCLIP
OPLOT, [1.,1.]*LonMax, [LatMin, LatMax], /NOCLIP

;****************Legend for Zonal********************

; Legend ( sideways from normal)
c_level=CU_level[1:21]
c_color=CU_color[1:21]

x0=.425
y0=.9
dx=.1
dy=.017
xb = [0, .1, .1, 0, 0]*dx
yb = [0, 0, 1.2, 1.2, 0]*dy

FOR icol = 0, N_ELEMENTS(c_level)-1 DO BEGIN
	POLYFILL, x0+xb, y0+yb + dy*(-1.*icol), COLOR=c_color[icol], /NORMAL
	PLOTS, x0+xb, y0+yb + dy*(-1.*icol), COLOR=0, /NORMAL
	IF icol NE -1 THEN $
	XYOUTS, x0+dx*.3, y0+dy*icol*(-1) +dy*.9, STRING(c_level[icol], FORMAT='(i3)'),$
		COLOR=0, /NORMAL, CHARSIZE=1
ENDFOR



; ***************Starting the V component.***********

; Read data.  File names
vfn = '/net/b/klotter/vwnd.mon.ltm.nc'

;same as from the U part above
;LatMin = -40. ; negative for South
;LatMax = 40.
;LonMin = -20. ; negative for West
;LonMax = 60.

;Pres = 850.  ; milibar level used
;Month = 9 ; October

map_limit = [LatMin, LonMin, LatMax, LonMax]; for maps

xx_level = INDGEN(21)*1 -10.
CU_level = xx_level*0.5
CU_level = [-9999., CU_level, 9999]
CU_color =  ['ffffff'x, '000088'x, '0000aa'x, '0000cc'x, '0000ff'x, $
  '0055ff'x, '0088ff'X, '00aaff'x, '00ccff'x, '00ffff'x, '11ffdd'x, '44ffaa'x, $
  '66ee44'x, '99ff88'x, 'ddff22'x, 'ffff00'x, 'ffcc00'x, 'ffaa00'X, $
  'ff8800'x, 'ff5500'x, 'ff2200'X, '555555'x, '111111'x ]
CU_color = REVERSE(CU_color)
CU_thick = FLTARR(N_ELEMENTS(CU_level)) +.001
y = WHERE(xx_level MOD 5 EQ 0, amt) + 1
CU_line = CU_level[y]
CU_style = INTARR(amt)+1
i = WHERE(CU_line EQ 0)
CU_style[i] = 2

; vwind data read
lun_vwnd = NCDF_OPEN(vfn)
; note that it is assumed that the files start at Jan 1948

; Read data.  Assume that the lon lat time and level are same for both U/V
NCDF_VARGET, lun_vwnd, "lon", lon
NCDF_VARGET, lun_vwnd, "lat", lat
NCDF_VARGET, lun_vwnd, "time", time
NCDF_VARGET, lun_vwnd, "level", level

NCDF_VARGET, lun_vwnd, "vwnd", vwind
; assuming that SF=1, add_offset=0
NCDF_CLOSE, lun_vwnd

; selection parameters used for plot.  limit to just those of interest here
; LONGITUDES.  Complicated due to wrap-around for Prime Meridian
lon_ind = [ WHERE(lon GE LonMin+360 AND lon LE LonMax+360), $
	WHERE(lon GE LonMin AND lon LE LonMax) ] ; index of longitudes used
check = WHERE(lon_ind GE 0, amt)
lon_ind = lon_ind[check]
lon_used = lon[lon_ind]
lon_ind2 = WHERE(lon_used GT 180, amt) ; wrap around longitudes to plot over PM
IF amt GT 0 THEN lon_used[lon_ind2] -= 360.  ; value of lon's used
n_lon = N_ELEMENTS(lon_used)

lat_ind = WHERE(lat GE LatMin AND lat LE LatMax, n_lat)
lat_used = lat[lat_ind]

prs_ind = WHERE(Pres EQ level, n_prs)

time_ind = Month

vwind_used = FLTARR(n_lon, n_lat, n_prs, 1)
vwind_used[*,*,*,*] = vwind[lon_ind[*], lat_ind[*], prs_ind[*], time_ind]
vwind_used = REFORM(vwind_used)

; test print
PRINT, 'ORIGINAL SIZE ', SIZE(vwind)
PRINT, 'NEW SIZE ', SIZE(vwind_used)

;DEVICE, DECOMPOSED=1
;WINDOW, 3, XSIZE=700, YSIZE=600
!P.background='ffffff'X ; white background, default black print
!P.color = 0
;ERASE, 'ffffff'X

MAP_SET, 0, 0, 0, /CYLINDRICAL,  $
		/ISOTROPIC, LIMIT=map_limit, POSITION=[.6, .6, .9, .9], $
		/NOERASE, /NOBORDER, TITLE='June Meridional Wind Plot 700 mb',CHARSIZE=1
CONTOUR, vwind_used, lon_used, lat_used, $
	/FILL, LEVELS=cu_level, $
	C_COLORS=cu_color, /OVERPLOT, COLOR=0
CONTOUR, vwind_used, lon_used, lat_used,  $
	LEVELS=cu_line, $
	 /OVERPLOT, COLOR=0, C_LINESTYLE=cu_style, C_LABELS=cu_style


MAP_CONTINENTS, /COASTS
;AXIS, XAXIS=0, XRANGE=[LonMin, LonMax], CHARSIZE=1.7, COLOR=0
;AXIS, YAXIS=0, YRANGE=[LatMin, LatMax], CHARSIZE=1.7, COLOR=0
AXIS, XAXIS=0, XRANGE=[LonMin, LonMax],XSTYLE=1, $
	XTITLE='Longitude in Degrees',CHARSIZE=1.2, COLOR=0
AXIS, YAXIS=0, YRANGE=[LatMin, LatMax],YSTYLE=1, $
	YTITLE='Latitude in Degrees',CHARSIZE=1.2, COLOR=0
OPLOT, [LonMin, LonMax], [.999,.999]*LatMax, /NOCLIP
OPLOT, [1.,1.]*LonMax, [LatMin, LatMax], /NOCLIP

;****************Meridional Legend********************

; Legend ( sideways from normal)
c_level=CU_level[1:21]
c_color=CU_color[1:21]

x0=.93
y0=.9
dx=.1
dy=.017
xb = [0, .1, .1, 0, 0]*dx
yb = [0, 0, 1.2, 1.2, 0]*dy

FOR icol = 0, N_ELEMENTS(c_level)-1 DO BEGIN
	POLYFILL, x0+xb, y0+yb + dy*(-1.*icol), COLOR=c_color[icol], /NORMAL
	PLOTS, x0+xb, y0+yb + dy*(-1.*icol), COLOR=0, /NORMAL
	IF icol NE -1 THEN $
	XYOUTS, x0+dx*.13, y0+dy*icol*(-1) +dy*.9, STRING(c_level[icol], FORMAT='(i3)'),$
		COLOR=0, /NORMAL, CHARSIZE=1
ENDFOR


;***************** VeloVect****************
u_grid=GRIDDATA(lon_used, lat_used, uwind_used)
v_grid=GRIDDATA(lon_used, lat_used, vwind_used)


VELOVECT, u_grid, v_grid, lon_used, lat_used, COLOR=0, /CYLINDRICAL, /ISOTROPIC,$
	DOTS=0, TITLE='Wind Vector with Magnitudes', XTITLE='Longitude in Degrees',$
	YTITLE='Latitude in Degrees', CHARSIZE=1,LENGTH=1,$
	POSITION=[.1,.1,.4,.4], /NOERASE,/NOBORDER

MAP_SET, 0, 0, 0, /CYLINDRICAL,  $
		/ISOTROPIC, LIMIT=map_limit, POSITION=[.1, .1, .4, .4], $
		/NOERASE;, /NOBORDER,CHARSIZE=1.5
MAP_CONTINENTS, /COASTS



END