;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;   Script to output ERA40 T2m for remove 9-year running  ;
;   mean with DJF varances.                               ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;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"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;
;---------------------------------------------------------------------------
; There are four region based indices used to monitor the tropical Pacific: 
; Nino 1+2 (0-10S, 90W-80W), Nino 3 (5N-5S, 150W-90W), 
; Nino 3.4/ONI (5N-5S, 170W-120W) and Nino 4 (5N-5S, 160E-150W)
;
; NOAA's operational definitions of El Nino and La Nina conditions are based
; upon the Oceanic Nino Index [ONI]. The ONI is defined as the 3-month running
; means of SST anomalies in the Nino 3.4 region [5N-5S, 120-170W]. The anomalies
; are derived from the 1971-2000 SST climatology.
;---------------------------------------------------------------------------
begin
;
  latS   = -5.0
  latN   =  5.0
  lonL   = 190.0
  lonR   = 240.0
  
  nrun   = 3             ; length of running average

  run1   = getenv("MODE")
  ymStrt = (stringtointeger(getenv("ISYY")))*100+1
  ymLast = (stringtointeger(getenv("IEYY")))*100+12

  model  = getenv("model")
  dir    = getenv("path")
;  dir    = "/archive/CMIP6/TaiESM1/"+run1+"/atmos/mon/r1i1p1f1/"
;  dir    = "/lfs/archive/taiesm_ppost/"+run1+"/v20210908/"
;---------------------------------------------------------------------
  ifn    = systemfunc("ls "+dir+"ts_Amon_*.nc")
  a      = addfiles(ifn,"r")
  data   = a[:]->ts;(:,{latS:latN},{lonL:lonR})
;---------------------------------------------------------------------
  yyyymm =cd_calendar( data&time,-1)
  delete(data&time)
  data&time = yyyymm
  wgt = (/1.,2.,1./)
  wgt = wgt/sum(wgt)

;*********************************
; Climatology and anomalies from base climatology
;*********************************

;  printVarSummary(xClm)

  tStrt = ind(yyyymm.eq.ymStrt)      ; indices of selected times
  tLast = ind(yyyymm.eq.ymLast)

;  print(date)
  x = data(tStrt:tLast,{latS:latN},{lonL:lonR})

  xClm     = clmMonTLL(x({ymStrt:ymLast},:,:))
  xAnom    = calcMonAnomTLL (x({ymStrt:ymLast},:,:),  xClm ) 
  xAnom@long_name = "SST Anomalies"
  printVarSummary(xAnom)

;*********************************
; Unweighted areal average anomalies (time series)
; Small latitudinal extent so no need to weight    
;*********************************

  xAnom_avg = wgt_areaave_Wrap(xAnom, 1.0, 1.0, 1)
  xAnom_avg@long_name = "areal avg anomalies"
  printVarSummary(xAnom_avg)

;*********************************
; Perform an unweighted 'nrun' month running average
;*********************************

   xAnom_avg = runave_n_Wrap (xAnom_avg, nrun, 1, 0)
   wAnom_avg = wgt_runave_n_Wrap (xAnom_avg, wgt, 1, 0)
   nino34_ndj = wAnom_avg(11:dimsizes(wAnom_avg)-13:12)   ; cannot count last 1yr as spatial composite uses +1yrs data betond NDJ..
   nino34_ndj!0 = "time"
   nino34_ndj = dtrend_msg(ispan(0,dimsizes(nino34_ndj&time)-1,1),nino34_ndj,True,False)
   nino34_ndj = dim_standardize(nino34_ndj,0)


   printVarSummary(nino34_ndj)
   print("time "+nino34_ndj&time+" :"+nino34_ndj)
;xit
   hicntr = 0
   locntr = 0
   hiyr = new(dimsizes(nino34_ndj&time),integer)
   loyr = hiyr

;           print(nino34_ndj)

           do hh = 0,dimsizes(nino34_ndj)-1
              if (nino34_ndj(hh).ge.1) then
                 ; print("High year = "+nino34_ndj&time(hh))
                 hiyr(hicntr) = nino34_ndj&time(hh)
                 hicntr = hicntr+1
              end if
              if (nino34_ndj(hh).le.-1) then
                 ; print("Low year = "+nino34_ndj&time(hh))
                 loyr(locntr) = nino34_ndj&time(hh)
                 locntr = locntr+1
              end if
           end do

   
    highyr = hiyr(:hicntr-1)
    delete(hiyr)
    lowyr = loyr(:locntr-1)
    delete([/loyr,hicntr,locntr/])

   print("highyr:"+highyr)
   print("lowyr:"+lowyr)
      
  ofn = model+"_"+run1+"_nino_case_"+ymStrt+"-"+ymLast+".nc"
  system("rm -rf "+ofn)
  odd = addfile(ofn,"c")
  odd->highyr = highyr
  odd->lowyr  = lowyr
exit

;do gg = 0,23
;   tt = gg/12.
;   print("gg="+gg+" ,tt="+tt)
;  dd= highyr+tt
;   print(dd)
;end do
;exit

;  print(xAnom_avg +" / "+ wAnom_avg)
;exit 
;*********************************
; output data
;*********************************
  ofn = "CAM5.3_"+ymStrt+"-"+ymLast+"_nino34_ONI.nc"
  system("rm -rf "+ofn)
  odd = addfile(ofn,"c")
  odd->oni = xAnom_avg
;*********************************
; plot graph
;*********************************
  yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0)

  wks = gsn_open_wks("pdf", "CAM5.3_nino34_"+ymStrt+"-"+ymLast)   

  res                  = True       
  res@gsnMaximize      = True
  res@gsnYRefLine      = 0.0            ; create a reference line 
  res@gsnAboveYRefLineColor = "red"              ; above ref line fill red
  res@gsnBelowYRefLineColor = "blue"             ; below ref line fill blue

  res@vpHeightF        = 0.4            ; change aspect ratio of plot
  res@vpWidthF         = 0.8            
  res@vpXF             = 0.1            ; start plot at x ndc coord 
  
  res@trYMinF          = -2.0           ; min value on y-axis
  res@trYMaxF          =  2.0           ; max value on y-axis

  res@tiMainString     = "CPC Ocean Nino Index"                            
  
  res@tiYAxisString    = "Anomalies (C)"    ; y-axis label      
  plot = gsn_csm_xy (wks,yrfrac,xAnom_avg,res)

end