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

GitHub Issue NOAA-EMC/GSI#604 Undefined values found in radar reflectivity direct DA #605

Merged

Conversation

shoyokota
Copy link
Collaborator

@shoyokota shoyokota commented Aug 4, 2023

Description

To prevent some undefined values found in the radar reflectivity direct DA (if_model_dbz=T and l_use_dbz_directDA=F), corresponding parts are fixed. It doesn't change the result except for the case of the execution with the debug option.

Fixes #604

Type of change

  • Bug fix (non-breaking change which fixes an issue)

How Has This Been Tested?

The radar reflectivity DA test with the RRFS setting was done on Orion. After this modification, EnVar was completed even with the debug option. This modification didn't change the result in the test without the debug option.

Checklist

  • My code follows the style guidelines of this project
  • I have performed a self-review of my own code
  • I have commented my code, particularly in hard-to-understand areas
  • New and existing tests pass with my changes
  • Any dependent changes have been merged and published

Due date for this PR is 9/15/2023. If this PR is not merged into develop by this date, the PR will be closed and returned to the developer.

if(l4dvar.or.l4densvar) then
timeb=real(mins_an-iwinbgn,r_kind)
else
timeb=zero
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why could we just use the " timeb=real(mins_an-iwinbgn,r_kind)" ? The obs relative to the analysis time could be non-zero even for 3dvar . When the fget is used for 3dvar, the non-zero timeb would also produce differences for timeb=0 case. So, I think if the general use of "timeb=real(mins_an-iwinbgn,r_kind)" does work, this should be unanimously used.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In my understanding, the radar reflectivity NetCDF file doesn't include the observation time and it is assumed that the observation time is the same as the analysis time. That is why I used "timeb=real(mins_an-iwinbgn,r_kind)" for 4DVar and 4DEnVar.

When the fget is used for 3dvar, the non-zero timeb would also produce differences for timeb=0 case.

Are you pointing out the case of FGAT (First Guess at Appropriate Time)? I forgot to think about FGAT. I will fix to use "timeb=real(mins_an-iwinbgn,r_kind)" also in the 3DVar case.

@shoyokota shoyokota force-pushed the feature/PR_NOAA-EMC_EnVar-DBZ2 branch from a4dc2cc to b96a464 Compare August 9, 2023 18:32
@@ -337,6 +338,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no

call w3fs21(iadate,mins_an) !mins_an -integer number of mins snce 01/01/1978
rmins_an=mins_an !convert to real number
timeb=real(mins_an-iwinbgn,r_kind)
Copy link
Contributor

@TingLei-NOAA TingLei-NOAA Aug 9, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure. this should work for various time-related situations.
An update: Please allow me share my clarified thoughts now:(.
I think here we should stick to the "assumption" that all dbz obs are at analysis time. Hence, the timeb should always be timeb=0. How will it be used in 4d da and fgat should be considered in other places including the specification of variables like t4df.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In other places, timeb is defined by (observation time) - (head of assimilation window) for 4D DA. In case (observation time) = (analysis time) is assumed, 'timeb=real(mins_an-iwinbgn,r_kind)' is required. In the case of 'timeb=zero', GSI regards that reflectivity was observed at the head of the assimilation window, but we don't expect that.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In read_dbz_nc.f90, timeb is commented as "observation time relative to analysis time" and is used as it as in

 cdata_all(7,iout) = timeb*r60inv                  ! obs time (analyis relative hour)

So, please stick to this definition of timeb. If this sub couldn't have any information about actual observation time, just let timeb=0 for both 4dvar/fgat and 3dvar.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I confirmed that reflectivity was assimilated at the head of the assimilation window (not analysis time) in the case of 'timeb=zero' in 4DEnVar. So, I think this comment is wrong. How should I do?

Copy link
Contributor

@TingLei-NOAA TingLei-NOAA Aug 9, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@shoyokota I think t4dv is also undefined as shown below. Is this an issue when you ran your case with debug mode gsi? That will be great if it is also taken care of in this PR. From my point of view, it think it could be defined as timeb+ (analysis time - the beginning time of the time window). You could see how it is set in other observation 's setup subroutine
image

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@TingLei-NOAA This t4dv is only used in thinning independently in each timeslot. However, it is assumed that all the inputted reflectivity observations were observed at the analysis time here. So, we don't have to input any specific times as t4dv here, and the result doesn't depend on t4dv as long as it is the same in all reflectivity observations. Here, I modified here to timedif=zero simply to prevent undefined values.

Copy link
Contributor

@TingLei-NOAA TingLei-NOAA Aug 16, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@shoyokota if you runs doesn't go through this part. We could let it be the current codes (with that comment) for future developers to take care . It should be analysis time - the begin of the time window. We don't need to enforce this part also work with that all obs are being at the analysis time . And, even with assumption/requirement that dbz obs being on the analysis time, the fgat/4dvar could still be invoked for other observations could be distributed on different time levels. So, I prefer leaving it to be updated in the future rather than making a change with unnecessary assumptions.

Copy link
Collaborator Author

@shoyokota shoyokota Aug 16, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This part is computed in the case of radar_no_thinning=F even now, so I think some kind of fix is needed. Is it better to add comments such as 'observation time (hour) should be input for 4D observations in 3D thinning, but 3D observations here' at the right of 'timedif=zero'?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@shoyokota Thanks for the clarification. I will add my thoughts on those specific codes respectively.

@shoyokota shoyokota force-pushed the feature/PR_NOAA-EMC_EnVar-DBZ2 branch from b96a464 to 664d072 Compare August 11, 2023 11:35
@@ -337,6 +338,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no

call w3fs21(iadate,mins_an) !mins_an -integer number of mins snce 01/01/1978
rmins_an=mins_an !convert to real number
timeb=real(mins_an-iwinbgn,r_kind)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I think now, we agree timeb should be the observation time relative to the begin of the time window.
You could change that wrong comment . And, could you confirm in "timeb=real(mins_an-iwinbgn,r_kind)" , iwindbgn is in the same time coordinates as mins_an?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I tested 'timeb=real(mins_an-iwinbgn,r_kind)' with 4DEnVar, and radar reflectivity was assimilated at the analysis time as we expected.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good. Thanks.

@@ -452,7 +454,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no


ntmp=ndata ! counting moved to map3gridS
timedif=abs(t4dv) !don't know about this
timedif=zero
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you add a comment that here, the obs are assumed on the analysis time?

Copy link
Contributor

@TingLei-NOAA TingLei-NOAA left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is great for this PR to sort out some obs time related issues in dbz obs processing.
Thanks to the developer!

Copy link
Contributor

@TingLei-NOAA TingLei-NOAA left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just , hopefully, the last question. For line 493 (read_dbz_nc.f90)

 cdata_all(7,iout) = timeb*r60inv                  ! obs time (analyis relative hour)

The comments should be changed to "!obs time relative to the beginning of the DA window".

@shoyokota shoyokota force-pushed the feature/PR_NOAA-EMC_EnVar-DBZ2 branch from 664d072 to 7b0d453 Compare August 16, 2023 16:48
@ShunLiu-NOAA
Copy link
Contributor

@chunhuazhou Could you review Sho's changes? Thanks.

@ShunLiu-NOAA
Copy link
Contributor

@hu5970 could you handle this PR during my absence?

@hu5970
Copy link
Collaborator

hu5970 commented Aug 30, 2023

I run regression tests on WCOSS2 (cactus):

1/9 Test #8: netcdf_fv3_regional ..............***Failed  544.25 sec
2/9 Test #5: hwrf_nmm_d3 ......................***Failed  553.58 sec
3/9 Test #7: rrfs_3denvar_glbens ..............   Passed  605.18 sec
4/9 Test #4: hwrf_nmm_d2 ......................   Passed  606.28 sec
5/9 Test #9: global_enkf ......................   Passed  609.68 sec
6/9 Test #6: rtma .............................   Passed  1087.63 sec
7/9 Test #3: global_4denvar ...................   Passed  1382.15 sec
8/9 Test #1: global_3dvar .....................   Passed  1501.75 sec
9/9 Test #2: global_4dvar .....................   Passed  1561.96 sec

78% tests passed, 2 tests failed out of 9

Total Test time (real) = 1561.96 sec

The following tests FAILED:
	  5 - hwrf_nmm_d3 (Failed)
	  8 - netcdf_fv3_regional (Failed)
Errors while running CTest

The netcdf_fv3_regional and hwrf_nmm_d3 failed because of the scalability test. Not the critical failure.

@hu5970 hu5970 self-assigned this Aug 30, 2023
@hongli-wang
Copy link
Collaborator

@shoyokota @hu5970
FYI. The below line in setupdbz.f90 crashes my GSI test run with DBZ DA when the code is compiled in debug model. Any thought?
call nc_diag_metadata("Station_Elevation", sngl(data(ielev,i)) )

@hu5970
Copy link
Collaborator

hu5970 commented Sep 7, 2023

@shoyokota @hu5970 FYI. The below line in setupdbz.f90 crashes my GSI test run with DBZ DA when the code is compiled in debug model. Any thought? call nc_diag_metadata("Station_Elevation", sngl(data(ielev,i)) )

The current GSI develop branch has some trouble to run under DEBUG mode.

@TingLei-NOAA
Copy link
Contributor

@hongli-wang Since such issue was not found in @shoyokota 's dbz DA using GSI at debug mode, could you do some further tracing on this issue? Missing values? a wrong array index? And you could make quick PR with the fix. Or if your case are on (or can be transfered to ) hera or orion, you can point me to your case and I could have a look into it. Thanks!

@shoyokota
Copy link
Collaborator Author

shoyokota commented Sep 7, 2023

@shoyokota @hu5970 FYI. The below line in setupdbz.f90 crashes my GSI test run with DBZ DA when the code is compiled in debug model. Any thought? call nc_diag_metadata("Station_Elevation", sngl(data(ielev,i)) )

@hongli-wang @hu5970 Thank you for letting me know the crash. Even after the current modification in this PR, it seems that some undefined values are still detected in the case of "DIAG_RADARDBZ=T." (My test was being done with "DIAG_RADARDBZ=F.") Now I'm working on additional modification to prevent these undefined values. I guess it can be prevented by small modification of read_dbz_nc.f90. Could you extend the due date of this PR until completing my additional modification and its review?

@hongli-wang
Copy link
Collaborator

@shoyokota @hu5970 FYI. The below line in setupdbz.f90 crashes my GSI test run with DBZ DA when the code is compiled in debug model. Any thought? call nc_diag_metadata("Station_Elevation", sngl(data(ielev,i)) )

@hongli-wang @hu5970 Thank you for letting me know the crash. Even after the current modification in this PR, it seems that some undefined values are still detected in the case of "DIAG_RADARDBZ=T." (My test was being done with "DIAG_RADARDBZ=F.") Now I'm working on additional modification to prevent these undefined values. I guess it can be prevented by small modification of read_dbz_nc.f90. Could you extend the due date of this PR until completing my additional modification and its review?

@shoyokota @TingLei-NOAA
diag_radardbz=.true. in my gsiparm.anl with your PR code. I guess that is the reason. Thanks for looking into this.

@hu5970
Copy link
Collaborator

hu5970 commented Sep 8, 2023

@shoyokota Please sync with develop and test the code again. This PR is straightforward and we should be able make it in time. Thanks, Ming

@shoyokota shoyokota force-pushed the feature/PR_NOAA-EMC_EnVar-DBZ2 branch 2 times, most recently from 4fea628 to be615b3 Compare September 9, 2023 03:14
src/gsi/read_dbz_nc.f90 Outdated Show resolved Hide resolved
@shoyokota
Copy link
Collaborator Author

@hu5970 @ShunLiu-NOAA The modification and the sync with develop have been done. Could you confirm if it is possible to be merged?

@ShunLiu-NOAA
Copy link
Contributor

@shoyokota I will run a regression test on WCOSS2. Could you post your ORION regression test result here?

@shoyokota
Copy link
Collaborator Author

@ShunLiu-NOAA I did sync with the current develop branch. Now the regression tests are running on Orion. Once they are finished, I'll show the results.

@ShunLiu-NOAA
Copy link
Contributor

WCOSS2 regression tests were completed on CACTUS. The result is posted here.
There are 4 failed tests. However, the results from all four test cases were reproducible. All four test cases had failed time scalability test. These are not fatal failures.

Test project /lfs/h2/emc/lam/noscrub/emc.lam/Shun.Liu/gsi_regression/GSI/build
Start 1: global_3dvar
1/9 Test #1: global_3dvar ..................... Passed 1671.04 sec
Start 2: global_4dvar
2/9 Test #2: global_4dvar ..................... Passed 1683.30 sec
Start 3: global_4denvar
3/9 Test #3: global_4denvar ...................***Failed 1616.01 sec
Start 4: hwrf_nmm_d2
4/9 Test #4: hwrf_nmm_d2 ...................... Passed 728.53 sec
Start 5: hwrf_nmm_d3
5/9 Test #5: hwrf_nmm_d3 ...................... Passed 493.23 sec
Start 6: rtma
6/9 Test #6: rtma ............................. Passed 1693.55 sec
Start 7: rrfs_3denvar_glbens
7/9 Test #7: rrfs_3denvar_glbens ..............***Failed 910.02 sec
Start 8: netcdf_fv3_regional
8/9 Test #8: netcdf_fv3_regional ..............***Failed 788.18 sec
Start 9: global_enkf
9/9 Test #9: global_enkf ......................***Failed 1776.11 sec

@shoyokota
Copy link
Collaborator Author

The ORION regression test was also completed. The test #4 exceeded maximum allowable threshold time and failed the scalability test. The test #8 failed only the scalability test. However, these are not fatal failures and all tests were reproducible.

Test project /home/syokota/da/gsi/ctest/PR_NOAA-EMC_EnVar-DBZ2/build
Start 1: global_3dvar
1/9 Test #1: global_3dvar ..................... Passed 1974.09 sec
Start 2: global_4dvar
2/9 Test #2: global_4dvar ..................... Passed 1863.71 sec
Start 3: global_4denvar
3/9 Test #3: global_4denvar ................... Passed 1682.52 sec
Start 4: hwrf_nmm_d2
4/9 Test #4: hwrf_nmm_d2 ......................***Failed 791.99 sec
Start 5: hwrf_nmm_d3
5/9 Test #5: hwrf_nmm_d3 ...................... Passed 498.67 sec
Start 6: rtma
6/9 Test #6: rtma ............................. Passed 1272.71 sec
Start 7: rrfs_3denvar_glbens
7/9 Test #7: rrfs_3denvar_glbens .............. Passed 670.20 sec
Start 8: netcdf_fv3_regional
8/9 Test #8: netcdf_fv3_regional ..............***Failed 424.61 sec
Start 9: global_enkf
9/9 Test #9: global_enkf ...................... Passed 835.90 sec

@hu5970
Copy link
Collaborator

hu5970 commented Sep 28, 2023

@ShunLiu-NOAA should one of us send a request to merge this PR tomorrow? Thanks, Ming

@ShunLiu-NOAA
Copy link
Contributor

@shoyokota could you please sync your branch with "develop"? Ming and I will run regression test again on WCOSS2 and HERA. Hopefully, we can merge this today.

@shoyokota
Copy link
Collaborator Author

@ShunLiu-NOAA @hu5970 I synced my branch with "develop" now.

@hu5970
Copy link
Collaborator

hu5970 commented Sep 29, 2023

@shoyokota Thanks, Shun and I are working on regression cases on Hera and WCOSS2.

@ShunLiu-NOAA
Copy link
Contributor

@shoyokota could you please sync your branch with develop again? there is another PR merged into develop. Thank you.

@shoyokota
Copy link
Collaborator Author

@ShunLiu-NOAA I synced it.

@hu5970
Copy link
Collaborator

hu5970 commented Sep 29, 2023

The regression test on Hera are all good:

[Ming.Hu@hfe02 build]$ ctest -j9
Test project /scratch1/BMC/wrfruc/mhu/gsi/sho/GSI/build
    Start 1: global_3dvar
    Start 2: global_4dvar
    Start 3: global_4denvar
    Start 4: hwrf_nmm_d2
    Start 5: hwrf_nmm_d3
    Start 6: rtma
    Start 7: rrfs_3denvar_glbens
    Start 8: netcdf_fv3_regional
    Start 9: global_enkf
1/9 Test #8: netcdf_fv3_regional ..............   Passed  486.33 sec
2/9 Test #4: hwrf_nmm_d2 ......................   Passed  492.15 sec
3/9 Test #5: hwrf_nmm_d3 ......................   Passed  504.07 sec
4/9 Test #7: rrfs_3denvar_glbens ..............   Passed  613.72 sec
5/9 Test #9: global_enkf ......................   Passed  1001.15 sec
6/9 Test #6: rtma .............................   Passed  1275.68 sec
7/9 Test #3: global_4denvar ...................   Passed  1612.58 sec
8/9 Test #2: global_4dvar .....................   Passed  1686.18 sec
9/9 Test #1: global_3dvar .....................   Passed  1910.31 sec

100% tests passed, 0 tests failed out of 9

Total Test time (real) = 1910.32 sec

@ShunLiu-NOAA
Copy link
Contributor

Test project /lfs/h2/emc/lam/noscrub/emc.lam/Shun.Liu/gsi_regression/GSI/build
Start 1: global_3dvar
1/9 Test #1: global_3dvar ..................... Passed 1736.38 sec
Start 2: global_4dvar
2/9 Test #2: global_4dvar ..................... Passed 1804.91 sec
Start 3: global_4denvar
3/9 Test #3: global_4denvar ................... Passed 1612.48 sec
Start 4: hwrf_nmm_d2
4/9 Test #4: hwrf_nmm_d2 ...................... Passed 726.46 sec
Start 5: hwrf_nmm_d3
5/9 Test #5: hwrf_nmm_d3 ...................... Passed 494.14 sec
Start 6: rtma
6/9 Test #6: rtma ............................. Passed 1628.39 sec
Start 7: rrfs_3denvar_glbens
7/9 Test #7: rrfs_3denvar_glbens .............. Passed 605.17 sec
Start 8: netcdf_fv3_regional
8/9 Test #8: netcdf_fv3_regional .............. Passed 483.08 sec
Start 9: global_enkf
9/9 Test #9: global_enkf ...................... Passed 1102.38 sec

100% tests passed, 0 tests failed out of 9

@ShunLiu-NOAA ShunLiu-NOAA merged commit c56d7bc into NOAA-EMC:develop Sep 30, 2023
4 checks passed
@TingLei-NOAA TingLei-NOAA mentioned this pull request Oct 17, 2023
6 tasks
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

Successfully merging this pull request may close these issues.

Undefined values found in radar reflectivity direct DA
6 participants