; ==============================================================
; eof_1.ncl
;
; Concepts illustrated:
;   - Calculating EOFs
;   - Drawing a time series plot
;   - Using coordinate subscripting to read a specified geographical region
;   - Rearranging longitude data to span -180 to 180
;   - Calculating symmetric contour intervals
;   - Drawing filled bars above and below a given reference line
;   - Drawing subtitles at the top of a plot
;   - Reordering an array
; ==============================================================
; NCL V6.4.0 has new functions eofunc_n_Wrap and 
; eofunc_ts_n_Wrap that allow you to calculate the EOFs without 
; first having to first reorder the data. See eof_1_640.ncl.
; ==============================================================
; Calculate EOFs of the Sea Level Pressure over the North Atlantic.
; ==============================================================
; The slp.mon.mean file can be downloaded from:
; http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.surface.html
; ==============================================================
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================
  latS   = -40.
  latN   =  40. 
  lonL   =  30.
  lonR   =  120.

  yrStrt = 1870
  yrLast = 2019

  season = "SON"    ; choose Dec-Jan-Feb seasonal mean

  neof   = 3        ; number of EOFs
  optEOF = True       
  optEOF@jopt = 0   ; This is the default; most commonly used; no need to specify.
;;optEOF@jopt = 1   ; **only** if the correlation EOF is desired

  optETS = False

; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
  f      = addfile("Hadsstyrmon_ano.nc", "r")
  g      = addfile("sstson-ano.nc", "r")
  t      = g->time	
  YYYY   = cd_calendar(t,-1)/100                 ; entire file
  jYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
  sst1   = g->sst(jYYYY,:,:)
  printVarSummary(sst1)
  
 

  TIME   = f->time
  YYYY   = cd_calendar(TIME,-1)/100                 ; entire file
  iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)

  sst2    = f->sst(iYYYY,:,:)
  printVarSummary(sst2)                              ; variable overview

; ==============================================================
; dataset longitudes span 0=>357.5
; Because EOFs of the North Atlantic Oscillation are desired
; use the "lonFlip" (contributed.ncl) to reorder 
; longitudes to span -180 to 177.5: facilitate coordinate subscripting
; ==============================================================
  sst2    = lonFlip( sst2)
  printVarSummary(sst2)                              ; note the longitude coord
; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
; =================================================================
  rad    = 4.*atan(1.)/180.
  clat   = f->latitude
  clat   = sqrt( cos(rad*clat) )                 ; gw for gaussian grid

; =================================================================
; weight all observations 
; =================================================================
  wSLP2   = sst2                                   ; copy meta data
  wSLP2   = sst2*conform(sst2, clat, 1)
  wSLP2@long_name = "Wgt: "+wSLP2@long_name
  wSLP2!0 = "time"
  wSLP2!1 = "latitude"
  wSLP2!2 = "longitude"
  printVarSummary(wSLP2)

; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
  xw     = wSLP2({latitude|latN:latS},{longitude|lonL:lonR},time|:)
  x      = wSLP2(time|:,{latitude|latN:latS},{longitude|lonL:lonR})

  eof      = eofunc_Wrap(xw, neof, optEOF)      
  eof_ts   = eofunc_ts_Wrap (xw, eof, optETS)

  printVarSummary( eof )                         ; examine EOF variables
  printVarSummary( eof_ts )

; =================================================================
; Normalize time series: Sum spatial weights over the area of used
; =================================================================
  dimxw  = dimsizes( xw )
  mln    = dimxw(1)
  sumWgt = mln*sum( clat({latitude|latS:latN}) )
  eof_ts = eof_ts/sumWgt
;==================================================================
	;Calculate DMIV and its correlation with PC
;==================================================================
  weio = wgt_areaave_Wrap(sst1(:, {-10:10}, {50:70}), 1.0, 1.0, 0)
  eeio = wgt_areaave_Wrap(sst1(:, {-10:0}, {90:110}), 1.0, 1.0, 0)
  printVarSummary(weio)
  printVarSummary(eeio)
  DI = weio-eeio
  printVarSummary(DI)
  
  cc0 = escorc(eof_ts(0,:),DI(:))
  cc1 = escorc(eof_ts(1,:),DI(:))
  cc2 = escorc(eof_ts(2,:),DI(:))
  print(cc0)
  print(cc1)
  print(cc2)
  
  n = dimsizes(DI)
  df = n-2
  
  tt = cc0*sqrt((n-2)/(1-cc0^2))      
  p    = student_t(tt, df)
  psig = 0.05                       ; test significance level                     
  print("tt="+tt+"  p="+p)            ; t=2.02755  p=0.0732238
  if (p.le.psig) then
     print("cc0="+cc0+" is significant at the 95% level")
  else
     print("cc0="+cc0+" is NOT significant at the 95% level")
  end if
  
  tt = cc1*sqrt((n-2)/(1-cc1^2))      
  p    = student_t(tt, df)
  psig = 0.05                       ; test significance level                     
  print("tt="+tt+"  p="+p)            ; t=2.02755  p=0.0732238
  if (p.le.psig) then
     print("cc1="+cc1+" is significant at the 95% level")
  else
     print("cc1="+cc1+" is NOT significant at the 95% level")
  end if
  
  tt = cc2*sqrt((n-2)/(1-cc2^2))      
  p    = student_t(tt, df)
  psig = 0.05                       ; test significance level                     
  print("tt="+tt+"  p="+p)            ; t=2.02755  p=0.0732238
  if (p.le.psig) then
     print("cc2="+cc2+" is significant at the 95% level")
  else
     print("cc2="+cc2+" is NOT significant at the 95% level")
  end if
  
end