<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>