[bug] rrfs-workflow [dev-sci] no observations DA doesn't reproduce 1h forecast

[delippi]

Current behavior (describe the bug)

Our recent retros have shown that JEDI performs slightly worse than GSI in the lower troposphere. While investigating this, @SamuelDegelia-NOAA and I discovered that cycling DA (no observations) produces vastly different results.

At present, JEDI does not directly write (that we know of) the D-grid u/v winds needed to restart FV3, so an external tool rdas_ua2u has been used.

We believe there is some issue with this tool to convert from A-grid winds (ua/va) to D-grid winds (u/v) using this tool causing the degraded performance in our current retros.

Steps to Reproduce (if applicable)

1. Setup a retro with DA_SYSTEM=GSI and another with DA_SYSTEM=JEDI.
2. Configure both GSI and JEDI to monitor/passivate all observations (i.e., run through the DA with no observations assimilated) to test all the machinery.
3. Compare any corresponding forecasts from any warm start cycle.

Expected behavior

We expect exactly reproducible results due to zero observations assimilated. There should be no (minimal) impacts from stopping and restarting the model.

Suggested Solution (if known)

Either

1. Correct the rdas_ua2u tool.
2. Learn how to use (if exists) or implement (if not) such variable transform directly in JEDI.

[delippi]

This is 1 hour forecast differences at roughly 850hPa between GSI and JEDI (both assimilating no observations) in the 01Z cycle in our simplified workflow (cold starts are twice per day at 00/12Z). The 01Z cycle 1h forecast would see smallest impacts which are pretty large (as large as 7 m/s). The forecast differences are much larger in the lower troposphere and much larger in the 11Z cycle (which would see greatest impact from this).

This corresponds exactly to what we see in the MATS verification. This verification is performed against RAOB for the 1 hour forecast of wind; therefore, it is verifying the 11Z 1h forecast valid at 12Z and the 23Z 1h forecast valid at 00Z (both of which would see the greatest impact from this bug). This result is over a 1 week retro period.

[delippi]

@xyzemc Thanks for the suggestion to plot the difference relative to the analysis file in order to know for sure if it is the rdas_ua2u tool. I think this also answers @TingLei-NOAA's suggestion to run the tool from a state with known correct ua/va (from GSI) and run the ua2u tool then compare the resulting u/v. This is exactly showing that, as you can see the ua option is grayed out denoting that when I clicked on it all diffs were 0.

[TingLei-NOAA]

@delippi Thanks for sharing. But the differences from the analysis might come from use of the different analysis variables (ua/va vs ud/vd) not necessarily from conversion between ua/va to ud/vd. But this difference could also be the reason for better or worse forecast.
My suggestion was to directly use ua/va and ud/vd in the background , both of which are generated from fv3 model (the truth) to test the conversion tool (in analysis, gsi won't update ua/va using the updated ud/vd).

[delippi]

@TingLei-NOAA I think you suggesting that u/v and ua/va aren't necessarily consistent in an analysis file but should be in a background file. I just tried again with a background file and it seems to give the same result.

Here are some clear steps for what I tried

cp ../det_c13_3dvar_noDA_gsi/2024050801/fcst_fv3lam/INPUT.gsi/bk_fv_core.res.tile1.nc bkg_truth.nc
ncks -O -x -v u,v bkg_truth.nc bkg_test_uava.nc
./rdas_ua2u.x ua_update_u --in_grid=fv3_grid_spec --in_file=bkg_test_uava.nc --out_file=bkg_test_out.nc
ncdiff bkg_truth.nc bkg_test_out.nc bkg_truth_diff.nc

[TingLei-NOAA]

@delippi Thanks. That is definitely pointing to some problems in the conversion.
Can we directly use hdas_ua2u.x to do the conversion to see where the problem is?

[delippi]

@TingLei-NOAA, yes. I still had that executable and it produced the exact same result. Which makes sense because when I transitioned the code, I had only tested that our tool replicated the hafs_datool and didn't check things otherwise.

Are there differences in the definitions of winds between RRFS and HAFS? Perhaps the tool is correct for them but is incorrect for us?

How does global handle restarting the forecast?

This may take some time to understand and find a solution for. This step is critical before we can run any new retros.

[TingLei-NOAA]

@yonghuiweng and @XuLu-NOAA what are your opinions on the issue with the the wind conversion found in rdasapp cases?
Have you compare the conversion with the model's own D-wind in hdas cases?
Thanks

[ShunLiu-NOAA]

I had a discussion with @delippi . We think we should interpolate the increment of A-grid wind to D-grid and then combine D-grid increment with D-grid background. It may reduce the difference from the interpolation of the full A-grid wind.

[XuLu-NOAA]

@yonghuiweng and @XuLu-NOAA what are your opinions on the issue with the the wind conversion found in rdasapp cases? Have you compare the conversion with the model's own D-wind in hdas cases? Thanks

I know there will be some info losses due to the grid conversion between u and ua. But I don't know if we have compared it with the model conversion. @yonghuiweng may have more info. My memory is that we grabbed the conversion from GSI, not from FV3 directly.

[TingLei-NOAA]

@XuLu-NOAA Thanks. I think I did some comparison with GSI 's conversion (forgot the details) and found the gsi conversion produced results the same with the model's conversion. but maybe, at that time, I focused on some special aspect where, then, fv3jedi had a problem. I will figure out a way to do the verification again for the conversion in gsi.

[TingLei-NOAA]

@delippi and @XuLu-NOAA I did some tweeks in gsi , to make gsi directly convert the analysis A wind (with 0 wind increment) to D wind( they first went through interpolation from analysis grid to D grid and then and then converted to D -wind and saved to the analysis file). The difference between such analysis file and background is diff.inc in ursa://home/Ting.Lei/dr-out/.

Ignore the large values along the lateral boundaries (related the treatment of there along the boundaries), I still see differences up to 2 m/s in u in the inner domain. So, the gsi's conversions method alongside the above conversion (from analysis grids to model grids) could yield the differences in the D wind of the similar magnitudes as you showed. Maybe more tweeks could give a more apple to apple comparison. But for being now, I can't rule out the possibilities that gsi 's conversion could cause the differences you showed here.
One way is as @ShunLiu-NOAA and you discussed , first convert the A increment to D wind increment and then add them on the background as GSI does.
Or dig more to see the conversion methods could be improved still using the model's version as the "truth"?

[delippi]

Edit: I should note these experiments were all done after implementing the rdas_ua2u tool on the analysis increments and then applying them on top of the background (as opposed to regridding the full A winds). The tests below just show there is more to it than just the rdas_ua2u tool. There will still likely be more differences from the tool when increments are > 0.

In a series of zero-observation experiments where DA is run in both systems, I consistently see a 0.5-1 m/s difference in the 1-hour forecast between JEDI and GSI. Diagnostics show that in GSI, the zero-obs DA cycle still produces very small but nonzero analysis minus background differences (order 1e-5 in temperature and 1e-7 in sphum), while JEDI’s zero-increment analysis does not introduce these offsets. To test whether the discrepancy comes from the DA step itself or from the forecast configuration, I ran a control experiment where I forced GSI’s background to be used as its analysis (i.e., background equals analysis, as zero-obs should ideally produce). In this configuration, the 1-hour forecasts from JEDI and GSI reproduce each other, confirming that the underlying forecast model setup is consistent, and that the discrepancy in the original zero-obs runs originates from how GSI constructs its analysis in a nominal zero-obs configuration, not from JEDI or the forecast integration.

1. Here are the 1h forecast differences where both systems assimilate zero-observations but GSI uses its analysis restarts.

Model top:

Somewhere in the middle of atmosphere. The rest of the layers looks similar but of different magnitudes.

2. For the 1h forecast differences for the case where I force GSI to use it's background restarts, it is sufficient to show the checksum that they are identical. I confirmed the result is the same here for JEDI with and without using rdas_ua2a with 0-increments.

diffing file: fcst_fv3lam/dynf001.nc
  GSI:  0af0ef740b582f705c7d53e61e1e41cfdcb71a7d  det_c13_3dvar_noDA_gsi/2024050601/fcst_fv3lam/dynf001.nc
  JEDI: 0af0ef740b582f705c7d53e61e1e41cfdcb71a7d  det_c13_3dvar_noDA_jedi/2024050601/fcst_fv3lam/dynf001.nc

3. Here are some analysis minus background differences for GSI (zero-obs analysis). The differences mainly occur in first 4 levels (top of model), but I did see some nonzero diffs for Q around model level 40 too.

T model level=3

Q model level=3

Q model level=40

[HuiLiu-NOAA]

Edit: I should note these experiments were all done after implementing the rdas_ua2u tool on the analysis increments and then applying them on top of the background (as opposed to regridding the full A winds). The tests below just show there is more to it than just the rdas_ua2u tool. There will still likely be more differences from the tool when increments are > 0.

In a series of zero-observation experiments where DA is run in both systems, I consistently see a 0.5-1 m/s difference in the 1-hour forecast between JEDI and GSI. Diagnostics show that in GSI, the zero-obs DA cycle still produces very small but nonzero analysis minus background differences (order 1e-5 in temperature and 1e-7 in sphum), while JEDI’s zero-increment analysis does not introduce these offsets. To test whether the discrepancy comes from the DA step itself or from the forecast configuration, I ran a control experiment where I forced GSI’s background to be used as its analysis (i.e., background equals analysis, as zero-obs should ideally produce). In this configuration, the 1-hour forecasts from JEDI and GSI reproduce each other, confirming that the underlying forecast model setup is consistent, and that the discrepancy in the original zero-obs runs originates from how GSI constructs its analysis in a nominal zero-obs configuration, not from JEDI or the forecast integration.

1. Here are the 1h forecast differences where both systems assimilate zero-observations but GSI uses its analysis restarts.

Model top:

Somewhere in the middle of atmosphere. The rest of the layers looks similar but of different magnitudes.

2. For the 1h forecast differences for the case where I force GSI to use it's background restarts, it is sufficient to show the checksum that they are identical. I confirmed the result is the same here for JEDI with and without using rdas_ua2a with 0-increments.

diffing file: fcst_fv3lam/dynf001.nc
  GSI:  0af0ef740b582f705c7d53e61e1e41cfdcb71a7d  det_c13_3dvar_noDA_gsi/2024050601/fcst_fv3lam/dynf001.nc
  JEDI: 0af0ef740b582f705c7d53e61e1e41cfdcb71a7d  det_c13_3dvar_noDA_jedi/2024050601/fcst_fv3lam/dynf001.nc

3. Here are some analysis minus background differences for GSI (zero-obs analysis). The differences mainly occur in first 4 levels (top of model), but I did see some nonzero diffs for Q around model level 40 too.

T model level=3

Q model level=3

Q model level=40

It would be interesting to show a couple vertical profiles of the differences for u/v, t, and Q at the max areas..

[delippi]

I found the culprit for why there is a difference of 1.e-07 in q for levels 1, 2, 3, and 4 (at model top).
In constants.f90:

  real(r_kind),parameter:: qmin   = 1.e-07_r_kind   ! lower bound on ges_q

I was suspecting this to be something intentional rather than a bug. Without this capability in JEDI, even zero-observation DA experiments can never have identical results.

Now, its a question of

1. Does this safety net already exist in JEDI and is just not currently turned on?
2. If not, do we want a similar function?
3. Where to implement it?

I think the answer is yes, we should add the capability if it doesn't already exist but make it configurable.