CAPE problem

Hi Dennis

      this is the program of CAPE by using the wrf_cape_3d function and the
result from this program shows in attachment, so do you think the problem in my
program or in the function because different models.

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/wrf/WRFUserARW.ncl"
 begin

class = 4

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

 f = addfile(diri+"evb_minclass_200noev_2095.nc","r")

  topo = f->terrain
  t = f->t(class-1,:,:,:,:)
  q = f->q(class-1,:,:,:,:)
  pp = f->pp(class-1,:,:,:,:)
  sigma = f->sigma_level
  ps = f->pstarcrs
  ptop = f->ptop
  lon2d= f->longicrs
  lat2d=f->latitcrs
  time = f->time

  dsize = dimsizes(q)

 pres = new(dsize,float)

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

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

new_rel = relhum(t,q,pres)

 ;to calculate the specific humidity
 spec_hum = mixhum_ptrh(pres*0.01,t,new_rel,2)

;calculate virtual temperature
 vir_temp = t*(1.+spec_hum*0.61)

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

;calculate the geopotential height
 gh = hydro(pres(time|:,i_cross|:,j_cross|:,sigma_level|::-1)*0.01, \
             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 reorder it from (bottom to top) to (top to bottom) like in q,t and pres
z = gh(time|:,sigma_level|::-1,i_cross|:,j_cross|:)

; to calculate surface pressure
  psur=conform(pp(:,0,:,:),ps,(/1,2/))+ptop+pp(:,dsize(1)-1,:,:)
  psur = psur*0.01

ter = conform(psur,topo,(/1,2/))

cinfo= wrf_cape_3d(pres,t,q,z,ter,psur,True)

cape = cinfo(0,:,:,:,:)

; to calculate the average of the CAPE
 avCape = dim_avg(dim_avg(dim_avg(cape(time|:,sigma_level|:,i_cross|:,j_cross|:))))

  wks = gsn_open_wks("x11","cape")
  

 res = True
 res@tiMainString = "Convective Available Potential Energy (class 4)"
 res@tiMainFontHeightF = 0.015
 res@tiXAxisString = "Time Steps"
 res@tiXAxisFontHeightF = 0.015
 res@tiYAxisFontHeightF = 0.015
 res@tiYAxisString = "CAPE (J/kg)"
 res@xyLineThicknessF = 1.5
 res@tmXBMinorValues = ispan(1,25,1)

plot= gsn_xy(wks,ispan(1,25,1),avCape,res)

 
end

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