;
; Reads in HADCM2 control run from 40-yr segment crunc files, extracts
; various windows and then computes different seasons from each window
; average.
;
nmon=12
yrstart=2260              ; 1860 for full record
yrend=3259                ; 3259 for full record
nyr=yrend-yrstart+1
seglen=40
nseg=nyr/seglen
fnstem='/cru/u2/f055/data/mod/hadcm2/con1000/nc/HCCI_temk_'
fntail='.mon.nc'
;
; Get HadCM2 land-sea mask
;
lsmask=linkmask('hadcm2')
ml=where(lsmask eq 0.,nsea)
if nsea gt 0 then lsmask(ml)=!values.f_nan
;
; Define arrays
;
nts=9
timey=fltarr(nyr)            ; Fill years in as we go
hccits=fltarr(nyr,nts)*!values.f_nan
moncet=fltarr(nmon,nyr)
;
; Do each segment in turn
;
for iseg = 0 , nseg-1 do begin
  ;
  ; Create input file name and open the file
  ;
  i1=iseg*seglen
  y1=yrstart+i1
  fn=fnstem+string([y1,y1+seglen-1],format='(I4,"-",I4)')+fntail
  print,fn
  ncid=ncdf_open(fn)
  ;
  ; Now read in the entire monthly fields
  ;
  fd40=crunc_rddata(ncid,grid=g,info=i)  
  ncdf_close,ncid 
  if seglen ne g.nt/nmon then message,'File contains wrong number of years'
  ;
  segyrs=reform(g.year,nmon,seglen)
  timey(i1:i1+seglen-1)=segyrs(0,*)
  ;
  ; The first time round locate the various windows required
  ;
  if iseg eq 0 then begin
    ; Locate Northern Hemisphere
    nhjkeep=where(g.y ge 0.)
    nhy=g.y(nhjkeep)
    ; Locate Jasper box
    jasikeep=where(g.x eq 243.75)
    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 63.75)
    uraikeep2=where(g.x eq 67.5)
    urajkeep=where(g.y eq 65.)
    ; Locate CET boxes (two)
    cetikeep1=where(g.x eq 0.)
    cetikeep2=where(g.x eq 356.25)
    cetjkeep=where(g.y eq 52.5)
  endif
  ;
  ; 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 segment of the series.
  ;
  fd1=fd40(*,nhjkeep,*)           ; extract northern hemisphere
  fd2=globalmean(fd1,nhy)         ; average NH, weighted by latitude
  fd2=reform(fd2,nmon,seglen)     ; separate months and years
  hccits(i1:i1+seglen-1,0)=mkseason(fd2,5,7)      ; JJA
  hccits(i1:i1+seglen-1,1)=mkseason(fd2,0,11)     ; Jan-Dec
  ;
  fd0=fd40                        ; set all sea points to missing
  for i = 0 , nmon*seglen-1 do fd0(*,*,i)=fd0(*,*,i)*lsmask(*,*)
  fd1=fd0(*,nhjkeep,*)            ; extract northern hemisphere
  fd2=globalmean(fd1,nhy)         ; average NH, weighted by latitude
  fd2=reform(fd2,nmon,seglen)     ; separate months and years
  hccits(i1:i1+seglen-1,2)=mkseason(fd2,3,8)      ; AMJJAS
  ;
  fd3=fd40(jasikeep,jasjkeep,*)
  fd3=reform(fd3,nmon,seglen)
  hccits(i1:i1+seglen-1,3)=mkseason(fd3,3,7)      ; AMJJA
  ;
  fd4=0.5*(fd40(uraikeep1,urajkeep,*)+fd40(uraikeep2,urajkeep,*))
  fd4=reform(fd4,nmon,seglen)
  hccits(i1:i1+seglen-1,4)=mkseason(fd4,4,8)      ; MJJAS
  ;
  fd5=fd40(fenikeep,fenjkeep,*)
  fd5=reform(fd5,nmon,seglen)
  hccits(i1:i1+seglen-1,5)=mkseason(fd5,3,7)      ; AMJJA
  ;
  fd6=0.5*(fd40(cetikeep1,cetjkeep,*)+fd40(cetikeep2,cetjkeep,*))
  fd6=reform(fd6,nmon,seglen)
  hccits(i1:i1+seglen-1,6)=mkseason(fd6,0,11)      ; Jan-Dec
  hccits(i1:i1+seglen-1,7)=mkseason(fd6,3,8)       ; AMJJAS
  ;
  ; Store monthly timeseries for doing the Oct-Mar average on all 1400yr at end
  ;
  moncet(*,i1:i1+seglen-1)=fd6(*,*)
  ;
  ; repeat for next segment
  ;
endfor
;
; Do the winter CET now
;
hccits(*,8)=mkseason(moncet,9,2)
;
; Now save all the series for trend and means analysis
;
save,filename='hcci_ts.idlsave',$
  hccits,timey,nts,nyr
;
end