interpolation to PBL height

Hi guys,

I want to interpolate u, v, and w to the PBL height from my wrfoutput. Since PBL change/may change on every grid point, i therefore used wrf_interp_1d to interpolate every grid point value of u, v and w to each grid point value of PBL. I applied double do loop to accomplish this. Apparently it runs fine also however, this double do loop take almost one hour to process a single frame of my domain that comprises of 600 x 600 grid points. This makes wrf_interp_1d impractical for use.

I would greatly appreciate if somebody could suggest a better and faster way to process these data. My ncl code is copied below.

Best regards,
Basit.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

 a = addfile("/project/k14/basit/wrf-sims/wrf_168/wrfout_d02_2006-05-31_06_60min.nc","r")

; We generate plots, but what kind do we prefer?
  type = "x11"
  wks = gsn_open_wks(type,"plt_HeightLevel")
  gsn_define_colormap(wks,"ViBlGrWhYeOrRe") ;WhiteYellowOrangeRed")

; Set some basic resources
  res = True
  res@MainTitle = "REAL-TIME WRF"
  res@Footer = False

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 times = wrf_user_getvar(a,"times",-1) ; get all times in the file
 ntimes = dimsizes(times) ; number of times in the file
 print(ntimes)
 lat2d = a->XLAT(0,:,:)
 lon2d = a->XLONG(0,:,:)

  pltres = True
  mpres = True
  mpres@gsnMaximize = True
  mpres@mpOutlineOn = True
  mpres@mpFillOn = False
  mpres@mpGridLineThicknessF = 0.5
  mpres@mpGridSpacingF = 3
  mpres@gsnDraw = False ; don't draw the plots
  mpres@gsnFrame = False ; don't advance the frame
  mpres@tmYROn = False
  mpres@tmXTOn = False
  mpres@mpGeophysicalLineColor = "Black"
  mpres@mpLimbLineColor = "Blue"
  mpres@mpGeophysicalLineThicknessF = 1.5
  mpres@mpDataBaseVersion = "HighRes" ; Use the high-res database
  mpres@mpDataResolution = "FinestResolution"
  mpres@pmTickMarkDisplayMode = "Always" ; Turn on map tickmarks.
  mpres@mpNationalLineColor = "Black"
  mpres@mpNationalLineThicknessF = 2.5
  mpres@mpGridLineColor = "Gray49" ; Grid line color
  mpres@mpGridLineDashPattern = 14.
  mpres@mpGridLineThicknessF = 0.8
  mpres@mpGridLatSpacingF = 1.
  mpres@mpGridLonSpacingF = 1.
  mpres@mpOutlineBoundarySets = "National"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

 do it = 10, ntimes-1, 1 ; TIME LOOP

print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need

ua = wrf_user_getvar(a,"ua",it) ; u averaged to mass points
va = wrf_user_getvar(a,"va",it) ; v averaged to mass points
wa = wrf_user_getvar(a,"wa",it)
z = wrf_user_getvar(a, "z",it) ; grid point height
pbl= wrf_user_getvar(a,"PBLH",it)

ua@lat2d = lat2d
ua@lon2d = lon2d
va@lat2d = lat2d
va@lon2d = lon2d
wa@lat2d = lat2d
wa@lon2d = lon2d
z@lat2d = lat2d
z@lon2d = lon2d
pbl@lat2d = lat2d
pbl@lon2d = lon2d

; ==================================================================================

dims_zZ = dimsizes(z)
xn = dims_zZ(2)
yn = dims_zZ(1)

print(dims_zZ)

ua_plane = new((/dims_zZ(1), dims_zZ(2)/), float)
va_plane = new((/dims_zZ(1), dims_zZ(2)/), float)
wa_plane = new((/dims_zZ(1), dims_zZ(2)/), float)

pbl_ASL = z(0,:,:) + pbl ; pbl is above surface while z is above sea level (ASL).
printVarSummary(pbl_with_ter)

ua_plane@_FillValue = -999.99
va_plane@_FillValue = -999.99
wa_plane@_FillValue = -999.99

do xi =0, xn-1, 1
do yi = 0, yn-1, 1
ua_plane(yi,xi) = wrf_interp_1d(ua(:,yi,xi), z(:,yi,xi), pbl_ASL(yi,xi))
va_plane(yi,xi) = wrf_interp_1d(va(:,yi,xi), z(:,yi,xi), pbl_ASL(yi,xi))
wa_plane(yi,xi) = wrf_interp_1d(wa(:,yi,xi), z(:,yi,xi), pbl_ASL(yi,xi))
print(ua_plane(yi,xi)+"----"+z(0,yi,xi)+"---"+pbl_ASL(yi,xi))
end do
end do

printVarSummary(ua_plane)
wa_plane@description = "Vertical Velocity at PBL height"
wa_plane@units = "~N~m s~S~-1~N~"
ua_plane@units = "~N~m s~S~-1~N~"
va_plane@units = "~N~m s~S~-1~N~"

; ==================================================================================

; Plotting options for Wa
opts = res
opts@cnFillOn = True
opts@ContourParameters = (/-1., 1., 0.1/)
cnWa = wrf_contour(a,wks,wa_plane,opts)
delete(opts)

; -------------------------------------------

; Plotting options for Wind Vectors

optsVc = res
optsVc@FieldTitle = "Wind vector at PBL height"
optsVc@vcRefAnnoOn = True ; to turn on/off vector/wbarb legend
optsVc@gsnDraw = False ; Don't draw individual plot
        optsVc@gsnFrame = False ; Don't advance frame.
optsVc@lbLabelStride = 2 ; every other color
optsVc@gsnSpreadColors = True ; use full range of color map
optsVc@vcRefMagnitudeF = 10.0 ; define vector ref mag
optsVc@vcRefLengthF = 0.030 ; define length of vec ref
optsVc@vcGlyphStyle = "CurlyVector" ;"CurlyVector" ; turn on curley vectors
optsVc@vcMinDistanceF = 0.04 ; thin out vectors
optsVc@vcVectorDrawOrder= "PostDraw"

vector = wrf_vector(a,wks,ua_plane, va_plane, optsVc)
delete(optsVc)

; MAKE PLOTS
plot = wrf_map_overlays(a,wks,(/cnWa,vector/),pltres,mpres)

; end do ; END OF LEVEL LOOP

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

delete(va_plane)
delete(ua_plane)
delete(wa_plane)
delete(pbl_with_ter)

 end do ; END OF TIME LOOP

end

________________________________

This message and its contents including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Dec 24 04:01:16 2012