Stage IV precip data: conversion error

Hi Everybody-
I am using NCEP Stage 4 data. I used ncl's .Grb to .nc converter.

I can view the data just fine using ncview, and it appears that I can print
the array values correctly in NCL, BUT I cannot plot the data. There is no
error message so I am struggling to understand what is wrong in order to fix
it.

Perhaps an NCL developer can take a look at the data and assist? I uploaded
the file via ftp: ST2ml2002092008.Grb

Any help is greatly appreciated.

Below is the code.

;************************************************
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/wrf/WRF_contributed.ncl"
begin
;************************************************
; open file and read in data
;************************************************

mon = (/"04","05","06","07","08","09"/)
hr =
(/"00","01","02","03","04","05","06","07","08","09","10","11","12","13","14"
,"15","16","17","18","19","20","21","22","23"/)
day =
(/"01","02","03","04","05","06","07","08","09","10","11","12","13","14","15"
,"16","17","18","19","20","21","22","23","24","25","26","27","28","29","30",
"31"/)

suma = new((/881,1121/),"float")
suma = 0.
suma@_FillValue = 1e20
sumh = new((/24,881,1121/),"float")
sumh = 0.
sumh@_FillValue = 1e20

  wks = gsn_open_wks("ps" ,"test") ; ps,pdf,x11,ncgm,eps

y1 = 2002
m1 = 5
if(m1 .eq. 5)then
nt = 30*24
end if

do yr=y1,y1
do month=m1,m1
if(month .eq. 0)then
nd = 29
end if
if(month .eq. 1)then
nd = 30
end if
if(month .eq. 2)then
nd = 29
end if
if(month .eq. 3)then
nd = 30
end if
if(month .eq. 4)then
nd = 30
end if
if(month .eq. 5)then
nd = 29
end if
nt = -1
do d=19,19
do h=8,8
nt = nt + 1

files = "ST2ml"+yr+mon(month)+day(d)+hr(h)+".nc"
  g = addfile(files, "r")
print(files+" ")
if(nt .eq. 0)then
xlat = g->gridlat_240
xlon = g->gridlon_240
end if

;need to do qc?

var = g->A_PCP_240_SFC_acc1h
;print(var)
;print(files+" "+max(var)+" "+min(var))
;printVarSummary(var)
suma = suma + var(:,:)
;print(suma)

sumh(h,:,:) = sumh(h,:,:) + var(:,:)

end do
end do
end do
end do

print("max: "+max(suma))
print("min: "+min(suma))

do h=0,23
print(h+" "+"max: "+max(sumh(h,:,:)))
print(h+" "+"min: "+min(sumh(h,:,:)))
end do

;************************************************
; create plots
;************************************************
  gsn_define_colormap(wks ,"BlAqGrYeOrReVi200"); choose colormap

  res = True ; plot mods desired
  res@gsnMaximize = True ; uncomment to maximize size
  res@gsnSpreadColors = True ; use full range of colormap
  res@cnFillOn = True ; color plot desired
  res@cnLinesOn = False ; turn off contour lines
  res@cnLineLabelsOn = False ; turn off contour labels
  res@lbLabelAutoStride = True ; let NCL figure lb stride
  res@cnLevelSelectionMode = "ExplicitLevels" ; explicit [unequal] cn
levels
  res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
  res@cnLevels = (/0,.254,1,8,16,32,48,64,80,96,128/)

; WRF_map_c(gg,res,0) ; set map resources
;************************************************
; set True for native mapping (faster plotting)
; set to False othewise
;************************************************
   res@tfDoNDCOverlay = False

;************************************************
; associate the 2-dimensional coordinates to variables for plotting
; only if res@tfDoNDCOverlay=False
;************************************************
  if (.not.res@tfDoNDCOverlay) then
      lat2d = xlat
      lon2d = xlon
  end if
lat2d@units = "degrees_north"
lon2d@units = "degrees_east"

; lata = gg->XLAT(0,:,:)
; lona = gg->XLONG(0,:,:)
;x2 = 143
;y2 = 118
;x1 = 10
;y1 = 10
;print("wrf")
; res@mpRightCornerLonF = lona(y2,x2)
; res@mpRightCornerLatF = lata(y2,x2)
; res@mpLeftCornerLonF = lona(y1,x1)
; res@mpLeftCornerLatF = lata(y1,x1)

;lets try the stereographic suggestion from the archives
res@mpProjection = "Stereographic"
res@mpLimitMode = "Corners"
res@mpCenterLonF = 255.
res@mpCenterLatF = 90.
res@cnInfoLabelOn = False
res@mpLeftCornerLatF = lat2d(0,0)
res@mpLeftCornerLonF = lon2d(0,0)
res@mpRightCornerLatF = lat2d(880,1120)
res@mpRightCornerLonF = lon2d(880,1120)
res@gsnAddCyclic = False
;res@mpGridSpacingF = 300.
;res@mpGridLatSpacingF = 3.0
;res@mpGridLonSpacingF = 3.0

suma!0 = "lat"
suma!1 = "lon"
;suma&lat2d = lat2d
;suma&lon2d = lon2d
;sumh&lat2d = lat2d
;sumh&lon2d = lon2d

;map it to the wrfgrid perhaps? Wow that takes entirely too long.
;nlat = fspan(20.,55.,360)
;nlon = fspan(-120.,-60.,611)
;print("here")
;suma2 = rcm2rgrid(xlat,xlon,suma,nlat,nlon,1)
;print("here2")
;suma1 = rgrid2rcm(nlat,nlon,suma2,lata,lona,1)

plot = gsn_csm_contour_map(wks,suma,res)

;delete(suma2)
;delete(suma1)
;suma2 = rcm2rgrid(lat2d,lon2d,sumh,nlat,nlon,1)
;suma1 = rgrid2rcm(nlat,nlon,suma2,lata,lona,1)

;do h=0,23
;plot = gsn_csm_contour_map(wks,sumh(h,:,:),res)
;end do

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Jun 14 17:17:49 2010