;+
; .. _statistics_3var_v1.pro:
;
; ======================
; statistics_3var_v1.pro
; ======================
;
; DESCRIPTION
; ===========
;
; SEE ALSO
; ========
;
; :ref:`project_profile.sh`
;
; :ref:`net_flux_validation_scatter_2000_2007.pro`
; :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`
; :ref:`lwr_validation_scatter_2000_2007.pro`
; :ref:`swr_validation_scatter_2000_2007.pro`
;
; EXAMPLES
; ========
;
; ::
;
;  IDL> tao=++
;  IDL> ++
;  IDL> statistics_3var_v1, tao, var1, var2, cor1, cor2, bias1, bias2, std1, std2, rmsd1, rmsd2
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20110411T142955Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
;-
pro statistics_3var_v1, tao,var1,var2, $
    cor1, cor2, bias1, bias2, std1, std2, rmsd1, rmsd2
;
@cm_project
x=tao
y=var1
z=var2
;
ind1=where(finite(x,/nan))
ind2=where(finite(y,/nan))
ind3=where(finite(z,/nan))
ind=inter(ind1,union(ind2,ind3))
;
if (ind ne -1.) then begin
   x(ind)=!Values.f_nan
   y(ind)=!Values.f_nan
   z(ind)=!Values.f_nan
endif
;
nsmooth=5   ; applying a 5 day smoothing
;
x=ts_smooth(x,nsmooth,/nan)
y=ts_smooth(y,nsmooth,/nan)
z=ts_smooth(z,nsmooth,/nan)
;
ind1=where(finite(x))
ind2=where(finite(y))
ind3=where(finite(z))
ind=inter(ind1,inter(ind2,ind3))
x=x(ind)
y=y(ind)
z=z(ind)
nn=n_elements(ind)
;
cor1=correlate(x,y)
cor2=correlate(x,z)
; cor1=correlation(tropflux,tao) cor2=correlation(tropflux,rec)
xxx=moment(x)
; mean1 and std1 --> for tropflux
mean1=xxx(0)
std1=sqrt(xxx(1))
xxx=moment(y)
; mean2 and std2 --> for tao
mean2=xxx(0)
std2=sqrt(xxx(1))
xxx=moment(z)
; mean3 and std3 --> for rec
mean3=xxx(0)
std3=sqrt(xxx(1))
;
bias1=mean2-mean1
st1=std2/std1
bias2=mean3-mean1
st2=std3/std1
;
std1=st1
std2=st2
rmsd1=sqrt(total((x-y)*(x-y))/n_elements(x))
rmsd2=sqrt(total((x-z)*(x-z))/n_elements(x))
;
end