[ncl-talk] Strange output with vinth2p when P0>1000hPa

[Wade Qiang Wei]

Hi all,

I am trying to analyze cam model results in a CCSM3.0 simulation of an exoplanet with P0=1379hPa. It kept giving me weird data after converting from hybrid coordinates to pressure coordinates with vinth2p. So I turned off extrapolation and noticed every output with pnew>=1000hPa is filled with missing values. The original data has levels across 3hpa to 1379hpa and therefore extrapolation should not be necessary from my point of view.

Could vinth2p deal with P0 greater than 1000hPa? Has anyone encountered similar problems?

Thanks a lot.


Best Rards,
Wade

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

begin

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

;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 at units     = "hPa" 
  P0mb = 1379 ;mb
  K = 0.286
  g=13.5


;===============================================
 hyam = a->hyam                           ; get a coefficiants
 hybm = a->hybm                             ; get b coefficiants
 PS   = a->PS                     ; get pressure in Pa

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

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

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

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

  plot = gsn_csm_contour(wks,Tout,res)

end 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150518/731ea8a3/attachment.html