Fix bug for ocean vector pairs

[KarinaAsmar-NOAA]

Description

This PR fixes a bug found in the ocnicepost tool affecting vector pairs. 2d vectors in the 'Cv' grid were coming out with zero values. The bug has been fixed and results confirmed by EMC. Please see note regarding bug fix.

Resolves #172 and #174
-->

Type of change

• Bug fix (fixes something broken)

Change characteristics

• Is this a breaking change (a change in existing functionality)? NO
• Does this change require a documentation update? NO

How has this been tested?

• clone build with bug fix in Ursa
• Ran ocnicepost executable on history netcdf file with SFS variables
• results in /scratch4/NCEPDEV/ovp/Karina.Asmar/ocnpost
• file with interpolated variables (SSV and tauy have been fixed) is ocean.1p00.nc
  -->

Checklist

• [x ] Any dependent changes have been merged and published
• [ x] My code follows the style guidelines of this project
• [ x] I have performed a self-review of my own code
• [ x] I have commented my code, particularly in hard-to-understand areas
• [ x] My changes generate no new warnings
• [ x] New and existing tests pass with my changes
• [ x] I have made corresponding changes to the documentation if necessary

[aerorahul]

Tag @DeniseWorthen for a review on this code as she is the author of this utility.

[DeniseWorthen]

This will take me some time to untangle.

[DeniseWorthen]

I don't understand the setup in /scratch4/NCEPDEV/ovp/Karina.Asmar/ocnpost. Your source file (ocean.nc) in that directory is 1/4 degree (1440x1080). Your srcdims in ocnicepost.nml is 360,320 (1-deg).

[KarinaAsmar-NOAA]

I don't understand the setup in /scratch4/NCEPDEV/ovp/Karina.Asmar/ocnpost. Your source file (ocean.nc) in that directory is 1/4 degree (1440x1080). Your srcdims in ocnicepost.nml is 360,320 (1-deg).

Thanks for noticing! I've run another test with the 1/4 deg source dimensions and 1deg output: /scratch4/NCEPDEV/ovp/Karina.Asmar/ocnpost_v2

[DeniseWorthen]

It appears this PR will revert the change to the 3d packing that was committed in #117. I think (?) I understand the git history to know what happened and why reversion of that fix is the right action.

The initial commit of this code contained a bug in the pack3d routine

gfs-utils/src/ocnicepost.fd/utils_mod.F90

Line 123 in 3457e53

do n = 1,dims(3)

Here, dims(3) is the number of vertical levels, which is 40 for the ocean model. This routine for the 3-d packing is only called by when doing ocean post (ice post is all 2-d).

The loop should have been over nflds. I believe this was a non-breaking bug because maxvars (the number of possible fields) was 50, so it didn't matter if you looped over nlevs because it was less than maxvars.

gfs-utils/src/ocnicepost.fd/init_mod.F90

Line 7 in 3457e53

real, parameter :: maxvars = 50 !< The maximum number of fields expected in a source file

As part of #117, the number of maxvars was increased

gfs-utils/src/ocnicepost.fd/init_mod.F90

Line 7 in 3240ee0

real, parameter :: maxvars = 70 !< The maximum number of fields expected in a source file

I'm guessing that now caused the original bug of looping (incorrectly) over 40 vertical levels to break the code. A modification was made to change the packing to resolve this, but for the wrong reason.

[KarinaAsmar-NOAA]

It appears this PR will revert the change to the 3d packing that was committed in #117. I think (?) I understand the git history to know what happened and why reversion of that fix is the right action.

The initial commit of this code contained a bug in the pack3d routine

gfs-utils/src/ocnicepost.fd/utils_mod.F90

Line 123 in 3457e53

do n = 1,dims(3)

Here, dims(3) is the number of vertical levels, which is 40 for the ocean model. This routine for the 3-d packing is only called by when doing ocean post (ice post is all 2-d).

The loop should have been over nflds. I believe this was a non-breaking bug because maxvars (the number of possible fields) was 50, so it didn't matter if you looped over nlevs because it was less than maxvars.

gfs-utils/src/ocnicepost.fd/init_mod.F90

Line 7 in 3457e53

real, parameter :: maxvars = 50 !< The maximum number of fields expected in a source file

As part of #117, the number of maxvars was increased

gfs-utils/src/ocnicepost.fd/init_mod.F90

Line 7 in 3240ee0

real, parameter :: maxvars = 70 !< The maximum number of fields expected in a source file

I'm guessing that now caused the original bug of looping (incorrectly) over 40 vertical levels to break the code. A modification was made to change the packing to resolve this, but for the wrong reason.

The bug fix from PR #117 was for the packarrays2d function. In the original line where the vecpairs were allocated, there was no corresponding deallocation. During initial testing, when there were 2 or more vector pairs of the same 2d remapping method (bilinear or conserve), this code crashed because the vecpair array was already allocated. By reverting back to the original code, the method works because we trimmed recently the ocean_sfs.csv file and now there is only one vector pair for each dimension and remapping method. In the future, if a csv file is given with more than one vector pair of the same 2d remapping method, this code could crash again. However, the current SFS netcdf ocean ouputs only produce the variables for the current ocean_sfs.csv file (one vector pair for each 2d remapping method), so I cannot test a more robust bug fix for that issue.

[DeniseWorthen]

@KarinaAsmar-NOAA Thanks. I realized I had a muddled interpretation of the fix #117. The problem is that apparently whatever triggered the error that prompted #117 has now been resolved because you have trimmed the variable lists.

I'm a little confused as to whether your original fix was required because of ocean fields or ice. In the ocean.csv that was committed with your #117, I don't see any variables in the ocean list that would have created more than 1 vector pair with the same mapping. I do see such variables in the ice.csv, where both the 2d ocean and atm velocites were being mapped bilinearly.

It should be possible to know how many vector pairs are going to be used for a given dimensionality and a given mapping type. The vectorpairs can then have an additional dimension to track each pair and loaded into the flds array.

Since this is clearly a bug which could arise in the future, I think an issue should be created that the fix PR you are making here will lead to failure when more than 1 vector pair for a given mapping type is requested.

[KarinaAsmar-NOAA]

@KarinaAsmar-NOAA Thanks. I realized I had a muddled interpretation of the fix #117. The problem is that apparently whatever triggered the error that prompted #117 has now been resolved because you have trimmed the variable lists.

I'm a little confused as to whether your original fix was required because of ocean fields or ice. In the ocean.csv that was committed with your #117, I don't see any variables in the ocean list that would have created more than 1 vector pair with the same mapping. I do see such variables in the ice.csv, where both the 2d ocean and atm velocites were being mapped bilinearly.

It should be possible to know how many vector pairs are going to be used for a given dimensionality and a given mapping type. The vectorpairs can then have an additional dimension to track each pair and loaded into the flds array.

Since this is clearly a bug which could arise in the future, I think an issue should be created that the fix PR you are making here will lead to failure when more than 1 vector pair for a given mapping type is requested.

@DeniseWorthen I agree. I'll create the issue in this repository and link it to this PR and #117 . If one of the models using the ocnicepost tool generates history files with more than 1 vector pair for a mapping type, that can be used to create a fix for this bug.

[DeniseWorthen]

@KarinaAsmar-NOAA Sounds good. I'm also going to create an issue to add some hardening to the code. For example, it is setup to write out the min/max of the original source grid fields, those fields packed, and those fields mapped. All of those min/max values should be w/in tolerances of each other (the packed field min/maxs should actually be identical to the original field min/max).

The code should fail if any of the steps produces non-sensical values. It shouldn't require eyes on the output; since this code is being used as a black box, it should simply fail if a change has been made such that SSV=0.0 is created. Right now, it relies on the user to read the ocean.post.log for example to see that there is an obvious problem.

[KarinaAsmar-NOAA]

@DeniseWorthen Please see #174. I was able to reproduce the bug with the ice_sfs.csv list. At this time, we are not using the ocnicepost tool for SFS ice products so this bug is not affecting global-workflow SFS runs.

[KarinaAsmar-NOAA]

@DeniseWorthen Would you please confirm whether the fix for the vector pair bug is acceptable? Test results are in /scratch4/NCEPDEV/ovp/Karina.Asmar/vecpairbug. Both ice and ocean vector pairs have non-zero values now. There are no errors now when processing more than one vector pair with the same remapping method. Please see Issue #174 for reference.

[DeniseWorthen]

@KarinaAsmar-NOAA Thanks. Have you confirmed that both vector pairs (with your current fix) are B4B with the case where you have only one vector pair with the change in #117 reverted?

In other words, we know if you revert your changes in #117, the mapping of a single vector pair is correct. If your new fix is indeed correct, I think each of those vector pairs should be B4B with what reversion of #117 gives you.

[KarinaAsmar-NOAA]

@KarinaAsmar-NOAA Thanks. Have you confirmed that both vector pairs (with your current fix) are B4B with the case where you have only one vector pair with the change in #117 reverted?

In other words, we know if you revert your changes in #117, the mapping of a single vector pair is correct. If your new fix is indeed correct, I think each of those vector pairs should be B4B with what reversion of #117 gives you.

What is B4B?

[DeniseWorthen]

B4B is bit-identical comparison. You can use nccmp for example

nccmp -d -S -q -f -B --Attribute=checksum --warn=format file1.nc file2.nc

If they are not identical, you will get identifying information for what is not identical.

[KarinaAsmar-NOAA]

nccmp -d -S -q -f -B --Attribute=checksum --warn=format file1.nc file2.nc

I ran this command and got no warnings:

[Karina.Asmar@ufe04 grads]$ nccmp -d -S -q -f -B --Attribute=checksum --warn=format ocean.1p00.nc fix172.ocean.1p00.nc
[Karina.Asmar@ufe04 grads]$

ocean.1p00.nc - bug fix from #174 (works for more than one vector pair)
fix172.ocean.1p00.nc - bug fix from #172 (crashes with more than one vector pair)

Visual comparisons are included here

[KarinaAsmar-NOAA]

nccmp -d -S -q -f -B --Attribute=checksum --warn=format file1.nc file2.nc

I ran this command and got no warnings:

[Karina.Asmar@ufe04 grads]$ nccmp -d -S -q -f -B --Attribute=checksum --warn=format ocean.1p00.nc fix172.ocean.1p00.nc [Karina.Asmar@ufe04 grads]$

ocean.1p00.nc - bug fix from #174 (works for more than one vector pair) fix172.ocean.1p00.nc - bug fix from #172 (crashes with more than one vector pair)

Visual comparisons are included here

Cleaner comparison:

[Karina.Asmar@ufe04 Karina.Asmar]$ nccmp -d -S -q -f -B --Attribute=checksum --warn=format vecpairbug/ocean.1p00.nc ocnpost_v2/ocean.1p00.nc
[Karina.Asmar@ufe04 Karina.Asmar]$

where ocnpost_v2/ocean.1p00.nc is the file with reverted #117 changes and vecpair/ocean.1p00.nc is the most recent file with the bug fix for more than one vector

[DeniseWorthen]

@KarinaAsmar-NOAA The fix for more than one vector pair with a given mapping type does seem to be working now. I'm going to approve. I still think the code needs hardening, but that is of course a separate issue.

[KarinaAsmar-NOAA]

@DavidHuber-NOAA @TravisElless-NOAA Can this PR be merged?

[DavidHuber-NOAA]

Changes look good to me, though I do not have the expertise to comment on the science.