;+
; .. _read_lh.pro:
;
; ===========
; read_lh.pro
; ===========
;
; DESCRIPTION
; ===========
;
; .. graphviz::
;
;    digraph read_lh {
;
;       qlat [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/qlat*_dy.cdf"];
;
;       lh [shape=tab,fontname=Courier,label="lh"]
;
;       read_lh [shape=box,
;       fontname=Courier,
;       color=blue,
;       URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/paper01/fig2/read_lh.pro",
;       label="${TROPFLUX}/src/paper01/fig2/read_lh.pro"];
;
;       {qlat} -> {read_lh} -> {lh}
;    }
;
; SEE ALSO
; ========
;
; :ref:`project_profile.sh`
; :ref:`project_init.pro`
; :ref:`cm_project.pro`
;
; :func:`time_lec`
;
; :ref:`fig2_timeline_diagram.pro`
;
; EXAMPLES
; ========
;
; ::
;
;  csite = '8s67e'
;  date1 = 19890101L
;  date2 = 20091231L
;  nsmooth = 1.
;  read_lh, csite, date1, date2, nsmooth, lh
;
; TODO
; ====
;
; missing input files
;
; coding rules
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20110412T093201Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * externalize function
;   * remove hard coding path
;   * add graphviz
;
; - fplod 20110411T142955Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
;-
pro read_lh, csite, date1, date2, nsmooth, lh
;
@cm_project
;
; 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=project_id_env+'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 else begin
  print, fi, ' not found'
endelse
;
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