GSI DA app errors with ncio/1.1.1 create_dataset

[RussTreadon-NOAA]

NOAA-EMC/GSI DA apps built using ncio/1.1.1 do not correctly copy all data from the input netcdf file to the output file via the create_dataset call with copy_vardata=.true.

grid_xt, grid_yt, lat, lon, pfull, and ``phalf are incorrect in the output file opened by create_dataset. For example, while the input file has

 grid_xt = 0, 1.875, 3.75, 5.625, 7.5, 9.375, 11.25, 13.125, 15, 16.875,
    18.75, 20.625, 22.5, 24.375, 26.25, 28.125, 30, 31.875, 33.75, 35.625,
    37.5, 39.375, 41.25, 43.125, 45, 46.875, 48.75, 50.625, 52.5, 54.375,
    56.25, 58.125, 60, 61.875, 63.75, 65.625, 67.5, 69.375, 71.25, 73.125,
    75, 76.875, 78.75, 80.625, 82.5, 84.375, 86.25, 88.125, 90, 91.875,
    93.75, 95.625, 97.5, 99.375, 101.25, 103.125, 105, 106.875, 108.75, ...

the output file created by create_dataset contains

 grid_xt = 0, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _ ;

This comment applies to both atmfXXX.nc and sfcfXXX.nc files created by the current operational gfs.v16.1.8 and atm and sfc nc files created by ufs_model tag Prototype-P8c

This issue is opened to report this issue and track its resolution.

[RussTreadon-NOAA]

Tagging @jswhit (ncio developer), @MichaelLueken-NOAA (GSI app CM), @WalterKolczynski-NOAA (g-w CM) , and @kgerheiser (nceplibs CM) for awareness.

Resolution of this issue impacts

• NCEPLIBS-ncio issues #58, #60
• hpc-stack issue #462
• g-w issue #711
• GSI issue #348

and perhaps other open issues.

[edwardhartnett]

Do we have a sample input file we can use for testing?

[RussTreadon-NOAA]

It seems we can demonstrate the behavior with any atmfXXX.nc or sfcfXXX.nc. I took source code for GSI DA app getsfcensmean.x and trimmed down to simply read/write a sfcfXXX.nc file. I did not create this as a self-contained test. I added test.x as an app built by the GSI cmake.

The test source code is on Orion in /work/noaa/da/Russ.Treadon/git/global_workflow/develop/sorc/gsi.fd/util/EnKF/gfs/src/test.fd

test.x was built by running /work/noaa/da/Russ.Treadon/git/global_workflow/develop/sorc/gsi.fd/ush/build.sh This builds all DA apps (overkill but easiest for me).

test.x is in the build directory, /work/noaa/da/Russ.Treadon/git/global_workflow/develop/sorc/gsi.fd/build/util/EnKF/gfs/src/test.fd. It was tested using gdas.t18.sfcf006.nc as input. gdas.t18z.sfcf006_output.nc is the output file. grid_xt and other coordinate fields are incorrect in the output file.

[edwardhartnett]

Can you submit a PR with these files added to the test directory?

[MichaelLueken]

@RussTreadon-NOAA @jswhit The issue is with write_vardata_code.f90. Similar to what was happening in read_vardata_code_1d.f90, when n==nd (which, for one dimensional arrays, is always the case), count is hardwired to 1. When I replaced:

     if (n == nd) then
        start(n) = ncount
        count(n) = 1
     else
        start(n) = 1
        count(n) = dset%variables(nvar)%dimlens(n)
        dimlens(ndim) = dset%variables(nvar)%dimlens(n)
        ndim = ndim + 1
     end if

with:

     start(n) = ncount
     count(n) = dset%variables(nvar)%dimlens(n)
     dimlens(ndim) = dset%variables(nvar)%dimlens(n)
     ndim = ndim + 1

The full analysis for global_fv3_4denvar_C192 gave:

grid_xt = 0, 0.46875, 0.9375, 1.40625, 1.875, 2.34375, 2.8125, 3.28125,
    3.75, 4.21875, 4.6875, 5.15625, 5.625, 6.09375, 6.5625, 7.03125, 7.5,
    7.96875, 8.4375, 8.90625, 9.375, 9.84375, 10.3125, 10.78125, 11.25,
    11.71875, 12.1875, 12.65625, 13.125, 13.59375, 14.0625, 14.53125, 15,
    15.46875, 15.9375, 16.40625, 16.875, 17.34375, 17.8125, 18.28125, 18.75,
    19.21875, 19.6875, 20.15625, 20.625, 21.09375, 21.5625, 22.03125, 22.5,
    22.96875, 23.4375, 23.90625, 24.375, 24.84375, 25.3125, 25.78125, 26.25, ...

rather than:

grid_xt = 0, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _,
    _, _ ;

@RussTreadon-NOAA Would you mind running a test using:

/gpfs/dell2/emc/modeling/noscrub/Michael.Lueken/gsi_common.lua

and

/gpfs/dell2/emc/modeling/noscrub/Michael.Lueken/gsi_wcoss_dell_p3.lua

[RussTreadon-NOAA]

I can open a PR but I do not have a branch to merge.

I opened this issue on behalf of others. Developers can not run cycled global parallels from the head of the NOAA-EMC/GSI and global-workflow develop branches due to this issue.

[MichaelLueken]

@edwardhartnett The file, dynf000_template.nc.in, that currently resides in NCEPLIBS-ncio/tests, contains grid_xt:

ncdump -v grid_xt dynf000_template.nc.in

 grid_xt = 0, 1.40625, 2.8125, 4.21875, 5.625, 7.03125, 8.4375, 9.84375,
    11.25, 12.65625, 14.0625, 15.46875, 16.875, 18.28125, 19.6875, 21.09375,
    22.5, 23.90625, 25.3125, 26.71875, 28.125, 29.53125, 30.9375, 32.34375,
    33.75, 35.15625, 36.5625, 37.96875, 39.375, 40.78125, 42.1875, 43.59375,
    45, 46.40625, 47.8125, 49.21875, 50.625, 52.03125, 53.4375, 54.84375,
    56.25, 57.65625, 59.0625, 60.46875, 61.875, 63.28125, 64.6875, 66.09375,
    67.5, 68.90625, 70.3125, 71.71875, 73.125, 74.53125, 75.9375, 77.34375,
    78.75, 80.15625, 81.5625, 82.96875, 84.375, 85.78125, 87.1875, 88.59375,
    90, 91.40625, 92.8125, 94.21875, 95.625, 97.03125, 98.4375, 99.84375,
    101.25, 102.65625, 104.0625, 105.46875, 106.875, 108.28125, 109.6875,
    111.09375, 112.5, 113.90625, 115.3125, 116.71875, 118.125, 119.53125,
    120.9375, 122.34375, 123.75, 125.15625, 126.5625, 127.96875, 129.375,
    130.78125, 132.1875, 133.59375, 135, 136.40625, 137.8125, 139.21875,
    140.625, 142.03125, 143.4375, 144.84375, 146.25, 147.65625, 149.0625,
    150.46875, 151.875, 153.28125, 154.6875, 156.09375, 157.5, 158.90625,
    160.3125, 161.71875, 163.125, 164.53125, 165.9375, 167.34375, 168.75,
    170.15625, 171.5625, 172.96875, 174.375, 175.78125, 177.1875, 178.59375,
    180, 181.40625, 182.8125, 184.21875, 185.625, 187.03125, 188.4375,
    189.84375, 191.25, 192.65625, 194.0625, 195.46875, 196.875, 198.28125,
    199.6875, 201.09375, 202.5, 203.90625, 205.3125, 206.71875, 208.125,
    209.53125, 210.9375, 212.34375, 213.75, 215.15625, 216.5625, 217.96875,
    219.375, 220.78125, 222.1875, 223.59375, 225, 226.40625, 227.8125,
    229.21875, 230.625, 232.03125, 233.4375, 234.84375, 236.25, 237.65625,
    239.0625, 240.46875, 241.875, 243.28125, 244.6875, 246.09375, 247.5,
    248.90625, 250.3125, 251.71875, 253.125, 254.53125, 255.9375, 257.34375,
    258.75, 260.15625, 261.5625, 262.96875, 264.375, 265.78125, 267.1875,
    268.59375, 270, 271.40625, 272.8125, 274.21875, 275.625, 277.03125,
    278.4375, 279.84375, 281.25, 282.65625, 284.0625, 285.46875, 286.875,
    288.28125, 289.6875, 291.09375, 292.5, 293.90625, 295.3125, 296.71875,
    298.125, 299.53125, 300.9375, 302.34375, 303.75, 305.15625, 306.5625,
    307.96875, 309.375, 310.78125, 312.1875, 313.59375, 315, 316.40625,
    317.8125, 319.21875, 320.625, 322.03125, 323.4375, 324.84375, 326.25,
    327.65625, 329.0625, 330.46875, 331.875, 333.28125, 334.6875, 336.09375,
    337.5, 338.90625, 340.3125, 341.71875, 343.125, 344.53125, 345.9375,
    347.34375, 348.75, 350.15625, 351.5625, 352.96875, 354.375, 355.78125,
    357.1875, 358.59375 ;

The issue appears to be one that the use of write_vardata isn't strenuously tested (at least, not for 1d arrays).

[RussTreadon-NOAA]

@MichaelLueken-NOAA , thank you for looking into this and developing a fix.

I used your lua files to rebuild gsi.x on Mars. The global_fv3_4denvar_C192 regression test was rerun with the test configured to write out the analysis not the increment. I find the same result as you. The coordinate fields in the resulting siganl are correct. For example, grid_yt contains

 grid_yt = 89.6416480725934, 89.1774327980165, 88.7104760813539,
    88.2428999560018, 87.7750888743169, 87.3071640967019, 86.8391757792909,
    86.3711483847176, 85.9030952554449, 85.435024284034, 84.966940436571,
    84.498846993417, 84.0307462082699, 83.5626396805404, 83.0945285766584,
    82.626413767268, 82.1582959153746, 81.6901755347503, 81.2220530296775, ...

'pfull` contains

 pfull = 0.01278146, 0.02033404, 0.03177342, 0.04878282, 0.07361853,
    0.1092587, 0.1595392, 0.2292877, 0.3244748, 0.4523215, 0.621393,
    0.8416426, 1.124391, 1.482229, 1.928879, 2.478976, 3.147755, 3.950706,
    4.903192, 6.020019, 7.315024, 8.800693, 10.48782, 12.38528, 14.49982,
    16.83605, 19.39651, 22.18178, 25.1909, 28.42169, 31.87125, 35.53667,
    39.41547, 43.5065, 47.81049, 52.33096, 57.07489, 62.0536, 67.28372, ...

I also reran an eobs job on Mars using the recompiled GSI DA apps. getsfcensmean.x and getsigensmeanp_smooth.x create ensemble mean files with correct values in coordinate variables. I followed this by using the resulting ensemble mean files in eobs. The gsi.x in eobs successfully ran to completion. No errors occurred when reading the ensemble mean atm and sfc nc files.

Your fix is working!

[jswhit]

I've created a PR (PR #63) to fix the problem in write_vardata_code.f90 (thanks @MichaelLueken-NOAA). I've also added a test for reading the grid_xt variable in tst_ncio.F90.

One thing that needs a little more thought is how unlimited dimensions are treated in write_vardata_code.f90