[ncl-talk] calculate the ensemble mean power spectrum

[Sri.durgesh Nandini-Weiss]

Dear ncl users,

I am trying to perform power spectral analysis of a time series of 
ensemble data (time*ens==>20years*100ensemble to compute over the year 
and each ensemble member and than calculate the ensemble mean of the 
spectra.

The best way would be to use specx_anal() and write it as a loop but 
this is where i am failing.

I would later on like to know how to add the Markov red noise and the 
upper and lower confidence levels?

would someone advice me on my script below?


begin

    f     = addfile ("anom_hist_zo_slr.nc", "r")
    hist_anom    = f->hist_anom
    printVarSummary(hist_anom)

    f1     = addfile ("anom_rcp45_zo_slr.nc", "r")
    rcp45_anom   = f1->rcp45_anom
    printVarSummary(rcp45_anom)

    dimx = dimsizes(hist_anom)
    ntim = dimx(0)           ; 240
    nens = dimx(1)           ; 100
    nlat = dimx(2)           ; 45
    nlon = dimx(3)           ; 90
    lat1=nlat
    lon1=nlon

    nmos = 12
    nyrs   = ntim/nmos       ; 20
    printVarSummary(nyrs)


;==================================================================
; print ("Data read in - start calculations")
; Region domain lat and lon, ensemble average
;==================================================================

   hist_anom1=hist_anom(:,:,{-7},{160}) ;(ntim, ens)
   printVarSummary(hist_anom1)
   aveX=dim_avg_n_Wrap(hist_anom1,(/1/))
   printVarSummary(aveX)     ; [time]

   rcp45_anom1=rcp45_anom(:,:,{-7},{160})                ;(ntim, ens)
   printVarSummary(rcp45_anom1)
   aveY=dim_avg_n_Wrap(rcp45_anom1,(/1/))
   printVarSummary(aveY)      ; [time]


   d   = 0                      ; detrending opt: 0=>remove mean 
1=>remove mean + detrend
   sm  = 7                    ; smooth periodogram:: should be at least 
3 and odd
                                    ; Average 7 periodogram estimates 
using modified Daniell
   pct = 0.10                 ; percent taper: (0.0 <= pct <= 1.0) 0.10 
common.

;==================================================================
   ;i can either use ndtooned and carry out specx_anal on this or make a 
loop.

   aveX = ndtooned(hist_anom1)
   printVarSummary(aveX)

   aveY = ndtooned(rcp45_anom1)
   printVarSummary(aveY)

   spec = specx_anal(aveX,d,sm,pct)
   printVarSummary(spec)

   spec1 = specx_anal(aveY,d,sm,pct)
   printVarSummary(spec1)


Thanx

Sri


-- 
Dr. Sri Nandini-Weiß
Research Scientist

Universität Hamburg
Center for Earth System Research and Sustainability (CEN)
Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)

Bundesstrasse 53, 20146 Hamburg
Tel: +49 (0) 40 42838 7472