nc = ncdf_open('/Users/rachel/idl_course_week4/zg_ukmo.nc')
	
	; Extract the variables of interest
	ncdf_varget, nc, 'lat', lat, count = 14, stride = 2
	ncdf_varget, nc, 'lon', lon, stride = 2
	ncdf_varget, nc, 'plev', level, count = 1, offset = 5
	ncdf_varget, nc, 'zg', height, count = [48, 14, 1, 600], offset = [0,0,5,0],  stride = [2,2,1,1]
	ncdf_close, nc
	

	height = reform(height, 48, 14, 600)
	
	dim = size(height)
	isz = dim(1)
	jsz = dim(2)
	ksz = dim(3)
	
	time_grid = reform(height, isz, jsz, 12, ksz/12)
	climatology = total(time_grid,4)/(total((finite(time_grid)),4))
	climate_grid = rebin(climatology, isz, jsz, 12, ksz/12)
	height_anomalies = time_grid - climate_grid
	height_anomalies = reform(height_anomalies, isz, jsz, ksz)
	
	djf_height = fltarr(isz, jsz, ksz/4)
	FOR i = 0, ksz/12-2 DO BEGIN		
		djf_height(*,*,(i*3)+2:(i*3)+4) = height_anomalies(*,*,(i*12)+11:(i*12)+13)
	ENDFOR

	weight = rebin(reform(sqrt(cos(lat*!DTOR)), 1, jsz), isz, jsz, ksz)

	DJF_height_weight = djf_height*weight
	
	DJF_height_weight = reform(DJF_height_weight, isz*jsz, ksz/4)
	

	SVDC, DJF_height_weight, W, U, V, /DOUBLE
	
	PC1 = REFORM(U(0,*)/STDDEV(U(0,*)))
	PC2 = REFORM(U(1,*)/STDDEV(U(1,*)))
	PC3 = REFORM(U(2,*)/STDDEV(U(2,*)))
		
	djf_height = REFORM(djf_height, isz*jsz, ksz/4)
	
	reg_grid_1 = DBLARR(isz*jsz)
	reg_grid_2 = DBLARR(isz*jsz)
	reg_grid_3 = DBLARR(isz*jsz)
		
	FOR i = 0, isz*jsz-1 DO BEGIN
		reg_grid_1(i) = REGRESS(PC1,REFORM(DJF_height(i,*)))
		reg_grid_2(i) = REGRESS(PC2,REFORM(DJF_height(i,*)))
		reg_grid_3(i) = REGRESS(PC3,REFORM(DJF_height(i,*)))
	ENDFOR
	
	grid_1 = reform(reg_grid_1, isz, jsz)
	grid_2 = reform(reg_grid_2, isz, jsz)
	grid_3 = reform(reg_grid_3, isz, jsz)
	
	lon = [lon, lon(0)]
	grid_1 = [grid_1, grid_1(0,*,*)]
	
	
	
	device,true=24,decompose=0,retain=2
	colors_anomalies
	
	!p.background = 255
	!p.color = 0
	cmax = 60
	cmin = -60
	step = (cmax - cmin) /247.0
	clevels = indgen(247)*step+cmin

	!p.font = 1
	map_set,/STEREO,-90,-60,/ISOTROPIC,/HORIZON,color=0,position=[.2,.2,.8,.8], $
		title = 'PC1', charsize = 1.5, limit = [-90, 0, -20, 360] 
	contour, grid_1, lon, lat, /overplot, /noerase, $
		c_colors = indgen(247)+1, levels = clevels, /cell_fill
	map_continents,/continents,mlinethick=2,/overplot
	cbar, bottom = 1, topp = 247, charsize = 1.5, divisions = 12, $
		position = [0.3, 0.1, 0.7, 0.14], range = [min(grid_1), max(grid_1)], $
		format = '(F8.1)'
	write_jpeg, 'map1.jpg', tvrd(true=3), quality=100, true=3	

	

	
END