source: trunk/src/read_basic_var.pro @ 152

;+
;
; .. _read_basic_var.pro:
;
; ==================
; read_basic_var.pro
; ==================
;
; DESCRIPTION
; ===========
;
; SEE ALSO
; ========
;
; :ref:`project_profile.sh`
;
; called by :ref:`basic_variable_evaluation_tpr_timeseries.pro`
;
; TODO
; ====
;
; supprimer ./paper01/fig3/read_basic_var.pro
; car ce module a l'air d'être le même
;
; make it work
;
; coding rules
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20111130T133857Z cratos (Linux)
;
;   * try to make it work on my account on cratos
;
;-
pro read_basic_var, csite,date1,date2,nsmooth, $
                     at, rh, sst, wu, wv, ws

@common
@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

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

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


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

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