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

Refine PM2.5 DA for the RRFS_SD model #609

Merged
merged 15 commits into from
Sep 29, 2023

Conversation

hongli-wang
Copy link
Collaborator

@hongli-wang hongli-wang commented Aug 11, 2023

Description

Refine the PM2.5 DA for the RRFS_SD model by use of veg_type, which is used to decide whether the obs is in urban area or not. Different thresholds for innovations outside/inside urban areas will be used. Add new namelist parameters, such as threshold for innovations, anowbufr type
Read in station terrain height, PM10 et al if extended BUFR format for anow air quality data is used

Fixes #606

Type of change

Please delete options that are not relevant.

  • Bug fix (non-breaking change which fixes an issue)
  • [X ] New feature (non-breaking change which adds functionality)
  • Breaking change (fix or feature that would cause existing functionality to not work as expected)
  • This change requires a documentation update

How Has This Been Tested?

Checklist

  • [X ] 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/22/2023. If this PR is not merged into develop by this date, the PR will be closed and returned to the developer.

@hu5970
Copy link
Collaborator

hu5970 commented Aug 25, 2023

@hongli-wang Could you suggest two reviewers for this PR?

@hongli-wang
Copy link
Collaborator Author

hongli-wang commented Aug 25, 2023

@hongli-wang Could you suggest two reviewers for this PR?

@hu5970 Previously, Ting and Cory reviewed the initial development of PM2.5 DA for RRFS_SD. They are the best to review this update.

@RussTreadon-NOAA
Copy link
Contributor

Cory is not available as a peer reviewer. Cory is a member of the Handling Review team. Handling Review team members process PRs. We need two peer reviewers not from the Handling Review team.

@hongli-wang
Copy link
Collaborator Author

Cory is not available as a peer reviewer. Cory is a member of the Handling Review team. Handling Review team members process PRs. We need two peer reviewers not from the Handling Review team.

@RussTreadon-NOAA Thanks for letting me know. Guoqing @guoqing-noaa can help to review the code. Thanks.

@CoryMartin-NOAA
Copy link
Contributor

@hongli-wang can you start by running the GSI regression tests on Hera? I will also run them on WCOSS2 once we know the all pass on Hera. Thanks!

@hongli-wang
Copy link
Collaborator Author

hongli-wang commented Aug 28, 2023

Cory is not available as a peer reviewer. Cory is a member of the Handling Review team. Handling Review team members process PRs. We need two peer reviewers not from the Handling Review team.

@RussTreadon-NOAA Thanks for letting me know. Guoqing @guoqing-noaa can help to review the code. Thanks.

@hongli-wang
Copy link
Collaborator Author

Cory is not available as a peer reviewer. Cory is a member of the Handling Review team. Handling Review team members process PRs. We need two peer reviewers not from the Handling Review team.

@RussTreadon-NOAA Thanks for letting me know. Guoqing @guoqing-noaa can help to review the code. Thanks.

Youhua ([email protected]) can also review the code. Thanks.

@@ -287,7 +292,9 @@ subroutine init_chem

berror_chem=.false.
berror_fv3_cmaq_regional=.false. ! Set .true. to use berror for fv3_cmaq_regional, whose cv has 10 characters
berror_fv3_sd_regional=.false. ! Set .true. to use berror for rrfs_sd model, whose cv has 10 characters
Copy link
Contributor

Choose a reason for hiding this comment

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

Would add an explanation for =.false.?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Would add an explanation for =.false.?

Thanks for the comments. I will revise soon.

oneobtest_chem=.false.
anowbufr_ext=.false.
Copy link
Contributor

Choose a reason for hiding this comment

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

documentation for this variable

if ( iret > 0 .and. (conc < conc_missing ) .and. &
(conc >= zero)) then
if(indata(nxob) > r360)cycle
if(indata(nyob) > 0.5_r_kind*r360)cycle
Copy link
Contributor

Choose a reason for hiding this comment

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

why don't directly use r180 (if r180 is not defined, it could be defined along with r360). The extra r180 will avoid the multiplication for 0.5*r360 for each obs.

@hongli-wang
Copy link
Collaborator Author

@hongli-wang can you start by running the GSI regression tests on Hera? I will also run them on WCOSS2 once we know the all pass on Hera. Thanks!

A bit busy this week. I will do it as soon as I can.

@CoryMartin-NOAA
Copy link
Contributor

@hongli-wang no big rush, we have a few weeks.

@@ -103,7 +108,7 @@ module chemmod
real(r_kind),parameter :: pm2_5_teom_max=900.0_r_kind !ug/m3
!some parameters need to be put here since convinfo file won't
!accomodate, stands for maximum realistic value of surface pm2.5
real(r_kind),parameter :: pm10_teom_max=150.0_r_kind !ug/m3
real(r_kind),parameter :: pm10_teom_max=3000.0_r_kind !ug/m3
Copy link
Contributor

Choose a reason for hiding this comment

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

This comment is not related to the code. But just curious whey bumping the max value from 150 to 3000?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

This comment is not related to the code. But just curious whey bumping the max value from 150 to 3000?

@guoqing-noaa The default value may not consider the heavy fire event. The value of 3000 is a reasonable PM10 value when a heavy fire event or dust storm takes place.

obstr='TPHR QCIND COPOPM ELV COPOPM10 COPOCO'
anowbufr=.true.
write(6,*)'READ_PM2_5: AIRNOW data type, subset=',subset
end if
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggest reverse the if and else block, i.e: if anowbufr_ext then; else ... (although they function the same, the latter reads more straightforward :) )

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Accepted. Thanks.

src/gsi/read_anowbufr.f90 Show resolved Hide resolved
obstr='TPHR QCIND COPOPM ELV COPOPM10 COPOCO'
anowbufr=.true.
write(6,*)'READ_PM10: AIRNOW data type, subset=',subset
end if
Copy link
Contributor

Choose a reason for hiding this comment

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

same comments as the above

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Accepted. Thanks.

if (trim(obstype)=='pm2_5') indata(ncopopm)=indata_b(3)
if (trim(obstype)=='pm10') indata(ncopopm)=indata_b(5)
site_elev = indata_b(4)
end if
Copy link
Contributor

Choose a reason for hiding this comment

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

remove two spaces to align the indentation

end if
if (pm25wc_ges(2) >= 1.0_r_kind) then
if (pm25wc_ges(2) >= 2.0_r_kind) then
Copy link
Contributor

Choose a reason for hiding this comment

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

Is the change from 1.0 to 2.0 a bug fix or a threshold change? For the latter, will it be better to use a parameter to control this?

src/gsi/setuppm2_5.f90 Outdated Show resolved Hide resolved
@@ -90,7 +93,9 @@ module chemmod

logical :: aero_ratios

logical :: oneobtest_chem,diag_incr,berror_chem,berror_fv3_cmaq_regional
logical :: oneobtest_chem,diag_incr,berror_chem
logical :: berror_fv3_cmaq_regional,berror_fv3_sd_regional
Copy link
Contributor

Choose a reason for hiding this comment

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

I searched this PR, it looks like we don't need separate berror_fv3_sd_regional and berror_fv3_cmaq_regional?

berror_fv3_sd_regional is only used at Line 460 of src/gsi/m_berror_stats_reg.f90 and there is no difference if we use berror_fv3_cmaq_regional instead. Is it possible to just use berror_fv3_cmaq_regional, or combine them to one logical variable, such as berror_fv3_aero_regional?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

The current ones make it more clear to know which B file will be used, and the previous retros runs can repeat without any changes of GSI namelist parameter. This is good to users.
Your suggestion will make the changes be general and good in code development.
So what do you think?

Copy link
Contributor

Choose a reason for hiding this comment

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

Do you plan to use different BEC for cmaq and sd? Or will they be the same as they are right now?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@guoqing-noaa
Different B files are used. The one for CMAQ has 70 aerosol species whereas
rrfs-ad only has 3 tracers. I think a separate parameter is good for users. I have no problem in combining then if that is the way you prefer. Thanks.

Copy link
Contributor

Choose a reason for hiding this comment

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

Got it. Thanks for the explanation!


character(len=max_varname_length), dimension(naero_smoke_fv3), parameter :: &
aeronames_smoke_fv3=[character(len=max_varname_length) :: 'smoke','dust' ]
aeronames_smoke_fv3=[character(len=max_varname_length) :: 'smoke','dust','coarsepm']
Copy link
Contributor

Choose a reason for hiding this comment

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

"coarsepm" is added here. However it is never used. And we have to modify Line 282 of src/gsi/setuppm2_5.f90
to remove the "coarsepm" contribution

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, you are correct. This was done by -1 in the do loop. For example,
do i=2,naero_smoke_fv3-1 ! remove contribution from coarsepm

Copy link
Contributor

Choose a reason for hiding this comment

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

Right, if we will remove coarsepm in src/gsi/setuppm2_5.f90 (line 281), why would we add it here?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@guoqing-noaa
This is because the dust PM10/AOD DA is not included in this PR. But we do consider EnVar to update the coarsepm, possibly in a next PR with FED DA.

@@ -76,10 +78,11 @@ module chemmod
logical :: luse_deepblue
character(len=max_varname_length) :: crtm_aerosol_model,crtm_aerosolcoeff_format,crtm_aerosolcoeff_file
integer(i_kind) :: aod_qa_limit ! qa >= aod_qa_limit will be retained
integer(i_kind) :: icvt_cmaq_fv3
integer(i_kind) :: icvt_cmaq_fv3,pm10_da_scheme
Copy link
Contributor

Choose a reason for hiding this comment

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

pm10_da_scheme is not used anywhere in this PR.

@hongli-wang
Copy link
Collaborator Author

@TingLei-NOAA @guoqing-noaa
Hi Ting and Guoqing,
Thanks for your comments. I have taken almost all of your suggestions and modified the codes. Please review the code again. Thanks. Hongli

@TingLei-NOAA
Copy link
Contributor

TingLei-NOAA commented Sep 8, 2023

@hongli-wang just maybe last two question: did you run your test for this PR with GSI at debug mode? :). and did you finish the regression tests?

@hongli-wang
Copy link
Collaborator Author

@hongli-wang just maybe last two question: did you run your test for this PR with GSI at debug mode? :). and did you finish the regression tests?

@TingLei-NOAA yes tested my code on debug mode. Not finish regression test yet, and I will do it.

@TingLei-NOAA
Copy link
Contributor

@hongli-wang Thanks. Then, the success of your regression tests are the only thing I need to know to give my approval for this PR. Thanks a lot for your efforts for keeping improving the codes.

@hongli-wang
Copy link
Collaborator Author

@hongli-wang can you update to the head of GSI develop? I think the regression test issues I'm seeing on both Hera and Orion are related to issues with changes from #571 not being included here

@CoryMartin-NOAA yes, just rebased my code, will have some tests and run a regression test.

@hongli-wang
Copy link
Collaborator Author

@hongli-wang can you update to the head of GSI develop? I think the regression test issues I'm seeing on both Hera and Orion are related to issues with changes from #571 not being included here

@CoryMartin-NOAA yes, just rebased my code, will have some tests and run a regression test.

@CoryMartin-NOAA @TingLei-NOAA

Regression tests failed due to GSI can't find correct intel lib, for example:
/scratch1/BMC/wrfruc/hwang/ctest/tmpreg_rrfs_3denvar_glbens/rrfs_3denvar_glbens_loproc_updat/gsi.x: error while loading shared libraries: libmkl_intel_lp64.so.2
Dir on Hera:
/scratch1/BMC/wrfruc/hwang/ctest

Any suggestions?

@CoryMartin-NOAA
Copy link
Contributor

@hongli-wang that's because the modulefiles changes in develop recently. Make sure your develop build is also up to date now that you've pulled in the recent updates.

I'm running the tests on Orion and Cactus now.

@hongli-wang
Copy link
Collaborator Author

@hongli-wang that's because the modulefiles changes in develop recently. Make sure your develop build is also up to date now that you've pulled in the recent updates.

I'm running the tests on Orion and Cactus now.

@CoryMartin-NOAA
Both develop and my branch are up to date and recompiled. Is there anything else I can do to keep the build to to date?
Thanks for performing regression tests on Orion and Cactus.

@CoryMartin-NOAA
Copy link
Contributor

Rerunning on Cactus:
All tests pass except hwrf_nmm_d3 which fails the scalability test and rtma which fails the runtime maximum for a couple of the cases.

On Orion:

The following tests FAILED:
	  2 - [=[global_4dvar]=] (Failed)
	  5 - [=[hwrf_nmm_d3]=] (Failed)
	  6 - [=[rtma]=] (Failed)
	  7 - [=[rrfs_3denvar_glbens]=] (Failed)
	  8 - [=[netcdf_fv3_regional]=] (Failed)

The 4DVar test failed the scalability test.
The HWRF test barely failed the scalability test.
The RTMA test also failed the scalability test.
The RRFS test exceeded the runtime for one test by 3 seconds.
Finally, the FV3 regional test also failed for scalability.

I think given these results, the tests "pass"

@hongli-wang
Copy link
Collaborator Author

Rerunning on Cactus: All tests pass except hwrf_nmm_d3 which fails the scalability test and rtma which fails the runtime maximum for a couple of the cases.

On Orion:

The following tests FAILED:
	  2 - [=[global_4dvar]=] (Failed)
	  5 - [=[hwrf_nmm_d3]=] (Failed)
	  6 - [=[rtma]=] (Failed)
	  7 - [=[rrfs_3denvar_glbens]=] (Failed)
	  8 - [=[netcdf_fv3_regional]=] (Failed)

The 4DVar test failed the scalability test. The HWRF test barely failed the scalability test. The RTMA test also failed the scalability test. The RRFS test exceeded the runtime for one test by 3 seconds. Finally, the FV3 regional test also failed for scalability.

I think given these results, the tests "pass"

@CoryMartin-NOAA
That is a good news. Thanks!
I have the issue that GSI can't find the lib to run on Hera. The same issue if I ran GSI on Jet new.

Do you think we can move forward, ask Ting and Guoqing to review and approve this PR?

Thank again!
Hongli

@CoryMartin-NOAA
Copy link
Contributor

@guoqing-noaa @TingLei-NOAA can you re-review this? We can double check the regression tests at the end but I think all we are seeing here is machine instability. Thanks @hongli-wang for your work on this so far!

@CoryMartin-NOAA
Copy link
Contributor

GitHub says 76 files changed, I think this needs updated to the head of develop again before it can be accurately reviewed.

@hongli-wang
Copy link
Collaborator Author

GitHub says 76 files changed, I think this needs updated to the head of develop again before it can be accurately reviewed.

@CoryMartin-NOAA
I noticed this too. It was caused my last rebase (not sure why it happened). I will rebase/merge the code again.

@hongli-wang
Copy link
Collaborator Author

GitHub says 76 files changed, I think this needs updated to the head of develop again before it can be accurately reviewed.

@CoryMartin-NOAA I noticed this too. It was caused my last rebase (not sure why it happened). I will rebase/merge the code again.

@CoryMartin-NOAA
Remerged code to EMC develop branch and delete a unused namelist variable. Could you please help to run regression tests if it is required? I can try if the reg tests work on Hera tomorrow. Thanks.

@hongli-wang
Copy link
Collaborator Author

hongli-wang commented Sep 26, 2023 via email

@hongli-wang
Copy link
Collaborator Author

hongli-wang commented Sep 26, 2023 via email

@CoryMartin-NOAA
Copy link
Contributor

@hongli-wang at this point, I'll re-run the regression tests on WCOSS once we have approvals from the reviewers.

Copy link
Contributor

@guoqing-noaa guoqing-noaa left a comment

Choose a reason for hiding this comment

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

LGTM. Thanks for the discussions and changes!

@hongli-wang
Copy link
Collaborator Author

@hongli-wang at this point, I'll re-run the regression tests on WCOSS once we have approvals from the reviewers.

@CoryMartin-NOAA
Ting and Guoqing has approved this PR. Could you please run the regression test again?

Thanks,
Hongli

@CoryMartin-NOAA
Copy link
Contributor

@hongli-wang already submitted them this morning on Cactus. I'll post when they are finished.

@CoryMartin-NOAA
Copy link
Contributor

Only 2 regression tests fail on Cactus, and they are due to scalability, which are nonfatal. I think this is good to go now.

@hongli-wang
Copy link
Collaborator Author

Only 2 regression tests fail on Cactus, and they are due to scalability, which are nonfatal. I think this is good to go now.

@CoryMartin-NOAA
Thank you!

@CoryMartin-NOAA CoryMartin-NOAA merged commit ba5a2ca into NOAA-EMC:develop Sep 29, 2023
4 checks passed
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.

Refine PM2.5 DA for RRFS-SD
6 participants