cam4="./bbS866C200000CPfb8.cam2.h0.1901-1999.avg.nc"
;---Read data
  a   = addfile(cam4,"r")
  T   = a->T
  
  printVarSummary(T)
  printMinMax(T,0)

;original T&lev=xxxxxxxxxxxx, 960, 1086., 1195, 1281, 1338, 1368

  pnew = (/3,10,20,30,50,70,100,150,200,250,300,400,500,600,700,850,1000,1150,1250,1330,1379/)
  pnew@units     = "hPa" 
 ;P0mb = 1379 ;mb
 
  P0 = a->P0
  P0mb = P0*0.01
  P0mb@units = "hPa"
  print(P0mb)


;===============================================
 hyam = a->hyam                           ; get a coefficiants
 hybm = a->hybm                             ; get b coefficiants
 PS   = a->PS                     ; get pressure in Pa
 printVarSummary(PS)
 printMinMax(PS,0)

;************************************************
   interp = 2
; is extrapolation desired if data is outside the range of PS
   extrap = True
;************************************************
 T2 = vinth2p(T,hyam,hybm,pnew,PS,interp,P0mb,1,extrap)
 
 printVarSummary(T2)
 printMinMax(T2,0)

 Tout=dim_avg_n_Wrap(T2 (0,:, :,:), 2)
 
 printVarSummary(T2)
 printMinMax(T2,0)
 
 
 printVarSummary(Tout)
 printMinMax(Tout,0)


  wks_type="x11"
  wks = gsn_open_wks(wks_type,"Zonal_Mean_T_Slice"+"_top")

; -- set resources
  res = True
  res@cnFillOn = True ; -- turn on color fill
  res@cnLineLabelsOn = False ; -- turns off contour line labels
  res@cnInfoLabelOn = False ; -- turns off contour info label
  res@cnLinesOn = False
  res@tiYAxisString = "Pressure"+" [hPa]" ; -- append units to y-axis label
  res@trYReverse = True ; -- reverses y-axis
  res@gsnMaximize=True

  plot = gsn_csm_contour(wks,Tout,res)