further conversation on vertical plots

Hi NCL users,

  I am working on building a curtain plot of water vapor (vertical contour
plot). To do this and get a correct physical representation, I will need
to plot water vapor as a function of height rather than level. Currently,
the water vapor data is in (time, k, j , i) space. I would like to replace
the k (level) values with height (have already calculated this) so that I
can then use the contour plotting function. At present, I am only able to
plot k (level) on the y-axis. I am wondering if this is a doable request,
and if so, where am I going wrong in my code (see below)?

Many thanks!
Kara

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;This script is for a vapor contour plot where the axis are
;distance from the left edge and height. The location of this
;vertical slice is determined by setting J.
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
;;;;;;;;;;;;;;
; Read in data
;;;;;;;;;;;;;;
  I = 10 ;this is a generic number
  J = 97
  dx = 25
  Ti = 17
  Tf = 120
  Out = "x11"
  Title = "base_unstable_64A"

  f =
addfile("/mnt/pvfs2/karay/Rath_SA/StabilitySound_xy4/"+Title+"/wrfout_d01_2006-02-23_12:00:
00.nc" , "r")
  qv = f->QVAPOR(Ti,:,J,:) ;(t,k,j,i) Water vapor
  PH = f->PH(Ti,:,J,I) ;(t,k,j,i) Perturbation geopotential
(m2 s-2)
  PHB = f->PHB(Ti,:,J,I) ;(t,k,j,i) Base-state geopotential
(m2 s-2)

;*********************
; height of each layer
;*********************
nlayer = 36 ;number of vertical layers
nlayerm1 = nlayer - 1 ;number of vertical layers minus one

    ; height above ground level at full level [=] m
        height1D= PH(:)
        do ilevel = 0,(nlayer)
           height1D(ilevel) = (PH(ilevel)+PHB(ilevel))/9.8
        end do

    ; height above ground level at mid level [=] m
        height_mid1D= qv(:,1)
        do ilevel = 0,(nlayerm1)
           height_mid1D(ilevel) = 0.5 *
(height1D(ilevel)+height1D(ilevel+1))
        end do

;************************
;Convert to distance
;************************
printVarSummary(qv)

e_we = 640
x = fspan(0,(dx * e_we),e_we)
z = height_mid1D(:)
print(z)

z!(0) = "z"
z&z = z
z@units= "m"

x!(0) = "x"
x&x = x
x@units= "m"

qv!0 = "z"
qv!1 = "x"
qv@z = z
qv@x = x

slice = qv(z|::,x|::)

print(slice)

;************************
;Plot Resources
;************************

wks = gsn_open_wks(Out,Title)

  gsn_define_colormap(wks,"hotres") ; choose a colormap used hotres
before
  res = True ; plot mods desired
  res@lbLabelBarOn = False ; turn off individual label bars
  res@cnLineLabelsOn = False
  res@cnInfoLabelOn = False
  res@gsnDraw = False ; don't draw
  res@gsnFrame = False ; don't adnvance frame
  res@cnFillOn = True ; turn on color
  res@gsnSpreadColors = True ; use full range of colormap
  res@cnLinesOn = False ; turn off contour lines
  res@lbLabelAutoStride = True ;
  res@lbLabelAngleF = 45 ; angle labels
  res@lbLabelFontHeightF =.022 ; make labels larger
  res@lbTitleOn = True ; turn on title
  res@lbTitleFontHeightF = .015 ; make title smaller
  res@pmLabelBarOrthogonalPosF = .10 ; move whole thing down

;************************
;Make Plot
;************************

  plot = gsn_csm_contour(wks,slice,res)
  draw(plot)
  frame(wks)

end

-- 
"Don't knock the weather; nine-tenths of the people couldn't start a
conversation if it didn't change once in a while."  ~Kin Hubbard

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Jun 20 12:42:25 2012