GEOS-ESM / swell

Workflow system for coupled data assimilation applications
https://geos-esm.github.io/swell/
Apache License 2.0
15 stars 4 forks source link

YAML Configuration | CRIS #65

Closed danholdaway closed 1 year ago

danholdaway commented 2 years ago

Agreement with GSI using GeoVaLs in UFO (Hofx no bias correction)

Agreement with GSI using GeoVaLs in UFO (Bias correction)

Agreement with GSI using GeoVaLs in UFO (Quality control)

Agreement with GSI using GeoVaLs in UFO (Hofx with bias correction)

Agreement with GSI using Swell

Pull requests in JEDI

Current issues:

gmao-wgu commented 2 years ago

For CrIS-fsr_NPP, the UFO has the good agreement with the GSI in terms of the number of observations that can pass the QC, OMF with bias correction and the effective observation errors. Here, the window channel of 962.500cm-1(10.39micron) is selected to show the agreement:

cris-fsr_npp Swell_Obs_Number O-F_BC_QC ch219 ps

cris-fsr_npp Swell_ObsValue-hofx ObsValue-GsiHofXBc_BC_QC ch219 ps

cris-fsr_npp Swell_EffectiveError GsiFinalObsError_QC ch219 ps

gmao-wgu commented 2 years ago

So is the CrIS-fsr_n20, the consistency can also be achieved between the GSI and UFO.

danholdaway commented 2 years ago

Great stuff @gmao-wgu. Thanks!

gmao-wgu commented 1 year ago

CRIS-FSR is being revisited with the aim of finalizing its YAML fil to ensure its compatibility with x0048 for usage in UFO. Changes were made in the YAML to reflect the reconfigured cloud detection and the updated gross check adopted in x0048.

gmao-wgu commented 1 year ago

@danholdaway @gmao-jjin3 The following plots indicate that the number of observations passing QC from GSI for channel 440 or 200 (using different subset indices in Eva and IDL) is 297. However, this count is not consistent with the number (298) reported by both ana_stats and diag for this channel. Should we investigate this issue further?

gsi_omb_jedi_omb_density_cris-fsr_n20_brightnessTemperature_440 cris-fsr_n20 UFO_Obs_Number O-F_BC_QC ch200 ps

gmao-wgu commented 1 year ago

@gmao-jjin3 The count for Cris_FSR_n20 channel 440 from GSI can be found on the line of 14025 in the file of x0048.ana_stats.log.20211212_00z.txt under /gpfsm/dnb05/projects/p139/dao_it/archive/x0048/etc/Y2021/M12

gmao-jjin3 commented 1 year ago

@gmao-wgu I just recall that that QC flags don't correctly show whether IR data are actively assimilated. "varinv" are set to be zero but new QC flags values are NOT assigned for many IR observations which fail QC in setuprad.f90. Now I suspect that "ana_stats" table is correct for this reason. GSI final observational error values are used to set GSI QC flags in both IODA and my IDL programs. I am not sure whether you checked those GSI final observational values in the binary file. If not, can you use those values to determine if IR data pass or fail QC?

gmao-wgu commented 1 year ago

@gmao-jjin3 @danholdaway I reviewed the notes I made for AIRS and qcmod for IR in GSI, it seemed to me the only filter being assigned a non zero QC flag when doing the error inflation is the short wave check for IR. so it should be reliable to utilize the QC flag as a determining factor for whether an observation is assimilated in GSI, except for AIRS short-wave channels. We already reached a conclusion utilizing observation error as an indicator is not feasible, particularly for conventional data, due to the presence of non-zero observation errors in GSI for passive data. Additionally, since the error value is a real number, it is not reliable to use "error = 0" as a condition.

gmao-jjin3 commented 1 year ago

@gmao-wgu I think errinv is a reliable QC index in GSI for all radiance and conventional observations. It is zero for an observation that fails in QC and is a negative value for passive observations.

gmao-wgu commented 1 year ago

@gmao-jjin3 It has been already tested that errinv is not a reliable index for conventional data from GSI, For conventional data from GSI, only GsiUseFlag is reliable.

gmao-wgu commented 1 year ago

@gmao-jjin3 For radiance data, GsiEffectiveQC should generate the identical number as in GSI for all except AIRS shortwave channels. For cirs-fsr, we do not have short-wave channels assimilated, the difference in counts we saw between GSI and GSI in IODA-convertor might be associated with "errinv" you use in the IODA-convetor.

gmao-wgu commented 1 year ago

@gmao-jjin3 The plots I posted yesterday already indicate errinv is not reliable to use for radiance usage from GSI.

gmao-jjin3 commented 1 year ago

Okay. "GsiEffectiveQC" is made according to "errinv" IODA converter. We definitely don't assimilate an observation if its errinv = 0 which means error R is a huge value and weight 1/R = 0. I really don't see a case in GSI that we assign a qc_flag but don't set errinv = 0. On the other hand, we actually assimilate an observation If its errinv is a positive value.

gmao-wgu commented 1 year ago

If that is the case, then why we see the difference in the count between the one directly from GSI and the one from the IODA converter? Honestly, I trust the one we got directly from both GSI diag binary and ana_stats for all except short-wave channels.

gmao-wgu commented 1 year ago

GsiEffectiveQc was used to store id_qc in the version 2 of IODA-converter, then somehow it was switched to hold "errinv" if I understood correctly, that is the reason when we plot GsiEffectiveQc, it only has the value of either 0 or 1. Apparently, that is not the case in the version 2 of IODA converter.

gmao-jjin3 commented 1 year ago

Wei, Here is how you check the GSI binary output: _when I read the binary file, I use "lqcpass = luse(ich) .and. datachan(ich)%qcmark .eq. 0 .and. datachan(ich)%errinv .gt. tiny_single" to determine whether the observation is assimilated or not_

In the IODA converter, GsiEffectiveQC is 1 (failed) if errinv < 9e-12. Can you check the binary file again after replacing the tiny_single by this value and use a simple condition: "datachan(ich)%errinv .gt. tiny_single"?

gmao-wgu commented 1 year ago

I think it is better to still use id_qc as in the version 2 of IODA converter, so that we then have GSI qc mark stored in the GsiEffectiveQc. With the version 3 of IODA converter, we are like in the dark room and do not have any idea of what qc or id_qc cause the difference between GSI and UFO.

gmao-jjin3 commented 1 year ago

Wei, it was the same in the IODA v2 converter as in the v3 converter for radiance data. See my file for the v2 converter:
/discover/nobackup/jjin3/jedi/jedi-work/jedi_bundle/iodaconv/src/gsi-ncdiag/gsi_ncdiag.py

1472 # create a GSI effective QC variable 1473 gsiqcname = value, 'GsiEffectiveQC' 1474 errname = value, 'GsiFinalObsError' 1475 gsiqc = np.zeros_like(outdata[varDict[value]['valKey']]) 1476 gsiqc[outdata[errname] > 1e8] = 1 1477 gsiqc[np.reshape(self.var('QC_Flag'), (nlocs, nchans)) < 0] = 1

Note: GsiFinalObsError is an inverse of "Errinv_Final".

gmao-wgu commented 1 year ago

It is not Jianjun. In the version2, when I plot the GsiEffectiveQc, I got values like 0, 3, 7, 52, 53, 54 etc. While in the version 3, GsiEffectiveQC is either 0 or 1. Obviously, it is not assigned as the same thing between these two versions.

gmao-jjin3 commented 1 year ago

Can you point me your plotting program again?

gmao-wgu commented 1 year ago

/discover/nobackup/wgu/jedi_build/jjin3_plot and please compare the following two files: ufo_gsi_cris-fsr_n20_qc_ufo.pro with gsieffectiveqc from the version 3(x0048) and ufo_gsi_cris-fsr_n20_qc_ufo_x44.pro with gsieffectiveqc from the version 2(x0044)

gmao-jjin3 commented 1 year ago

@gmao-wgu Finally, I understand why UFO cannot reproduce GSI Cris-FSR_n20 results. Two observations are skipped by UFO because it thinks they are out of the time slot. There are 8311 profiles in UFO inputs, but there are 8309 profiles in the UFO output file. One of the two data passing QC in GSI. (Sign) The observation time is within one second after 21h, 12/11/2021 for those two data. However, UFO uses long integer values for time. As a result, the time for those two is exactly at 21:00:00 in UFO. JEDI folks think observations make exactly at the beginning of analysis cycle should be in the previous cycle and those data are masked out in OOPS: https://github.com/JCSDA-internal/oops/blob/develop/src/oops/base/GetValues.h#L327

This makes sense for conventional data, but not really for satellite observations though its impact is small. However, that makes a second shift in time slots between UFO and GSI, and a minor obstacle to reproduce GSI by UFO. Anyway, you can get a sense after running my IDL program /gpfsm/dnb34/jjin3/jedi/othters/Wei_Gu/x0048/gsi_nc4_rd.pro

As to QC_Flags, we did use PreQC values, which are GSI QC flags, in my IDL programs at the beginning. Later we abandoned the PreQC flag after we found GSI QC flags are not assigned for some IR data which don't pass QC in GSI. Therefore, they have to be used with a grain of salt. I tap @gmao-yzhu too.

gmao-wgu commented 1 year ago

@gmao-jjin3 @danholdaway @gmao-yzhu The 'id_qc' info is important for tracing the cause of the difference in observation usage between GSI and UFO, it is better to retain this info in the UFO until we have completely phased out GSI. It appears a bug for some short-wave observations from AIRS passing QC but being assigned a non-zero value. This issue should be addressed and fixed in GSI .

gmao-wgu commented 1 year ago

I have made code changes to the surface temperature check in UFO, ensuring that it uses the same algorithm as in GSI, but these modifications have had no impact on the CrIS-FSR QC, which has also been confirmed in GSI.

gmao-wgu commented 1 year ago

@gmao-jjin3 Based on the plots I have posted, it is obvious that more observations can pass QC in UFO although there are two observations fall outside the time window, indicating some QC filters are not functioning identically between GSI and UFO.

gmao-jjin3 commented 1 year ago

@gmao-wgu Thanks. Right. A few more observations should be tossed out by UFO. You can run ufo_gsi_cris-fsr_n20_qc_ufo.pro and see GSI flags for those data after copying it and plot_ufo_gsi_comp.pro from /discover/nobackup/jjin3/jedi/othters/Wei_Gu/x0048/

gmao-wgu commented 1 year ago

@gmao-jjin3 I was confused. In your directory, I could not spot any plots generated that include GSI flags("id_qc")

gmao-jjin3 commented 1 year ago

Wei, I didn't make that plot because we had better to checkout its printout. Anyway, I changed the program ufo_gsi_cris-fsr_n20_qc_ufo.pro, and now it make a figure named cris-fsr_n20.UFO_EffectiveQC.PreQC.ch200.ps.png. PreQC is the qc_flag in GSI. You can see a list of mismatched data after you run the program.

gmao-wgu commented 1 year ago

Do you mean that the testing files you provided last week already include the qc flag from GSI? It would have been helpful to clarify this at the beginning.

gmao-jjin3 commented 1 year ago

GSI qc flags are always saved as "PreQC" in UFO output files. We have to be careful about it. They may be misleading because of they are not assigned/re-assigned for some data tossed out in GSI.

gmao-wgu commented 1 year ago

what do you mean? as far as I knew, it is not reliable to use "errinv" to count the number of observations that pass QC from GSI, which is different from the number we get in both ana_stats file and diag binary file for CrIS-fsr_n20. Is "PreQC" assigned from the "id_qc" in setuprad.f90?

gmao-jjin3 commented 1 year ago

Wei, "errinv" is always assigned to zero for active channel observations to be tossed out. However, QC flag is not always assigned a non-zero value in setuprad.f90. Take a look at https://github.com/GEOS-ESM/GEOSana_GridComp/edit/develop/GEOSaana_GridComp/GSI_GridComp/setuprad.f90#L1377

As to the difference between data passing QC in UFO and GSI, which you checked in your binary file and the ana_state table, that is because of one-second time-slot shift and two observations were not even processed by UFO, as what I said in https://github.com/GEOS-ESM/swell/issues/65#issuecomment-1618825143.

gmao-wgu commented 1 year ago

@gmao-jjin3 Thanks for providing IDL program to read GSI diag output. However, I still have no clue why we are having trouble reading the PreQC from JEDI output using Eva.

gmao-wgu commented 1 year ago

With extensive debugging, it has been observed that incorrect sensor central wave numbers in the testing files are causing QC to behave differently between UFO and GSI, which are shown in the following:

init_ch= 22 water_frac= 1 solar= 61.62 wavenumber= 68500 init_ch= 23 water_frac= 1 solar= 61.62 wavenumber= 68562.5 init_ch= 24 water_frac= 1 solar= 61.62 wavenumber= 68625 init_ch= 25 water_frac= 1 solar= 61.62 wavenumber= 68687.5 init_ch= 26 water_frac= 1 solar= 61.62 wavenumber= 68750 init_ch= 27 water_frac= 1 solar= 61.62 wavenumber= 68812.5 init_ch= 28 water_frac= 1 solar= 61.62 wavenumber= 68875 init_ch= 29 water_frac= 1 solar= 61.62 wavenumber= 68937.5 init_ch= 30 water_frac= 1 solar= 61.62 wavenumber= 69000 init_ch= 31 water_frac= 1 solar= 61.62 wavenumber= 69062.5 init_ch= 32 water_frac= 1 solar= 61.62 wavenumber= 69125 init_ch= 33 water_frac= 1 solar= 61.62 wavenumber= 69187.5 init_ch= 34 water_frac= 1 solar= 61.62 wavenumber= 69250 init_ch= 35 water_frac= 1 solar= 61.62 wavenumber= 69312.5 init_ch= 36 water_frac= 1 solar= 61.62 wavenumber= 69375 init_ch= 37 water_frac= 1 solar= 61.62 wavenumber= 69437.5 init_ch= 38 water_frac= 1 solar= 61.62 wavenumber= 69500 init_ch= 39 water_frac= 1 solar= 61.62 wavenumber= 69562.5 init_ch= 40 water_frac= 1 solar= 61.62 wavenumber= 69625 init_ch= 41 water_frac= 1 solar= 61.62 wavenumber= 69687.5 init_ch= 43 water_frac= 1 solar= 61.62 wavenumber= 69812.5 init_ch= 44 water_frac= 1 solar= 61.62 wavenumber= 69875 init_ch= 46 water_frac= 1 solar= 61.62 wavenumber= 70000 init_ch= 47 water_frac= 1 solar= 61.62 wavenumber= 70062.5 init_ch= 49 water_frac= 1 solar= 61.62 wavenumber= 70187.5 init_ch= 51 water_frac= 1 solar= 61.62 wavenumber= 70312.5 init_ch= 52 water_frac= 1 solar= 61.62 wavenumber= 70375 init_ch= 53 water_frac= 1 solar= 61.62 wavenumber= 70437.5 init_ch= 54 water_frac= 1 solar= 61.62 wavenumber= 70500 init_ch= 55 water_frac= 1 solar= 61.62 wavenumber= 70562.5 init_ch= 56 water_frac= 1 solar= 61.62 wavenumber= 70625 init_ch= 57 water_frac= 1 solar= 61.62 wavenumber= 70687.5 init_ch= 58 water_frac= 1 solar= 61.62 wavenumber= 70750 init_ch= 59 water_frac= 1 solar= 61.62 wavenumber= 70812.5 init_ch= 60 water_frac= 1 solar= 61.62 wavenumber= 70875 init_ch= 61 water_frac= 1 solar= 61.62 wavenumber= 70937.5 init_ch= 62 water_frac= 1 solar= 61.62 wavenumber= 71000 init_ch= 64 water_frac= 1 solar= 61.62 wavenumber= 71125 init_ch= 65 water_frac= 1 solar= 61.62 wavenumber= 71187.5 init_ch= 66 water_frac= 1 solar= 61.62 wavenumber= 71250 init_ch= 67 water_frac= 1 solar= 61.62 wavenumber= 71312.5 init_ch= 68 water_frac= 1 solar= 61.62 wavenumber= 71375 init_ch= 69 water_frac= 1 solar= 61.62 wavenumber= 71437.5 init_ch= 71 water_frac= 1 solar= 61.62 wavenumber= 71562.5 init_ch= 72 water_frac= 1 solar= 61.62 wavenumber= 71625 init_ch= 73 water_frac= 1 solar= 61.62 wavenumber= 71687.5 init_ch= 74 water_frac= 1 solar= 61.62 wavenumber= 71750 init_ch= 76 water_frac= 1 solar= 61.62 wavenumber= 71875 init_ch= 77 water_frac= 1 solar= 61.62 wavenumber= 71937.5 init_ch= 78 water_frac= 1 solar= 61.62 wavenumber= 72000 init_ch= 79 water_frac= 1 solar= 61.62 wavenumber= 72062.5 init_ch= 85 water_frac= 1 solar= 61.62 wavenumber= 72437.5 init_ch= 86 water_frac= 1 solar= 61.62 wavenumber= 72500 init_ch= 87 water_frac= 1 solar= 61.62 wavenumber= 72562.5 init_ch= 88 water_frac= 1 solar= 61.62 wavenumber= 72625 init_ch= 89 water_frac= 1 solar= 61.62 wavenumber= 72687.5 init_ch= 91 water_frac= 1 solar= 61.62 wavenumber= 72812.5 init_ch= 92 water_frac= 1 solar= 61.62 wavenumber= 72875 init_ch= 93 water_frac= 1 solar= 61.62 wavenumber= 72937.5 init_ch= 95 water_frac= 1 solar= 61.62 wavenumber= 73062.5 init_ch= 101 water_frac= 1 solar= 61.62 wavenumber= 73437.5 init_ch= 106 water_frac= 1 solar= 61.62 wavenumber= 73750 init_ch= 111 water_frac= 1 solar= 61.62 wavenumber= 74062.5 init_ch= 117 water_frac= 1 solar= 61.62 wavenumber= 74437.5 init_ch= 118 water_frac= 1 solar= 61.62 wavenumber= 74500 init_ch= 123 water_frac= 1 solar= 61.62 wavenumber= 74812.5 init_ch= 124 water_frac= 1 solar= 61.62 wavenumber= 74875 init_ch= 126 water_frac= 1 solar= 61.62 wavenumber= 75000 init_ch= 128 water_frac= 1 solar= 61.62 wavenumber= 75125 init_ch= 129 water_frac= 1 solar= 61.62 wavenumber= 75187.5 init_ch= 130 water_frac= 1 solar= 61.62 wavenumber= 75250 init_ch= 131 water_frac= 1 solar= 61.62 wavenumber= 75312.5 init_ch= 132 water_frac= 1 solar= 61.62 wavenumber= 75375 init_ch= 133 water_frac= 1 solar= 61.62 wavenumber= 75437.5 init_ch= 134 water_frac= 1 solar= 61.62 wavenumber= 75500 init_ch= 135 water_frac= 1 solar= 61.62 wavenumber= 75562.5 init_ch= 136 water_frac= 1 solar= 61.62 wavenumber= 75625 init_ch= 137 water_frac= 1 solar= 61.62 wavenumber= 75687.5 init_ch= 138 water_frac= 1 solar= 61.62 wavenumber= 75750 init_ch= 139 water_frac= 1 solar= 61.62 wavenumber= 75812.5 init_ch= 191 water_frac= 1 solar= 61.62 wavenumber= 86312.5 init_ch= 192 water_frac= 1 solar= 61.62 wavenumber= 89250 init_ch= 193 water_frac= 1 solar= 61.62 wavenumber= 89937.5 init_ch= 194 water_frac= 1 solar= 61.62 wavenumber= 90062.5 init_ch= 195 water_frac= 1 solar= 61.62 wavenumber= 90187.5 init_ch= 196 water_frac= 1 solar= 61.62 wavenumber= 90312.5 init_ch= 197 water_frac= 1 solar= 61.62 wavenumber= 90562.5 init_ch= 198 water_frac= 1 solar= 61.62 wavenumber= 91625 init_ch= 200 water_frac= 1 solar= 61.62 wavenumber= 92437.5 init_ch= 213 water_frac= 1 solar= 61.62 wavenumber= 95312.5 init_ch= 214 water_frac= 1 solar= 61.62 wavenumber= 95375 init_ch= 215 water_frac= 1 solar= 61.62 wavenumber= 95562.5 init_ch= 217 water_frac= 1 solar= 61.62 wavenumber= 95937.5 init_ch= 218 water_frac= 1 solar= 61.62 wavenumber= 96125 init_ch= 219 water_frac= 1 solar= 61.62 wavenumber= 96250 init_ch= 236 water_frac= 1 solar= 61.62 wavenumber= 102188 init_ch= 247 water_frac= 1 solar= 61.62 wavenumber= 104062 init_ch= 250 water_frac= 1 solar= 61.62 wavenumber= 105312 init_ch= 253 water_frac= 1 solar= 61.62 wavenumber= 106125 init_ch= 265 water_frac= 1 solar= 61.62 wavenumber= 121250 init_ch= 282 water_frac= 1 solar= 61.62 wavenumber= 131000 init_ch= 283 water_frac= 1 solar= 61.62 wavenumber= 131500 init_ch= 290 water_frac= 1 solar= 61.62 wavenumber= 134938 init_ch= 297 water_frac= 1 solar= 61.62 wavenumber= 138562 init_ch= 301 water_frac= 1 solar= 61.62 wavenumber= 139375 init_ch= 308 water_frac= 1 solar= 61.62 wavenumber= 140250 init_ch= 325 water_frac= 1 solar= 61.62 wavenumber= 142500

gmao-jjin3 commented 1 year ago

Here is the reason: Original wavenumber values are multiplied by 100 in the IODA converter: https://github.com/JCSDA-internal/ioda-converters/blob/develop/src/gsi_ncdiag/gsi_ncdiag.py#L1653 Wavenumber is defined as 1/wave_length, and typically, in centimeter^-1. However, JEDI folks changed it to a SI unit, which is meter^-1. That's why those values are so larger now. I have to say new values are misleading simply because it is cm^-1 in all remote sensing literature and data documents. Anyway, we have to ask @BenjaminRuston whether JEDI wants to stick to this SI unit no one uses. Then, we'll decide how to solve this issue.

BenjaminRuston commented 1 year ago

thanks @gmao-jjin3 need to double check with Andrew Collard on this, as these conventions are a joint decision