Re: Height/Pressure versus time plot

Hi Helen,

First, gsn_csm_pres_hgt requires that the 1st dimension of the input
array be in units of Pa, hPa, or mb as the error message indicates. The
units attribute for the pressure coordinate variable must also be set
correctly for the function to work. I am guessing that your input array
isn't on pressure levels, and is actually on CAM hybrid_sigma levels,
which would cause the error messages you are reporting.

So, you do need to interpolate to pressure levels if you want to use
gsn_csm_pres_hgt. However, if you just want to plot the data on the
current hybrid-sigma levels, then you can skip this step by simply using
  gsn_csm_contour. Take a look at the examples on the Contours w/o Maps
Applications page here:
http://www.ncl.ucar.edu/Applications/conwomap.shtml

If you do want/need to interpolate to pressure levels, then the easiest
way to shorten the interpolation time is to reduce the number of
timesteps and/or the number of pressure levels that are interpolated to,
neither of which you likely want to do.

Perhaps others will have suggestions on how to speed up the interpolation..
Adam

Helen Parish wrote:
> I am aiming to plot a quantity as a function of pressure/height
> versus time (in this case zonal average zonal velocity as a function
> of pressure/height and time).
>
> In order follow the changes as a function of time I need to read in
> multiple netcdf files. In this case I am reading in 45 files (the
> maximum I can read in without running into local computer memory
> limits).
>
> I was thinking that I can make use of routine "gsn_csm_pres_hgt" to
> do the plot, provided that I reorder the subscripts from (time, lev,
> lat) to (lev,time,time), since routine "gsn_csm_pres_hgt" expects a
> height type variable in the leftmost position.
>
> However, since I am producing results for Venus, I believe I need to
> interpolate onto a new set of defined pressure levels for each file
> (the upper and lower pressure limits are very different from those on
> the Earth) . I am not sure how to do this within the context of my
> current script.
>
> I am currently getting the following error when I try to run my script :
>
> "(0) gsn_csm_pres_hgt: Fatal: The coordinate array for the first
> dimension of the input data must be in Pascals, Hecto-pascals, or
> millibars
> (0) and it must contain the attribute 'units' set to one of the
> following strings (depending on your units):
> (0) 'mb' 'Mb' 'MB' 'millibar' 'millibars' 'MILLIBARS'
> 'hybrid_sigma_pressure' 'Pa' 'pa' 'PA' 'Pascals' 'pascals' 'PASCALS'
> 'hpa' 'hPa' 'HPA' 'hecto-pascals' 'HECTO-PASCALS'
> (0) Cannot create plot.
> fatal:Illegal right-hand side type for assignment
> fatal:Execute: Error occurred at or near line 54 in file timely.ncl "
>
> Below, I include my current script. Following this, I include the
> first few lines of another script I have been using to interpolate
> onto a new set of pressure levels appropriate for Venus. I am not
> sure which lines (if any) I need to add from the second script into
> the first script, and specifically how exactly to do this, so that it
> will understand my pressure levels correctly.
>
> I am also hoping to avoid the whole script becoming unmanageable, if
> I add interpolation for every one of the files to my current script.
> It currently takes about 10 to 15 minutes for the script to run, and
> takes so much processing power that nothing else can be done on the
> computer whilst the script is running, when I generate a simple
> periodogram with this number of files. I know that it also takes a
> while to do the interpolation for a single file, which will be
> multiplied by 45 if I have to do this for all 45 files.
>
> Can anyone please advise me on what I need to do to my script to get
> a height / pressure level versus time plot to work ?.
>
> Thanks,
> Helen.
>
> 1) This is the script I am currently using :
>
> ;***********************
> ; timely.ncl
> ;***********************
> 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"
>
> begin
>
> ; Keep all files open when have a lot of files :
> setfileoption("nc","SuppressClose",False)
>
> diri = "./"
>
> fili = systemfunc("cd "+diri+" ; ls rand09125*.cam2.h0*nc")
>
> f = addfiles (fili,"r")
> U = addfiles_GetVar(f,fili,"U")
> lat = addfiles_GetVar(f,fili,"lat")
> lon = addfiles_GetVar(f,fili,"lon")
> lev = addfiles_GetVar(f,fili,"lev")
> time = addfiles_GetVar(f,fili,"time")
> printVarSummary( U )
>
> dimu = dimsizes( U )
> ntim = dimu(0)
> klvl = dimu(1)
> nlon = dimu(2)
> mlon = dimu(3)
>
> uavg = dim_avg_Wrap(U)
>
> uavg_reverse = uavg(lev | :, time | :, lat | :)
>
> ;***********************
> ; Create Plot
> ;***********************
> wks = gsn_open_wks ("pdf", "presstimetest" ) ; 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
>
> nl = 96
>
> res_at_gsnCenterString = "zeusorigspongej381121"
> ;; plot = gsn_csm_pres_hgt(wks, uavg(nt,:,:), res ) ; (lev,lat)
> plot = gsn_csm_pres_hgt(wks, uavg_reverse(:,:,nl), res ) ;
> (time,lev,lat)
>
> res_at_trYReverse = True
>
> end
>
>
> 2) These are the first few lines of another script which I have used
> to interpolate onto the new pressure levels, for a single file :
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> begin
> ; input file
> diri = "./"
> fili = "zeusorigspongej.cam2.h0.0038-11-21-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.,
> 14., 7., 3.5, 1.9, 0.95, 0.4, 0.1, 0.03 /)
>
> 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
> ;***********************
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
--------------------------------------------------------------
Adam Phillips			             asphilli_at_ucar.edu
National Center for Atmospheric Research   tel: (303) 497-1726
ESSL/CGD/CAS                               fax: (303) 497-1333
P.O. Box 3000				
Boulder, CO 80307-3000	  http://www.cgd.ucar.edu/cas/asphilli
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk