vertical interpolation

Hi,
   I am trying to interpolate OH at pressure levels to a different set of
pressure levels (from another netCDF file). I used the following command
line
 oh_out =
int2p_Wrap(Pin(time|:,lat|:,lon|:,level|:),OH(time|:,lat|:,lon|:,level|:),Pout(time|:,lat|:,lon|:,lev|:),
linlog)
 An error message " Number of dimensions in parameter (2) of (int2p_Wrap) is
(4), (1) dimensions were expected". In fact, this command line is very
similar to example 5
http://www.ncl.ucar.edu/Document/Functions/Contributed/int2p_Wrap.shtml.
  Any idea on what's wrong with the code? The whole code is copied below.
Thanks,
  Lei

~~~~~~~~~~~~~~~~~~~~~~~~~~
;*************************************************
; regrid_1.ncl
; this program does vertical interpolation

;
;************************************************
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;************************************************
begin
;************************************************
; read in netCDF file
;************************************************
  oh1 = addfile("newOH_hinterp.nc","r"); newOH_hinterp is output of
horizontal interpolation
  oh2 = addfile("oh_2004.nc","r")

; read in OH and other parameters
 OH = oh1->OH ; select variable to ave
 time=oh1->TIME
 lat=oh1->lat
 lon=oh1->lon
 A = oh1->A ; get a coefficiants
 B = oh1->B ; get b coefficiants
 level=oh1->LEVEL
 PSURF = oh1->PSURF ; get pressure
 LEVELP1 =oh1->LEVELP1

; read in pressure field from oh2
 PS=oh2->PS
 hyam=oh2->hyam
 hybm=oh2->hybm
 P0=oh2->P0
 lev=oh2->lev

   OH!3 = "lon"
   OH!2 = "lat"
   OH!1 = "level"
   OH!0 = "time"
   OH&lat=lat
   OH&lon=lon
   OH&level=level
   OH&time=time

 P_edge = new((/12,61,64,128/),float)
 Pin = new((/12,60,64,128/),float) ; P at middle point

   P_edge!3 = "lon"
   P_edge!2 = "lat"
   P_edge!1 = "levelp"
   P_edge!0 = "time"
   P_edge&lat=lat
   P_edge&lon=lon
   P_edge&levelp=LEVELP1
   P_edge&time=time

   Pin!3 = "lon"
   Pin!2 = "lat"
   Pin!1 = "level"
   Pin!0 = "time"
   Pin&lat=lat
   Pin&lon=lon
   Pin&level=level
   Pin&time=time

 do m=0,11
   do l=0,60
     P_edge(m,l,:,:) = A(l) + B(l)*PSURF(m,:,:)
   end do
 end do
 do l=0,59
   Pin(:,l,:,:) = (P_edge(:,l,:,:)+P_edge(:,l+1,:,:))/2.0 ; P at middle
point
 end do

   Pout = new((/12,28,64,128/),float)
   Pout!3 = "lon"
   Pout!2 = "lat"
   Pout!1 = "lev"
   Pout!0 = "time"
   Pout&lat=lat
   Pout&lon=lon
   Pout&lev=lev
   Pout&time=time
 do m=0,11
    do i=0,27
      Pout(m,i,:,:)=hyam(27-i)*P0+hybm(27-i)*PS(m,:,:)
    end do
 end do

oh_out = new((/12,28,64,128/),float)
linlog = 1

;do m=0,11
; do k=0, 63
; do j=0,127
; oh_out(m,:,k,j)
=int2p_Wrap(Pin(m,:,k,j),OH(m,:,k,j),Pout(m,:,k,j),linlog)
; end do
; end do
;end do

oh_out =
int2p_Wrap(Pin(time|:,lat|:,lon|:,level|:),OH(time|:,lat|:,lon|:,level|:),Pout(time|:,lat|:,lon|:,lev|:),
linlog)

 filo = "newOH_hint_vint.nc"
 system ("rm "+filo) ; remove any pre-existing file
 fo = addfile(filo , "c") ; open output file
 fo_at_title = "OH"
 fo->OH = oh_out; write OH to a file
 fo->PSURF =PSURF;
 fo->A =A
 fo->B =B
end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Sep 29 2009 - 13:28:15 MDT