Error in potential temperature values for cold-air-outbreak quality control in all-sky.

[jianjunj]

Atmospheric stability is used to check cold-air-outbreak in all-sky assimilations. It is the difference between background potential temperature near 700 hPa and those at surface. These potential temperature values near 700 hPa are incorrectly calculated in crtm_interface.f90.

Atmospheric pressure and temperature at the nearest model level to 700 hPa are used to calculate potential temperature in crtm_interface.f90. See line 2204 - 2208
https://github.com/NOAA-EMC/GSI/blob/develop/src/gsi/crtm_interface.f90#L2204
and
https://github.com/NOAA-EMC/GSI/blob/develop/src/gsi/crtm_interface.f90#L2214

Pressure differences "dprs" and "dprs_min" are used to search the nearest model grid. They are assigned by "real_kind" values. https://github.com/NOAA-EMC/GSI/blob/develop/src/gsi/crtm_interface.f90#L2067
However, they are defined as "i_kind" values by the line "integer(i_kind):: idx700,dprs,dprs_min":
https://github.com/NOAA-EMC/GSI/blob/develop/src/gsi/crtm_interface.f90#L1164
This creates wrong values, which are mostly zeros, for "dprs" and "dprs_min" and yields wrong model level id (idx700) for the calculation of potential temperature.

This bug potentially has a minor impact in some (AMSU-A and ATMS) microwave radiance data selection in all-sky conditions.

[emilyhcliu]

So, the dprs and dprs_min should be r_kind
The definition of i_kind rounds down the dprs and dprs_min

[emilyhcliu]

The top plot below is the stability calculation from potential temperature using the current GSI code. It looks like the stability calculation (< 12) catches the unstable regions (top plot). I think the impact of the bug in real and integer definitions for dprs and dprs_min may be small, I hope....!
Let's compare the CAO screening flags (left panel of the bottom plot) with and without the bug to check the impact (hopefully, it is small).

Thanks for finding and fixing the bug.

[jianjunj]

@emilyhcliu Yes, they should be "r_kind". I also had thought "i_kind" was fine. :)

[RussTreadon-NOAA]

WCOSS2 test

Clone GSI develop at 27c03e8 on Cactus. Make the following change to src/gsi/crtm_interface.f90

@@ -1161,7 +1161,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   integer(i_kind):: iqs,iozs,icfs
   integer(i_kind):: inis,inrs   
   integer(i_kind):: error_status_clr
-  integer(i_kind):: idx700,dprs,dprs_min  
+  integer(i_kind):: idx700
   integer(i_kind),dimension(8)::obs_time,anal_time
   integer(i_kind),dimension(msig) :: klevel
 ! ****************************** 
@@ -1178,6 +1178,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   real(r_kind):: hwp_total,theta_700,theta_sfc,hs
   real(r_kind):: dlon,dlat,dxx,dyy,yy,zz,garea
   real(r_kind):: radiance, radiance_overcast, radiance_ratio
+  real(r_kind):: dprs,dprs_min
   real(r_kind),dimension(0:3):: wgtavg
   real(r_kind),dimension(nsig,nchanl):: omix
   real(r_kind),dimension(nsig,nchanl,n_aerosols_jac):: jaero

Compile modified code along with unaltered copy of develop at 27c03e8. The only regression test that runs with cao_check=.true. is global_4denvar. Run global_4denvar ctest with the following result

Test project /lfs/h2/emc/da/noscrub/russ.treadon/git/gsi/test/build
    Start 1: global_4denvar
1/1 Test #1: global_4denvar ...................   Passed  1803.61 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) = 1803.63 sec

The global_4denvar regression case is a low resolution (C48L127) analysis for 20240223 00Z. This case may not contain condition which trigger the CAO check.

Run operational 20250110 00Z gdas case using unaltered develop and test code. This case shows a difference between the original and modified gsi.x. The initial radiance penalties differ

< radiance                     1.3409973133088651E+06
---
> radiance                     1.3409919627781105E+06

< is the original develop. > is the modified copy of develop.

The following satellite_sensor differ for o-g 01

< o-g 01 rad   n18          amsua           1594335       126936        83250    40153.       40153.      0.48231      0.48231    
> o-g 01 rad   n18          amsua           1594335       126936        83244    40152.       40152.      0.48235      0.48235    
< o-g 01 rad   n19          amsua           1636620       129862        91878    36386.       36386.      0.39602      0.39602    
> o-g 01 rad   n19          amsua           1636620       129862        91864    36383.       36383.      0.39606      0.39606    
< o-g 01 rad   n20          atms           15399934       236135       166446    29661.       29661.      0.17820      0.17820    
> o-g 01 rad   n20          atms           15399934       236135       166422    29660.       29660.      0.17822      0.17822    
< o-g 01 rad   metop-c      amsua           2220420       136052       106853    60008.       60008.      0.56159      0.56159    
> o-g 01 rad   metop-c      amsua           2220420       136052       106846    60006.       60006.      0.56161      0.56161    
< o-g 01 rad   n21          atms           15400308       242487       173890    28334.       28334.      0.16294      0.16294    
> o-g 01 rad   n21          atms           15400308       242487       173863    28333.       28333.      0.16296      0.16296    

For each case the modified gsi.x is assimilating fewer observations with a corresponding slightly smaller initial penalty.

Here are relevant directories for this test

• original develop: /lfs/h2/emc/da/noscrub/russ.treadon/git/gsi/develop
  • 20250110 00Z run directory: /lfs/h2/emc/ptmp/russ.treadon/tmp766/develop.gdas.2025011000
• modified develop: /lfs/h2/emc/da/noscrub/russ.treadon/git/gsi/test
  • 20250110 00Z run directory: /lfs/h2/emc/ptmp/russ.treadon/tmp766/test.gdas.2025011000

Note: While the GSI runs were configured to generate o-g and o-a diagnostics, the job script did not concatenate the pe specific diagnostic files into radstat files. The above run directories remain on disk so others can assemble radstat files if they want to more closely examine differences.

[emilyhcliu]

@RussTreadon-NOAA Thanks so much for doing these tests.
My understanding from your tests (original code vs. modified code) is the following:

• The low-resolution (regression test) global_4denvar tests do not show the difference
• The operational-resolution global_4denvar tests show differences - slightly more observations are filtered out with the modified code.

Do I understand correctly? If so, the impact may be small.
We will plot the CAO QC flags to verify and confirm.

[RussTreadon-NOAA]

Yes, @emilyhcliu , your understanding is correct.

[jianjunj]

Thank you very much, Russ. The difference is tiny. I could not match the potential temperature at 700 hPa made by GSI and by my UFO function. Eventual traced the difference to the data types.

Also, @emilyhcliu is right that "The definition of i_kind rounds down the dprs and dprs_min". I made a mistake to print out values when I debug it yesterday.

[jianjunj]

I had one simple run using the data for UFO tests and got 7 more observations tossed out by the cold-air-outbreak (qc_flag=12) but 7 fewer observations tossed out by another QC check. In other words, I got the same number of data passing quality control. In addition, there is zero difference in OmF and Obs Error after QC.

qc_flag= -3  Obs # (wrong idx700)= 11226  Obs # (correct idx700) = 11226  diff = 0
qc_flag= 0  Obs # (wrong idx700)= 59663  Obs # (correct idx700) = 59663  diff = 0
qc_flag= 3  Obs # (wrong idx700)= 7731  Obs # (correct idx700) = 7731  diff = 0
qc_flag= 4  Obs # (wrong idx700)= 1691  Obs # (correct idx700) = 1691  diff = 0
qc_flag= 7  Obs # (wrong idx700)= 210  Obs # (correct idx700) = 210  diff = 0
qc_flag= 8  Obs # (wrong idx700)= 396  Obs # (correct idx700) = 396  diff = 0
qc_flag= 12  Obs # (wrong idx700)= 1358  Obs # (correct idx700) = 1365  diff = 7
qc_flag= 50  Obs # (wrong idx700)= 17150  Obs # (correct idx700) = 17150  diff = 0
qc_flag= 51  Obs # (wrong idx700)= 2640  Obs # (correct idx700) = 2640  diff = 0
qc_flag= 53  Obs # (wrong idx700)= 33040  Obs # (correct idx700) = 33033  diff = -7
range of differences in omf after qc: 0.0  ~  0.0
range of differences in 1/error after qc: 0.0  ~  0.0

The difference in the cold-air-outbreak is because model pressure levels have the unit of "hPa" and the nearest grid level to 700 hPa can be incorrectly identified for the calculation of potential temperature at 700 hPa when the difference in pressure values is round off. Here is an example that the level is identified to level 87 instead of level 88 which is closer to 700 hPa.

84  air_pressure= 64492.26247414376 , |P/100-700|= 55 ,k for 700hPa= 84
85  air_pressure= 66092.30684694136 , |P/100-700|= 39 ,k for 700hPa= 85
86  air_pressure= 67664.65752080704 , |P/100-700|= 23 ,k for 700hPa= 86
87  air_pressure= 69206.93084346283 , |P/100-700|= 7 ,k for 700hPa= 87
88  air_pressure= 70716.95240472407 , |P/100-700|= 7 ,k for 700hPa= 87
89  air_pressure= 72192.77110642493 , |P/100-700|= 21 ,k for 700hPa= 87

In conclusion, this issue is negligible though the minor correction is needed.