[ncl-talk] Maybe I found a bug of the function "dpres_plevel"

[Bilbo Qiu]

Dear everyone,

It's my first time to send the email. I am trying to describe the problem
briefly as possible.

Version:  NCL 6.3.0

System: Ubuntu

No errors.

Warning: dpres_plevel: At one or more grid points the sum of the layer
thicknesses is not equal to (psfc-ptop). Are units of plev, psfc and ptop
matching?

Problems: I tried to calculate the pressure layer thickness "dp", using the
function "dpres_plevel". The input file is "ps_qv.nc", which includes two
variables, ps (surface pressure) and qv(specific humidity). I exported the
result to a file named "dp.nc". A warning displayed as the above. I checked
all of the units. It is "Pa". I visualized "dp.nc", using ncview and found
dp is about 3000 Pa on Tibet Plateau on the level of 1000 hPa where the
surface pressure is about 550 hPa. There should be no value on Tibet
Plateau where the surface pressure is much below 1000 hPa. In addition,
there exists extreme value about 30000 pa on the top level of 700 hPa on
the boundary of Tibet Plateau. The value is too big for dp. Below is my
script:
--------------------------------
begin
   diri = "./"
   fili = "ps_qv.nc"
   fi   = addfile(diri+fili, "r")

   q  = fi->qv(0,0:12,:,:) ; (levels, latitude, longitude), it will provide
coordinate information for dp

   psfc  = fi->ps(0,:,:)   ; (lat,lon), the unit of the variable "ps" in
the input file is Pa
   psfc at units = "Pa"
   printVarSummary(psfc)
   printMinMax(psfc, 0)

   ;isobaric levels
    plev = fi->levels(0:12) ;from 1000 to 700 hPa
  plev = plev*100. ;the unit of plev is hPa, convert it to Pa by *100
    plev at units = "Pa"
    printVarSummary(plev)
    printMinMax(plev, 0)

    ptop= 70000
  ptop at units = "Pa"

    dp  = dpres_plevel(plev, psfc, ptop, 0) ; the key function
  copy_VarCoords(q, dp)
    printVarSummary(dp)
    printMinMax(dp, 0)
;output dp to a netCDF file
system("/bin/rm -f dp.nc")   ; remove any pre-existing file
ncdf = addfile("dp.nc" ,"c")  ; open output netCDF file
fAtt               = True            ; assign file attributes
fAtt at title         = "NCL Simple Approach to netCDF Creation"
fAtt at source_file   =  "original-file.nc"
fAtt at Conventions   = "None"
fAtt at creation_date = systemfunc ("date")
fileattdef( ncdf, fAtt )            ; copy file attributes
filedimdef(ncdf,"time",-1,True)
ncdf->dp  = dp
end
-----------------------------------

I couldn't find the source code of the function "dpres_pleve" and had no
idea where the errors are.

I am looking forwards to your help.

Additions: I set the surface pressure to a constant, 101800 Pa. The result
makes sense. I am wondering if the function is wrong when the surface
pressure is a 2-D or 3-D array.

Bilbo Qiu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151014/567677ab/attachment.html