;
; Reads in observed gridded temperature and CET, extracts
; various windows and then computes different seasons from each window
; average.
;
nmon=12
yrstart=1900
yrend=1999     ; 2000 wasn't complete at the time
nyr=yrend-yrstart+1
fn1='/cru/u2/f055/data/obs/grid/surface/ist_ipcc_18561999.mon.nc'
fn2='/cru/u2/f055/data/obs/grid/surface/lat_jones_18512000.mon.nc'
fn3='cet.dat'
;
; Define arrays
;
nts=9
timey=findgen(nyr)+yrstart
hccits=fltarr(nyr,nts)*!values.f_nan
;
; First of all, read in the gridded land and marine temperatures
;
ncid=ncdf_open(fn1)
fd40=crunc_rddata(ncid,tst=[1900,0],tend=[1999,11],grid=g,info=i)  
ncdf_close,ncid
;
; Locate the various windows required
;
; Locate Northern Hemisphere
nhjkeep=where(g.y ge 0.)
nhy=g.y(nhjkeep)
; Locate Jasper box
jasikeep=where(g.x eq -117.5)
jasjkeep=where(g.y eq 52.5)
; Locate N. Fennoscandia box
fenikeep=where(g.x eq 22.5)
fenjkeep=where(g.y eq 67.5)
; Locate N. Urals boxes (two)
uraikeep1=where(g.x eq 62.5)
uraikeep2=where(g.x eq 67.5)
urajkeep=where(g.y eq 67.5)
uraikeep1=uraikeep1(0) & uraikeep2=uraikeep2(0) & urajkeep=urajkeep(0)
;
; For each required series, extract the appropriate window, average (with
; area-weighting if it spans latitudes) if necessary (i.e., not a single
; box), then make the appropriate seasonal average for that series.
; Store this series.
;
fd1=fd40(*,nhjkeep,*)           ; extract northern hemisphere
fd2=globalmean(fd1,nhy)         ; average NH, weighted by latitude
fd2=reform(fd2,nmon,nyr)     ; separate months and years
hccits(*,0)=mkseason(fd2,5,7)      ; JJA
hccits(*,1)=mkseason(fd2,0,11)     ; Jan-Dec
;
fd3=fd40(jasikeep,jasjkeep,*)
fd3=reform(fd3,nmon,nyr)
hccits(*,3)=mkseason(fd3,3,7)      ; AMJJA
;
fd4 = reform(fd40(uraikeep1:uraikeep2,urajkeep,*))
fd4 = total(fd4,1,/nan) / float(total(finite(fd4),1))
help,fd4
fd4=reform(fd4,nmon,nyr)
hccits(*,4)=mkseason(fd4,4,8)      ; MJJAS
;
fd5=fd40(fenikeep,fenjkeep,*)
fd5=reform(fd5,nmon,nyr)
hccits(*,5)=mkseason(fd5,3,7)      ; AMJJA
;
; Now do the CET series
;
openr,1,fn3
readf,1,ncet
rawdat=fltarr(18,ncet)
readf,1,rawdat,format='(I4,12F5.1,5F6.2)'
close,1
cetyr=reform(rawdat(0,*))
kl=where(cetyr ge yrstart,nk)
if nk ne nyr then message,'Ooooops!!'
moncet=rawdat(1:12,kl)
ts=mkseason(moncet,0,11)
hccits(*,6)=ts(*)
ts=mkseason(moncet,3,8)
hccits(*,7)=ts(*)
ts=mkseason(moncet,9,2)
hccits(*,8)=ts(*)
;
; Compute the >20N land instrumental temperature timseries and filter
;
print,'Reading temperatures'
ncid=ncdf_open(fn2)
tsmon=crunc_rddata(ncid,tst=[1900,0],tend=[1999,11],grid=gtemp)
ncdf_close,ncid
;
; Compute the northern hemisphere >20N land series
;
; First extract >20N rows
kl=where(gtemp.y gt 20.)
ylat=gtemp.y(kl)
tsnorth=tsmon(*,kl,*)
; Compute latitude-weighted mean
nhmon=globalmean(tsnorth,ylat)
; Compute Apr-Sep mean
nhmon=reform(nhmon,nmon,nyr)
hccits(*,2)=mkseason(nhmon,3,8)
;
; Now save all the series for trend and means analysis
;
save,filename='obs_ts.idlsave',$
  hccits,timey,nts,nyr
;
end