<div dir="ltr"><div>This type of question was recently answered on ncl-talk.</div><div><br></div><div>To avoid repeated emails, I extracted the pertinent section.</div><div>PLEASE look-at and think-about what is being done!</div><div>Add your own print statements to examine any issues.</div><div><br></div><div>This has not been tested on your data. There may be errors.</div><div>It is your responsibility to fix them. <br></div><div>+++++++++++++++++++++++++++++++++++++++++++++++++<br></div><div>hist_anom1=hist_anom(:,:,{-7},{160}) Â  Â  ;(ntim, ens) ; (240,100)</div><div>=====</div><div>  nLen = ntim/2<br></div><div>  spec Â = new ( nLen, typeof(hist_anom1))<br>  spec Â = 0.0<br>  z1 Â  Â = 0.0<br>  xtsVar= 0.0<br><br>  iopt Â = 0<br>  jave Â = 1 Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Â Â  ; periodigram for each segment<br>  pct Â  = 0.1<br>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â ; loop over Â every ensemble member<br>  do ns=0,nens-1 Â  Â  Â  Â  Â  Â  Â  Â  ; spectra for each ensemble member<br>  Â  Â sdof Â  = specx_anal(hist_anom1(:.ns),iopt,jave,pct)<br>  Â  Â spec Â  = spec + sdof@spcx      ; sum spectra<br>  Â  Â xtsVar = xtsVar + sdof@xvari   ; sum ensemble member variances</div><div>                                                      ; sum z-traansformed lag-1 correlations<br></div><div>  Â Â  z1 Â  Â  = z1 Â  + 0.5*log((1+sdof@xlag1)/(1-sdof@xlag1)) ; Fischer z-transform<br>  end do<br><br>  spec Â = spec/nens Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Â  ; average spectra ==> ensemble mean spectra<br>  z1 Â  Â = z1/nens Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; ensemble mean z=transformed lag1<br>  xvari = xtsVar/nens Â  Â  Â  Â  Â  Â  Â  Â  Â Â  ; mean variance of ensemble members<br><br>  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; create a new variable<br>  ss Â  Â  Â  Â Â Â  = (/ sdof /)<br>  ss@dof Â  = 2*nens Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  ; # degrees freedom<br>  ss@bw Â  Â = sdof@bw<br>  ss@spcx Â = spec<br>  ss@xvari = xvari<br>  ss@frq Â  Â  = sdof@frq<br>  ss@xlag1 = (exp(2*z1)-1)/(exp(2*z1)+1) Â ; match specx_anal <br></div><div><br></div><div>  printVarSummary(ss)</div><div><br></div><div>;************************************************<br>; plotting parameters<br>;************************************************<br>  pltDir Â = "./"<br>  pltName = "ens_mean"<br>  pltPath = pltDir+pltName<br>  pltType = "png"<br><br>  wks Â = gsn_open_wks(pltType,pltPath) Â  <br><br>  r  = True Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â Â  Â  ; plot mods desired<br>  r@gsnMaximize Â  = True Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â  Â ; blow up plot<br>  r@tiMainString Â = "???" Â  Â  Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â  ; title<br>  r@tiXAxisString = "Frequency (cycles/day)" Â  Â ; xaxis<br>  r@tiYAxisString = "Variance" Â  Â  Â  Â  Â  Â  Â  Â  Â ; yaxis<br>  splt = specx_ci(ss, 0.05, 0.95) Â  Â  Â  Â  Â  Â  Â Â  ; calc confidence interval<br>  r@xyLineThicknesses Â  = (/2.,1.,1.,1./) Â  ; Define line thicknesses <br>  r@xyDashPatterns Â  Â  Â = (/0,0,1,1/) Â  Â  Â  Â Â  ; Dash patterns <br>  plot = gsn_csm_xy(wks,ss@frq, splt,r)<br></div><div><br></div><div>     <br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Mar 3, 2020 at 7:08 AM Sri.durgesh Nandini-Weiss via ncl-talk <<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Dear ncl users,<br>
<br>
I am trying to perform power spectral analysis of a time series of <br>
ensemble data (time*ens==>20years*100ensemble to compute over the year <br>
and each ensemble member and than calculate the ensemble mean of the <br>
spectra.<br>
<br>
The best way would be to use specx_anal() and write it as a loop but <br>
this is where i am failing.<br>
<br>
I would later on like to know how to add the Markov red noise and the <br>
upper and lower confidence levels?<br>
<br>
would someone advice me on my script below?<br>
<br>
<br>
begin<br>
<br>
    f     = addfile ("<a href="http://anom_hist_zo_slr.nc" rel="noreferrer" target="_blank">anom_hist_zo_slr.nc</a>", "r")<br>
    hist_anom    = f->hist_anom<br>
    printVarSummary(hist_anom)<br>
<br>
    f1     = addfile ("<a href="http://anom_rcp45_zo_slr.nc" rel="noreferrer" target="_blank">anom_rcp45_zo_slr.nc</a>", "r")<br>
    rcp45_anom   = f1->rcp45_anom<br>
    printVarSummary(rcp45_anom)<br>
<br>
    dimx = dimsizes(hist_anom)<br>
    ntim = dimx(0)           ; 240<br>
    nens = dimx(1)           ; 100<br>
    nlat = dimx(2)           ; 45<br>
    nlon = dimx(3)           ; 90<br>
    lat1=nlat<br>
    lon1=nlon<br>
<br>
    nmos = 12<br>
    nyrs   = ntim/nmos       ; 20<br>
    printVarSummary(nyrs)<br>
<br>
<br>
;==================================================================<br>
; print ("Data read in - start calculations")<br>
; Region domain lat and lon, ensemble average<br>
;==================================================================<br>
<br>
   hist_anom1=hist_anom(:,:,{-7},{160}) ;(ntim, ens)<br>
   printVarSummary(hist_anom1)<br>
   aveX=dim_avg_n_Wrap(hist_anom1,(/1/))<br>
   printVarSummary(aveX) Â  Â  ; [time]<br>
<br>
   rcp45_anom1=rcp45_anom(:,:,{-7},{160})                ;(ntim, ens)<br>
   printVarSummary(rcp45_anom1)<br>
   aveY=dim_avg_n_Wrap(rcp45_anom1,(/1/))<br>
   printVarSummary(aveY) Â  Â Â  ; [time]<br>
<br>
<br>
   d   = 0                      ; detrending opt: 0=>remove mean <br>
1=>remove mean + detrend<br>
   sm  = 7                    ; smooth periodogram:: should be at least <br>
3 and odd<br>
                             Â Â  Â  Â  ; Average 7 periodogram estimates <br>
using modified Daniell<br>
   pct = 0.10                 ; percent taper: (0.0 <= pct <= 1.0) 0.10 <br>
common.<br>
<br>
;==================================================================<br>
   ;i can either use ndtooned and carry out specx_anal on this or make a <br>
loop.<br>
<br>
   aveX = ndtooned(hist_anom1)<br>
   printVarSummary(aveX)<br>
<br>
   aveY = ndtooned(rcp45_anom1)<br>
   printVarSummary(aveY)<br>
<br>
   spec = specx_anal(aveX,d,sm,pct)<br>
   printVarSummary(spec)<br>
<br>
   spec1 = specx_anal(aveY,d,sm,pct)<br>
   printVarSummary(spec1)<br>
<br>
<br>
Thanx<br>
<br>
Sri<br>
<br>
<br>
-- <br>
Dr. Sri Nandini-Weiß<br>
Research Scientist<br>
<br>
Universität Hamburg<br>
Center for Earth System Research and Sustainability (CEN)<br>
Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)<br>
<br>
Bundesstrasse 53, 20146 Hamburg<br>
Tel: +49 (0) 40 42838 7472<br>
<br>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote></div>