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
Received on Tue Feb 24 2009 - 00:09:22 MST
This archive was generated by hypermail 2.2.0 : Tue Mar 03 2009 - 09:53:57 MST