Satwind thinning in JEDI with Gaussian Thinning function

[BrettHoover-NOAA]

This issue wraps the thinning work that is being applied to AHI/Himawari (#741) and SEVIRI/METEOSAT (#758) satwinds.

In all testing, the Gaussian Thinning function appears to be more severe than GSI's thinning function when applying similar thinning meshes. The goal here is not to perfectly reproduce GSI's thinning, but rather to find a qualitatively similar thinning mesh for JEDI for proper comparisons to GSI.

[BrettHoover-NOAA]

First observation: The Gaussian Thinning function's impact on assimilated ob density is different depending on when the function is called in the QC YAML. Using the same QC filter:

  - filter: Gaussian Thinning
    horizontal_mesh: 200
    vertical_mesh: 10000

I tried putting this in the obs pre filters section and then again in the obs post filters section. Results on AHI/Himawari WVCT satwinds:

pre-filter: 2507 obs assimilated, 45771 obs thinned, 1589 obs rejected, 49867 total
post-filter: 3241 obs assimilated, 31950 obs thinned, 14676 obs rejected, 49867 total

By placing the same Gaussian Thinning filter in the post rather than pre filter section, 29% more satwinds are assimilated, 43% fewer satwinds are thinned, and nearly 10 times more satwinds are rejected. I believe this drastic difference in outcome is a function of the sequencing between thinning and QC. As a pre-filter, the sequence is thin-->QC, and as a post-filter (at the end of the list of post-filters), the sequence is QC-->thin. So as a pre-filter thinning is applied to all observations prior to QC and then QC is applied to what's left, which is why so many obs are thinned and so few are rejected. As a post-filter, QC is applied first to all observations and then thinning is applied to what's left. The post-filter version is more true to GSI's thinning methodology, where thinning is applied just before exiting setupw.f90 following all QC filters. However, it is worth noting that with these mesh-grids, which are identical to the meshes applied in GSI, we are still getting far fewer observations assimilated and too many thinned:

GSI: 6450 obs assimilated, 28630 obs thinned, 3619 obs rejected

It's also worth recognizing that the thinned ob-count is not unique in either GSI or UFO. In GSI, thinned observations appear to be dropped from the diag file entirely, so the number of thinned winds needs to be inferred by comparing a thinned-GSI experiment and a non-thinned-GSI experiment. In UFO, thinning as a pre-filter will apply an EffectiveQC value of 16 to the thinned obs and then QC will not be further applied, which is why the number of rejected obs changes based on where the thinning filter is applied - you end up thinning otherwise-rejected obs in a pre-filter and only thin otherwise-accepted obs in a post-filter.

[BrettHoover-NOAA]

Some updates to this work: I have what I believe is an acceptable thinning-setting framework for JEDI to qualitatively match GSI, at least in horizontal and vertical dimensions (i.e. not the time-dimension).

Some observations:

1. GSI performs thinning on some but not all non-GOES GEO AMVs. All three JMA types (Type 252/LWIR, 242/VIS, 250/WVCT) are thinned, as are all three EUMETSAT types from METEOSAT-10 (Type 253/LWIR, 243/VIS, 254/WVCT). EUMETSAT AMVs from METEOSAT-9 do not thin Type 254/WVCT, per global_convinfo.txt. There are very few Type 243/VIS winds that pass QC (<10 in the test data), so these are excluded from analysis here.

2. GSI provides very little priority for determining the observation that is retained in thinning. The value of crit1 determines thinning priority. This variable is initially assigned to zero when thin4d=.true., and no time-based prioritization is assigned. For EUMETSAT AMV types, a crit1 value is computed based on the quality indicator. A crit1 value is also assigned to GOES winds based on QI and EE scores, although no thinning of GOES AMVs is performed per global_convinfo.txt. All other thinned AMVs are treated as equal priority with crit1=0..

3. GSI performs vertical thinning on a reduced-resolution grid in model-level space. It appears that the retained levels are pressure-based with pressure differences defined by pmesh, but in testing it appears that GSI retains more observations in the upper troposphere than JEDI when using the same value for vertical_mesh, which may be a byproduct of the GSI operating on terrain-following coordinate values.

In testing, I find that the histogram of retained observation pressure values closely matches between JEDI and GSI when the vertical_mesh value is set to 10000. (Pa), equivalent to the pmesh value of 100. (hPa) for all thinned types. If I reduce the vertical_mesh value to retain more winds, it also pushes the distribution of retained observation pressure values to lower pressures, with peaks appearing in JEDI at discretely higher altitude / lower pressure than GSI. As an example, comparing GSI (blue) to JEDI (orange) histograms of retained pressure values among all thinned wind types when vertical_mesh is reduced to 9000. for observations below 50000. Pa produces these shifted distributions:

While peaks are in the right places in the lower troposphere, the peaks are shifted to lower pressures in the upper troposphere. The vertical position of AMVs can be very important for impact testing, as the observations can be sensitive to their height assignment, and as a result I consider it a high priority to keep the JEDI pressure distributions as close to the GSI distributions as possible for comparison purposes.

In general, JEDI's horizontal_mesh and vertical_mesh filter functions leave sparser observations in the upper troposphere than they do in the lower troposphere. As an example, a test with horizontal_mesh set to 100. (km) throughout the column with a constant vertical_mesh of 10000. (Pa) shows that the JEDI thinning is too severe in the upper troposphere:

Since altering the value of vertical_mesh will change JEDI's thinning to prioritize observations at lower pressure and create differences in the pressure distribution relative to GSI, it was decided that this over-thinning in the upper-troposphere should be addressed by allowing for a finer horizontal mesh in the upper troposphere. This is done with the use of where statements and 2 applications of the Gaussian Thinning filter, one with settings for observations in the lower troposphere and the other with settings for observations in the upper troposphere, delineated by pressure values at 50000. Pa. An example:

  - filter: Gaussian Thinning
    horizontal_mesh: 100
    vertical_mesh: 10000
    where:
    - variable:
        name: MetaData/pressure
      minvalue: 50001.
  - filter: Gaussian Thinning
    horizontal_mesh: 90
    vertical_mesh: 10000
    where:
    - variable:
        name: MetaData/pressure
      maxvalue: 50000.

This 2-filter approach with an adaptive horizontal mesh appears to work well to retain both the vertical distribution as well as total number of retained satwinds when comparing JEDI and GSI. With these filters (modified for use of priority_variable=MetaData/qiWithoutForecast for EUMETSAT winds and excluding thinning of METEOSAT-9 Type 254 winds), vertical distributions look close in both the upper/lower troposphere:

Coverage of AMVs across all thinned types is very close as well, differing by fewer than 500 winds in either the upper/lower troposphere:

However, the fit is not equal across all types. There is type-based variation in the number of retained winds and where there is under/over-sampling in the vertical. In this test, I find the following type-based dependencies:
AHI H9
Type 252 (LWIR): GSI retains 7401, JEDI retains 5151 (-30% change)
Type 242 (VIS): GSI retains 3190, JEDI retains 3034 (-5% change)
Type 250 (WVCT): GSI retains 7246, JEDI retains 9741 (+34% change)
METEOSAT-9
Type 253 (LWIR): GSI retains 5069, JEDI retains 5927 (+17% change)
Type 243 (VIS): GSI retains 617, JEDI retains 1117 (+81% change)
METEOSAT-10
Type 253 (LWIR): GSI retains 6673, JEDI retains 5142 (-23% change)
Type 254 (WVCT): GSI retains 5650, JEDI retains 6601 (+17% change)

In general, LWIR winds are undersampled in the upper-troposphere in H9 and M10, and undersampled in the lower-troposphere in M9, and WVCT winds are oversampled in the upper-troposphere in H9 and M10. The large percentage of oversampling of M9 Type 243/VIS winds is only among a base of 617 observations, so the effect-size is still small. These type differences tend to balance each other out with oversampling in one type and undersampling in another in a way that makes further modification of the thinning filters an undesirable way to approach reaching parity.

The thinning filters being applied here are as follows:
AHI HIMAWARI-9

  - filter: Gaussian Thinning
    horizontal_mesh: 100
    vertical_mesh: 10000
    where:
    - variable:
        name: MetaData/pressure
      minvalue: 50001.
  - filter: Gaussian Thinning
    horizontal_mesh: 90
    vertical_mesh: 10000
    where:
    - variable:
        name: MetaData/pressure
      maxvalue: 50000.

SEVIRI METEOSAT-9

  - filter: Gaussian Thinning
    horizontal_mesh: 100
    vertical_mesh: 10000
    priority_variable:
      name: MetaData/qiWithoutForecast
    where:
    - variable: ObsType/windEastward
      is_in: 243,253
    - variable:
        name: MetaData/pressure
      minvalue: 50001.
  - filter: Gaussian Thinning
    horizontal_mesh: 90
    vertical_mesh: 10000
    priority_variable:
      name: MetaData/qiWithoutForecast
    where:
    - variable: ObsType/windEastward
      is_in: 243,253
    - variable:
        name: MetaData/pressure
      maxvalue: 50000. 

SEVIRI METEOSAT-10

  - filter: Gaussian Thinning
    horizontal_mesh: 100
    vertical_mesh: 10000
    priority_variable:
      name: MetaData/qiWithoutForecast
    where:
    - variable:
        name: MetaData/pressure
      minvalue: 50001.
  - filter: Gaussian Thinning
    horizontal_mesh: 90
    vertical_mesh: 10000
    priority_variable:
      name: MetaData/qiWithoutForecast
    where:
    - variable:
        name: MetaData/pressure
      maxvalue: 50000.

[BrettHoover-NOAA]

I have submitted PRs in NOAA-EMC/jcb-gdas#25 and #1263 to enroll AHI-H9, SEVIRI-M9/10 satwinds and add these thinning filters to their corresponding yaml.j2 files.