ncl 原来的函数有问题，改成现在的：
function wrf_user_intrp3d( var3d:numeric, z_in:numeric, \
plot_type:string, \
loc_param:numeric, angle:numeric, isInUnitPressure:logical,opts:logical )
其中，多加了一个isInUnitPressure:logical,用于指示是否使用的Ｐ坐标。（注意，单位是Ｐａ，不是ｈＰａ）

原来的函数，使用z_in（0,0,0,0,0）来决断是否使用的P坐标。如果>500，则认为是P坐标，这显然是不合理的，海拔大于500m的地方很多！
zuohj, NingXiaCC, 2011.4.28

 
 
 undef("wrf_user_intrp3d")
 functionwrf_user_intrp3d( var3d:numeric, z_in:numeric, \
 plot_type:string, \
 loc_param:numeric, angle:numeric, isInUnitPressure:logical,opts:logical )
 
 ; var3d - 3d field to interpolate (all input fields must be unstaggered)
 ; z_in - interpolate to this field (either p/z)
 ; plot_type - interpolate horizontally "h", or vertically "v"
 ; loc_param - level(s) for horizontal plots (eg. 500hPa ; 3000m - scalar), 
 ; plane for vertical plots (2 values representing an xy point 
 ; on the model domain through which the vertical plane will pass
 ; OR 4 values specifying start and end values
 ; angle - 0.0 for horizontal plots, and 
 ; an angle for vertical plots - 90 represent a WE cross section
 ; opts Used IF opts is TRUE, else use loc_param and angle to determine crosssection
 
 begin
 
 if(plot_type .eq. "h" ) then; horizontal cross section needed
 
 dimL = dimsizes(loc_param)
 
 dims = dimsizes(var3d)
 nd = dimsizes(dims)
 
 dimX = dims(nd-1)
 dimY = dims(nd-2)
 dimZ = dims(nd-3)
 dim4 = 1
 dim5 = 1
 if ( nd .eq. 4 ) then
 dim4 = dims(nd-4)
 endif
 if ( nd .eq. 5 ) then
 dim4 = dims(nd-4)
 dim5 = dims(nd-5)
 endif
 
 var3 = new ( (/ dim5, dim4, dimZ, dimY, dimX /) , typeof(var3d) )
 z = new ( (/ dim5, dim4, dimZ, dimY, dimX /) , typeof(var3d) )
 var2d = new ( (/ dim5, dim4, dimL, dimY, dimX /) , typeof(var3d) )
 
 if ( nd .eq. 5 ) then
 var3 = var3d
 z = z_in
 endif
 if ( nd .eq. 4 ) then
 var3(0,:,:,:,:) = var3d(:,:,:,:)
 z(0,:,:,:,:) = z_in(:,:,:,:)
 endif
 if ( nd .eq. 3 ) then
 var3(0,0,:,:,:) = var3d(:,:,:)
 z(0,0,:,:,:) = z_in(:,:,:)
 endif
 
 
 if ( isInUnitPressure ) then
 z = z * 0.01
 loc_param = loc_param * 0.01
 endif
 
 do il = 0,dimL-1
 var = wrf_interp_3d_z(var3,z,loc_param(il))
 var2d(:,:,il,:,:) = var(:,:,:,:)
 enddo
 
 copy_VarAtts(var3d,var3)
 if(isatt(var3,"description")) then
 delete_VarAtts(var3,(/"description"/))
 endif
 if(isatt(var3,"units")) then
 delete_VarAtts(var3,(/"units"/))
 endif
 if(isatt(var3,"MemoryOrder")) then
 delete_VarAtts(var3,(/"MemoryOrder"/))
 endif
 if(isatt(var3,"_FillValue")) then
 delete_VarAtts(var3,(/"_FillValue"/))
 endif
 copy_VarAtts(var3,var2d)
 
 nn = nd-2
 var2d!nn = "plevs"
 
 if ( dimL .gt. 1 ) then
 if ( nd .eq. 5 ) then
 return( var2d )
 endif
 if ( nd .eq. 4 ) then
 return( var2d(0,:,:,:,:) )
 endif
 if ( nd .eq. 3 ) then
 return( var2d(0,0,:,:,:) )
 endif
 else
 if ( z(0,0,0,0,0) .gt. 500.) then
 var2d@PlotLevelID = loc_param + " hPa"
 else
 var2d@PlotLevelID = .001*loc_param + " km"
 endif
 if ( nd .eq. 5 ) then
 return( var2d(:,:,0,:,:) )
 endif
 if ( nd .eq. 4 ) then
 return( var2d(0,:,0,:,:) )
 endif
 if ( nd .eq. 3 ) then
 return( var2d(0,0,0,:,:) )
 endif
 endif
 
 
 endif
 
 
 
 
 if(plot_type .eq. "v" ) then; vertical cross section needed
 
 dims = dimsizes(var3d)
 
 if ( dimsizes(dims) .eq. 4 ) then
 if ( z_in(0,0,0,0) .gt. 500.) then
  ; We must be interpolating to pressure
  ; This routine needs input field and level in hPa - lets make sure of this
 if ( z_in(0,0,0,0) .gt. 2000. ) then
  ; looks like we have Pa as input - make this hPa
 z_in = z_in * 0.01
 endif
 endif
 z = z_in(0,:,:,:)
 else
 if ( isInUnitPressure) then
 z_in = z_in * 0.01
 endif
 z = z_in
 endif
 
 
 ; set vertical cross section
 if (opts) then
 xy = wrf_user_set_xy( z, loc_param(0)-1, loc_param(1)-1, \ ; the -1 is for NCL dimensions
 loc_param(2)-1, loc_param(3)-1, \
 angle, opts )
 else
 xy = wrf_user_set_xy( z, loc_param(0), loc_param(1), \
 0.0, 0.0, angle, opts )
 endif
 xp = dimsizes(xy)
 
 
 ; first we interp z
 var2dz = wrf_interp_2d_xy( z, xy)
 
 ; interp to constant z grid
 if(var2dz(0,0) .gt. var2dz(1,0) ) then; monotonically decreasing coordinate
 z_max = floor(max(z)/10)*10; bottom value
 z_min = ceil(min(z)/10)*10; top value
 dz = 10
 nlevels = tointeger( (z_max-z_min)/dz)
 z_var2d = new( (/nlevels/), typeof(z))
 z_var2d(0) = z_max
 dz = -dz
 else
 z_max = max(z)
 z_min = 0.
 dz = 0.01 * z_max
 nlevels = tointeger( z_max/dz )
 z_var2d = new( (/nlevels/), typeof(z))
 z_var2d(0) = z_min
 endif
 
 do i=1, nlevels-1
 z_var2d(i) = z_var2d(0)+i*dz
 enddo
 
 
 ; interp the variable
 if ( dimsizes(dims) .eq. 4 ) then
 var2d = new( (/dims(0), nlevels, xp(0)/), typeof(var2dz))
 do it = 0,dims(0)-1
 var2dtmp = wrf_interp_2d_xy( var3d(it,:,:,:), xy)
 do i=0,xp(0)-1
 var2d(it,:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d)
 enddo
 enddo
 var2d!0 = var3d!0
 var2d!1 = "Vertical"
 var2d!2 = "Horizontal"
 else
 var2d = new( (/nlevels, xp(0)/), typeof(var2dz))
 var2dtmp = wrf_interp_2d_xy( var3d, xy)
 do i=0,xp(0)-1
 var2d(:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d)
 enddo
 var2d!0 = "Vertical"
 var2d!1 = "Horizontal"
 endif
 
 
 
 st_x = tointeger(xy(0,0)) + 1
 st_y = tointeger(xy(0,1)) + 1
 ed_x = tointeger(xy(xp(0)-1,0)) + 1
 ed_y = tointeger(xy(xp(0)-1,1)) + 1
 if (opts) then
 var2d@Orientation = "Cross-Sesion: (" + \
 st_x + "," + st_y + ") to (" + \
 ed_x + "," + ed_y + ")"
 else
 var2d@Orientation = "Cross-Sesion: (" + \
 st_x + "," + st_y + ") to (" + \
 ed_x + "," + ed_y + ") ; center=(" + \
 loc_param(0) + "," + loc_param(1) + \
 ") ; angle=" + angle
 endif
 
 return(var2d)
 endif
 
 
 
 end