;-------------------------------------------------------------------------------------------------
pro read_lh, csite,date1,date2,nsmooth, $
                      lh

;
;-------------------------------------------------------------------------------------------------
@common

dir='/Users/pkb/data/TPR/'
;
;; DEFINE THE OUTPUT TIME AXIS
;
jda1=date2jul(date1)
jda2=date2jul(date2)
jpt=(jda2-jda1+1l)
time=jda1+dindgen(jpt)

;
;; FIRST READ ALL DATA FROM SITE
;;   - if file exists, extract correct time axis
;;   - if not, fill variable with missing values
;


;;LHF
fi=dir+'qlat'+csite+'_dy.cdf'
f=file_test(fi)
lh=fltarr(jpt)+!values.f_nan
lh_q=fltarr(jpt)+!values.f_nan
if (f) then begin
  tt0=time_lec(fi)
  lh0=ncdf_lec(fi,var='QL_137')
  lh_q0=ncdf_lec(fi,var='QQL_5137')
  tt0=tt0-time(0)
  ind=where((tt0 ge -1e-5) and (tt0 le jpt-1+1e-5))
  if (ind(0) ne -1) then begin
  lh(tt0(ind))=lh0(ind)
  lh_q(tt0(ind))=lh_q0(ind)
  endif
endif

ind=where(lh_q ne 1 and lh_q ne 2)


if (ind(0) ne -1) then begin
  lh(ind)=!Values.f_nan
endif
;
;; Replace missing values by "NaN"
;
tsvars=['lh'] 
vars=[tsvars]
nn=n_elements(vars)
for n=0,nn-1 do begin
  var=vars(n)
  com='ind=where('+var+' ge 1.e20) & if (ind(0) ne -1) then '+var+'(ind)=!values.f_nan'
  r=execute(com)
endfor

nsmooth=nsmooth
lh=smooth(lh,nsmooth,/nan) 

;-------------------------------------------------------------------------------------------------
end
;-------------------------------------------------------------------------------------------------


;-------------------------------------------------------------------------------------------------
function time_lec, fi
tt=ncdf_lec(fi,var='time')
fid=ncdf_open(fi) & vid=ncdf_varid(fid,'time')
ncdf_attget, fid,vid,'units',orig
orig=string(orig)
yy=long(strmid(orig,11,4))
mm=long(strmid(orig,16,2))
dd=long(strmid(orig,19,2))
tt=julday(mm,dd,yy,12,00)+tt
return, tt
end

;-------------------------------------------------------------------------------------------------