[ncl-talk] Divide by zero error using "pot_vort_isobaric" function

[Daniel Swain]

Hello,

I'm currently attempting to use the built-in NCL function to calculate
potential vorticity ("pot_vort_isobaric"). Unfortunately, I am repeatedly
getting a divide by zero error:

"
fatal:divide: Division by 0, Can't continue
fatal:Div: operator failed, can't continue
fatal:["Execute.c":8575]:Execute: Error occurred at or near line 3892 in
file $NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl
fatal:["Execute.c":8575]:Execute: Error occurred at or near line 3909 in
file $NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl
"

I have confirmed that my original units are in hPa (thus the *100 hPa-> Pa
conversion should be correct), and that the original latitudes were
north->south (so the lat flip should also be correct). I'm assuming this
occurring during the static stability calculation (as mentioned in previous
ncl-talk discussions), but I'm not sure whether there has yet been a
resolution.

I've attempted this with data from different years, but eventually a zero
division error always crops up. A short sample script is below, and
representative data files have been uploaded to the FTP site with the same
names as below.

I'm currently using NCL version 6.3.0.

Thanks very much for any advice!

~Daniel

"
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

infile1 = addfile ("air.2013.nc", "r")
infile2 = addfile ("uwnd.2013.nc", "r")
infile3 = addfile ("vwnd.2013.nc", "r")

     U   = infile2->uwnd   ; (time,lev,lat,lon)
     V   = infile3->vwnd
     T   = infile1->air    ; K

     lev = T&level  ; hPa
     lev = lev*100 ; convert units
     lev at units = "Pa"

     U   = U(:,:,::-1,:)    ; reorder to South -> North
     V   = V(:,:,::-1,:)
     T   = T(:,:,::-1,:)

     lat = T&lat
     gridType  = 0          ; Gaussian grid
     opt = 0
     pvort = pot_vort_isobaric(lev,U,V,T,lat, gridType, opt)

outfile1=addfile("pvort.2013.nc","c")
outfile1->pvort = pvort

end

"
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160309/c5500b2b/attachment.html