Incorporate Sonde DA into JEDI for RRFS

[delippi]

Description

The UFO in JEDI is the component that not only computes model-simulated observations, but also houses filters and methods for observation QC, ob errors, and bias correction. The GSI observer is the equivalent. Many forward operators for various observations have been developed for the UFO. These operators can be utilized in RDAS. However these operators still must be tested for RDAS. The steps for transition by observation platform are as follows:

1. Establish IODA processing
2. Establish operator in UFO
3. Establish YAML input file
4. Validate with test assimilation experiments
5. Create ctest, if necessary

Sonde data is bufr obtype=*20

Requirements

To create yaml configurations for assimilating sonde data.

Acceptance Criteria (Definition of Done)

• Phase 1 validation (hofx validation using GSI-IODA and GSI-geovals as necessary; GSI vs FV3-JEDI)

• Phase 2 validation (QC validation; no reliance on GSI except to be used as the baseline; GSI vs FV3-JEDI)

• Phase 3 validation (FV3-JEDI and GSI vs MPAS-JEDI)

• Link any relevant pull requests here:

  • Phase 2 validation: sonde (adpupa 120/220) #245
  • PR for Phase 3 (to be added later)

Dependencies

• Phase 2 validation: updates for msonet yamls #188 is useful to be used as a building block for this work.
• This is one of the more important ob types to start cycling MPAS-JEDI so we need this ready ASAP.

[delippi]

Phase 2 update
I did Phase 1 validation awhile ago for these types (unofficially; not documented). They use the exact same operator as mesonet, so I feel confident about it. I will likely still do Phase 1 again just to have it documented.

Here are some preliminary results with some sonde yamls mostly based on the mesonet yamls in PR #188. Ob counts are definitely incorrect. Even when removing all QC the ob counts are badt (for example, only 96 airTemp obs used with/without any QC). This has been noticed before; need to investigate. Increments look very different as well; probably due to mismatch in ob errors? Omb look pretty good (which should mean hofx in Phase 1 is fine).

airTemperature

specificHumidity

winds

[delippi]

Phase 2 update
The discrepancy in observation counts between JEDI and GSI for sounding observations was due to two factors:

1. GSI was using ext_sonde which isn't something we ought to be using. This namelist option interpolates the observations such that there is at least one pseudo ob per model level.
2. The time window of assimilation is different. GSI uses time_window_max=1.5, nhr_assimilation=3, and nhr_obsbin=3 where JEDI seems to only have a time window of 1 hour. For some reason I had to change the GSI time_window_max=0.99 since setting it to 1.0 still did not work. We need to figure out how to extend the window for JEDI to match the 1.5 window for obs. If I recall correctly, this is important to get higher latency obs that might not be otherwise available at the expense of assimilating some obs more than once.

Now observation counts match exactly (or within 4 for winds) which is the expected result for just a handful of sounding data.

The current issue is the observation errors do not match close enough. In fact, I've checked that they do not get inflated at all in JEDI (I checked this by prescribing a constant ob error in JEDI and there was no change). From phase 1 testing, I seem to recall that the ObsErrorFactorPressureCheck really only made a difference for surface obs. I've already tried ObsErrorFactorConventional, but I haven't been successful getting any movement in the errors with that (I had to test it with the offline domain check; it seems this function also doesn't correctly use the qcflags from the online domain check). I need to dig into the GSI code to find out what it is doing with these obs. Then I can more easily determine what I need to do in JEDI.

airTemperature

specificHumidity

winds

[SamuelDegelia-NOAA]

@delippi Glad we are understanding the time window better! Do you mean that we should disable ext_sonde for our GSI runs? I am setting up the GSI run directory for the updated FV3-LAM case and can make this change if needed.

Also, regarding the obserrs, I recall that GSI has a routine that inflates the obserrs for profiles such as sondes based on their density. I just checked and that inflation occurs in the errormod function in qcmod.f90. Not sure if that is the difference here but it could be worth checking.

[delippi]

@delippi Glad we are understanding the time window better! Do you mean that we should disable ext_sonde for our GSI runs? I am setting up the GSI run directory for the updated FV3-LAM case and can make this change if needed.

Also, regarding the obserrs, I recall that GSI has a routine that inflates the obserrs for profiles such as sondes based on their density. I just checked and that inflation occurs in the errormod function in qcmod.f90. Not sure if that is the difference here but it could be worth checking.

Our messages crossed paths. Yes I think we need to disable ext_sonde for anything RRFSv2. Shun said it is probably too late for v1.

[delippi]

Also, regarding the obserrs, I recall that GSI has a routine that inflates the obserrs for profiles such as sondes based on their density. I just checked and that inflation occurs in the errormod function in qcmod.f90. Not sure if that is the difference here but it could be worth checking.

@SamuelDegelia-NOAA, thanks for sharing this info. I've confirmed that errormod in qcmod.f90 is being used to inflate my obs in the GSI case. The ObsErrorFactorConventional routine was designed to mimic the GSI observer code (i.e., subroutine errormod in qcmod.f90). So I need to keep trying to use that function and see if I can get any movement in the ob errors.

[delippi]

I'm able to get ObsErrorFactorConventional work now. You must use obsgrouping as follows:

         obsdatain:
           engine:
             type: H5File
             obsfile: "@OBSFILE@"
           obsgrouping:
             group variables: ["stationIdentification", "dateTime"]
             sort variable: "pressure"
             sort order: "descending"

An example of ObsErrorFactorConventional is

         # Error inflation based on errormod (qcmod.f90)
         - filter: Perform Action
           filter variables:
           - name: windEastward
           where:
           - variable: ObsType/windEastward
             is_in: 220
           action:
             name: inflate error
             inflation variable:
               name: ObsFunction/ObsErrorFactorConventional
               options:
                 inflate variables: [windEastward]
                 pressure: MetaData/pressure

The ob errors look so much better, but I'm not sure why my offline domain check has more observations... I need to create an issue to add qcflag checks in the ObsErrorFactorConventional much like I recently did for ObsErrorFactorPressureCheck so we can just use the online domain check.

[delippi]

Phase 2 update
The difference in observation counts in the previous update was due to user error having modified the ioda data for debugging and forgot to revert--all time information was set to be the analysis time so all observations were accept.

There was also another user error in the configuration used in the previous configuration settings. This is the correct usage (there should NOT be "dateTime" in the obsgrouping.group variables):

         obsdatain:
           engine:
             type: H5File
             obsfile: "@OBSFILE@"
           obsgrouping:
             group variables: ["stationIdentification"]
             sort variable: "pressure"
             sort order: "descending"

Here are the phase2 results. Note this in a 19Z cycle so the impact is small. I'm including the analysis increments anyway since they show min/max increment values. These results are also using the online domain check obs. I think we're ready to submit a PR for phase2 sonde data:

airTemperature

specificHumidity

winds

[JingCheng-NOAA]

In my case, for airTemperature/specificHumidity data, the observation number looks similar to GSI after applying all the QC filters, however, the output QCflag doesn't show anything related to out of domain obs, which is strange because in wind data it shows the number out of domain.

[delippi]

Phase 3 update

Now seemed like a good time to start phase 3 since I wanted to start digging into some errors I was seeing before in MPAS-JEDI and just make sure they are ready for cycling.

1. MPAS-JEDI isn't currently able to get saturation_water_vapor_mixing_ratio_wrt_moist_air in by using request_saturation_specific_humidity_geovals: true in ObsFunction/ObsErrorFactorPressureCheck. A temporary solution is to just comment out that entire filter in the specificHumidity yamls. ObErrors just won't be inflated correctly. Hence the ob errors for specific humidity in all MPAS cases (for now) will look uninflated.
2. Obs counts for specificHumidity are strangely low in this case; they were exact in the old 19Z FV3-JEDI case. I found this can be fixed by removing the initial Reject all obs with QualityMarker > 3. I need to double check what we should be doing. I have a feeling the QM isn't as robust as I'm thinking and we could be more lenient. (see bonus figure at the end). The more relaxed QM handling also doesn't seem to change the final increment results. Not sure if it would be better to keep it more stringent to reduce number of obs since there isn't as much of an impact?
3. ObErrors for winds don't seem to follow GSI as expected.
4. Increments, for the most part, look similar spatially but GSI is more smooth. That shouldn't be related to the use of observations.

airTemperature

specificHumidity

winds

Bonus figures: specificHumidity without Reject all obs with QualityMarker > 3

[delippi]

Another Phase 3 update
This also includes some tuning associated with Phase 2-type updates. I updated all my testing to the new case which more accurately reflects the RRFSv1 configurations in errtable, convinfo, gsiparm.anl, etc which also needs to be reflected in the JEDI yamls. With Phase 3, we introduce MPAS-JEDI testing; although FV3-JEDI and MPAS-JEDI tend to be fairly similar. The most important thing to compare are the omb and oberrors. Differences in increments can be affected by other parts of the assimilation, of course. There were a few differences in omb and oberrors that I haven't accounted for, but might be okay:

1. specificHumidity: ob counts for JEDI are about half of that of GSI.
2. winds: ob errors for JEDI are much more flat compared to GSI.

From a use of observation validation perspective, I think sonde looks really good now. I'm not sure if additional tuning would provide much benefit over GSI at this point.

[delippi]

Another Phase 3 update
Turned off SDL/VDL in GSI for comparisons. Updated gross error settings to match RRFSv1. Fixed a bug with double inflating windEastward.