;+
; .. _icoads_ws_stats_paper.pro:
;
; =========================
; icoads_ws_stats_paper.pro
; =========================
;
; DESCRIPTION
; ===========
;
; .. graphviz::
;
;    digraph icoads_ws_stats_paper {
;
;       ws [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/ws_monthly_1960_2010_oafluxgrid.nc"];
;       ws_tropflux [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/ws_TropFlux_global_monthly_1989_2009.nc"];
;       ws_oaflux [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/ws_OAFlux_global_monthly_1989_2009.nc"];
;       ws_erai [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/ws_ERAI_global_monthly_1989_2009.nc"];
;       ws_ncep2 [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/ws_NCEP2_global_monthly_1989_2009.nc"];
;       ws_ncep1 [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/ws_NCEP1_global_monthly_1989_2009.nc"];
;       ws_noc [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/nocv2_ws_19890101_20091231_oafluxgrid.nc"];
;
;       figure [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/icoads_ws_stats_paper.ps"];
;
;       icoads_ws_stats_paper [shape=box,
;       fontname=Courier,
;       color=blue,
;       URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/paper01/fig6/icoads_ws_stats_paper.pro",
;       label="${TROPFLUX}/src/paper01/fig6/icoads_ws_stats_paper.pro"];
;
;       {ws ws_tropflux ws_oaflux ws_erai ws_ncep2 ws_ncep1 ws_noc} -> {icoads_ws_stats_paper} -> {figure}
;    }
;
; SEE ALSO
; ========
;
; :ref:`project_profile.sh`
; :ref:`project_init.pro`
;
; :ref:`statistics.pro`
;
; EXAMPLES
; ========
;
; ::
;
;  IDL> icoads_ws_stats_paper
;
;
; TODO
; ====
;
; make it work (missing data, box definition)
;
; coding rules
;
; complete description
;
; handle IO error
;
; EVOLUTIONS
; ==========
;
; $Id$
;
; $URL$
;
; - fplod 20110420T092454Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * remove hard coding path
;   * add graphviz
;
; - fplod 20110411T142955Z aedon.locean-ipsl.upmc.fr (Darwin)
;
;   * minimal header
;
;-
pro icoads_ws_stats_paper
@cm_general
@cm_project
reinitplt, /z,/invert
key_portrait = 1
coefpalit=.9

openps, FILENAME = project_od_env+"icoads_ws_stats_paper.ps"
; partie a changer
marge=[-2,-2, -4,2]
st=19890101 & en=20091231
domdef, box

file=project_id_env+'ws_monthly_1960_2010_oafluxgrid.nc"
initncdf, file
icoads=read_ncdf("ws", st-1, en , file=file,/nostr, box=box)
help, icoads

file=project_id_env+"ws_TropFlux_global_monthly_1989_2009.nc"
initncdf, file
trop=read_ncdf("ws", st, en , file=file,/nostr, box=box)
help, trop

file=project_id_env+"ws_OAFlux_global_monthly_1989_2009.nc"
initncdf, file
oaf=read_ncdf("ws", st, en , file=file,/nostr, box=box)
help, oaf

file=project_id_env+"ws_ERAI_global_monthly_1989_2009.nc"
initncdf, file
erai=read_ncdf("ws", st, en , file=file,/nostr, box=box)
help, erai

file=project_id_env+"ws_NCEP2_global_monthly_1989_2009.nc"
initncdf, file
ncep2=read_ncdf("wsm", st, en , file=file,/nostr, box=box)
help, ncep2

file=project_id_env+"ws_NCEP1_global_monthly_1989_2009.nc"
initncdf, file
ncep1=read_ncdf("ws", st, en , file=file,/nostr, box=box)
help, ncep1

file=project_id_env+'nocv2_ws_19890101_20091231_oafluxgrid.nc'
initncdf, file
noc=read_ncdf("ws", st, en , file=file,/nostr, box=box)
help, noc

si=size(trop) & nt=si(3)

lat_rama=[67, 55, 55, 80.5, 80.5, 80.5, 80.5, 90, 90, 90, 90, 90, 90, 90, 95]
lon_rama=[-8, -12, -8,-8, -1.5, 0, 1.5, -1.5, 0, 1.5, 4, 8, 12, 15, -5]

lat_pirata=[0, 350, 337, 325, 350, 337, 322, 328, 322, 326, 322, 337, 337, 322, 350, 322, 330]
lon_pirata=[0, 0, 0, 0, -10, 12, 12, -14, 15, -19, 20, 21, 4, 4, -6, 8, -8]

lat_tao=[165, 180, 205, 235, 250, 265, 156, 165, 180, 190, 205, 220, 235, 250, 265, 156, 165, 180, 190, 205, 220, 235, 250, 265, 147, 156, $
         165, 180, 190, 205, 220, 235, 250, 265, 147, 156, 165, 180, 190, 205, 220, 235, 250, 265, 147, 156, 165, 190, 205, 220, 235, 250, 265, $
         156, 165, 180, 190, 220, 235, 250, 265]

lon_tao=[-8, -8, -8, -8, -8, -8, -5, -5, -5, -5, -5, -5, -5, -5, -5, -2, -2, -2, -2, -2, -2, -2, -2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, $
          2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 5, 5, 5, 5, 5, 5, 5, 5, 5, 8, 8, 8, 8, 9, 8, 8, 8]

lat=[lat_rama, lat_pirata, lat_tao]
lon=[lon_rama, lon_pirata, lon_tao]

lat_tao=lat & lon_tao=lon
nn=n_elements(lat)

list=0
for n=0, nn-1 do begin
  x=lat_tao(n)
  y=lon_tao(n)
  dx=abs(reform(glamt-x)) & dy=abs(reform(gphit)-y)
  ind=where((dx le 0.5) and (dy le 0.5))
  if (total(ind) ge 0.) then begin
    list=[list,ind]
  endif
endfor

msk=replicate(1,nxt,nyt)
msk(list)=!Values.f_nan

mask=trop*0
for jt=0,jpt-1 do begin
  t=reform(msk(*,*))
  mask(*,*,jt)=t
endfor

param=icoads*mask
trop=trop*mask & oaf=oaf*mask & erai=erai*mask
ncep2=ncep2*mask & ncep1=ncep1*mask
noc=noc*mask

ind=where(finite(param) and finite(noc) and finite(erai) and finite(trop) and finite(oaf) and finite(ncep2) and finite(ncep1))
x=param(ind) & y=trop(ind)
param=param(ind) & trop=trop(ind) & erai=erai(ind) & oaf=oaf(ind) & ncep2=ncep2(ind) & ncep1=ncep1(ind)
noc=noc(ind)

;
statistics, param, trop, $
                 cor, bias, std, rmsd
print, cor, bias, std,  rmsd
cstat_trop=string(cor, bias, std, rmsd, format='(f4.2,2x,f6.3,2x,f4.2,1x,f6.2)')
abs_diff=total(abs(trop-param))/n_elements(trop)
abs_trop=string(abs_diff, format='(f5.2)')
;
statistics, param, oaf, $
                 cor, bias, std, rmsd
print, cor, bias, std,  rmsd
cstat_oaf=string(cor, bias, std, rmsd, format='(f4.2,2x,f6.3,2x,f4.2,1x,f6.2)')
abs_diff=total(abs(oaf-param))/n_elements(oaf)
abs_oaf=string(abs_diff, format='(f5.2)')
;
statistics, param, erai, $
                 cor, bias, std, rmsd
print, cor, bias, std,  rmsd
cstat_erai=string(cor, bias, std, rmsd, format='(f4.2,2x,f6.3,2x,f4.2,1x,f6.2)')
abs_diff=total(abs(erai-param))/n_elements(erai)
abs_erai=string(abs_diff, format='(f5.2)')

;
statistics, param, ncep2, $
                 cor, bias, std, rmsd
print, cor, bias, std,  rmsd
cstat_ncep2=string(cor, bias, std, rmsd, format='(f4.2,2x,f6.3,2x,f4.2,1x,f6.2)')
abs_diff=total(abs(ncep2-param))/n_elements(ncep2)
abs_ncep2=string(abs_diff, format='(f5.2)')
;
statistics, param, ncep1, $
                 cor, bias, std, rmsd
print, cor, bias, std,  rmsd
cstat_ncep1=string(cor, bias, std, rmsd, format='(f4.2,2x,f6.3,2x,f4.2,1x,f6.2)')
abs_diff=total(abs(ncep1-param))/n_elements(ncep1)
abs_ncep1=string(abs_diff, format='(f5.2)')
;
statistics, param, noc, $
                 cor, bias, std, rmsd
print, cor, bias, std,  rmsd
cstat_noc=string(cor, bias, std, rmsd, format='(f4.2,2x,f6.3,2x,f4.2,1x,f6.2)')
abs_diff=total(abs(noc-param))/n_elements(noc)
abs_noc=string(abs_diff, format='(f5.2)')
;
mio=0 & mao=15 & inx=1 & nx=(mao-mio)/inx+1l & xx=mio+indgen(nx)*inx
mis=0 & mas=15 & iny=1 & ny=(mas-mis)/iny+1l & yy=mis+indgen(ny)*iny
xp=4 & yp=13
x1=0 & x2=15 & y1=x1 & y2=x2
;
pdf, param,trop,mio,mao,inx,mis,mas,iny, $
    pdf,xx,yy
scontour, pdf,xx,yy ,/noer, charsize=1, nlevels=30,/fill, $
    title="ws - ICOADS/TropFlux", small=[2,3,1], $
    xrange=[x1, x2], yrange=[y1, y2]
xyouts, xp, yp, cstat_trop, charsize=1
xyouts, xp, yp-1, abs_trop, charsize=1

;
pdf, param,oaf,mio,mao,inx,mis,mas,iny, $
    pdf,xx,yy
scontour, pdf,xx,yy ,/noer, charsize=1, nlevels=30,/fill, $
    title="ws - ICOADS/OAFlux", small=[2,3,2], $
     xrange=[x1, x2], yrange=[y1, y2]
xyouts, xp, yp, cstat_oaf, charsize=1
xyouts, xp, yp-1, abs_oaf, charsize=1

;
pdf, param,erai,mio,mao,inx,mis,mas,iny, $
    pdf,xx,yy
scontour, pdf,xx,yy ,/noer, charsize=1, nlevels=30,/fill, $
    title="ws - ICOADS/ERAI", small=[2,3,3], $
    xrange=[x1, x2], yrange=[y1, y2]
xyouts, xp, yp, cstat_erai, charsize=1
xyouts, xp, yp-1, abs_erai, charsize=1

;
pdf, param,ncep2,mio,mao,inx,mis,mas,iny, $
    pdf,xx,yy
scontour, pdf,xx,yy ,/noer, charsize=1, nlevels=30,/fill, $
    title="ws - ICOADS/NCEP2", small=[2,3,4], $
    xrange=[x1, x2], yrange=[y1, y2]
xyouts, xp, yp, cstat_ncep2, charsize=1
xyouts, xp, yp-1, abs_ncep2, charsize=1

;
pdf, param,ncep1,mio,mao,inx,mis,mas,iny, $
    pdf,xx,yy
scontour, pdf,xx,yy ,/noer, charsize=1, nlevels=30,/fill, $
    title="ws - ICOADS/NCEP1", small=[2,3,5], $
    xrange=[x1, x2], yrange=[y1, y2]
xyouts, xp, yp, cstat_ncep1, charsize=1
xyouts, xp, yp-1, abs_ncep1, charsize=1
;
pdf, param,noc,mio,mao,inx,mis,mas,iny, $
    pdf,xx,yy
scontour, pdf,xx,yy ,/noer, charsize=1, nlevels=30,/fill, $
    title="ws - ICOADS/NOC", small=[2,3,6], $
    xrange=[x1, x2], yrange=[y1, y2]
xyouts, xp, yp, cstat_noc, charsize=1
xyouts, xp, yp-1, abs_noc, charsize=1

closeps
end