[ncl-talk] Inconsistency between wgt_arearmse and wgt_volrmse... possible bug in wgt_volrmse?

[Ben Green - NOAA Affiliate]

Hello,

I believe that the NCL function wgt_volrmse has a bug -- namely, that the
formula is wrong. I am using NCL V6.3.0, but downloaded the source code of
V6.4.0 and confirmed that the relevant source code is unchanged. I also
could not find any reports on a bug in wgt_volrmse on the NCL website.

I first became suspicious of a bug when my code using wgt_volrmse were an
order of magnitude different than results that a colleague of mine obtained
using raw Fortran code.
So I downloaded the source code for versions 6.3.0 and 6.4.0 to look at
wgt_volrmse. The relevant fortran file is under
./ncl_ncarg-6.4.0/ni/src/lib/nfpfort/volRmse.f

I looked at this code and it appears that a square is being taken twice:
First, SUMTQ = (T - Q)**2
Then, VARTQ = (SUMTQ/SUMWZ)**2
A square root is taken at the end. So that "negates" only one of the **2
operations.


I then looked at the code for wgt_arearmse. In theory, these two functions
should provide identical answers if the third-from-rightmost dimension is
of size 1 and no weighting is applied to this dimension in wgt_volrmse.
The relevant source code for wgt_arearmse is under
./ncl_ncarg-6.4.0/ni/src/lib/nfpfort/areaRmse.f

This code takes a square once (SUMD) then takes a square root at the end.
So I have no reason to believe that wgt_arearmse has a bug.


As I said earlier, in theory wgt_volrmse and wgt_arearmse should provide
identical answers if the third-from-rightmost dimension is of size 1 and no
weighting is applied to this dimension in wgt_volrmse.
I tested this out with a script, and with data samples. I upload a .tar
file (demo_RMSEbug.tar, *filesize 530K*) via ftp. The contents of this .tar
file are:
fcst.nc     (a file of size (1, 181 lat, 360 lon) of forecast 2-m
temperatures): *filesize 260K*
obs.nc     (a file of size (1, 181 lat, 360 lon) of observed 2-m
temperatures): *filesize 261K*
demo_bug.ncl    (a small sample script): *filesize 1.5K*

In my test script, I calculate RMSE using both wgt_arearmse and
wgt_volrmse, and print the output. As might be expected with the extra **2
operation in wgt_volrmse, the RMSE I get with this function is much higher
than the RMSE I get when using wgt_arearmse. Values from both are printed
out within the code.

I believe that this conclusively demonstrates a bug in wgt_volrmse (or at
least an inconsistency with wgt_arearmse). If not, I apologize.

Other relevant information:
ncl -V returns 6.3.0
uname -a returns   Linux tfe10 3.10.0-693.11.1.el7.x86_64 #1 SMP Fri Oct 27
05:39:05 EDT 2017 x86_64 x86_64 x86_64 GNU/Linux


Please let me know if you have any further questions.

Thank you very much for all of the work you do for NCL -- it is a great
tool.

Sincerely,

Benjamin Green

-- 
Benjamin Green
CIRES Research Associate at NOAA/ESRL/GSD
David Skaggs Research Center, Room 2B-501
Phone: 303-497-5132
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180110/47e980ed/attachment-0001.html>