[ncl-talk] ensemble SOI index

[Sri nandini]

Hello ncl-users,

I had a question on computing indices on model output data which also 
has a large ensemble members.

My script and the attached plot is below. I do not get any errors but 
was hoping someone could tell me if the following is correct. The idea 
is to first remove the ensemble mean, and compute the index, in this 
case for the Nino 3 region.

Would be grateful!

Sri

   latS  = -5.0
   latN  =  5.0
   lonL  = -150.0        ;original data is longititude [-178.. 178]

   lonR  = -90.0
   
   nrun  = 5             ; length of running average

   yrStrt= 1870
   yrLast= 2005

   clStrt= 1971          ; climatology start
   clLast= 2000          ;             last
  

;-------------------- End User Input ---------------------------------------
   ymStrt = yrStrt*100 +  1
   ymLast = yrLast*100 + 12

   clStrt = clStrt*100 +  1              ; redefine
   clLast = clLast*100 + 12

;=======================================
; The index code below assumes that each year has 12 months
; model data is in format [time, 100 ensemble members, lat, lon]
;=======================================

   f        = addfile ("sst_1850-2005_ens_1-100.nc", "r")
   X      = mon_fullyear( f->sst(:,:,0,{latS:latN},{lonL:lonR}), 0)
   printVarSummary(X)
   YYYYMM = cd_calendar(X&time, -1)      ; ALL dates assciated with X
   tStrt  = ind(YYYYMM.eq.ymStrt)        ; indices of selected times
   tLast  = ind(YYYYMM.eq.ymLast)
   delete(YYYYMM)
   x      = X(tStrt:tLast,:,:,:)           ; subset to desired time interval
   yyyymm = cd_calendar(x&time, -1)
   dimx   = dimsizes(x)
   ntim   = dimx(0)
   printVarSummary(x)
   delete(X)                             ; no longer needed
   hist_anom1=x-273.15
   
   ;=======================================
    ;An index of the ensemble mean means nothing since these are supposed to characterise variability.
    ;So first remove the ensemble mean when calculating the indices such that there is no impact from a changing baseline.
  ;=======================================
    dimx = dimsizes(hist_anom1)
    ntim = dimx(0)          ; 240
    nens = dimx(1)

   hist_ensavg=dim_avg_n(hist_anom1,1)
   printVarSummary(hist_ensavg)
   hist = hist_anom1
   do n=0,nens-1
   hist(:,n,:,:) =  hist_anom1(:,n,:,:) - hist_ensavg
   end do
   printVarSummary(hist)
   printMinMax(hist,True)

;=======================================
; time indices for base climatology
;=======================================

   iClmStrt = ind(yyyymm.eq.clStrt)
   iClmLast = ind(yyyymm.eq.clLast)
   
;=======================================
; Climatology and anomalies from base climatology AFTER removing ensmean
;=======================================

   xClm     = clmMonTLLL(hist(iClmStrt:iClmLast,:,:,:))
   printVarSummary(xClm)

   xAnom    = calcMonAnomTLLL (hist,  xClm )
   xAnom at long_name = "SST Anomalies"
   printVarSummary(xAnom)
; The ensemble dimension is averaged out.
   xAnom1=dim_avg_n(xAnom,1)

;=======================================
; Unweighted areal averages & anomalies (time series)
; Small latitudinal extent so no need to weight
;=======================================

   xAnom_avg = wgt_areaave_Wrap(xAnom1, 1.0, 1.0, 1)
   xAnom_avg at long_name = "areal means of SST anomalies in the Nino 3 region"
   printVarSummary(xAnom_avg)
   
;=======================================
; Compute standardized anomalies; use clm period
;=======================================

   xAnom_std = xAnom_avg
   xAnom_std = xAnom_avg/stddev(xAnom_avg(iClmStrt:iClmLast))
   xAnom_avg at long_name = "areal avg standardized  anomalies"
   printVarSummary(xAnom_std)

; Perform an unweighted nrun-month running average on the index
; 2 months lost at start & end if endopt=0  ... reflective if endopt=1

   endopt    = 1
   ii = ind(.not.ismissing(xAnom_std))
   xAnom_std(ii) = runave_n_Wrap (xAnom_std(ii), nrun, endopt, 0)
   printVarSummary(xAnom_std)
;=======================================
   yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0)

   wks = gsn_open_wks("png","Hist_nino3_index_ensmean_rv")
   
   res                       = True
   res at gsnDraw               = False
   res at gsnFrame              = False
   res at vpHeightF             = 0.4            ; change aspect ratio of plot
   res at vpWidthF              = 0.8
   res at vpXF                  = 0.1            ; start plot at x ndc coord
   res at gsnYRefLine           = 0.0            ; create a reference line
   res at gsnAboveYRefLineColor = "red"              ; above ref line fill red
   res at gsnBelowYRefLineColor = "blue"             ; below ref line fill blue
   
   res at trYMinF               = -4.0           ; min value on y-axis
   res at trYMaxF               =  4.0           ; max value on y-axis
   res at gsnMaximize           = True
   
   res at gsnCenterString= "a) Hist: Nino3: 1870-2005 Base:1971-2000"

   res at tiYAxisString    = "SST Standardized Anomalies"    ; y-axis label
   plot = gsn_csm_xy (wks,yrfrac,xAnom_std,res)
   draw(plot)
   frame (wks)


-------------- next part --------------
A non-text attachment was scrubbed...
Name: hist_soi.png
Type: image/png
Size: 64056 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210707/7a3570ed/attachment.png>