Height/Pressure versus time plot

From: Helen Parish <hparish_at_nyahnyahspammersnyahnyah>
Date: Mon, 23 Feb 2009 23:09:22 -0800

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

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
(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 ?.


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"


; Keep all files open when have a lot of files :

   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
   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 ) ;

      res_at_trYReverse = True


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"
                                        ; 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
   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


ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
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