pb2nc results not consistently reproducible

[gsketefian]

Hi,

While running METplus's PB2NC tool from within the SRW App, I found that I get varying results between different runs of PB2NC with identical inputs. With @michelleharrold's help, I was able to reproduce the issue using a direct call to pb2nc without going through the SRW or METplus.

I've attached a shell script (test_pb2nc.sh.txt; I added the .txt extension is so I can upload it) that can be run on Hera that reproduces the problem. The script calls pb2nc 10 times and saves the output of each call in a subdirectory (call001, call002, etc). Each subdirectory contains the NetCDF output from pb2nc as well as a log file containing the output of the pb2nc command to stdout. After each call is complete, the script prints out to screen the line from the log file that shows the number of "retained or derived observations".

Below is the output from a sample run of the script. You can see that the number of retained or derived observations takes on one of two numbers (apparently randomly): 127785 and 127787. You'll see further differences if you compare the log file from a 127785 run to one from a 127787 run. Specifically, the 127787 runs have lines like

DEBUG 4: Not saved the computed PBL (30859 from 31090) because of the MAX PBL 10000  (20190615_000000 ADPUPA 74455)

where the computed PBL height of around 30000 (feet or meters, not sure) is at least remotely physical, while the 127785 runs have lines like

DEBUG 4: Not saved because the computed PBL (-3.381817831e+19) is less than the station elevation (231). 20190615_000000 ADPUPA 74455

where the computed PBL height of -3e+19 is completely unphysical (maybe because the code is accessing a random memory location).

Note that the netcdf files from the 127785 runs are all identical, and the ones from the 127787 runs are identical (I verified this using nccmp). Also, the number of "retained or derived observations" does not take on any value other than 127785 and 127787.

This seems like a bug to me. Could you please take a look and confirm, or let me know why it isn't a bug (maybe my modules are inconsistent)?

Thanks.

Sample output from test_pb2nc.sh on Hera:

$ ./test_pb2nc.sh 

Inactive Modules:
  1) ncview

The following have been reloaded with a version change:
  1) intel/2023.2.0 => intel/2022.1.2


Currently Loaded Modules:
  1) hpss/hpss    3) ncl/6.6.2      5) nco/5.1.6     7) imagemagick/7.1.1-11   9) rocoto/1.3.7    11) anaconda/2.3.0
  2) gnu/13.2.0   4) netcdf/4.9.2   6) nccmp/1.9.1   8) xxdiff/5.1            10) intel/2022.1.2  12) met/11.1.0

Inactive Modules:
  1) ncview




Calling pb2nc for 001-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call001/pb2nc_call001.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10454:DEBUG 2: Observations retained or derived	= 127787

Calling pb2nc for 002-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call002/pb2nc_call002.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10456:DEBUG 2: Observations retained or derived	= 127785

Calling pb2nc for 003-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call003/pb2nc_call003.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10456:DEBUG 2: Observations retained or derived	= 127785

Calling pb2nc for 004-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call004/pb2nc_call004.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10456:DEBUG 2: Observations retained or derived	= 127785

Calling pb2nc for 005-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call005/pb2nc_call005.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10454:DEBUG 2: Observations retained or derived	= 127787

Calling pb2nc for 006-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call006/pb2nc_call006.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10456:DEBUG 2: Observations retained or derived	= 127785

Calling pb2nc for 007-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call007/pb2nc_call007.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10454:DEBUG 2: Observations retained or derived	= 127787

Calling pb2nc for 008-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call008/pb2nc_call008.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10454:DEBUG 2: Observations retained or derived	= 127787

Calling pb2nc for 009-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call009/pb2nc_call009.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10454:DEBUG 2: Observations retained or derived	= 127787

Calling pb2nc for 010-th time...
Log file is at:  /scratch2/BMC/det/Gerard.Ketefian/UFS_CAM/tmp3/test_PB2NC/call010/pb2nc_call010.log
Done.
Line from log file indicating number of observations retained or derived (obtained via grep; includes line number):
10454:DEBUG 2: Observations retained or derived	= 127787

test_pb2nc.sh.txt

[JohnHalleyGotway]

@gsketefian thanks for pointing out this issue and providing test data to demonstrate it.

I strongly suspect that it's related to the derivation of PBL height in this FORTRAN code ported over from NOAA/EMC:
https://github.com/dtcenter/MET/blob/develop/src/tools/other/pb2nc/calpbl.f

My first guess is that it's an issue with precision, perhaps in this comparison.

In the MET C++ code, we use an is_eq() function to test for equality within a certain tolerance:

int is_eq(double a, double b, double tol)

And perhaps something similar is needed here?

I'll try to replicate the issue locally and test to see if that type of solution resolves it.

[gsketefian]

@JohnHalleyGotway Thanks for your quick response. Yes, it seems reasonable that FORTRAN code is involved, but it's not clear to me what is changing from one run to the next to make the output change. Could that code be accessing an uninitialized variable that sometimes happens to have a reasonable value and sometimes not? From my FORTRAN experience (admittedly many years ago now), this kind of random change in output usually happens when there is a memory bug, e.g. accessing uninitialized variables or array elements. Might be worth recompiling that code in debug mode if that's easy to do.

Thanks for your help.

[JohnHalleyGotway]

@gsketefian, @hsoh-u and I just finished testing this out and summarized it in this issue:
dtcenter/MET#2902

The good news is that porting over the fix from develop into main_v11.1 is a one-liner. I'll do that now and include the fix in the MET-11.1.1 bugfix release.

[JohnHalleyGotway]

FYI, here's the corresponding dtcenter/MET#2903 pull request to fix this in the main_v11.1 branch for inclusion in the MET-11.1.1 bugfix release.