Interpolation Anomaly

I wanted to make a height versus time type plot for the zonally
averaged zonal wind. Data for the 44 different times is stored in
separate netcdf files (44 total).

In the first case I simply plotted the zonally averaged zonal wind as
a function of time in years on the x-axis, and model hybrid level on
the y-axis. (Note that this is not the Earth - there are 50 levels
total, and the pressure is 92 bars at the surface). (Also note that
the y-axis label is produced automatically). Since the highest
altitude level is level 0, and level 49 is at the surface, this
simple plot is effectively upside down, with level 0, representing
the highest altitude, at the bottom of the plot. The y-axis is also
not directly scaled according to pressure.

In the second case, I have tried to plot exactly the same results, as
a function of time in years on the x-axis, but with the y-axis values
given in terms of pressure / height. In order to do this I needed to
interpolate the data from hybrid coordinates onto pressure
coordinates, using NCL function vinth2p, for all 44 files. In this
case the plot is the right way up, with the surface at the bottom of
the plot.

The two plots below show the non-interpolated, and the interpolated
results respectively. The non-interpolated results look relatively
smooth, but the interpolated results show jagged structures, which do
not appear to have been in the original data. Also, I am not sure if
corresponding features are still at the same times as they were in
the non-interpolated version. This could either be due to something
that happens in the interpolation process, or perhaps because I have
done something incorrect in doing the interpolation.

If it is a problem with how the interpolation works, perhaps there
are some parameters I can adjust that may affect the way the results
look ?. I am sure that these jagged features should not be there.

I also include the two versions of the scripts I used. The script
"zon.ncl" is the non-interpolated version. The script "zoninterp.ncl"
it the interpolated version.

Can anyone see what is going wrong here ?.

Thanks,
Helen.

a). NON-INTERPOLATED PLOT :


b). INTERPOLATED PLOT :

c). SCRIPT TO MAKE THE NON-INTERPOLATED PLOT :

;***********************
; zon.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 zeusrand09125*.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", "zon" ) ; 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 ; equator

      res_at_gsnLeftString = ""
      res_at_gsnCenterString = "Zon. av. Zonal Wind, 0 deg lat, years
155-198"
      res_at_gsnRightString = ""

      res_at_tiXAxisString = "Time (years)"

      plot = gsn_csm_contour(wks, uavg_reverse(:,:,nl), res ) ;
(time,lev,lat)

      res_at_trYReverse = True

end

d). SCRIPT TO MAKE INTERPOLATED PLOT :

;***********************
; zoninterp.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 zeusrand09125*.cam2.h0*nc")
   print(fili)
   nfili = dimsizes(fili)

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

; ***********************************************
; Read in multiple files, interpolate, and write to temporary output
files
; **********************************************
   do nf = 0, nfili-1

   fi = addfile (diri+fili(nf), "r")

   diro = "./"
   filo = "uthelen"+nf+".nc"
   print(filo)
   system ("/bin/rm -f "+diro+filo) ; remove any pre-exist file
   fo = addfile (diro+filo, "c")

      hyam = fi->hyam ; read hybrid info
      hybm = fi->hybm
      hyai = fi->hyai ; read hybrid info
      hybi = fi->hybi
      PS = fi->PS
      P0mb = 0.01*fi->P0
      U = fi->U
      lat = fi->lat

   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"

      Up = vinth2p (U, hyam, hybm, lev_p ,PS, 1, P0mb, 2, False)
      copy_VarAtts(U, Up)
      fo->U = Up ; write to netCDF file

      end do
; **************************
; calculate zonal average of zonal wind versus time for multiple
interpolated files
; *************************
   diri = "./"
   fili2 = systemfunc("cd "+diri+" ; ls uthelen*nc")
   print(fili2)
   nfili = dimsizes(fili2)

   f2 = addfiles (fili2,"r")
   U2= addfiles_GetVar(f2,fili2,"U")

   printVarSummary( U2 ) ; (time, lev_p, lat,lon)

   dimt2 = dimsizes( U2 )
   ntim2 = dimt2(0)
   klvl2 = dimt2(1)
   nlat2 = dimt2(2)
   mlon2 = dimt2(3)

; take zonal average

  uavg = dim_avg_Wrap(U2) ;
(time,lev_p,lat)
   printVarSummary( uavg )

; create variable "work" of appropriate size and with appropriate
variable names

  ulatav = dim_avg_Wrap(uavg) ; (time,lev_p)
   printVarSummary( ulatav )
  ulatav_reorder = ulatav(lev_p | :, time | :) ; (lev_p,time)
   printVarSummary( ulatav_reorder )
  work = ulatav_reorder
   printVarSummary( work )

; reorder data to plot with gsn_csm_pres_hgt

  uavg_reorder = uavg(lev_p | :, time | :, lat | :) ;
(lev_p,lat,time)
   printVarSummary( uavg_reorder )

;***********************
; Create Plot
;***********************
  wks = gsn_open_wks ("pdf", "zoninterp" ) ; 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 ; equator

     work(lev_p|:,time|:) = uavg_reorder(:,:,nl) ;
(lev_p,time)
     printVarSummary( work )

      res_at_gsnLeftString = ""
      res_at_gsnCenterString = "Zon av. ZW, 0 deg lat, years 155-198"
      res_at_gsnRightString = ""

      res_at_tiXAxisString = "Time (years)"

      plot = gsn_csm_pres_hgt(wks, work, res ) ; (lev,time)

      res_at_trYReverse = True

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Apr 15 2009 - 23:46:14 MDT