;==============================================================
; Match the AAO pattern at: 
;  JISAO
; 08/28/2023  from NOAA/PSL
;
; To run
;   ncl -Q -n newaao.20crv3.mon.ncl > aao.20crv3.long.data
; Use method: 
;==============================================================
; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================
  iyst   = 1836
  iyend  = 2015
  latS   = -90.
  latN   = -20.

  yrStrt = 1979
  yrLast = 2001
  
  cl1    = 1979
  cl2    = 2001

  var    = "prmsl"
  title  = str_upper(var)+" "
  f      = addfile ("/Datasets/20thC_ReanV3/Monthlies/miscMO/prmsl.mon.mean.nc", "r")

  ymStrt = yrStrt*100 +  1
  ymLast = yrLast*100 + 12

  neof   = 1                                  ; Leading EOF only
  optEOF = True       
  optETS = False

; ==============================================================
; Open the file: Read only the user specified period and level
; ==============================================================

  TIME   = f->time
  YYYY   = cd_calendar(TIME,-1)/100                 ; entire file
  iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
  prmsl  = ( f->$var$(iYYYY,{latS:latN},:) )
  xall   = ( f->$var$(:,{latS:latN},:) )


; ==============================================================
; compute climatology and Anomalies
; ==============================================================
  xClm   = clmMonTLL(prmsl)                             ; (12,lat,lon)
  xAnom  = calcMonAnomTLL ( prmsl, xClm)                ; (time, lat,lon)
  xAnomall=calcMonAnomTLL ( xall, xClm)                ; (time, lat,lon)

;  printMinMax(xAnom, True)

; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
; =================================================================
  rad    = 4.*atan(1.)/180.
  clat   = xAnom&lat            
  clat   = sqrt( cos(rad*clat) )                 ; gw for gaussian grid

  clat(35)=0

; =================================================================
; weight all observations 
; =================================================================
  xw     = xAnom*conform(xAnom, clat, 1)
  xwall  = xAnomall*conform(xAnomall, clat, 1)
  copy_VarMeta(prmsl, xw)
  copy_VarMeta(xall, xwall)
  xw@long_name = "Wgt: "+prmsl@long_name
  xwall@long_name = "Wgt: "+xall@long_name

; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Compute EOFs & Standardize time series
; =================================================================
  wx     = xw(lat|:,lon|:,time|:)                ; convenience, cleaner code
  wxall     = xwall(lat|:,lon|:,time|:)                ; convenience, cleaner code
  delete(xw)
  delete(xwall)

  eof    = eofunc_Wrap(wx, neof, optEOF)      
  eof    = -1*eof                                ; *special* match sign of CPC

  file_eof = "eof.aao.v3.nc"
  if(fileexists(file_eof))then
	system("rm "+file_eof)
  end if
  
  ncdf     = addfile(file_eof ,"c")  ; open output netCDF file
  ncdf->eof  = eof  
  eof_ts = eofunc_ts_Wrap (wx, eof, optETS)
  eof_tsall = eofunc_ts_Wrap (wxall, eof, optETS)
  eof_ts = dim_standardize_n( eof_ts, 0, 1)      ; normalize

;  standardized here using 1979-2001 values

k1=(cl1-iyst)*12
k2=-1+((cl2+1)-iyst)*12
 

 std=stddev(eof_tsall(:,k1:k2))
 mean=avg(eof_tsall(:,k1:k2))
; print(mean+" "+std+"k1 k2 "+k1+" "+k2)
 eof_tsall(:,:)=-1*(eof_tsall(:,:)-mean)/std

; =================================================================
; Regress
; =================================================================

  eof_regres = eof                               ; create an array w meta data
  eof_regresall = eof                               ; create an array w meta data
  do ne=0,neof-1
     eof_regres(ne,:,:) = (/ regCoef(eof_ts(ne,:), xAnom(lat|:,lon|:,time|:)) /)
     eof_regresall(ne,:,:) = (/ regCoef(eof_tsall(ne,:), xAnomall(lat|:,lon|:,time|:)) /)
  end do

  neof=1
  iflag=0
  n=0
  eof_regres(0,:,:)=-1*eof_regres(0,:,:)
  eof_tsall(0,:)=-1*eof_tsall(0,:)

; =================================================================
; Extract the YYYYMM from the time coordinate 
; associated with eof_ts [same as x&time] 
; =================================================================

  yyyymm = cd_calendar(eof_ts&time,-1)  
  yyyymmall=cd_calendar(eof_tsall&time,-1)  
  yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here
  yrfracall = yyyymm_to_yyyyfrac(yyyymmall, 0.0); not used here

;============================================================
; PLOTS
;============================================================
  wks = gsn_open_wks("ps","eof")
  gsn_define_colormap(wks,"BlueWhiteOrangeRed")                    
                                          ; only needed if paneling
; EOF patterns

  res = True                                     ; plot mods desired

  res@cnLinesOn            = True
  res@tiMainString         = "20CR V3 AAO Based on "+yrStrt+" to "+yrLast +" SLP" ; plot title
  res@gsnCenterString      = "AAO Spatial pattern"         ; center title

  res@cnFillOn             = True               ; turns on the color
  res@mpFillOn             = False              ; turns off continent gray
  res@cnLinesOn            = False              ; turn off contour lines
  res@cnFillOn             = True     ; turn on color fill
  res@gsnPolar             = "SH"               ; specify the hemisphere
  res@mpMinLatF            = -20                 ; specify min lat
  res@cnLineLabelsOn       = False        ; True is default
  res@lbLabelBarOn         = True        ; turn off individual lb's
  res@gsnAddCyclic=True
  res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levelsNH"
  res@cnLevels             =   (/ -450.,-400.,-350.,-300.,-250.,-200,-150.,-100.,-50.,50.,100.,150.,200.,250.,300.,350.,400.,450./)	; set levels
  res@cnFillColors         =  (/   8,  20,    32,   44,   56,   68,  80,   92,   104,  125, 153,  166, 177, 189  , 201, 213,  225,  237,250 /)

;*******************************************
	; first plot
;*******************************************

; get std monthly timeseries

 stddev_f = stddev(eof_tsall(0,:))

 n = 0

 res@gsnLeftString  = "EOF "+(n+1)
 res@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"

 contour = gsn_csm_contour_map_polar (wks, eof_regres(0,:,:),res)
  print(iyst+" "+iyend)
  do n = iyst ,iyend
    nn = (n - iyst)*12
    print(n +" "+sprintf("%7.3f",eof_tsall(0,nn))+" "+sprintf("%7.3f",eof_tsall(0,nn+1))+" "+sprintf("%7.3f",eof_tsall(0,nn+2))+" "+sprintf("%7.3f",eof_tsall(0,nn+3))+" "+sprintf("%7.3f",eof_tsall(0,nn+4))+" "+sprintf("%7.3f",eof_tsall(0,nn+5))+" "+sprintf("%7.3f",eof_tsall(0,nn+6))+" "+sprintf("%7.3f",eof_tsall(0,nn+7))+" "+sprintf("%7.3f",eof_tsall(0,nn+8))+" "+sprintf("%7.3f",eof_tsall(0,nn+9))+" "+sprintf("%7.3f",eof_tsall(0,nn+10))+" "+sprintf("%7.3f",eof_tsall(0,nn+11)))
  end do
  print(" -999.0")
  print ("Monthly AAO using 1979-2001 climo and standardization period. See ")
  print("https://psl.noaa.gov/data/20thC_Rean/data/timeseries/")
  print("created at NOAA/PSL")
  print("SLP is used to create EOF.")