FUNCTION get_lag,doy,yr,hours,lagheader=lagheader
;opens and reads the lags for a day
;read lag file unless nolag is set
;set - if can not find lag - use median of all day
  IF NOT keyword_set(yr) THEN year=co2_const('cur_year') ELSE year = yr
  year2 = year MOD 100
  yrstr = string(year2,format='(i2.2)')
  doystr = string(doy,format='(i3.3)')
  IF NOT keyword_set(hours) THEN hours = indgen(48) /2.0
  ntimes = n_elements(hours)
  fhead = co2_const('file_head')
  lagfile = co2_basedir('lag',yr=year)+fhead+yrstr+'.lag.xdf'
  lags = make_array(3,ntimes)
  lags[0,*] = hours
  lags[1:2,*] = -1.0 ;default c and q lags
  IF file_test(lagfile,/read) THEN BEGIN
    print,'  Reading LAGS from ',lagfile
    lagfp = assoc_xdf(lagfile,header=lagheader) ;doy,hhmm,ww_lag,wt_lag,wc_lag,wq_lag
    lagdata = line_xdf(lagfp,(doy-1)*48,numlines=48)
    close_xdf,lagfp
    goodlag = where(isnotnan(lagdata[4,*]) OR isnotnan(lagdata[5,*]),ngood,complement=badlag)
    IF ngood GT 0 THEN BEGIN    ;some or all good data - can use
      IF badlag[0] NE -1 THEN BEGIN ;some bad data
        medc = median(lagdata[4,*])
        medq = median(lagdata[5,*])
        lagdata[4,badlag] = medc
        lagdata[5,badlag] = medq
      ENDIF
      timelocs = fix(hours * 2)
      lags[1,*] = lagdata[4,timelocs]
      lags[2,*] = lagdata[5,timelocs]
    ENDIF ELSE print,'    No good LAG data for the day, using -1'
  ENDIF ELSE BEGIN
    print,'  Creating LAG file ',lagfile
    print,'  Using -1 for current lags'
    lagarray = make_array(6,17568,/float,value=nan())
    lagarray[0,*] = float(fix(findgen(17568) / 48) + 1)
    lagarray[1,*] = float(fjday_to_hhmm((findgen(17568) MOD 48) / 48.0))
    lagheader = ['doy','hhmm','ww_lag','wt_lag','wc_lag','wq_lag']
    create_dir,co2_basedir('lag',yr=year)
    write_xdf,lagfile,lagarray,header=lagheader
    lagarray = 0b
  ENDELSE
  return,lags
END

PRO calc_lag,lagwindow,dur,dvr,dwr,dt,dc,dq,lag=lag,cov=cov,time=time,nolag=nolag,noplot=noplot,plottitle=plottitle,hz=hz
  IF NOT keyword_set(nolag) THEN begin
    maxWin =lagwindow(0,0)
    minWin =lagwindow(0,1)
    cmaxWin=lagwindow(1,0) & cminWin =lagwindow(1,1)
    qmaxWin=lagwindow(2,0) & qminWin =lagwindow(2,1)
    IF n_elements(hz) EQ 0 THEN hz = 10.0
    lagcovar_window,dwr,dwr,hz,maxWin,minWin,lcov=wwcov,tlag=wwtime,tpeak=wwlag
    lagcovar_window,dwr,dt,hz,maxWin,minWin,lcov=wtcov,tlag=wttime,tpeak=wtlag
    lagcovar_window,dwr,dc,hz,cmaxWin,cminWin,lcov=wccov,tlag=wctime,tpeak=wclag
    lagcovar_window,dwr,dq,hz,qmaxWin,qminWin,lcov=wqcov,tlag=wqtime,tpeak=wqlag
    ;cov=[transpose(wwcov),transpose(wtcov),transpose(wccov),transpose(wqcov)]
    ;time=[transpose(wwtime),transpose(wttime),transpose(wctime),transpose(wqtime)]
    lag=[wwlag,wtlag,wclag,wqlag]
    print,'    Lags (ww,wt,wc,wq) ',lag

    IF NOT keyword_set(noplot) THEN BEGIN
      oldmult = !p.multi
      !p.multi = [0,2,2]
      qq = where(isnan(wwcov),numww)
      if numww ne n_elements(wwcov) then plot,wwtime,wwcov,chars=1.3,thick=2,charthick=2,title='ww Lag',xtitle='Time (s)',ytitle='ww covariance',xrange=[-50,50]
      qq = where(isnan(wtcov),numwt)
      if numwt ne n_elements(wtcov) then plot,wttime,wtcov,chars=1.3,thick=2,charthick=2,title='wt Lag',xtitle='Time (s)',ytitle='wt covariance',xrange=[-50,50]
      qq = where(isnan(wccov),numwc)
      if numwc ne n_elements(wccov) then plot,wctime,wccov,chars=1.3,thick=2,charthick=2,title='wc Lag',xtitle='Time (s)',ytitle='wc covariance',xrange=[-50,50]
      qq = where(isnan(wqcov),numwq)
      if numwq ne n_elements(wqcov) then plot,wqtime,wqcov,chars=1.3,thick=2,charthick=2,title='wq Lag',xtitle='Time (s)',ytitle='wq covariance',xrange=[-50,50]
      IF keyword_set(plottitle) THEN BEGIN
        xyouts,0.5,1.0,plottitle+' LAGS',/normal,alignment=0.5,chars=1.5
      ENDIF
      !P.multi = oldmult
    ENDIF
  ENDIF ELSE BEGIN
    print,'    Skipping lag calculation (nolag is set)'
    cov = nan() & time = nan() & lag = [nan(),nan(),nan(),nan()]
  ENDELSE
END

PRO save_lag,lagarr,doy,yr
  IF NOT keyword_set(yr) THEN year=co2_const('cur_year') ELSE year = yr
  year2 = year MOD 100
  fhead = co2_const('file_head')
  yrstr = string(Year2,format='(i2.2)')
  lagfile = co2_basedir('lag',yr=year)+fhead+yrstr+'.lag.xdf'
  IF file_test(lagfile,/read) THEN BEGIN
    lagfl = read_xdf(lagfile,header=lagheader)
  ENDIF ELSE BEGIN
    print,'  Creating LAG file ',lagfile
    lagfl = make_array(6,17568,/float,value=nan())
    lagheader = ['doy','hhmm','ww_lag','wt_lag','wc_lag','wq_lag']
    lagfl[0,*] =  float(fix(findgen(17568) / 48) + 1)
    lagfl[1,*] =  float(fjday_to_hhmm((findgen(17568) MOD 48) / 48.0))
  ENDELSE
  locs = ((doy-1)*48) + hhmm_to_fjday(lagarr[0,*]) * 48
  lagfl[2:5,locs] = lagarr[1:4,*]
  print,'  Writing LAG file ',lagfile
  create_dir,co2_basedir('lag',yr=year)
  write_xdf,lagfile,lagfl,header=lagheader
END

PRO fill_lag,yr,nosave=nosave
;fill missing data with 20 day median
  IF NOT keyword_set(yr) THEN year=co2_const('cur_year') ELSE year = yr
  year2 = year MOD 100
  fhead = co2_const('file_head')
  yrstr = string(Year2,format='(i2.2)')
  lagfile = co2_basedir('lag',yr=year)+fhead+yrstr+'.lag.xdf'
  IF file_test(lagfile,/read) THEN BEGIN
    print,'Reading lagfile ',lagfile
    lagdata = read_xdf(lagfile,header=lagheader)
  ENDIF ELSE BEGIN
    print,'Cannot find lagfile ',lagfile
    GOTO,leave
  ENDELSE

;fill c holes
  badlag = where(isnan(lagdata[4,*]),nbadlag)
  IF nbadlag GT 0 THEN BEGIN
    print,'  Filling C lag holes n: ',nbadlag
    startloc = badlag - 960
    endloc = badlag + 960
    badstart = where(startloc LT 0,nbadstart)
    IF nbadstart GT 0 THEN BEGIN
      endloc[badstart] = endloc[badstart] - startloc[badstart]
      startloc[badstart] = 0
    ENDIF
    endpt = days_in_year(year)*48
    badend = where(endloc GE endpt,nbadend)
    IF nbadend GT 0 THEN BEGIN
      startloc[badend] = startloc[badend] - (endloc[badend] - endpt - 1)
      endloc[badend] = endpt - 1
    ENDIF
    medianval = make_array(nbadlag,/float,value=nan())
    FOR i = 0,nbadlag-1 DO medianval[i] = median(lagdata[4,startloc[i]:endloc[i]])
    badmedian = where(isnan(medianval),nbadmedian)
    IF nbadmedian GT 0 THEN BEGIN
      doys = lagdata[0,badlag[badmedian]]
      doys = doys[uniq(doys,sort(doys))]
      print,'  Cannot fill holes for all days (DOYS,nvals) ',doys,nbadmedian
    ENDIF
    lagdata[4,badlag] = medianval
    choles = 1b
  ENDIF ELSE BEGIN
    print,'  No C holes to fill'
    choles = 0b
  ENDELSE


;fill q holes
  badlag = where(isnan(lagdata[5,*]),nbadlag)
  IF nbadlag GT 0 THEN BEGIN
    print,'  Filling Q lag holes n: ',nbadlag
    startloc = badlag - 960
    endloc = badlag + 960
    badstart = where(startloc LT 0,nbadstart)
    IF nbadstart GT 0 THEN BEGIN
      endloc[badstart] = endloc[badstart] - startloc[badstart]
      startloc[badstart] = 0
    ENDIF
    endpt = days_in_year(year)*48
    badend = where(endloc GE endpt,nbadend)
    IF nbadend GT 0 THEN BEGIN
      startloc[badend] = startloc[badend] - (endloc[badend] - endpt - 1)
      endloc[badend] = endpt - 1
    ENDIF
    medianval = make_array(nbadlag,/float,value=nan())
    FOR i = 0,nbadlag-1 DO medianval[i] = median(lagdata[4,startloc[i]:endloc[i]])
    badmedian = where(isnan(medianval),nbadmedian)
    IF nbadmedian GT 0 THEN BEGIN
      doys = lagdata[0,badlag[badmedian]]
      doys = doys[uniq(doys,sort(doys))]
      print,'  Cannot fill holes for all days (DOYS,nvals) ',doys,nbadmedian
    ENDIF
    lagdata[5,badlag] = medianval
    qholes = 1b
  ENDIF ELSE BEGIN
    print,'  No Q holes to fill'
    qholes = 0b
  ENDELSE

  IF ((choles EQ 1) OR (qholes EQ 1)) AND (NOT keyword_set(nosave)) THEN BEGIN
    print,'  Writing cleaned lag file'
    write_xdf,lagfile,lagdata,header=lagheader
  ENDIF

  IF NOT keyword_set(nosave) THEN update_version,'fill_lag',1,year

leave:
  print,'Lag filling complete.'
END