Spectral analysis problems

I am trying to do spectral analysis on a netcdf file.

a. Firstly I want to find out what the major periodicities are in the
data (these may vary from months to several years). I am assuming
they are periodic in longitude. I would like to plot the range of
periods which are present, as a function of pressure and / or height
at a given latitude. I would also like to plot the period present as
a function of latitude at a given height or pressure. There are a lot
of possible routines available, and I am not sure if I am using the
right one(s) (see below).

b. I would also like to be able to find a way to follow a specific
wave as it travels upwards (or downwards), and to plot this. I am not
sure what the best routine to use would be, and how to plot this?.

in order to do a. above, I have tried to incorporate specx_anal into
a script, but I am not sure how I should be inputting the arguments,
and whether this will do what I want if I succeed ?. The error I am
getting is below, and underneath that is the script I have been
trying to use (please bear in mind that my pressure levels apply to
Venus, and the magnitudes are correct).

I am not sure how to do b. above/

Thanks,
Helen.

Here is the error message:

(0) T: interpolated and written to netCDF
(0) Q: interpolated and written to netCDF
(0) U: interpolated and written to netCDF
(0) V: interpolated and written to netCDF
fatal:syntax error: line 76 in file four.ncl before or near time
   spec = specx_anal(u(lev
-------------------------------------------------------^

fatal:Error in subscript, named subscripting is being used, make sure
each subscript has a name
fatal:syntax error: function specx_anal expects 4 arguments, got 0
fatal:error at line 76 in file four.ncl

fatal:Syntax Error in block, block not executed
fatal:error at line 114 in file four.ncl

Here is the script I have been trying to use:

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
                                        ; input file
   diri = "./"
   fili = "rand09125k.cam2.h0.0026-01-17-00000.nc"
   fi = addfile (diri+fili, "r")
                                        ; output file
   diro = "./"
   filo = "helen.nc"
   system ("/bin/rm -f "+diro+filo) ; remove any pre-exist file
   fo = addfile (diro+filo, "c")
                                        ; add any file attributes
   fo_at_title = fi_at_title + ": Selected Variables at pressure levels"
   fo_at_history = systemfunc ("date")
   fo_at_source = fi_at_source
   fo_at_case = fi_at_case
   fo_at_Conventions = fi_at_Conventions

   Var = (/ "T" , "Q", "U", "V"/) ; select variable to be
interpolated.
   nVar = dimsizes (Var)
                                        ; desired output levels
;lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100.,
  50.,25., 20., 15.,10., 8., 5.,3.,1. /)
;lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100.,
  50.,25., 14., 7., 3.5, 1.9, 0.95, 0.4, 0.1, 0.01 /)
lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100.,
50.,25., 14., 7., 3.5, 1.9, 0.95, 0.4, 0.1, 0.03 /)
;lev_p = (/ 91500., 90000., 87500., 80000., 67500., 50000., 30000.,
20000., 10000., 5000., 3000., 1500., 750., 390., 200., 100.,
  50.,25., 14., 7., 3.5, 2.0 /)
   lev_p!0 = "lev_p" ; variable and dimension name
the same
   lev_p&lev_p = lev_p ; create coordinate variable
   lev_p_at_long_name = "pressure" ; attach some attributes
   lev_p_at_units = "hPa"
   lev_p_at_positive = "down"

   hyam = fi->hyam ; read hybrid info
   hybm = fi->hybm
   PS = fi->PS
   P0mb = 0.01*fi->P0

   do n=0,nVar-1 ; loop over the variables
      X = fi->$Var(n)$
      Xp = vinth2p (X, hyam, hybm, lev_p ,PS, 1, P0mb, 2, True)
      copy_VarAtts(X, Xp)
      fo->$Var(n)$ = Xp ; write to netCDF file
      print (Var(n)+": interpolated and written to netCDF")
   end do

end
;***********************
; zonal.ncl
;***********************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;***********************
begin
   f = addfile ("helen.nc","r")
   print(f)
   u = f->U
   printVarSummary( u )

   dimu = dimsizes( u )
   ntim = dimu(0)
   klvl = dimu(1)
   nlon = dimu(2)
   mlon = dimu(3)

  uavg = dim_avg_Wrap(u)

;***********************
; Create Plot
;***********************
   d = 0
   sm = 3
   pct = 0.1

    kl = 10
    nl = 95
    ml = 40
   spec = specx_anal(u(lev|kl:kl,lat|nl:nl,lon|ml:ml,time:),d,sm,pct)

   wks = gsn_open_wks ("pdf", "rand09125k260117spec" ) ;
open workstation
   gsn_define_colormap(wks,"rainbow") ; choose colormap

   res = True
   res_at_cnFillOn = True
   res_at_lbLabelAutoStride = True
   res_at_gsnMaximize = True ; if [ps, eps, pdf] make large
   res_at_gsnSpreadColors = True ; span color map

   do nt=0,ntim-1

      res_at_gsnCenterString = "rand09125k260117spec"
      res_at_tiXaxisString = "Frequency"
      res_at_tiYaxisString = "Variance"

    plot = gsn_csm_xy(wks,spec_at_frq,spec_at_spcx,res)

; plot = gsn_csm_pres_hgt(wks, u(nt,:,:,40), res ) ; (lev,lat)
;; plot = gsn_csm_pres_hgt(wks, uavg(nt,:,:), res ) ; (lev,lat)
; plot = gsn_csm_contour (wks, uavg(t,:,:), res ) ; (lev,lat)

      res_at_trYReverse = True
; plot = gsn_csm_contour (wks, u(nt,:,:,0), res ) ; (lev,lat)
; plot = gsn_csm_contour (wks, u(nt,:,0,:), res ) ; (lev,lon)

      kl = 5
      res_at_mpFillOn = False
      res_at_mpGridAndLimbOn = True
      res_at_mpGridLineDashPattern = 2
      res_at_mpOutlineBoundarySets = "NoBoundaries"
      res_at_mpCenterLonF = 180.
      res_at_gsnCenterString = res_at_gsnCenterString+" p="+u&lev_p(kl)
; plot = gsn_csm_contour_map_ce (wks, u(nt,kl,:,:), res )

   end do

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Feb 02 2009 - 22:41:57 MST