Fix to ensure the initialization step will be performed for uninitialized channels #439

emilyhcliu · 2022-07-18T03:02:29Z

Issue

A very slow and less optimal initialization of radiance bias correction was observed in the gfsv16.3.0 parallel experiment (v163, blue line). The first analysis is 2021101600, and the initial bias estimation occurs in 2021101606. The issue described here did NOT occur in the low-resolution v16.x parallel experiments on HERA because older version of GSI was used in these experiments.

Example: AMSU-A NOAA-15 Channel 1
Vertical axis is the Number of Observation used (passed QC) in the analysis.

Vertical axis is Bias

The v163 (blue line) is observed from gfsv16.3.0 parallel experiment. The v163t (red line) is an experimental run with the fix proposed in this PR.

Please see Issue #438 for background, diagnostics, and the fix (more details).

bias correction coefficients. This modification ensures that the satellite bias correction initialization procedure will be performed.

ADCollard

Changes look good to me. I would suggest that @jderber-NOAA should give his approval before merging.

emilyhcliu · 2022-07-18T15:34:22Z

@ADCollard Thanks for your review. Let's wait for John's (@jderber-NOAA ) review and approval.

jderber-NOAA

I am a bit confused by this change. As I understand it, the maximum value for the bias correction variance is r10. It is not clear to me what the values for the variance is being set to for channels that are new.

It makes sense that the old code fails and does not work properly because it would assume any channel with a variance not equal to zero was a old already used channel.

However, with the new code, it seems that one of the varx values would have to be > r10 to be satisfied to make it an old channel. It would make more sense to me that if the change was to say if any(varx<r10) ... This would then say if any of the bias correction coefficients has a reasonable value, we should use the input file values and inew_rad(j) == .false.. Also, I think we have to be careful with varx values = r10. The code makes this the maximum value. So when starting a new bias correction calculation, the varx values should be > r10 so that the condition is not accidentally satisfied.

emilyhcliu · 2022-07-18T15:38:02Z

@ADCollard The current v16.3.0 parallel started from 20211015 18Z and is not on 20211029. Do you think we need to rewind the run with the fix?

emilyhcliu · 2022-07-18T16:42:53Z

@jderber-NOAA
Using amsua noaa15 channel 1 as an example:
The bias coefficient background error variance is set to zero at initial time (20211015 18Z):

 1 amsua_n15           1   0.000000E+00   0.000000E+00     0
      0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000   0.000000    0.000000    0.000000    
      0.000000    0.000000

On the first analysis date (20211016 00Z):
The iuse_flag will be set to -1 and amsua_n15 diagnostic file will be generated.
The output bias coefficient background error variance is set to 10 for all predictors

 1 amsua_n15                 1  0.0000000E+00
      0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  
      0.1000000E+02  0.1000000E+02

On the second analysis date 20211016 06Z::
In radinfo_read, it reads in satbias_pc from previous cycle and checks if any of the predictor variance (varx) is greater than 10. If so, then this channel would be determined as an old channel and no initialization required. If not, then this channel would be initialized based on samples from quasi mode.

      if ((any(varx>r10) .and. iuse_rad(j)>-2) .or. iuse_rad(j)==4) &
           inew_rad(j)=.false.
      cycle read3

Since in the previous cycle, for this channel, all varx values are 10, so inew_rad is still true (new channel), the channel will be initialized based on quasi mode in init_predx.
Here is the output satbias_pc for this cycle.

    1 amsua_n15                1  0.3744000E+04
      0.9712486E-03  0.1000000E+02  0.1000000E+02  0.4372668E-01  0.2901824E-01  0.1000000E+02  0.1000000E+02  0.1132186E-05  0.8866306E-01  0.3538301E-01
      0.1332834E-01  0.4459499E-02

emilyhcliu · 2022-07-18T16:44:57Z

I am a bit confused by this change. As I understand it, the maximum value for the bias correction variance is r10. It is not clear to me what the values for the variance is being set to for channels that are new.

@jderber-NOAA It is set to 10 for new channel (please see the example in the comment above).
In the old code, it is set to 10000.

It makes sense that the old code fails and does not work properly because it would assume any channel with a variance not equal to zero was a old already used channel.

However, with the new code, it seems that one of the varx values would have to be > r10 to be satisfied to make it an old channel. It would make more sense to me that if the change was to say if any(varx<r10) ... This would then say if any of the bias correction coefficients has a reasonable value, we should use the input file values and inew_rad(j) == .false.. Also, I think we have to be careful with varx values = r10. The code makes this the maximum value. So when starting a new bias correction calculation, the varx values should be > r10 so that the condition is not accidentally satisfied.

@jderber-NOAA The initial bias variance assigned to new channel is the maximum value 10 (defined in set_predictors_var (berror.f90)

jderber-NOAA · 2022-07-18T17:28:50Z

Using amsua noaa15 channel 1 as an example:
The bias coefficient background error variance is set to zero at initial time (20211015 18Z):

  1 amsua_n15                1   0.000000E+00   0.000000E+00     0
    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000
    0.000000    0.000000

@emilyhcliu The above values are 0., not 10. This is what does not make sense to me. 0. means the values have no error. Should be 10. or another larger number.

On the first analysis date (20211016 00Z):
The iuse_flag will be set to -1 and amsua_n15 diagnostic file will be generated.
The output bias coefficient background error variance is set to 10 for all predictors

 1 amsua_n15                 1  0.0000000E+00
      0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02  0.1000000E+02
      0.1000000E+02  0.1000000E+02

On the second analysis date 20211016 06Z::
In radinfo_read, it reads in satbias_pc from previous cycle and checks if any of the predictor variance (varx) is greater than 10.

@emilyhcliu So all input variances are equal to 10. Since 10 is not greater than 10, I would assume this makes this fail the test.
I don't know how an input value for varx would ever be greater than 10.

If so, then this channel would be determined as an old channel and no initialization required. If not, then this channel would be initialized based on samples from quasi mode.

   if ((any(varx>r10) .and. iuse_rad(j)>-2) .or. iuse_rad(j)==4) &
      inew_rad(j)=.false.
   cycle read3

Since in the previous cycle, for this channel, all varx values are 10, so inew_rad is still true (new channel), the channel will be initialized based on quasi mode in init_predx.
Here is the output satbias_pc for this cycle.

1 amsua_n15                1  0.3744000E+04
  0.9712486E-03  0.1000000E+02  0.1000000E+02  0.4372668E-01  0.2901824E-01  0.1000000E+02  0.1000000E+02  0.1132186E-05  0.8866306E-01  0.3538301E-01
  0.1332834E-01  0.4459499E-02

@emilyhcliu So with this output the maximum value is still 10. Wouldn't that make it fail the test again? Then it will reinitialize the bias correction? Wouldn't it do this for all channels and instruments?

emilyhcliu · 2022-07-18T22:32:32Z

Please see my comment below

emilyhcliu · 2022-07-18T22:33:16Z

@jderber-NOAA Your comment totally make sense to me. The code I modified (see below) only works for initial time when all varx values are 10. For the following cycle, the old channels will still be set as new channel because the valid varx values are always less than 10.

I did not look careful enough to realize that 10 is the maximum rather than a minumum value. The varx is not a physical variable like wind speed or temperature, so 10 for varx could be minimum or maximum depending on how the preconditioning is designed together with all other control variables.

Definitely, I should change "any(varx>r10)" to "any(varx<r10)".

if ((any(varx>r10) .and. iuse_rad(j)>-2) .or. iuse_rad(j)==4) &
inew_rad(j)=.false.
cycle read3

But, still, I have to answer why my test (v163t) looks correct.

The blue line (v163) is the original (no quasi-mode initialization)
The red line (v163t) is the test with my fix (with quai-mode initialization)
Here is the link to the RadMon for v163 and v163t comparison

The answers are the following:

As I mentioned before, the satbias and satbias_pc were zeroed out in the ICs, and the AMSU-A, ATMS, MHS and CrIS NPP diagnostic files were removed. There is no bias initialization for the first cycle (2021101518) since there are no diagnostic files. However, at the end of this first cycle, the all zeros in the input satbias_pc file will be replaced by 10 (in set_predictors_var).
With the diagnostic files generated from the first cycle, the second cycle (2021101600) should initiate the bias correction estimation using quasi model for new channels. My modification (even though with flaw) worked in this case: the input satbias_pc has all varx filled with 10. So, it will not trigger the condition to set a new channnel into an old one. The new channels are initialized using quasi mode as expected.
But, you are right, the following cycle (the third cycle, 2021101606), new channels will not be set to old channels since the valid varx will be all smaller than 10.
Fortunately, GSI code has its own idiot-proof mechanism to deal with people like me. Following the inew_chan logic in the read section for satbias_pc, the read section for satbias_in actually checks for inew_chan again (see below):

          if(  .not. cold_start_seviri .or. isis(1:6) /= 'seviri' .or. .not. bias_zero_start ) then
             do j =1,jpch_rad
                if(trim(isis) == trim(nusis(j)) .and. ichan == nuchan(j))then
                   cfound = .true.
                   nfound(j) = .true.
                   do i=1,npred
                      if (iuse_rad(j)==4) then
                         predx(i,j)=zero
                      else
                         predx(i,j)=predr(i)
                      end if
                   end do
                   if (adp_anglebc) then
                      tlapmean(j)=tlapm
                      tsum_tlapmean(j)=tsum
                      count_tlapmean(j)=ntlapupdate
                      if (ntlapupdate > ntlapthresh) update_tlapmean(j)=.false.
                   end if
                   if (any(predr/=zero) .or. iuse_rad(j)==4) inew_rad(j)=.false.
                   cycle read4
                end if
             end do
          endif

The inew_chan logic is checked again based on the predr (bias correction coefficients) from satbias_in file. So, even though my modification failed to detect old channels in the previous section, the inew_chan is checked again in the read section for satbias.

I corrected my mistake and is running a short experiment again. I will post result soon.
@jderber-NOAA and @ADCollard Thanks for your review and discussion.
@jderber-NOAA Thanks for spotting my mistake.

… satbias_pc file. Modify varx condition to detect old channel.

emilyhcliu · 2022-07-19T15:59:33Z

@jderber-NOAA and @ADCollard
I pushed code changes to the branch. Please let me know if you have any comments. I am running a short experiment with this code change. I will post the result here.

jderber-NOAA

Looks good now!

emilyhcliu · 2022-07-21T20:57:55Z

@jderber-NOAA Thanks for your review and discussion. I merge the bug fix into gfsv16.3.0 release.

Modify the threshold in checking the minminum backgound variance for

69e83f6

bias correction coefficients. This modification ensures that the satellite bias correction initialization procedure will be performed.

emilyhcliu added the bug Something isn't working label Jul 18, 2022

emilyhcliu requested review from jderber-NOAA, ADCollard, RussTreadon-NOAA and CatherineThomas-NOAA July 18, 2022 03:06

emilyhcliu self-assigned this Jul 18, 2022

emilyhcliu requested a review from dtkleist July 18, 2022 05:09

ADCollard approved these changes Jul 18, 2022

View reviewed changes

jderber-NOAA reviewed Jul 18, 2022

View reviewed changes

Add handling for channel initialized with zero background variance in…

87aac5f

… satbias_pc file. Modify varx condition to detect old channel.

jderber-NOAA approved these changes Jul 20, 2022

View reviewed changes

emilyhcliu merged commit 82620df into release/gfsda.v16.3.0 Jul 23, 2022

emilyhcliu deleted the release/gfsda.v16.3.0_emily branch July 27, 2022 20:38

This was referenced Aug 2, 2022

GFS v16.3 retro parallel for implementation NOAA-EMC/global-workflow#951

Closed

v16.3 DA pre-implementation parallel NOAA-EMC/global-workflow#776

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix to ensure the initialization step will be performed for uninitialized channels #439

Fix to ensure the initialization step will be performed for uninitialized channels #439

emilyhcliu commented Jul 18, 2022 •

edited

Loading

ADCollard left a comment

emilyhcliu commented Jul 18, 2022

jderber-NOAA left a comment

emilyhcliu commented Jul 18, 2022

emilyhcliu commented Jul 18, 2022 •

edited

Loading

emilyhcliu commented Jul 18, 2022 •

edited

Loading

jderber-NOAA commented Jul 18, 2022 •

edited by emilyhcliu

Loading

emilyhcliu commented Jul 18, 2022 •

edited

Loading

emilyhcliu commented Jul 18, 2022 •

edited

Loading

emilyhcliu commented Jul 19, 2022

jderber-NOAA left a comment

emilyhcliu commented Jul 21, 2022

Fix to ensure the initialization step will be performed for uninitialized channels #439

Fix to ensure the initialization step will be performed for uninitialized channels #439

Conversation

emilyhcliu commented Jul 18, 2022 • edited Loading

Issue

Please see Issue #438 for background, diagnostics, and the fix (more details).

ADCollard left a comment

Choose a reason for hiding this comment

emilyhcliu commented Jul 18, 2022

jderber-NOAA left a comment

Choose a reason for hiding this comment

emilyhcliu commented Jul 18, 2022

emilyhcliu commented Jul 18, 2022 • edited Loading

emilyhcliu commented Jul 18, 2022 • edited Loading

jderber-NOAA commented Jul 18, 2022 • edited by emilyhcliu Loading

emilyhcliu commented Jul 18, 2022 • edited Loading

emilyhcliu commented Jul 18, 2022 • edited Loading

emilyhcliu commented Jul 19, 2022

jderber-NOAA left a comment

Choose a reason for hiding this comment

emilyhcliu commented Jul 21, 2022

emilyhcliu commented Jul 18, 2022 •

edited

Loading

emilyhcliu commented Jul 18, 2022 •

edited

Loading

emilyhcliu commented Jul 18, 2022 •

edited

Loading

jderber-NOAA commented Jul 18, 2022 •

edited by emilyhcliu

Loading

emilyhcliu commented Jul 18, 2022 •

edited

Loading

emilyhcliu commented Jul 18, 2022 •

edited

Loading