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

[Dennis Shea]

Thank you for sending the file. There is a bug.

If you use ptop=0 or ptop=min(levels)=200, there is no bug. The [
psfc-ptop ] differences match the sum of the layer thickness.

However, if you use a ptop > min(levels) [for example, ptop=500] a
small subset [0.7%] of the grid points  do not sum correctly.

I am out tomorrow. I will look on Wednesday.

Regards
Dennis Shea




On Tue, Oct 13, 2015 at 11:56 PM, Bilbo Qiu <biloisbrave at gmail.com> wrote:
> 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
>
>
>
>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>