Re: About Wrf_cape_3d.shtml function

The documentation says the fortran codes used are from RIP.
Regardless of the data source, if you input in the correct
variables *with the correct units*, the function should
calculate CAPE.

The CAPE code used by the skewT script is specific
to the skewT script. It can not be used independently.

--
Kerry Emanuel [MIT] makes his CAPE fortran code available at
http://wind.mit.edu/~emanuel/ftpacc.html
Perhaps, you can call that from NCL
http://www.ncl.ucar.edu/Document/Manuals/Getting_Started/beyondbasics.shtml#IncorporateYourOwn
Good luck
z3303149@unsw.edu.au wrote:
> Hello Dennis 
> 
>   i used these functions before but i got wrong result and when i asked the wrf
> help about the problem, they told me this functions was developed for wrf model
> and there is a large change that this will not work for other model output. i
> work in MM5 model.
> 
> Thanks 
> 
> Ali
> 
> 
> Hi group
> 
>     i'm trying to calculate the CAPE from skew T digram , all my variables 4D
> (time, level,lat,lon) i attachment the figure, it work just in (eps)  and i
> can't open it as (x11) or (pdf),but i got the figure without the red line of the
> parcel lapse rate like in example no. 2 in skew examples because i want to
> calculate the CAPE, if anyone can help me please.
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 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 = 6
> 
> 
> 
> ;1999
>  diri = "/srv/science/maths1/z3236814/data02/swap/co2/fc/"
> 
> 
> 
>   f = addfile(diri+"evb_minclass_200noev_2095.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
> 
> ----- End forwarded message -----
> 
> 
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk