anyone can help me to plot the skew T

Hi group

    i'm trying to plot skew T digram , i got the plot but without the
pseudoadiabates curve (red line) like in example 2 in skewt examples , i have
attachment the figure i got it if anyone can help me.

thanks

Ali

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/skewt_func.ncl"

begin

colmin = 40
colmax = 47
rowmin = 33
rowmax = 40

class = 4

 diri = "/srv/science/maths1/z3236814/data02/swap/co2/fc/"

  f = addfile(diri+"evb_minclass_200noev_2000.nc","r") ; open data netcdf file

  topo = f->terrain
  w = f->w(class-1,:,:,:,:)
  u = f->u(class-1,:,:,:,:)
  v = f->v(class-1,:,:,:,:)
  t = f->t(class-1,:,:,:,:)
  q = f->q(class-1,:,:,:,:)
  pp = f->pp(class-1,:,:,:,:)
  sigma = f->sigma_level
  ps = f->pstarcrs
  time = f->time
  lat = f->latitcrs
  lon = f->longicrs
  ptop = f->ptop

  dsize = dimsizes(q)
  pres = new(dsize,float)

;to calculate the pressure for each grid point
 pres = pp + ptop +conform(pp,ps,(/2,3/))*conform(pp,sigma,1)

rel = relhum(t,q,pres)

 ;to calculates the dewpoint temperature
  dewpoint= dewtemp_trh(t,rel)

dewpoint = dewpoint -273 ;returned depoint temperature in (C)

 ; to define the dewpoint
 dewpoint!0 = "time"
 dewpoint!1 = "sigma_level"
 dewpoint!2 = "i_cross"
 dewpoint!3 = "j_cross"

 pres = pres*0.01 ; to converte it to (hPa)

;calculate virtual temperature
 vir_temp = new(dsize,float)
 vir_temp = t*(1.+q*0.61)

 ; to define the pres
 pres!0 = "time"
 pres!1 = "sigma_level"
 pres!2 = "i_cross"
 pres!3 = "j_cross"

 ; to define the virtual temperature
 vir_temp!0 = "time"
 vir_temp!1 = "sigma_level"
 vir_temp!2 = "i_cross"
 vir_temp!3 = "j_cross"

t = t-273 ; to converte it to (C)

 gh = hydro(pres(time|:,i_cross|:,j_cross|:,sigma_level|::-1), \
             vir_temp(time|:,i_cross|:,j_cross|:,sigma_level|::-1), \
             conform(vir_temp(:,0,:,:),topo,(/1,2/)))

 ; to define the geopotential height
 gh!0 = "time"
 gh!1 = "i_cross"
 gh!2 = "j_cross"
 gh!3 = "sigma_level"

;to calculates the meteorological wind direction
r2d = 45.0/atan(1.0) ;to conversion from radians to degrees because the return
value in radians
wdir = atan2(u,v)* r2d + 180

 wdir!0 = "time"
 wdir!1 = "sigma_level"
 wdir!2 = "i_cross"
 wdir!3 = "j_cross"

  wks = gsn_open_wks ("eps", "skewclass62095")

pres1=
dim_avg(dim_avg(pres(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
  z =
dim_avg(dim_avg(gh(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
wind =
dim_avg(dim_avg(v(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
windd=
dim_avg(dim_avg(wdir(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
 temp =
dim_avg(dim_avg(t(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))
dew1 =
dim_avg(dim_avg(dewpoint(time|:,sigma_level|:,i_cross|rowmin:rowmax,j_cross|colmin:colmax)))

  skewtOpts = True
  skewtOpts@DrawFahrenheit = False
  skewtOpts@DrawWind = True
   skewtOpts@DrawHeightScale = True

    dataOpts = True ; options describing data and ploting
    dataOpts@linePatternDewPt = 1
    dataOpts@colDewPt = "red"

;add legend
  legend = create "Legend" legendClass wks
    "vpXF" : 0.15
    "vpYF" : 0.88
    "vpWidthF" : 0.24 ; width
    "vpHeightF" : 0.09 ; height
    "lgPerimOn" : True ; perimeter
    "lgItemCount" : 2 ; how many
    "lgLabelStrings" : (/"Temperature","Dew point Temperature"/)
    "lgLineThicknessF" : 3.0
    "lgBoxMinorExtentF" : 0.2
  end create

  skewt_bkgd = skewT_BackGround (wks, skewtOpts)

 skewt_data=skewT_PlotData(wks,skewt_bkgd,pres1(15,:),temp(15,:),dew1(15,:),z(15,:),wind(15,:),windd(15,:),dataOpts)

   draw (skewt_data)
   draw (skewt_bkgd)

   draw(legend)
   frame(wks)

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk