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.

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

YAML Configuration | Surface ship #77

Closed
3 tasks
danholdaway opened this issue May 12, 2022 · 54 comments
Closed
3 tasks

YAML Configuration | Surface ship #77

danholdaway opened this issue May 12, 2022 · 54 comments
Assignees

Comments

@danholdaway
Copy link
Contributor

danholdaway commented May 12, 2022

Agreement with GSI using GeoVaLs in UFO (Hofx)

  • Done

Agreement with GSI using GeoVaLs in UFO (Quality control)

  • Done

Agreement with GSI using Swell

  • Done

Current issues:

  • None

Notes

@gmao-jjin3 gmao-jjin3 self-assigned this Dec 28, 2022
@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Dec 28, 2022

corscat_omb_vs_gsiombbc

Ship data only.

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Dec 28, 2022

density_qc_omb_vs_gsiombbc

Ship data only.

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Dec 28, 2022

corscat_eferr_x_GsiDup_vs_finaler

Ship data only.

@gmao-jjin3
Copy link
Contributor

corscat_omb_vs_gsiombbc

Ship data + drifting buoys

@gmao-jjin3
Copy link
Contributor

density_qc_omb_vs_gsiombbc

Ship data + drifting buoys

@gmao-jjin3
Copy link
Contributor

corscat_eferr_x_GsiDup_vs_finaler

Ship data + drifting buoys

@danholdaway
Copy link
Contributor Author

Looking good @gmao-jjin3. Thanks for sharing the plots

@gmao-jjin3
Copy link
Contributor

density_qc_omb_vs_gsiombbc b
density_qc_omb_vs_gsiombbc
corscat_eferr_x_GsiDup_vs_finaler

@gmao-jjin3
Copy link
Contributor

density_qc_omb_vs_gsiombbc b
density_qc_omb_vs_gsiombbc
corscat_eferr_x_GsiDup_vs_finaler

@gmao-jjin3
Copy link
Contributor

Observational errors for airTemperature and virtaulTemperature in UFO don't match those values in GSI even after inflated errors made by UFO functions are replaced by GSI adjusted obs errors. This is because UFO does not have the additional error inflation as in GSI:
https://github.com/CoryMartin-NOAA/GSI/blob/ea04c69c4a3c29fc2109b6d7cdd21dfd049b1cf0/src/gsi/setupt.f90#L819
See JEDI team's explanation in #78

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Feb 13, 2023

density_qc_omb_vs_gsiombbc b
density_qc_omb_vs_gsiombbc
corscat_eferr_vs_finaler
corscat_eferr_x_GsiDup_vs_finaler

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Feb 13, 2023

density_qc_omb_vs_gsiombbc b
density_qc_omb_vs_gsiombbc
corscat_eferr_vs_finaler

corscat_eferr_x_GsiDup_vs_finaler

@gmao-jjin3
Copy link
Contributor

UFO air temperature and virtual temperature results match GSI values after using (1) updated pressure-check inflation function written by @gmao-wgu , and (2) GSI adjust errors which are configured in YAML files.

@gmao-wgu
Copy link
Contributor

@gmao-jjin3 great. I will remove gsiadjusterr from my code and use gsiadjusterr in the yaml instead. Thanks!

@gmao-jjin3
Copy link
Contributor

@gmao-wgu Not sure why I don't see your last post about updating testing environment:

You can copy, edit and run this bash script to change names after you build your JEDI checkout:
$NOBACKUP/jedi/jedi-work2/Wei_Gu/ioda_v3_temperature/rename_all_name_map_files.sh

The new testing data are in /discover/nobackup/jjin3/jedi/For_others/x0044.jj_20230201.iodav3/

BTW, I have updated your PressureCheck codes and you copy them over, and edit them before making a PR:
/discover/nobackup/jjin3/jedi/jedi-work2/Wei_Gu/ioda_v3_temperature/ObsErrorFactorPressureCheck.*

@gmao-wgu
Copy link
Contributor

gmao-wgu commented Feb 13, 2023

@gmao-jjin3 Thanks. I already did that, but I did not realize we have a script ready to convert the yaml into the one with the new naming convention though.

@gmao-wgu
Copy link
Contributor

@gmao-jjin3 By the way, can you cut your yaml to include only those necessary so that we can see clearly about how to handle gsiadjusterror in the yaml?

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Feb 13, 2023

@gmao-wgu Here it is in (/discover/nobackup/jjin3/jedi/jedi-work2/Wei_Gu/ioda_v3_temperature/sondes_geos_qc.yaml)


310   ## a tempporary solution: Replace error by GsiAdjustObsError
311   ## overwite above error assignments and inflation.
312   - filter: Perform Action
313     filter variables:
314     - name: airTemperature
315     action:
316       name: assign error
317       error function: GsiAdjustObsError/airTemperature
318     where:
319     - variable:
320         name: GsiAdjustObsError/airTemperature
321       is_defined:
322     defer to post: true

@gmao-wgu
Copy link
Contributor

@gmao-jjin3 Thanks. So I only need to replace airTemperature with windEastward, it that correct?

@gmao-jjin3
Copy link
Contributor

Yes, it is.

@gmao-jjin3
Copy link
Contributor

corscat_eferr_vs_finaler

@gmao-jjin3
Copy link
Contributor

corscat_eferr_vs_finaler

@gmao-jjin3
Copy link
Contributor

Duplicate factors are applied for above two figures. Need further check.

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Mar 1, 2023

The reason for why there are a few larger observational errors in UFO after duplication than GSI final errors is that location (longitudes values) in UFO are different from those in GSI because of machine errors. For example: Duplicate factor = 1 in UFO at (0N, 220.1E) :

lat[i],lon[i], pressure[i], datetime[i],ObsValue[i], GsiAdjustObsError[i],GsiFinalObsError[i],EffectiveError[i],                 
int(PreQC[i]),int(PreQCps[i]), Dupobs_Factor[i], i
0.0 220.1 100970.0 1607986800 297.35 1.75 1.7533125 2.4797208 2 2 1.0 11204
0.0 220.1 100970.0 1607986800 297.35 1.75 1.7533147 2.4797208 2 2 1.0 11471

However, they are at (0N, 220.1E) and (0N, 220.10001E) in GSI:

jjdebug1,mype,k,a.late=         329          17  0.000000000000000E+000
 jjdebug1,mype,k,a.ilat=         329          17   181.000000000000     
 jjdebug1,mype,k,b.lone=         329          17   **220.100000000000**     
 jjdebug1,mype,k,b.ilon=         329          17   353.160000000000     
 jjdebug1,mype,k,c.pre=         329          17   100970.000000000     
 jjdebug1,mype,k,d.time=         329          17   2.00000000000000     
 jjdebug1,mype,k,e.iqt=         329          17           1
 jjdebug1,mype,k,f.obs=         329          17   297.350000000000     
  jjdebug1,mype,k,a.late=         329          19  0.000000000000000E+000
 jjdebug1,mype,k,a.ilat=         329          19   181.000000000000     
 jjdebug1,mype,k,b.lone=         329          19   **220.100010000000**     
 jjdebug1,mype,k,b.ilon=         329          19   353.160016000000     
 jjdebug1,mype,k,c.pre=         329          19   100970.000000000     
 jjdebug1,mype,k,d.time=         329          19   2.00000000000000     
 jjdebug1,mype,k,e.iqt=         329          19           1
 jjdebug1,mype,k,f.obs=         329          19   297.350000000000     

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Mar 1, 2023

The reason for under-estimate of observation errors for Air Temperature in UFO than in GSI at quite a few locations is that some data which have the same locations, observation time, and values are treated as one observation in IODA and UFO. However, they are treated as different observations in GSI. As a result, there are larger duplicate factors in GSI than in UFO. Here is an example:
UFO sees tow observations at (30.71N, 271.96E) with a pressure level 102010 Pa.

lat[i],lon[i], pressure[i], datetime[i],ObsValue[i], GsiAdjustObsError[i],GsiFinalObsError[i],EffectiveError[i],                 int(PreQC[i]),int(PreQCps[i]), Dupobs_Factor[i], i, tmpXdup-GSIErr
30.71 271.96 102010.0 1607983200 283.55 1.75 7.5306387 6.148573 2 2 **1.7319707** 25470 -3.0994415e-05   
30.71 271.96 102010.0 1607983560 283.45 1.75 7.470293 6.099872 2 2 **1.7318904** 25471 0.0003209114 

However, there are three obervations in GSI (duplicate factor is about 1.73 which is sqrt(3):

jjdebug2,mype,k,station_id=         444        2797  OBLA1   
jjdebug2,mype,k,a.late=         444        2797   30.7100000000000     
jjdebug2,mype,k,a.ilat=         444        2797   242.420000000000     
jjdebug2,mype,k,b.lone=         444        2797   271.960000000000     
jjdebug2,mype,k,b.ilon=         444        2797   436.136000000000     
jjdebug2,mype,k,c.pre=         444        2797   102010.000000000     
jjdebug2,mype,k,d.time=         444        2797   1.00000000000000     
jjdebug2,mype,k,e.iqt=         444        2797           1
jjdebug2,mype,k,f.obs=         444        2797   283.550000000000     

jjdebug2,mype,k,station_id=         444        2798  OBLA1   
jjdebug2,mype,k,a.late=         444        2798   30.7100000000000     
jjdebug2,mype,k,a.ilat=         444        2798   242.420000000000     
jjdebug2,mype,k,b.lone=         444        2798   271.960000000000     
jjdebug2,mype,k,b.ilon=         444        2798   436.136000000000     
jjdebug2,mype,k,c.pre=         444        2798   102010.000000000     
jjdebug2,mype,k,d.time=         444        2798   1.10000002384186     
jjdebug2,mype,k,e.iqt=         444        2798           1
jjdebug2,mype,k,f.obs=         444        2798   283.450000000000     

jjdebug2,mype,k,station_id=         444        3204  OBLA1   
jjdebug2,mype,k,a.late=         444        3204   30.7100000000000     
jjdebug2,mype,k,a.ilat=         444        3204   242.420000000000     
jjdebug2,mype,k,b.lone=         444        3204   271.960000000000     
jjdebug2,mype,k,b.ilon=         444        3204   436.136000000000     
jjdebug2,mype,k,c.pre=         444        3204   102010.000000000     
jjdebug2,mype,k,d.time=         444        3204   1.00000000000000     
jjdebug2,mype,k,e.iqt=         444        3204           1
jjdebug2,mype,k,f.obs=         444        3204   283.550000000000     

@gmao-jjin3
Copy link
Contributor

corscat_hofx_vs_gsihofx

It is strange that Hofx values without any QC don't match.

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Mar 6, 2023

It is strange that Hofx values without any QC don't match.

Okay, it is because that GSI hofx "qges" and saturated moisture "qsges" are values at the 1st model level for marine surface specific humidity! The grid number "dpres" is reset to 1
https://github.com/GEOS-ESM/GEOSana_GridComp/blob/develop/GEOSaana_GridComp/GSI_GridComp/setupq.f90#L528
As a result, these values are not interpolated from the reported pressure levels.
CVS version history shows this change was added in NCEP's August 2005 update. @gmao-msienkie and @gmao-yzhu , do you know the reason for it? Thanks.

@gmao-msienkie
Copy link
Contributor

I don't know why they did it, but it seems to me to make sense to place the surface marine observations at the surface level of the model. The observations are taken at the surface of the ocean (though I guess for ships it could be slightly higher than the surface) and there is no issue with differing model topography vs real world topography.
We used to carry files with a log of changes to the early GSI in a 'Changes' directory but have dropped that in recent versions. There is a copy in /home/dao_ops/GEOSadas-5_12_4_SLES12/GEOSadas/src/GEOSgcs_GridComp/GEOSana_GridComp/GEOSaana_GridComp/GSI_GridComp/Changes but I didn't see anything that looked like it related to this unless it was related to adding a 2D-var surface analysis option (20050601.txt).

@gmao-jjin3
Copy link
Contributor

Thank you @gmao-msienkie . It makes sense to me too, but have to think about how to do that in UFO.

@CoryMartin-NOAA
Copy link

It's possible that this could be related to type 183 see: https://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/table_2.htm which is "SURFACE MARINE (SHIP, BUOY, C-MAN, TIDE GAUGE) OR LAND [SYNOPTIC (fixed and mobile), METAR] WITH MISSING STATION PRESSURE - Tv, q, Pstn, sst"

@gmao-jjin3
Copy link
Contributor

Right. It looks like that. Thanks @CoryMartin-NOAA !

@gmao-jjin3
Copy link
Contributor

corscat_hofx_vs_gsihofx
@CoryMartin-NOAA It looks that I get values at the top of the model which are extremely small. Here is my configuration:

  obs operator:
    name: Composite
    components:
     - name: VertInterp
       observation alias file: /gpfsm/dnb31/jjin3/jedi/JediUfoTests_iodav3/Config/obsop_name_map.yaml
       apply near surface wind scaling: true
       variables:
       - name: airTemperature
       - name: virtualTemperature
       - name: windEastward
       - name: windNorthward
     - name: SfcPCorrected
       variables:
       - name: stationPressure
       da_psfc_scheme: GSI
       geovar_geomz: geopotential_height
       geovar_sfc_geomz: surface_geometric_height
       station_altitude: height
     - name: Identity
       observation alias file: /gpfsm/dnb31/jjin3/jedi/JediUfoTests_iodav3/Config/obsop_name_map.yaml
       variables:
       - name: specificHumidity

Any idea?

@CoryMartin-NOAA
Copy link

@gmao-jjin3 you may need this: https://github.com/JCSDA-internal/ufo/blob/785cc133070e37860425e96e8605b26e3e562c1c/src/ufo/operators/identity/ObsIdentityParameters.h#L37
This decides if identity grabs the first or last vertical value for the interpolation.

@gmao-jjin3
Copy link
Contributor

A PR to make GeoVals top_down during reading the geoval file if the levels are not top_down in it.
https://github.com/JCSDA-internal/ufo/pull/2708

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented Mar 22, 2023

corscat_eferr_vs_finaler
corscat_eferr_vs_finaler

Air temperature and virtual temperature errors in UFO match those in GSI after additional changes in GSI: 1) discard observations which have the same location, pressure, and time; 2) determine duplicated observations by single precision location an pressure values; 3) not to consider co-located air temperature and virtual temperature as duplicated observations.

@asewnath
Copy link
Collaborator

asewnath commented May 2, 2023

For surface marine moisture, the obserr plot was not perfectly matching, and the conclusion was that we were having challenges with the saturated specific humidity calculation. To prove this, we ran a test where we provided UFO with the saturated specific humidity value from the GSI for a given obs location. The obserr plots were perfectly matched. Then, we provided the saturated specific humidity profiles in the geoval file such that the only calculation needed for UFO's saturated specific humidity at a given obs location was an interpolation. The resulting obserr plot shown below indicates a discrepancy with the interpolation.
Screenshot 2023-05-02 at 10 15 53 AM
The error difference (UFO error - GSI error) is illustrated in the plot below:
Screenshot 2023-05-02 at 10 16 04 AM
Zoomed in
Screenshot 2023-05-02 at 10 27 48 AM

The number of error difference values that are less than -1e-7 is 39.

This seems to corresponds to the calculation of saturated specific humidity. There are 39 locations where the interpolated saturated specific humidity is not equal to the one provided at the obs location.

@asewnath
Copy link
Collaborator

asewnath commented May 4, 2023

On further inspection of the 39 locations where the GSI-provided saturated specific humidity for a given location did not equal the UFO interpolated saturated specific humidity, we've found that the pressure at the obs location is between the first two model levels. The other locations where the two values for saturated specific humidity are equal have the pressure at obs location outside of the pressure model levels.

The following is a print out of the first two levels of the saturated specific humidity profile, the GSI-provided saturated specific humidity for a given location, and the UFO interpolated saturated specific humidity. These values are printed out when a difference between the two saturated specific humidity values occur:

q_profile[0]: 0.00446344
q_profile[1]: 0.00415493
GSI provided sat_specific_humidity[iloc]: 0.00446344
interpolated sat specific humidity: 0.004214

q_profile[0]: 0.00446462
q_profile[1]: 0.00415794
GSI provided sat_specific_humidity[iloc]: 0.00446462
interpolated sat specific humidity: 0.0042152

q_profile[0]: 0.00446581
q_profile[1]: 0.00416095
GSI provided sat_specific_humidity[iloc]: 0.00446581
interpolated sat specific humidity: 0.00421848

In all these locations, the GSI seems to pick the first saturated specific humidity level as the saturated specific humidity for the given location, whereas UFO tries to interpolate it. If saturated specific humidity is set to the first level of the saturated specific humidity profile, then the obserror perfectly matches the GSI for sfcship.

Screenshot 2023-05-04 at 11 05 29 AM

@gmao-jjin3
Copy link
Contributor

gmao-jjin3 commented May 12, 2023

There is a good agreement between the UFO effective errors and GSI final errors for marine ships wind data when the same "obs" pressure values are used for duplicate factors. In GSI, model (BKG) pressure values at the center time within a time window are used for duplicate factor.
corscat_eferr_vs_finaler

@gmao-jjin3
Copy link
Contributor

It turns out that GSI HofX for marine ships wind data cannot be reproduced by UFO even they are interpolated in a geometric coordinate in UFO:


       apply near surface wind scaling: true
       variables:
       - name: windEastward
       - name: windNorthward
       interpolation method: linear
       vertical coordinate: geometric_height
       observation vertical coordinate: height


The likely reason is that observation heights for some of wind data are manipulated between interpolation is conducted. See code between lines 588 - 600
https://github.com/GEOS-ESM/GEOSana_GridComp/blob/develop/GEOSaana_GridComp/GSI_GridComp/setupw.f90#L600
Not sure there is a way to make this kind of twist to observation heights before HofX is calculated in UFO. @CoryMartin-NOAA Did you or anyone else met this issue before? Thanks.

corscat_hofx_vs_gsihofx

@asewnath
Copy link
Collaborator

GSI HofX is reproduced in UFO when the GSISfcModel obs operator is used for windEastward and windNorthward

@CoryMartin-NOAA
Copy link

@asewnath that operator does not yet have a TL/AD included in UFO for it, and I was planning on removing it... Is this a case where the VertInterp can reproduce GSI if you provide it the 10m wind reduction factor? Or does that need to be computed by UFO?

@danholdaway
Copy link
Contributor Author

Thanks for flagging that @CoryMartin-NOAA. That would indeed by an issue. I don't think the wind scaling factor is sufficient. Have you been able to reproduce the winds from surface marine using the regular vertical interpolation operator?

@CoryMartin-NOAA
Copy link

I thought that I had done so back in the fall with just a combination of 10m wind reduction factor and the VertInterp. Perhaps we use GSI slightly differently for these obs? I would still think it would be a good idea to remove that GsiSfcModel operator, so perhaps we can code up a variable transform/obs function to do what we need to scale the winds?

@danholdaway
Copy link
Contributor Author

With adding the surface wind scaling to the Identity operator in UFO (https://github.com/JCSDA-internal/ufo/pull/2853) we can reproduce the wind h(x).

corscat_hofx_vs_gsihofx

density_qc_omb_vs_gsiomb

corscat_eferr_vs_finaler

Even though the counts are the same it looks like different observations might be passing QC as the stats in the density plots reveal very small differences.

@asewnath @gmao-jjin3

@danholdaway
Copy link
Contributor Author

Never mind, I misread the plots, there is a count difference of 10. If you want to check the YAML to see if all filters are active it is here: /discover/nobackup/drholdaw/JediUfoTests/Config/sfcship_geos_qc.yaml

@danholdaway
Copy link
Contributor Author

corscat_hofxdiff_vs_gsi
corscat_hofxdiffnorm_vs_gsihofx

@CoryMartin-NOAA
Copy link

Looks like a TIE fighter 😄 :
image

@danholdaway
Copy link
Contributor Author

@CoryMartin-NOAA this means this code is ready for JEDI©!

@gmao-jjin3
Copy link
Contributor

@danholdaway @asewnath @gmao-wgu
In regards to these 10 mismatched observations, they failed gross-error check in GSI: sqrt(U_omf^2 + V_omf^2) / observational_error > cgross_benchmark which is 6.
Here is the line in the GSI code
https://github.com/GEOS-ESM/GEOSana_GridComp/blob/develop/GEOSaana_GridComp/GSI_GridComp/setupw.f90#L1113

Here is my search to these 10 data in the GSI ncdiag output file
/discover/nobackup/jjin3/jedi/For_others/conv_obs_check/ioda_satwind/gsi_use_flag.py.log
You can see the above ration_final is larger than qcgross_final (=6):
14:qcgross_final, ratio_final= 6.0 7.2893376
25:qcgross_final, ratio_final= 6.0 6.6631603
36:qcgross_final, ratio_final= 6.0 6.8705525
47:qcgross_final, ratio_final= 6.0 10.145198
58:qcgross_final, ratio_final= 6.0 7.6063986
69:qcgross_final, ratio_final= 6.0 6.6306386
80:qcgross_final, ratio_final= 6.0 6.5806003
91:qcgross_final, ratio_final= 6.0 6.0000467
102:qcgross_final, ratio_final= 6.0 6.6451516
113:qcgross_final, ratio_final= 6.0 7.1501427

I don't see this kind of gross error check in the
/discover/nobackup/drholdaw/JediUfoTests/Config/sfcship_geos_qc.yaml

@asewnath
Copy link
Collaborator

Thanks for looking into this Jianjun! Adding Wei's gross check to the sfcship yaml took care of the remaining 10 obs.

@gmao-jjin3
Copy link
Contributor

Validation against GEOS x0048 results (/discover/nobackup/jjin3/jedi/For_others/x0048.jj_20230818/conv/combine_conv/obs/sfcship_obs_2021121200.nc4).

density_qc_omb_vs_gsiomb
corscat_qc_omb_vs_hofx-GsiHofXBc
corscat_EffectiveQC_vs_GsiEffectiveQC
corscat_eferr_x_GsiDup_vs_finaler
corscat_eferr_vs_finalerr

@gmao-jjin3
Copy link
Contributor

density_qc_omb_vs_gsiomb
corscat_qc_omb_vs_hofx-GsiHofXBc
corscat_EffectiveQC_vs_GsiEffectiveQC
corscat_eferr_x_GsiDup_vs_finaler
corscat_eferr_vs_finalerr

@Arcayou
Copy link

Arcayou commented Apr 18, 2024

It's possible that this could be related to type 183 see: https://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/table_2.htm which is "SURFACE MARINE (SHIP, BUOY, C-MAN, TIDE GAUGE) OR LAND [SYNOPTIC (fixed and mobile), METAR] WITH MISSING STATION PRESSURE - Tv, q, Pstn, sst"

This website cannot be opend. Is there any new website viewing this table ?

@CoryMartin-NOAA
Copy link

We had a data center meltdown and non-operational things such as the EMC website are currently down. I hope its back sometime next week.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

7 participants