I am trying to create a panel plot using the script below. Within the
loop, the script reads in two files, a wrfout file and a gridded PRISM
file, and calculates the difference, xdiff. I am saving the difference in
a variable I created using new called PrWrfdiff. However, I'm getting an
error when I try to save xdiff in PrWrfdiff that says "Subscript out of
range, error in subscript #0" (line 95). The dimensions of xdiff are 192
x 132 and the dimensions of PrWrfdiff are (6,192,132), six time steps one
for each month in the six-month simulation.
Thanks,
Liz.
Script:
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/esmf/ESMF_regridding.ncl"
;======================================================================
; The main code
;======================================================================
begin
;---Specify interpolation to be used
method = "bilinear"
wgtFileDir = "/glade/p/work/burakows/plot/NCL/ESMF_regridding/"
wgtFileName = "PRISM_4kmM2_to_WRF."+method+"_wgt.nc"
wgtFilePath = wgtFileDir+wgtFileName
print("wgtFilePath="+wgtFilePath)
;---PRISM files [input]
varp = "tmax" ; PRISM name
dirp = "/glade/p/work/burakows/prism/data/"
filp = systemfunc("cd "+dirp+" ; ls PRISM_"+varp+"*bil.nc")
nfilp = dimsizes(filp)
print(filp(0:nfilp-1)) ; 1st 5 file
;---WRF directory
varw = "T2" ; WRF name
dirw =
"/glade/p/work/burakows/plot/albedo_CLMV351current_200811-200904/"
varm = "T2max" ; WRF monthly name
prtFlg = True
;---Create new ncl variable to hold panel plot contents
plot_diff = new(6,graphic) ; from p.80 of NCL manual, right
before looping (panel1e.ncl)
WrfPrdiff = new((/6,192,132/),float,1e20) ; create new variable
WrfPr(time,y,x)
print(dimsizes(WrfPrdiff))
;----------------------------------------------------------------------
; Loop through six PRISM and WRF monthly files to calculate difference
;----------------------------------------------------------------------
do nfp=0,nfilp-1
yyyymm = toint( str_get_field(filp(nfp),5,"_.") ) ; parse
PRISM name for time
yyyy = yyyymm/100 ; current
year as integer
mm = yyyymm-(yyyy*100) ; current month
fp = addfile(dirp+filp(nfp), "r")
XP = fp->z ; original
PRISM variable (all 'z')
xp = ESMF_regrid_with_weights(XP,wgtFilePath,False); create new
variable on WRF grid
printVarSummary(XP)
printVarSummary(xp)
; matching WRF file
filw = "avg_"+varm+"_"+yyyy+"-"+sprinti("%0.2i",mm)+".nc"
if (isfilepresent(dirw+filw)) then
fw = addfile(dirw+filw, "r")
xw = fw->$varm$ ; WRF
monthly variable
;xw@lat2d = fw->XLAT(0,:,:)
;xw@lon2d = fw->XLONG(0,:,:)
printVarSummary(xw)
;----------------------------------------------------------------------
; Subtract PRISM from WRF
;----------------------------------------------------------------------
fwo = addfile(dirw+"wrfout_d03_2008-11-01_00:00:00.nc","r"); Open any
wrfout daily file to obtain
; lat and lon
xwo = fwo->T2 ; read T2 and preserve metadata
xdiff = xw(0,:,:)-xp ; Avg monthly WRF minus Prism for month i
print(dimsizes(xdiff))
xdiff@long_name = "difference: "+yyyymm ; Assign long name
to difference
xdiff@units = "Degrees C" ; Assign units to difference
printVarSummary(xdiff)
xdiff@lat2d = fwo->XLAT(0,:,:) ; Assign lat to xdiff
xdiff@lon2d = fwo->XLONG(0,:,:) ; Assign lon to xdiff
;print("xdiff: min="+min(xdiff(nfp,:,:))+" max="+max(xdiff(nfp,:,:)))
minlat = min((/min(xdiff@lat2d),min(fp&y)/))
maxlat = max((/max(xdiff@lat2d),max(fp&y)/))
minlon = min((/min(xdiff@lon2d),min(fp&x)/))
maxlon = min((/max(xdiff@lon2d),max(fp&x)/))
bndadd = 0.10 ; spacing around the plot edges
WrfPrdiff(nfp,:,:)=xdiff ; assign xdiff to WrfPr variable
plot_diff(nfp) = gsn_csm_contour_map(wks,WrfPrdiff(nfp),res)
; create plot of month nfp in plot_diff
if (prtFlg) then ; first time only
printVarSummary(XP)
printVarSummary(xp)
printVarSummary(xw)
prtFlg = False
end if ; nfp
end if ; isfilepresent
end do
printVarSummary(WrfPrdiff)
;----------------------------------------------------------------------
; Plotting options section
;----------------------------------------------------------------------
pltType = "ps" ; plot type
pltDir = "./" ; plot directory
pltName = "Diff_WRF_PRISM_2009" ; plot name (ps file)
pltPath = pltDir+pltName ; plot path
wks = gsn_open_wks(pltType,pltPath) ; create
workstation for ps file
gsn_define_colormap(wks,"temp_diff_18lev") ; define color table
res = True
res@gsnDraw = False
res@gsnFrame = False
res@cnFillOn = True ; color plot desired
res@cnLinesOn = False ; turn off contour lines
res@cnLineLabelsOn = False ; turn off contour labels
res@cnFillMode = "RasterFill" ; turn raster on
res@lbLabelBarOn = False ; Will turn on in panel
later
res@lbOrientation = "Horizontal" ; Horizontal label bar
res@mpFillOn = False
res@mpOutlineOn = True
res@mpOutlineBoundarySets = "AllBoundaries"
res@mpProjection = "CylindricalEquidistant"
res@mpLimitMode = "LatLon" ; required
res@mpMinLatF = minlat-bndadd
res@mpMaxLatF = maxlat+bndadd
res@mpMinLonF = minlon-bndadd
res@mpMaxLonF = maxlon+bndadd
res@mpCenterLonF = (minlon + maxlon)*0.5
res@mpCenterLatF = (minlat + maxlat)*0.5
res@gsnMajorLatSpacing = 2
;res@gsnMajorLonSpacing = 2
res@gsnMinorLonSpacing = 2
res@cnLinesOn = False
res@cnLineLabelsOn = False
res@cnFillOn = True
res@cnFillMode = "RasterFill"
res@gsnAddCyclic = False ; regional grid (changes central meridian)
;----------------------------------------------------------------------
; Panel Plotting options section
;----------------------------------------------------------------------
pres = True
pres@gsnMaximize = True
pres@gsnPanelLabelBar = True
pres@gsnCenterString = "WRF minus PRISM" + yyyymm
pres@txString = "Monthly averaged WRF T2 minus PRISM"
gsn_panel(wks,plot_diff,(/2,3/),pres)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Feb 11 11:18:07 2014
This archive was generated by hypermail 2.1.8 : Wed Feb 12 2014 - 09:35:52 MST