GPSRO Validation

[ADCollard]

Review GSI GPSRO QC code and existing options within JEDI UFO to determine whether we want to a) Reproduce the latest GPSRO QC in JEDI b) Take one of the exisiting GPSRO QC schemes For the latter we need to determine whether the committed YAML is acceptable based on statistical considerations.

[XuanliLi-NOAA]

The flow chart of GSI GPSRO QC steps are listed below:

[XuanliLi-NOAA]

4 GNSSRO forward operators are available in JEDI:

• NBAM
• BendMetOffice
• BndROPP1D
• BndROPP2D

NBAM is the same as GSI. Two QC steps are included in the NBAM operator:

1. Model top and bottom Check
2. Super refraction check (2 options available: NBAM and ROPP)

As illustrated below:

Statistic QC check is included in the script in the operators directory ufo/src/ufo/operators/gnssro/QC/BackgroundCheckRONBAM.cc

Obs error assignment has 4 options: NRL, ECMWF, NBAM, and latitude which are provided in ufo/src/ufo/operators/gnssro/utils/operators/gnssro/utils/gnssro_mod_obserror.F90 ufo/src/ufo/operators/gnssro/QC/ufo_roobserror_mod.F90 and ufo_roobserror_utils_mod.F90

RO Obs error inflation for superobs is in ufo/src/ufo/operators/gnssro/QC/action/ROobserrInflation.cc.

As illustrated below:

[XuanliLi-NOAA]

In the 4 forward operators for GPSRO data, NBAM, MetOffice, and ROPP1D are 1D bending angle. ROPP2D is 2D ray tracing operator. The operators are compared below:

[XuanliLi-NOAA]

A 1DVAR check is also available in ufo/src/ufo/filters/gnssroonedvarcheck directory. The flowchart is shown below:

[XuanliLi-NOAA]

Emily modified the iodoconverter (feature/geoval-fix branch) and created the test input files. The diag file looks fine. But the geoval file is empty and many information are missing from the obs file. The following variables need to be included in the obs file: impactHeightRO, impactParameterRO, airTemperature, specificHumidity, airPressure, pccf, earthRadiusCurvature, geoidUndulation, satelliteAscendingFlag, satelliteConstellationRO, sensorAzimuthAngle, qfro, sequenceNumber

The following variables need to be included in the geoval file: geoPotentialHeight, geoPotentialHeight_levels, heightOfSurface

Variables have been added to obs file: impactHeightRO, airTemperature, specificHumidity, airPressure

[XuanliLi-NOAA]

All variables have been added into the obs and geoval files. The scatterplots for JEDI HofX vs. GSI HofX indicated some inconsistency. Here are the plots without any QCs:

Out of the total of 496413 data, GsiEffectiveQC count 0 393458 1 102955

Counts of JEDI HofX - GSI HofX > 0.0000001: 7725/496413 = 1.5%

All of the 7725 data points have GSIEffectiveQC=1 and all of the data are below 500 hPa with most of them are below 850 hPa.

A few data have very large JEDI HofX (>0.05 in the first scatterplot), here is the list: pressure lat lon hofx gsihofx GsiEffectiveQC, EffectiveQC 90913.14 -3.64896 112.98019 0.06438975 0.02737023 1 0 90603.78 4.5712 343.8958 0.15309183 0.02130947 1 0 88755.63 2.8237 340.59164 0.07307973 0.02280927 1 0 93614.54 -23.07203 2.85219 0.053725474 0.020253355 1 0 92376.3 -8.54056 278.83075 0.14327128 0.02553164 1 0 94499.516 26.72656 336.0147 0.07363514 0.03510589 1 0 93843.76 13.49242 62.64037 0.07347244 0.03469352 1 0 97365.99 21.41083 242.72414 0.11198245 0.01606566 1 0 89258.7 -24.82843 357.80927 0.112184495 0.01445685 1 0 95190.625 -32.61547 281.68127 0.050192155 0.0194675 1 0 97310.08 20.68423 267.4164 0.17397434 0.02748836 1 0

[XuanliLi-NOAA]

Looking into the JEDI code and compared with GSI. There are a few differences (might have more):

1. All of the observations in the lowest two model levels are rejected in GSI. Line 549-554 of setupbend.f90 These observations are kept in JEDI, unless they are below the lowest model level. Line 270 of ufo_gnssro_bndnbam_mod.F90

2. The super-refraction were handled differently in JEDI from GSI. This check affects mostly the low troposphere observations. The super-refraction were considered for all obs below 5 km height in GSI. Line 521-538, 588-605 of setupbend.f90. The data within 5 levels above the top of super-refraction are rejected in GSI. This also affect the computation of dN/dx_s levels (Line 726-732). Super-refraction were consider for obs below 6 km in JEDI. The criteria in JEDI are also somewhat different from GSI. The data within 2 levels above the top of super-refraction are rejected. Line 292-326 of ufo_gnssro_bndnbam_mod.F90

3. HofX computation:
   In GSI, if the refractivity derivative dN/dx_s at any grids > 0, the observation is rejected. Line 754-761 of setupbend.f90 In JEDI, dN/dx_s is defined to be max(0, abs(dN/dx_s)). So if this value is positive at a certain level, it is set as 0, but the observation is not rejected. Line 97 -98 of ufo_gnssro_bndnbam_util_mod.F90 (The bending angle is a function of the integral of the derivative dN/dx_s over 80 equally-spaced vertical grids.) Line 105-111 of ufo_gnssro_bndnbam_util_mod.F90

[XuanliLi-NOAA]

The scatterplots for HofX with QCs in the swell yaml file are shown below:

[XuanliLi-NOAA]

A few tests were done.

1. Reject data below model level 2. This helps to remove the very large bending angles in JEDI, but the count with HofX difference < 0.0000001 did not change much. In the new plots below, only 1 (previously 11) value is larger than 0.05.
2. Changes in dN/dx_s did not help much.

[XuanliLi-NOAA]

A few tests were done.

1. Changed the limit to model level 3 to be consistent with GSI, line 549-554 of setupbend.f90.
2. Impact height: In GSI, the height (impact height - surface altitude, Line 520 of setupbend.f90) is used for SR check. While in JEDI, the impact height is used, line 259 of ufo_gnssro_bndnbam_mod.F90.
   These changes didn't change much about the total count with HofX difference < 0.0000001, but when the missing values were excluded, there are 1594 data left (previously 3290). Although the counts go down, but these counts are within 0.0001 to 0.0000001, so it doesn't help much with the plots. The large differences still exist.

[XuanliLi-NOAA]

A couple of changes were made to GSI.

1. Fixed a bug in the dimension of hges (originally nsig, should be nsig+1), line 224 of setupbend.f90. This did not cause any problems in our system though.
2. The HofX in GSI was computed as a local variable dbend. But when ddnj > 0, this variable was NOT updated (see line 754-761 of setupbend.f90), and the updates of rdiagbuf(17) and rdiagbuf(5) were also skipped. This caused the dbend for the previous observation data being reused and incorrect rdiagbuf(17) and rdiagbuf(5) were written to the gsi_diag file. This wouldn't cause problems in GSI since these data would be tossed due to 0 values of obs error, but it lead to the inconsistent HofX between JEDI and GSI. A fix was made by assigning the rdiagbuf(17) and rdiagbuf(5) to be missing for ddnj>0.

A change was made to UFO. In GSI, there is a large bending angle check. If HofX > 0.05, the data was tossed. This check is now added in the GNSSRO NBAM operator ufo_gnssro_bndnbam_mod.F90.

Here is the new plots for HofX:

[XuanliLi-NOAA]

Hailing has made several changes to the UFO NBAM operator, her PR is https://github.com/JCSDA-internal/ufo/pull/2891.

Now QC numbers in JEDI are consistent with GSI except for the commercial data (wait for GMAO changes, their PR is https://github.com/JCSDA-internal/ufo/pull/2872).

In the following plots, SR1 is flag for super refraction condition 1 (refractivity derivative > 0.75 critical_value) and SR2 is super refraction condition 2 (bending angle >0.03 and refractivity derivative > 0.5 critical_value)

Kompsat-5 (SAID=825):

There are 1 more data with JEDI SR2 than GSI SR2. This data has JEDIQC=19 and GSIPreQC=10, so it is already included in cutoffs.

COSMIC-2 SAID=750:

There are 100 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=19 and GSIPreQC=10, so these are already included in cutoffs.

COSMIC-2 SAID=751:

151 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=19 and GSIPreQC=10, so these are already included in cutoffs.

COSMIC-2 SAID=752: There are 109 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=19 and GSIPreQC=10, so these are already included in cutoffs.

COSMIC-2 SAID=753: There are 133 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=19 and GSIPreQC=10, so these are already included in cutoffs.

COSMIC-2 SAID=754: There are 69 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=19 and GSIPreQC=10, so these are already included in cutoffs.

COSMIC-2 SAID=755: GSI tossed 3 more data than JEDI due to SR2 (JEDI SR2 and GSIPreQC=8).
51 data with JEDI SR2 and GSIPreQC=10, so these are included in cutoffs.

MetOp SAID=3: There are 149 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=12 and GSIPreQC=12 due to MetOp 8km check.

MetOp SAID=4: There are 42 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=12 and GSIPreQC=12 due to MetOp 8km check.

MetOp SAID=5: There are 104 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=12 and GSIPreQC=12 due to MetOp 8km check.

PAZ SAID=44: There are 4 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=19 and GSIPreQC=10, so these are already included in cutoffs.

TerraSAR SAR SAID=42: All numbers are consistent.

GeoOptics SAID=265: Need to rerun the number for cutoffs when the cutoff is updated with GMAO changes. GPSTOP: 9459 data with GSIQC=2. Out of these data, 3124 data have JEDIQC=13 (domain check) and 6335 data have JEDIQC=12 (bounds check) There are 24 more data with JEDI SR2 than GSI SR2. These data have JEDIQC=19 and GSIPreQC=10, so these are already included in cutoffs.

[XuanliLi-NOAA]

Worked with Hailing to solve the inconsistent observation error issue. When we compared the initial error (before inflation, saved as GSI_adjust_obs_error in the diag file), JEDI and GSI are already different.

As a result, the final obs errors are different as well.

A bug was found in ufo_roobserror_mod.F90 (https://github.com/JCSDA-internal/ufo/blob/develop/src/ufo/operators/gnssro/QC/ufo_roobserror_mod.F90#L206C12-L206C12). Line 206, obsImpA should be obsImpH. This solved the inconsistency in the initial obs error.

Still needs more work on the obs error inflation.

[XuanliLi-NOAA]

Hailing fixed a bug in ufo_roobserror_mod.F90 in which obsImpA should be used, rather than ogsImpH. This fixed the inconsistency between GSI adjust error and JEDI obs error before inflation. Another bug was found in ROobserrInflation.cc where the record number was not read properly. With both fixes, the JEDI effective error is consistent with GSI final obs error.

[XuanliLi-NOAA]

There remains one issue. The QC flags in JEDI is slightly different with GSI for commercial data and COSMIC-2.
Total tossed numbers: 14043 in JEDI and 14025 in GSI which means 18 more data being tossed either due to cutoffs or SR2.

[XuanliLi-NOAA]

The inconsistency in the QC flag arises from various reasons. One reason that we found is the truncate difference between JEDI and GSI when doing the computation of cutoff values. There is one data out of the whole dataset having this issue:   JEDIQC=19 GSIQC=0 satid=265 imph=40294.50 Seq#=1785 obsValue=6.5000e-02 JEDIhofx=6.8384e-02 GSIhofx= 6.8384e-02 JEDIobserr=-3.3688e+38 GSIobserr=4.1317e-06 According to my discussion with Hailing, it appears that this is ok.

As for the remaining the cases, the reason for the disagreement is due to super refraction 2. Conditions for super refraction 2:

1. gradRef > critical_value * 0.5
2. obsValue >= 0.03 Rad
3. ImpP <= ImpP(max_obs_profile) (All data below the maximum obs within the profile will be tossed) GSI uses all three conditions. Hailing is adding an option in UFO to bypass the last one.  

The differences in super refraction 2 between JEDI and GSI cause two situations:

1) GSIQC = 0 while JEDIQC ≠ 0. There are 3 data falling into this category. JEDIQC=15 SR_flag_JEDI=2 GSIeffQC=0 GSIQC=0 satid=750 imph=2519.50 Seq#=704 obs=3.2821e+01 JEDIhofx=-inf GSIhofx=3.6405e+01 JEDIobserr=-3.3688e+38 GSIobserr=1.8419e-03 JJEDIQC=15 SR_flag_JEDI=2 GSIeffQC=0 GSIQC=0 satid=750 imph=2838.00 Seq#=833 obs=3.2939e+01 JEDIhofx=-inf GSIhofx=3.4840e+01 JEDIobserr=-3.3688e+38 GSIobserr=1.7704e-03 JEDIQC=15 SR_flag_JEDI=2 GSIeffQC=0 GSIQC=0 satid=750 imph=2439.50 Seq#=1478 obs=3.1054e+01 JEDIhofx=-inf GSIhofx=2.9244e+01 JEDIobserr=-3.3688e+38 GSIobserr=1.8467e-03 The reason for this discrepancy is that when GSI encounters ddnj>0, it discards the data, and this data is not considered when searching for the max obs value of the entire profile. On the other hand, JEDI includes this data when determining the max obs value. Since all data below this height of the max obs is discarded, this difference results in 3 additional data being tossed in JEDI.     2) GSIQC=8 while JEDIQC=0. There are 200+ data falling into this category. Several adjustments have been made to ufo_gnssro_bndnbam_mod.F90 to determine the accurate maximum obs value of the profile. These changes include removing 'super(k)>0' from the list of the conditions for SR2, separately searching max obs and its height, etc. These changes have resulted in reducing the number of inconsistency to 43 data points. 

[XuanliLi-NOAA]

A few more changes have been made to ufo_gnssro_bndnbam_mod.F90.  In JEDI, the SR2 QC searches within the entire profile for the maximum obs value with SR_flag = 0. All of the obs below the maximum will be rejected. In contrast, in GSI, the SR2 QC searches for the maximum obs that satisfies the criteria listed above regardless of SR_flag value. These differences caused inconsistencies. With the new changes in ufo_gnssro_bndnbam_mod.F90 , we now have consistent QC flags except for one data point (see the plot for sid=265). This discrepancy exists because JEDI and GSI use different methods to apply the cutoffs. In JEDI, the comparison is between O-HofX and cutoff*obs. While in GSI, it compares (O-HofX)/O with the cutoff value. For this specific data point, both O and the cutoff value are small which has resulted in the difference in truncation error.

The following are the plots of QC flags for each satellite:

[XuanliLi-NOAA]

Code changes are merged into ufo develop branch. More details about the changes in the operator, super refraction QC filter, and obs error assignment can be found in PR 2960 and PR 3040.