;==============================================================
; Match the AO pattern at: 
;   http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/ao/ao.loading.shtml
;
; Use method: 
;   http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/ao/ao.shtml
;==============================================================
; Data source was Reanalysis-2 geopotential height.

;==============================================================

; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================
  islp=1
  iraw=1
  ncdf     = addfile("eof.20crv3.nc" ,"r")  ; open output netCDF file
  eof=ncdf->eof
  eof=eof/100.
  latS   = 20.
  latN   = 90.

  fclimo    = addfile ("/Datasets/20thC_ReanV3/Derived/Dailies/miscSI/prmsl.day.ltm.nc", "r")
  var    = "prmsl"
 
; GET DAILY CLIMOS. HANDLE LEAP YEARS

  xclimo   = ( fclimo->$var$(:,{latS:latN},:) )
  xclimon= array_append_record (xclimo(0:58,:,:), xclimo(58:58,:,:), 0)
  xclimom= array_append_record (xclimon, xclimo(59:364,:,:), 0)
  nn=366*(2015-1835)
  xout= new( nn, float)

  kk=0
  kindex=0
  k1=0
  k2=0

do iy=1836,2015
    if(islp.eq.1)then
     fmt="%4.0i"
     title  = str_upper(var)+" "
     if(iy.lt.1981)then
     filein="/Datasets/20thC_ReanV3/Dailies/miscSI/prmsl."+sprinti(fmt,iy)+".nc"
     else
     filein="/Datasets/20thC_ReanV3/Dailies/miscMO/prmsl."+sprinti(fmt,iy)+".nc"
end if
     fmean    = addfile (filein, "r")
    end if

; ==============================================================
; Open the file: Read only the user specified period and level
; ==============================================================
  
  if(islp.eq.1)then
    xall   = ( fmean->$var$(:,{latS:latN},:) )
    if((mod(iy,4).eq.0).and.(iy.ne.1900))then
      xall=xall/100.-xclimom/100.
      kindex=kindex+366
    else
      xall=xall/100.-xclimo/100.
      kindex=kindex+365
    end if
    if(iy.eq.1981)then
     k1=kindex-365
    end if
    if(iy.eq.2010)then
     k2=kindex
    end if
;    print("iy and k1,k2 "+iy+" "+k1+" "+k2)
  end if

; ==============================================================
; compute climatology and Anomalies
; ==============================================================
  rad    = 4.*atan(1.)/180.
  clat   = xall&lat            
  clat   = 1
  clat(35)=1

; =================================================================
; weight all observations 
; =================================================================

  xwall  = xall*conform(xall, clat, 1)
  copy_VarMeta(xall, xwall)
  xwall@long_name = "Wgt: "+xall@long_name

; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; eof_tsall needs the reordering.
; =================================================================

  wxall     = xwall(lat|:,lon|:,time|:)                ; convenience, cleaner code

  delete(xwall)
  delete(xall)
  optETS = False

  eof_tsall = eofunc_ts_Wrap (wxall, eof, optETS)
  eof_tsall=eof_tsall*-1.0	
  delete(wxall)

  k=0
  do im=1,12
    xmm=0.
    dim=days_in_month(iy,im)
    do id=1,dim
	xout(kk)=eof_tsall(0,k)
	k=k+1
	kk=kk+1
    end do
    end do
   delete(eof_tsall)
end do
klast=kk
;print("after years")
std=stddev(xout(:))
std2=stddev(xout(k1:k2))
;print("std and std1 "+std+" "+std2)

kk=0

; note; there is a slope in the AO

yDtrend = dtrend(xout(:klast-1),True)
;print("slope "+yDtrend@slope)

; now print the values out, standardize by the daily std

 print(" AO from 20CRV3 daily")
 do iy=1836,2015
  do im=1,12
    dim=days_in_month(iy,im)
    do id=1,dim
	print(iy+" "+im+" "+id+" "+xout(kk)/std2)
	k=k+1
	kk=kk+1
    end do
    end do
end do