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

Add NSSL 2-Moment Microphysics #478

Closed
wants to merge 20 commits into from

Conversation

tsupinie
Copy link

This pull request adds NSSL microphysics.

Major contributions:

  1. Add NSSL 2-moment microphysics code to CCPP
  2. Add NSSL reflectivity to the hourly-max diagnostics

Minor contributions:

  1. Update .gitignore to clean up the git status listings

@grantfirl
Copy link
Collaborator

Associated PRs:
NOAA-EMC/fv3atm#144

@grantfirl
Copy link
Collaborator

@tsupinie Did you do any work to get this to work in the SCM? No problem if not, just curious. If so, though, a PR there would be useful.

@@ -1,2 +1,11 @@
# Auto-generated include files in Fortran code
*.inc
*_cap.F90
Copy link
Collaborator

Choose a reason for hiding this comment

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

This is helpful, thanks!

@@ -232,7 +232,7 @@ subroutine GFS_MP_generic_post_run(im, levs, kdt, nrcm, ncld, nncl, ntcw, ntrac,
!
! HCHUANG: use new precipitation type to decide snow flag for LSM snow accumulation

if (imp_physics /= imp_physics_gfdl .and. imp_physics /= imp_physics_thompson) then
if (imp_physics /= imp_physics_gfdl .and. imp_physics /= imp_physics_thompson .and. imp_physics /= imp_physics_nssl) then
Copy link
Collaborator

Choose a reason for hiding this comment

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

@tsupinie This isn't a problem with your PR or even a question for you necessarily, but I'm just flagging this for the future, perhaps. Should this IF statement be posed in terms of the positive -- that is, what MP schemes should this run for? Only those schemes without explicit frozen hydrometeors?

@@ -82,7 +82,7 @@ end subroutine GFS_PBL_generic_pre_finalize
subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, &
ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, &
ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef, trans_aero, ntchs, ntchm, &
imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, &
imp_physics, imp_physics_gfdl, imp_physics_nssl, imp_physics_thompson, imp_physics_wsm6, &
Copy link
Collaborator

Choose a reason for hiding this comment

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

This new variable doesn't appear to be used in this scheme.

@@ -209,6 +209,14 @@
type = integer
intent = in
optional = F
[imp_physics_nssl]
Copy link
Collaborator

Choose a reason for hiding this comment

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

Same comment as above. Please remove if this variable is not actually used.

@@ -864,7 +864,8 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input
! clouds, cldsa, mtopa, mbota, de_lgth) ! --- outputs
endif

elseif(Model%imp_physics == 6 .or. Model%imp_physics == 15) then
elseif(Model%imp_physics == 6 .or. Model%imp_physics == 15 .or. &
Model%imp_physics == 17) then
Copy link
Collaborator

Choose a reason for hiding this comment

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

@tsupinie Again, not your problem to clean this up, but just flagging for a future cleanup. These tests should use index variables rather than hardcoded integers. (This problem obviously existed before this PR!)

@@ -630,7 +630,7 @@ end subroutine GFS_suite_interstitial_4_finalize
!! \htmlinclude GFS_suite_interstitial_4_run.html
!!
subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_total, ntrac, ntcw, ntiw, ntclamt, &
ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, &
ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_nssl, imp_physics_thompson, &
Copy link
Collaborator

Choose a reason for hiding this comment

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

imp_physics_nssl doesn't appear to be used in GFS_suite_interstitial_4_run.

@@ -1657,6 +1665,14 @@
type = integer
intent = in
optional = F
[imp_physics_nssl]
Copy link
Collaborator

Choose a reason for hiding this comment

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

This variable doesn't appear to be used in this scheme.

Copy link
Author

Choose a reason for hiding this comment

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

The unused variables were removed from the interstitial schemes. If my interpretation is correct, these interstitials are used for feeding the deep and shallow convective schemes the correct microphysical variables. We weren't using deep or shallow convective schemes for our purposes, so we didn't need to use these, but maybe someone in the future will want to?

!
! constants
!
real, parameter :: cp608 = 0.608 ! constant used in conversion of T to Tv
Copy link
Collaborator

Choose a reason for hiding this comment

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

Not a huge deal, but it is best practice for some of these constants (if they exist in the CCPP from host models) to come through the argument list to maintain consistency with other physics. Off the top of my head, cp608, pi, and gr should be available from the CCPP/host.

Copy link
Author

Choose a reason for hiding this comment

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

Pulled some constants from CCPP.

real, parameter :: cbwbolton = 29.65 ! constants for Bolton formulation
real, parameter :: cawbolton = 17.67

real, parameter :: tfr = 273.15, tfrh = 233.15
Copy link
Collaborator

Choose a reason for hiding this comment

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

Some of these constants as well (e.g., tfr, cp, rd, rw, cv?)

ENDIF
ELSEIF ( ipconc >= 6 ) THEN
write(0,*) 'NSSL microphysics has not been compiled for 3-moment. Sorry.'
STOP
Copy link
Collaborator

Choose a reason for hiding this comment

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

CCPP schemes are not supposed to STOP the model. Ideally, you should pass in the CCPP error message and error flag to whichever subroutines are currently stopping the model. Set the CCPP error message to what is currently being written out, set the error flag to something other than zero, and return. In the calling routine, check for the error flag different than zero and return all the way up the CCPP entry point, which should also return. The CCPP framework then writes out the error message and signalas the host to stop gracefully. I searched for STOP in this PR and found something like 20 instances. All should be handled like I describe for CCPP-compliancy.

Copy link
Author

Choose a reason for hiding this comment

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

I added the errflag and errmsg arguments to the NSSL routines to fix this.

ngscnt = ngscnt + 1
igs(ngscnt) = ix
kgs(ngscnt) = kz
if ( ngscnt .eq. ngs ) goto 1100
Copy link
Collaborator

@grantfirl grantfirl Jul 21, 2020

Choose a reason for hiding this comment

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

Not a huge deal, but CCPP schemes shouldn't rely on GOTOs. A search showed several apparently active GOTO statements.

Copy link
Author

Choose a reason for hiding this comment

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

I've removed some of the GOTO statements. Mostly the ones that involve breaking out of loops early. Currently testing to make sure I didn't mess anything up.

I didn't remove all of them because there's a section that's kind of a rat's nest of GOTOs, and I'm not confident I can do it without messing up the logic. I'd say that's a project for the original author of the scheme.

mixphase = 0
ihvol = 1

nssl_params = 0.
Copy link
Collaborator

Choose a reason for hiding this comment

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

These could be set via namelist parameters if you want. Let me know if you would like help setting that up.

Copy link
Author

Choose a reason for hiding this comment

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

A couple considerations here:

  1. mixphase is hard-coded to be 0 in the original WRF implementation, and I'm not sure what other valid values are.
  2. ipconc basically controls the number of moments (the WRF implementation has a 1-moment option), and ihvol controls whether or not hail is present. The other options in nssl_params provide options, most of which are only used if ipconc == 0 (meaning we're running the 1-moment version). WRF hides all this from the user, and uses different values of mp_physics to trigger the different settings. I could do something similar here, but it would take more testing.


! Convert specific humidity/moist mixing ratios to dry mixing ratios
qv_mp = spechum/(1.0_kind_phys-spechum)
qc_mp = qc/(1.0_kind_phys-spechum)
Copy link
Collaborator

Choose a reason for hiding this comment

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

You may want to double-check whether this conversion is necessary. I thought that water vapor is carried as specific humidity, but that other species are mixing ratios. There may be inconsistencies in other CCPP schemes vis-a-vis specific humidity vs mixing ratios too.

Copy link
Author

Choose a reason for hiding this comment

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

Yeah, I'm not 100% sure on that, either. I'm IIRC, we pulled this from the Thompson implementation, which (still) does translate between moist and dry hydrometeor mixing ratios. After digging around in the model code for a bit, I didn't find anything that definitively said one way or the other whether hydrometeors masses are relative to moist or dry air.

qg = qg_mp/(1.0_kind_phys+qv_mp)
qh = qh_mp/(1.0_kind_phys+qv_mp)

prcp = max(0.0, delta_rain_mp/1000.0_kind_phys)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Dividing by water density? Might want to comment just for clarity.

Copy link
Author

Choose a reason for hiding this comment

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

Actually, converting mm to m. Added a comment to clarify.

@@ -351,6 +351,8 @@ subroutine cld_init &
print *,' --- WSM6 cloud microphysics'
elseif (imp_physics == 10) then
print *,' --- MG cloud microphysics'
elseif (imp_physics == 17) then
Copy link
Collaborator

Choose a reason for hiding this comment

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

Same comment WRT hard-coded integers, but you're obviously following precedent here...

Copy link
Collaborator

@grantfirl grantfirl left a comment

Choose a reason for hiding this comment

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

@tsupinie This is excellent work, especially since it was done independently! I left some comments about a few things that would be good to fix for "true" CCPP-compliancy, especially error handling internal to the scheme. There were also a few unnecessary variable additions. Writing another CCPP-compliant driver layer on top of the existing code isn't necessarily ideal, but totally understandable given 1) the complexity of this particular scheme and 2) the underlying scheme code is likely used as-is in other models (e.g., WRF). Should this scheme become operational, it might be advantageous computationally to make the code in module_mp_nssl_2mom.F90 CCPP-compliant instead, but that certainly isn't necessary at the moment.

Other than that, successful regression testing is required before approval to master.

@tsupinie
Copy link
Author

tsupinie commented Aug 8, 2020

@grantfirl, no, I haven't tried in the SCM.

Copy link
Collaborator

@climbfuji climbfuji left a comment

Choose a reason for hiding this comment

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

@tsupinie @grantfirl where are we with this PR? We should probably add a regression test for this new microphysics scheme. We also need accompanying fv3atm and ufs-weather-model PRs. I can take care of those and set up the regression test once we are satisfied with the ccpp-physics PR.

real(kind_phys), intent(in ) :: con_g
real(kind_phys), intent(in ) :: con_rd
! Hydrometeors
real(kind_phys), intent(inout) :: spechum(1:ncol,1:nlev)
Copy link
Collaborator

Choose a reason for hiding this comment

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

We are currently in the process to declare all interface variable arrays (i.e. those that come in via the argument list) as assumed-size arrays to enable bounds checking by the compiler. We should do this right away for this PR, too, and use

real(kind_phys),           intent(inout) :: spechum(:,:)

instead of

real(kind_phys),           intent(inout) :: spechum(1:ncol,1:nlev)

and similar for all other arrays that come in via the argument list.

@grantfirl
Copy link
Collaborator

@tsupinie @grantfirl where are we with this PR? We should probably add a regression test for this new microphysics scheme. We also need accompanying fv3atm and ufs-weather-model PRs. I can take care of those and set up the regression test once we are satisfied with the ccpp-physics PR.

@climbfuji From my perspective, it really just needs testing of some sort. I was planning on at least testing in the SCM since that is easy for me to do, but someone (DTC?) also needs to run FV3 RTs as you pointed out.

@tsupinie Thanks for addressing some of my comments. One thing that I noticed with dc275b2 was that you passed in physical constants with the physcons module. We're trying to get away from this. Instead, these same constants should go through the scheme's argument list and propagated down to child subroutines if necessary. For FV3, the constants that are available through the CCPP are found at the end of GFS_typedefs.meta. For this model, they are being imported from the physcons module, but other host models that might use this scheme will not have a physcons module, so if they get passed in via the argument list using CCPP metadata, this allows for other host models to handle their constants in their own way.

@MicroTed
Copy link
Contributor

MicroTed commented Sep 3, 2020

@tsupinie I wish I had known about this earlier! I've spent some time recently adapting NSSL scheme for CCPP and SCM (both versions 3 and 4). I have github forks for these, but need to push my working branches. Let's compare notes on the implementation. Between the two of us we seem to cover SCM and FV3? I also worked from the Thompson example, so there is probably a good bit of congruence.

@grantfirl
I have made semi-successful tests with SCM 3 and 4. I say semi because I need to understand the forcing better for the twpice case. It seems to keep adding cloud droplet concentration, so I get tiny droplets that never make rain. (And resulting droplet mixing ratios of 30 g/kg!) In SCM v3 it runs to conclusion, but SCM v4 I get a segfault (memory access) on the next-to-last output time. There are some code variable details I don't fully understand yet, and I have an email drafted for a support ticket I opened last March (right at v4 release time, of course).

Your complaints about the code details can be directed to me. As Tim noted, that section with goto statements is an example needing more effort to update 1980s code than is probably worthwhile.

@tsupinie
Copy link
Author

@grantfirl Okay, I've modified the implementation to use the constants from CCPP.

@climbfuji
Copy link
Collaborator

@tsupinie @grantfirl this PR is aging off quickly. Shall we close it and open a new one when the changes are ready?

@climbfuji
Copy link
Collaborator

Not having heard of or seen any updates since September, I am closing this PR. Please reopen it or open a new PR when the code changes are ready.

@climbfuji climbfuji closed this Feb 18, 2021
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.

4 participants