;************************************************ ; These files are loaded by default in NCL V6.2.0 and newer ;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" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl" load "./psdfunc.ncl" begin ;---------------------------------------------------------------------- ; Time and Data Call ;---------------------------------------------------------------------- yrStrt = 1951 yrLast = 2015 season = "OND" latS = -10 ; Latitude of interest latN = 10 ; to make the average lonW = 50. lonE = 120. f = addfile("ersst.mnmean.nc","r") ;TIME = f->time ;YYYY = cd_calendar(TIME,-1)/100 ;iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast) ;sst = f->sst(iYYYY,{latS:latN},{lonW:lonE}) timeUnits = f->time@units ti = cd_inv_calendar(1951,01,01,00,0,0,timeUnits,0) tf = cd_inv_calendar(2015,12,31,00,0,0,timeUnits,0) sst = f->sst({ti:tf},{latS:latN},{lonW:lonE}) lat=f->lat({latS:latN}) lon=f->lon({lonW:lonE}) ; ============================================================== ; compute desired global seasonal mean: month_to_season (contributed.ncl) ; ============================================================== ssta_OND = month_to_season (sst, season) ; Extracts only the average of MAM. ;************************************************ ; detrend sst ;************************************************ ssta_OND = dtrend_msg_n(ispan(0,dimsizes(ssta_OND&time)-1,1),ssta_OND,True,False,0) ; Average in the latitude dimension (1) between latS and latN. ; With this, only the dimensions of [time|66]x[lon|36] remain. sstM = dim_avg_n_Wrap(ssta_OND,1) timeUnits = "days since 1800-01-01 00:00:00" sstM&time = ut_convert(sstM&time,timeUnits) var=sstM(lon|:,time|:) ;printVarSummary(var) ;printMinMax(var,0) ;printVarSummary(sstM&time) ;exit ntaper=0.1 var=dtrend(var,False) ; in my case I rewrote the var to be 2D (dim0 =points in grid space, dim1 = time); var=taper(var,ntaper,0) ntim = dimsizes(var&time) nsmooth=9 ; width of Daniell window ibase2 = floattointeger(log(ntim)/log(2.0)) + 1 npad = floattointeger(2.0^ibase2) - ntim ; padding primarily to increase frequency resolution ; any way I forgot to pad the data in the following -- take care of that by yourselves dan_win= Daniell_Window(nsmooth) sclfactor = 2.0/sum(dan_win^2) df = Calculate_Dof(dan_win,ntim,npad,ntaper) ; calculate Degrees of freedom,after accounting for effects of ;spectral smoothing, padding and tapering (see the notes from 'dmeko' above) cf = ezfftf(var)*dble2flt(sclfactor) psd = cf(0,:,:)^2 + cf(1,:,:)^2 ;- sum of power of real and imaginary parts nfrq = dimsizes(psd(0,:)) psd(:,nfrq-1) = psd(:,nfrq-1)/2.0 ; why? i forgot why ;() frq=ispan(0,nfrq,1)/int2flt(ntim) freq=frq(1:) ;---calculate red spectra rspc=rspec(var,freq) add_dimensions(rspc, (/"lon", "frq"/)) add_dimensions(psd,(/"lon","frq"/)) psd&frq=freq rscale = dim_sum(psd) rspc=rspc*conform(rspc,rscale,0) xHigh = chiinv(0.95, df)/df ; estimating 95% siglev rspc=rspc*xHigh ; peaks higher than this are significant at 95% sig lev rspc!0 = "lon" rspc!1 = "frq" rspc&lon=sst&lon ; -- smooth spectra with dan_win psd = wgt_runave(psd,dan_win,0) psd!0 = "lon" psd!1 = "frq" psd&lon=sst&lon ;psd&frq=psd&frq printVarSummary(psd) printMinMax(psd, 0) printMinMax(freq, 0) ;printVarSummary(rspc) ;printVarSummary(var) ; ================================================================= ; Plot ; ================================================================= wks = gsn_open_wks("png","hovmoller_psd") gsn_define_colormap(wks,"WhBlGrYeRe") ; choose a color map ;gsn_reverse_colormap(wks) res = True res@gsnMaximize = True res@cnFillOn = True res@cnLineLabelsOn = False res@cnInfoLabelOn = False res@cnLinesOn = False res@gsnSpreadColors =True ;res@gsnAddCyclic =True res@trXMaxF = 0.5 ; Maximum value of x. res@trXMinF = 0.0 ; Minimum value of x. res@tmXBMode = "Manual" ; X axis customization my way. res@tmXBTickStartF = res@trXMinF res@tmXBTickSpacingF = 0.1 res@tmXBTickEndF = res@trXMaxF res@tmYLMinorOn = True ; Disable minor tickmark on left y-axis (YL). res@tmXBMinorOn = True ; Disable minor tickmark on bottom x-axis (XB). ;contour line attributes res@cnLevelSelectionMode = "ExplicitLevels" ;res@cnLevels = (/60,80,100,120/) ;title res@gsnLeftString = "" res@gsnCenterString = "" res@gsnRightString = "" res@tiMainString = "" res@tiYAxisString = "" ;reshape color bars ;res@pmLabelBarWidthF =0.05 ;res@pmLabelBarHeightF =0.4 res@lbLabelBarOn =True res@lbOrientation ="vertical" res@pmLabelBarOrthogonalPosF=-0.04 ;res@lbOrientation ="horizontal" ;res@pmLabelBarOrthogonalPosF=0.1 res@lbLabelStride =2 res@lbBoxMinorExtentF =0.2 ;resize res@vpWidthF =0.60 res@vpHeightF =0.45 ;adjust X-AXis ;res@tmYLTickSpacingF =20 ;adjust Y-axis ;res@gsnXAxisIrregular2Linear =True ;res@trYReverse = True ;font size res@tiMainFontHeightF =0.022 res@txFontHeightF =0.02 res@tmXBLabelFontHeightF =0.018 res@tmYLLabelFontHeightF =0.018 ;fontsize of tickmark Labels res@tiXAxisFontHeightF =0.02 res@tiYAxisFontHeightF =0.02 ; Fixes special features for time shaft . ;resTick = True ;resTick@ttmFormat = "%Y" ;resTick@ttmAxis = "YL" ;resTick@ttmMajorStride = 10 ;time_axis_labels(psd&lon, res, resTick) plot = gsn_csm_hov(wks,psd,res) end