source: trunk/src/paper01/fig3/read_variables_v2.pro @ 70

;+
;
; .. _read_variables_v2.pro:
;
; =====================
; read_variables_v2.pro
; =====================
;
; DESCRIPTION
; ===========
;
; .. graphviz::
;
;    digraph read_variables_v2 {
;       graph [
;       rankdir="LR",
;       ]
;       rad [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/rad*_dy.cdf"];
;       qlat[shape=ellipse,fontname=Courier,label="${PROJECT_ID}/qlat*_dy.cdf"];
;       pos[shape=ellipse,fontname=Courier,label="${PROJECT_ID}/pos*_dy.cdf"];
;       met[shape=ellipse,fontname=Courier,label="${PROJECT_ID}/met*_dy.cdf"];
;
;       at [shape=tab,fontname=Courier,label="at"]
;       sw [shape=tab,fontname=Courier,label="sw"]
;       rh [shape=tab,fontname=Courier,label="rh"]
;       sst [shape=tab,fontname=Courier,label="sst"]
;       wu [shape=tab,fontname=Courier,label="wu"]
;       wv [shape=tab,fontname=Courier,label="wv"]
;       ws [shape=tab,fontname=Courier,label="ws"]
;       lh [shape=tab,fontname=Courier,label="lh"]
;
;       read_variables_v2 [shape=box,
;       fontname=Courier,
;       color=blue,
;       URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/paper01/fig2/read_variables_v2.pro",
;       label="${TROPFLUX}/src/paper01/fig2/read_variables_v2.pro"];
;
;       {qlat rad pos met} -> {read_variables_v2} -> {at sw rh sst wu wv ws lh}
;    }
;
; SEE ALSO
; ========
;
; :ref:`project_profile.sh`
;
; :ref:`air_validation_scatter_2000_2009_v50.pro`
; :ref:`q2m_validation_scatter_2000_2009_v50.pro`
; :ref:`sst_validation_scatter_2000_2009_v50.pro`
; :ref:`ws_validation_scatter_2000_2009_v50.pro`
; :ref:`air_validation_scatter_2000_2009_basin.pro`
; :ref:`q2m_validation_scatter_2000_2009_basin.pro`
; :ref:`sst_validation_scatter_2000_2009_basin.pro`
; :ref:`ws_validation_scatter_2000_2009_basin.pro`
; :ref:`lhf_validation_scatter_2000_2009.pro`
; :ref:`shf_validation_scatter_2000_2009.pro`
;
; EXAMPLES
; ========
;
; ::
;
;  IDL> csite=++
;  IDL> ++
;  IDL> read_variables_v2, csite, date1, date2, nsmooth, at, sw, rh, sst, wu, wv, ws, lh
;
; TODO
; ====
;
; externalize function
;
; make it work (missing data)
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; - fplod 20110412T110534Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * remove hard coding path
;   * add graphviz
;
; - fplod 20110411T142955Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
;-
pro read_variables_v2, csite, date1, date2, nsmooth, $
                     at, sw, rh, sst, wu, wv, ws, 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
;
;;SHORTWAVE
fi=project_id_env+'rad'+csite+'_dy.cdf'
f=file_test(fi)
sw=fltarr(jpt)+!values.f_nan
sw_q=fltarr(jpt)+!values.f_nan
;tt=fltarr(jpt)+!values.f_nan
if (f) then begin
  tt0=time_lec(fi)
  sw0=ncdf_lec(fi,var='RD_495')
  sw_q0=ncdf_lec(fi,var='QSW_5495')
  tt0=tt0-time(0)
  ind=where((tt0 ge -1e-5) and (tt0 le jpt-1+1e-5))
  if (ind(0) ne -1) then begin
     sw(tt0(ind))=sw0(ind)
     sw_q(tt0(ind))=sw_q0(ind)
  endif
endif else begin
  print, fi, ' not found'
endelse

ind_sw=where(sw_q ne 1. and sw_q ne 2.)

;;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

;;POSITION
fi=project_id_env+'pos'+csite+'_dy.cdf'
f=file_test(fi)
lat=fltarr(jpt)+!values.f_nan
lon=fltarr(jpt)+!values.f_nan
if (f) then begin
  tt0=time_lec(fi)
  lat0=ncdf_lec(fi,var='LAT_500')
  lon0=ncdf_lec(fi,var='LON_502')
  tt0=tt0-time(0)
  ind=where((tt0 ge -1e-5) and (tt0 le jpt-1+1e-5))
  if (ind(0) ne -1) then begin
    lat(tt0(ind))=lat0(ind)
    lon(tt0(ind))=lon0(ind)
  endif
endif else begin
  print, fi, ' not found'
endelse


;;MET PARAMETERS
fi=project_id_env+'met'+csite+'_dy.cdf'
f=file_test(fi)
wu=fltarr(jpt)+!values.f_nan
wv=fltarr(jpt)+!values.f_nan
ws=fltarr(jpt)+!values.f_nan
ws_q=fltarr(jpt)+!values.f_nan
at=fltarr(jpt)+!values.f_nan
at_q=fltarr(jpt)+!values.f_nan
rh=fltarr(jpt)+!values.f_nan
rh_q=fltarr(jpt)+!values.f_nan
sst=fltarr(jpt)+!values.f_nan
sst_q=fltarr(jpt)+!values.f_nan

if (f) then begin
  tt0=time_lec(fi)
  wu0=ncdf_lec(fi,var='WU_422')
  wv0=ncdf_lec(fi,var='WV_423')
  ws0=ncdf_lec(fi,var='WS_401')
  ws_q0=ncdf_lec(fi,var='QWS_5401')
  at0=ncdf_lec(fi,var='AT_21')
  at_q0=ncdf_lec(fi,var='QAT_5021')
  rh0=ncdf_lec(fi,var='RH_910')
  rh_q0=ncdf_lec(fi,var='QRH_5910')
  sst0=ncdf_lec(fi,var='T_25')
  sst_q0=ncdf_lec(fi,var='QT_5025')
  tt0=tt0-time(0)
  ind=where((tt0 ge -1e-5) and (tt0 le jpt-1+1e-5))
  if (ind(0) ne -1) then begin
    wu(tt0(ind))=wu0(ind)
    wv(tt0(ind))=wv0(ind)
    ws(tt0(ind))=ws0(ind)
    ws_q(tt0(ind))=ws_q0(ind)
    at(tt0(ind))=at0(ind)
    at_q(tt0(ind))=at_q0(ind)
    rh(tt0(ind))=rh0(ind)
    rh_q(tt0(ind))=rh_q0(ind)
    sst(tt0(ind))=sst0(ind)
    sst_q(tt0(ind))=sst_q0(ind)
  endif
endif else begin
  print, fi, ' not found'
endelse

ind_at=where(at_q ne 1 and at_q ne 2)   &  ind_ws=where(ws_q ne 1 and ws_q ne 2)
ind_rh=where(rh_q ne 1 and rh_q ne 2)   &  ind_sst=where(sst_q ne 1 and sst_q ne 2)
ind_lh=where(lh_q ne 1 and lh_q ne 2)

ind=union(ind_at, union(ind_rh, union(ind_ws, union(ind_lh, ind_sst))))

if (ind(0) ne -1) then begin
  sw(ind)=!Values.f_nan   &   at(ind)=!Values.f_nan
  rh(ind)=!Values.f_nan   &   ws(ind)=!Values.f_nan
  sst(ind)=!Values.f_nan  &   lh(ind)=!Values.f_nan
  wu(ind)=!Values.f_nan  &    wv(ind)=!Values.f_nan
endif
;
;; Replace missing values by "NaN"
;
tsvars=['at','sw','rh','sst','wu','wv','ws','lat','lon']
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
at=smooth(at,nsmooth,/nan) & sst=smooth(sst,nsmooth,/nan)
lh=smooth(lh,nsmooth,/nan) & rh=smooth(rh,nsmooth,/nan)
wu=smooth(wu,nsmooth,/nan) & wv=smooth(wv,nsmooth,/nan)
ws=smooth(ws,nsmooth,/nan) & sw=smooth(sw,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