Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Incorporate Sondes DA into JEDI for HAFS #16

Open
JingCheng-NOAA opened this issue Sep 9, 2024 · 31 comments
Open

Incorporate Sondes DA into JEDI for HAFS #16

JingCheng-NOAA opened this issue Sep 9, 2024 · 31 comments

Comments

@JingCheng-NOAA
Copy link
Collaborator

JingCheng-NOAA commented Sep 9, 2024 •
edited
Loading

As the single ob DA is done within JEDI for HAFS, this issue is opened for more sounding obs DA in JEDI and comparison with GSI.
This issue will record the process for IODA file generation, JEDI yaml file for upper air DA, validation of Hx.
@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Sep 9, 2024 •
edited
Loading

Start with the airTemperature with otype=120. Using the 2024 02L Beryl test case.

First split ADPUPA from gfs prepbufr file and edit the gsiparm.anl to take in OBS t only.
Then convert the GSI diag file into IODA NC file to provide observations for JEDI analysis.
The conversion produced two observation files:
sondes_tsen_obs_2024063012.nc4
sondes_tv_obs_2024063012.nc4

Run JEDI analysis with observation tsen, using VertInterp.

       obs operator:
         name: VertInterp
         variables:
         - name: airTemperature
       obs error:
         covariance model: diagonal

fv3jedi_vs_gsi_validation_120_airTemperature
fv3jedi_vs_gsi_OMB_120_t

The hofx of the two are highly agree to each other. However the observation error in JEDI are smaller comparing with GSI

@JingCheng-NOAA
Copy link
Collaborator Author

I've assigned the error as what HAFS GSI does with errtable, but the results of final ob error still doesn't look close to GSI.

           action:
             name: assign error
             error function:
               name: ObsFunction/ObsErrorModelStepwiseLinear
               options:
                 xvar:
                   name: MetaData/pressure
                 xvals: [110000, 105000, 100000, 95000, 90000, 85000, 80000, 75000, 70000, 65000, 60000, 55000, 50000, 45000, 40000, 35000, 30000, 25000, 20000, 15000, 10000, 7500, 5000, 4000, 3000, 2000, 1000, 500, 400, 300, 200, 100, 0]
                 errors: [1.2696, 1.3282, 1.3932, 1.439, 1.4354, 1.3669, 1.2552, 1.1362, 1.0397, 0.98016, 0.94757, 0.92875, 0.91714, 0.91625, 0.93506, 0.98447, 1.0699, 1.1816, 1.2862, 1.3542, 1.3831, 1.3832, 1.3534, 1.3109, 1.3018, 1.3414, 1.4017, 1.4471, 1.4744, 1.4892, 1.4964, 1.4988, 1.4962]

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Sep 19, 2024 •
edited
Loading

Adding observation inflation using ObsFunction/ObsErrorFactorPressureCheck, but the final observation error of JEDI is unchanged as not use. Guess need to find a better way to inflate the observation error.
added section:

         - filter: Perform Action
           filter variables:
           - name: airTemperature
           where:
           - variable: ObsType/airTemperature
             is_in: 120
           action:
             name: inflate error
             inflation variable:
               name: ObsFunction/ObsErrorFactorPressureCheck
               options:
                 variable: airTemperature
                 inflation factor: 8.0
                 #geovar_sfc_geomz: surface_geometric_height
           defer to post: true

@JingCheng-NOAA
Copy link
Collaborator Author

Test on specificHumidity
fv3jedi_vs_gsi_validation_120_specificHumidity

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Sep 24, 2024 •
edited
Loading

Test with uv. Here for the wind observations, introducing hofx scaling into the observation operator.
Use GSI geovals will get identical hofx.
GSIgeovals_hofx_220_windEastward

if not use GSI geovals, then there are slight difference in JEDI and GSI hofx. The hofx scaling also did some correction here.
jedi_gsi_hofx_220_windEastward

HOFX after scaling - HOFX no scaling
fv3jedi_hofxscaling_vert_validation_220_windEastward

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Sep 30, 2024 •
edited
Loading

Rerun the GSI with grid_ratio_fv3_regional=1, but still find the difference in hofx is large.
Plot out the difference of hofx and found the large mismatch on the top of model level.

hofxhist_windEastward
hofxdiff_vs_pressure_windEastward
hofxhist_windNorthward
hofxdiff_vs_pressure_windNorthward

@JingCheng-NOAA
Copy link
Collaborator Author

Though the IODA NC file is created from GSI diag file, but it tends out several points are out of JEDI domain boundary. Using the RRFS domain filter tool to further filter out points outside or near the boundary.
domain_check_FV3

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Oct 2, 2024 •
edited
Loading

After applying the domain check, the hofx difference between jedi and gsi are smaller. But the large difference near top of the model remains.
domainfilter_fv3jedi_vs_gsi_validation_220_windEastward
domainfilter_fv3jedi_vs_gsi_validation_220_windNorthward
domainfilter_hofxhist_windEastward
domainfilter_hofxdiff_vs_pressure_windEastward
domainfilter_hofxhist_windNorthward
domainfilter_hofxdiff_vs_pressure_windNorthward

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Oct 3, 2024 •
edited
Loading

Select the point where has max hofx difference, the location is lat=4.75 lon=285.80
FROM JEDI DIAG:
obs index 1431, location: 4.75 285.80,
pressure: 8419.999694824219,
obsvalue:-10.199999809265137
jhofx: -12.158164978027344, ghofx: -13.30904769897461

FROM GSI DIAG:
matched obs index 587
location: 4.754079818725586 285.79833984375,
eastward_wind -13.183977127075195
pressure level: 8431.2216796875
u_obs:-10.199999809265137

FROM ncdiag file converted IODA nc:
location: 4.754079818725586 285.79833984375
obs value -10.199999809265137
FROM bufr2ioda converted IODA nc:
location: 4.754079818725586 -74.20166015625
obs value -10.199999809265137

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Oct 7, 2024 •
edited
Loading

Check the gsi and jedi geoval files to see the "air_pressure" distribution difference. GSI geofiles contains lat and lon information, so it's easier to locate the obs location that located near "4.75 285.80".
For JEDI, all the files do not contain any lat/lon information, the only way to locate is using the index number the same as we found in GSI geofiles, the results are quite different too

geoval_airpressure

From GSI
surface pressure 75540.2734375
air pressure [ 325. 600. 957.5685 1440.7946 2017.909 2646.323
3325.437 4054.6511 4833.365 5660.979 6536.893 7460.5073
8431.222 9448.436 10494.844 11540.1 12575.047 13607.0205
14641.574 15682.952 16734.43 17798.547 18877.28 19972.168
21084.408 22214.926 23364.422 24533.428 25722.33 26931.395
28160.795 29410.621 30680.898 31971.592 33282.613 34613.848
35965.13 37336.27 38727.06 40137.258 41566.61 43013.66
44474.168 45942.13 47411.387 48875.785 50329.164 51765.37
53178.25 54561.9 55910.97 57220.7 58486.99 59706.445
60876.363 61994.676 63059.98 64071.48 65028.94 65932.64
66783.3 67581.99 68330.1 69029.31 69681.47 70288.586
70852.86 71376.42 71861.51 72310.33 72725.13 73108.086
73461.21 73786.61 74086.26 74361.95 74615.46 74848.445
75062.484 75259.05 75446.484 ]

For JEDI,
surface pressure [96645.22]
air pressure [ 319.10495 595.48126 952.16327 1434.4716 2012.5248 2641.5
3320.988 4050.47 4829.386 5657.159 6533.2026 7456.9243
8427.73 9445.024 10505.213 11602.284 12734.048 13901.352
15104.716 16344.424 17620.58 18933.156 20282.023 21666.967
23087.717 24543.947 26035.299 27561.373 29121.754 30715.994
32343.637 34004.207 35697.22 37422.17 39178.56 40965.87
42783.586 44631.184 46508.13 48413.887 50347.92 52308.094
54288.496 56280.83 58276.555 60267.156 62244.094 64198.84
66122.875 68008.016 69846.836 71632.73 73359.984 75023.88
76620.63 78147.336 79602. 80983.48 82291.38 83526.05
84688.414 85779.914 86802.414 87758.18 88649.71 89479.75
90251.266 90967.18 91630.51 92244.29 92811.58 93335.336
93818.32 94263.39 94673.27 95050.41 95397.19 95715.91
96008.71 96277.62 96525.76 ]

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Oct 8, 2024 •
edited
Loading

The index match method may not correct. So check the 1st point in both GSI and JEDI GeoVals.
From GSI for idx 0
surface pressure 101772.7421875
From JEDI for idx 0
surface pressure 101786.336
geoval_airpressure_0point

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Oct 9, 2024 •
edited
Loading

Found a bug in finding the matched point index in JEDI GEOVaLs.
And @delippi mentioned about the difference in air_pressure calculation recipes in GSI and JEDI https://github.com/NOAA-EMC/RDASApp/issues/54#issuecomment-2150339798.
So modified the code to use the same middle point of pressure level and it did change the first 15 levels of air_pressure values.
geoval_airpressure
geoval_GEOVAL_uwind
geoval_GEOVAL_vwind

From gsi geoval file nc:
idx 587
location 4.754079818725586 285.79833984375
surface pressure 75540.2734375
air pressure [ 325. 600. 957.5685 1440.7946 2017.909 2646.323
3325.437 4054.6511 4833.365 5660.979 6536.893 7460.5073
8431.222 9448.436 10494.844 11540.1 12575.047 13607.0205
14641.574 15682.952 16734.43 17798.547 18877.28 19972.168
21084.408 22214.926 23364.422 24533.428 25722.33 26931.395
28160.795 29410.621 30680.898 31971.592 33282.613 34613.848
35965.13 37336.27 38727.06 40137.258 41566.61 43013.66
44474.168 45942.13 47411.387 48875.785 50329.164 51765.37
53178.25 54561.9 55910.97 57220.7 58486.99 59706.445
60876.363 61994.676 63059.98 64071.48 65028.94 65932.64
66783.3 67581.99 68330.1 69029.31 69681.47 70288.586
70852.86 71376.42 71861.51 72310.33 72725.13 73108.086
73461.21 73786.61 74086.26 74361.95 74615.46 74848.445
75062.484 75259.05 75446.484 ]
From jedi geoval file nc:
idx 354
location 28.25101089477539 279.1451416015625
surface pressure [75602.24]
air pressure [ 325. 600. 957.5685 1440.7946 2017.909 2646.323
3325.437 4054.6511 4833.365 5660.979 6536.893 7460.5073
8431.222 9448.436 10494.884 11540.291 12575.522 13607.895
14642.944 15684.905 16737.045 17801.893 18881.422 19977.166
21090.314 22221.791 23372.297 24542.355 25732.355 26942.557
28173.13 29424.17 30695.695 31987.672 33300.01 34632.594
35985.254 37357.81 38750.035 40161.7 41592.543 43041.113
44503.16 45972.67 47443.49 48909.44 50364.37 51802.113
53216.504 54601.64 55952.16 57263.297 58530.95 59751.727
60922.906 62042.426 63108.887 64121.484 65079.98 65984.66
66836.234 67635.79 68384.71 69084.69 69737.55 70345.33
70910.21 71434.34 71919.95 72369.26 72784.516 73167.88
73521.4 73847.15 74147.12 74423.12 74676.9 74910.14
75124.414 75321.195 75508.805 ]

At least we know the difference in hofx shouldn't be coming from air_pressure.
hofxdiff_vs_pressure_windEastward_modifiedReceipe

@delippi
Copy link

delippi commented Oct 9, 2024

@JingCheng-NOAA, Thanks for testing this. It looks like you are using GSI-IODA but not using GSI-geovals and just letting JEDI compute its own geovals, is that correct? Did you say that you get identical results when using the GSI-geovals? Can you add the the same pressure vs hofx difference plot for when using GSI-geovals?

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Oct 9, 2024 •
edited
Loading

@JingCheng-NOAA, Thanks for testing this. It looks like you are using GSI-IODA but not using GSI-geovals and just letting JEDI compute its own geovals, is that correct? Did you say that you get identical results when using the GSI-geovals? Can you add the the same pressure vs hofx difference plot for when using GSI-geovals?

Yes the results of using same GSI GeoVaLs are showed in https://github.com/NOAA-EMC/HDASApp/issues/16#issuecomment-2372193384
GSIgeovals_hofx_220_windEastward

And if plot the difference in vertical distribution:
GSIGeoVals_hofxdiff_vs_pressure_windEastward

@delippi
Copy link

delippi commented Oct 9, 2024

@JingCheng-NOAA, interesting. I would expect that if you get identical results with GSI-geovals then you would have identical results at levels 1-15 and still have differences elsewhere. But you might be using gph as vertical coordinate or something else instead of pressure. Check what is used in GSI.

@JingCheng-NOAA
Copy link
Collaborator Author

@JingCheng-NOAA, interesting. I would expect that if you get identical results with GSI-geovals then you would have identical results at levels 1-15 and still have differences elsewhere. But you might be using gph as vertical coordinate or something else instead of pressure. Check what is used in GSI.

I checked the setupw.f90 in GSI and for this type of observations (otype=220), the GSI is using pressure levels.

@JingCheng-NOAA
Copy link
Collaborator Author

An update was made in comparison GSI test, block the otype 232 and only keep the otype 220.
Without domain filter, the results are updated here
fv3jedi_vs_gsi_validation_220_windEastward
fv3jedi_vs_gsi_validation_220_windNorthward

@JingCheng-NOAA
Copy link
Collaborator Author

fv3jedi_vs_gsi_validation_120_specificHumidity

@JingCheng-NOAA
Copy link
Collaborator Author

fv3jedi_vs_gsi_validation_120_stationPressure

@JingCheng-NOAA
Copy link
Collaborator Author

fv3jedi_vs_gsi_validation_120_airTemperature

@JingCheng-NOAA
Copy link
Collaborator Author

fv3jedi_vs_gsi_validation_120_virtualTemperature

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Oct 19, 2024 •
edited
Loading

I've noticed this discussion in GDASApp, which compared the difference in Cube and Gaussian grids. This may explain our issue found in the difference of geovals in GSI and JEDI.
https://github.com/NOAA-EMC/GDASApp/issues/750#issuecomment-2079617045
ps, similar difference in sounding data was found in GDAS test too
https://github.com/NOAA-EMC/GDASApp/issues/750#issuecomment-2077646953

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Oct 30, 2024 •
edited
Loading

Move to phase 2 for sounding data.
Start with prepareing the IODA nc using the bufr2ioda in GDASApp for wind observations (otype=220).
Meanwhile, block other observation types of data and only keep otype=220 in GSI parm file.
Observation count difference in GSI and JEDI diag file. Need to find out which data were rejected by GSI but accepted by JEDI
GSI: 1456
JEDI: 1614

@JingCheng-NOAA
Copy link
Collaborator Author

Modified the bufr2ioda python script to add in MetaData/timeOffset. This variable will be used to QC the observation that fall out of time window

@JingCheng-NOAA
Copy link
Collaborator Author

Regenerated IODA nc file with updated bufr2ioda yaml file.
Added in the on-line domain filter with updated HDASApp submodules.
The number used by JEDI is reduced a lot comparing with GSI.
obsnumber

Here are some statistics:
*JEDI without any obs filter:

  • Number used in analysis: 1150
  • Missing value prevent use of observation: 30510
  • hofx computation failed: 37399

*JEDI with obs filter applied:

  • Number used in analysis: 1138
  • Missing value prevent use of observation: 30510
  • Observation value out of bounds: 35327
  • Observation not within domain: 878
  • Observation removed due to thinning: 1206

@JingCheng-NOAA
Copy link
Collaborator Author

The observation distribution plot illustrated where the missing observations in JEDI are.
image

@delippi delippi mentioned this issue Nov 22, 2024
Closed
@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Dec 4, 2024 •
edited
Loading

Updated IODA yaml file with correct time and recreated observation file shows correct observation numbers that similar to GSI.
Also fixed a bug in python script to correctly calculate the wind increment of GSI.
diag_wind_220_hofxobs_on_domain

@JingCheng-NOAA
Copy link
Collaborator Author

Check for airTemperature and specificHumidity
diag_t_120_hofxobs_on_domain
diag_q_120_hofxobs_on_domain

@JingCheng-NOAA
Copy link
Collaborator Author

Thanks to @delippi 's plotting scripts, I've checked the stats after QC for each variable. Some of the observation error still has difference.
windstats
airTemperatureStats
specificHumidityStats

@JingCheng-NOAA
Copy link
Collaborator Author

An update was made in Wind yaml file, including adjusting the gross error and a bug fix in ObsErrorFactorConventional. Meanwhile, rerun the GSI with assimilate wind 220 only (previously was assimilating every type of observations).

image

@JingCheng-NOAA
Copy link
Collaborator Author

JingCheng-NOAA commented Dec 16, 2024 •
edited
Loading

For specific Humidity, adjusted the gross error min/max in JEDI yaml; And fixed a bug in the initial error assignment in JEDI yaml (the value grabbed from GSI errtable should be divided by 10 for q)
image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@delippi @JingCheng-NOAA