Option to roll back Thompson MP to WRFv3.8.1 (RAPv5/HRRRv4), add stochastic perturbations code

[climbfuji]

This PR implements an option to roll back the version of the Thompson microphysics scheme ( module_mp_thompson.F90) to match WRFv3.8.1, which is also the version used in RAPv5/HRRRv4. The switch is currently made by manually activating a preprocessor definition at the top of module_mp_thompson.F90:

! DH* 2020-06-05
! Use the following preprocessor directive to roll back
! to the WRFv3.8.1, used in RAPv5/HRRRv4 for more reasonable
! representation of mesoscale storms and reflectivity values
!#define WRF381

For the purpose of merging this PR, this option will remain commented out for running the regression tests (to make sure that the modifications do not change the answer). Manual tests will be run with the WRFv3.8.1 version to make sure that the code runs successfully when using the older version of the scheme.

In addition to the roll back, this PR adds the stochastic physics perturbations to the graupel intercept parameter, the cloud water shape parameter, and the number concentration of nucleated aerosols. This code is taken from newer versions of module_mp_thompson.F90 and retrofitted here. This capability will remain untested for the time being, and a guard has been put in place to prevent its use without testing.

Associated PRs:
#34
NOAA-GSL/fv3atm#35
NOAA-GSL/ufs-weather-model#28

For regression testing information, see NOAA-GSL/ufs-weather-model#28.

[climbfuji]

@johnbrownsbody @joeolson42 @tanyasmirnova @evankalina @jaymes-kenyon please see this PR for adding (but not yet making default) the WRF v3.8.1 flavor of Thompson MP, as well as the stochastic perturbations. My plan is to merge this with the WRF v4 flavor default (so that all the regression tests give b4b identical results), and then give GSL developers time to test how the WRFv3.8.1 flavor performs in this setting, compared to the WRF v4. We can make the switch in a follow-up commit.

[climbfuji]

Note: a newer version of the Thompson MP scheme that I received from Greg in January this year has a slightly updated version of the random perturbation code, but most of the changes compared to the code that @evankalina ported over are commented out (rand4). Since we are not using those random perturbations for the moment, I suggest to keep the version that Evan ported over and update to a more mature version later if needed.

[tanyasmirnova]

Dom,
I think that on Line 5306 the limit should 4.99e-6 for re_qc and 9.99e-6 on line 5321 for re_qi.

[climbfuji]

Dom,
I think that on Line 5306 the limit should 4.99e-6 for re_qc and 9.99e-6 on line 5321 for re_qi.

WRFV381 doesn't have any initial value here (a bug?) the initial value is for WRFV4 (see #ifndef WRFv381 not ifdef), as in the original code. Attached a screenshot of the diff between the RAP/HRRR version I got from Jaymes (left) and the combined version of this PR.

[tanyasmirnova]

Also, I remember that Greg recommended to use 1000. in the lines below for global applications:
!tgs - fix from Greg
! if (rg(kts).gt.R1 * 10.) &
if (rg(kts).gt.R1 * 1000.) &
This code is from my implementation of Thompson MP in FIM.

[tanyasmirnova]

Sorry, I did not notice that it is "ifndef" and not "ifdef".
Maybe it is OK, because we set initial values before the call to calc_effectRad.
do k = kts, kte
re_qc1d(k) = 2.49E-6
re_qi1d(k) = 4.99E-6
re_qs1d(k) = 9.99E-6
enddo
But they should be consistent with the actual values used in the computations of re_qc, re_qi and re_qs.

[johnbrownsbody]

Also, I remember that Greg recommended to use 1000. in the lines below for global applications:
!tgs - fix from Greg
! if (rg(kts).gt.R1 * 10.) &
if (rg(kts).gt.R1 * 1000.) &
This code is from my implementation of Thompson MP in FIM.

Thanks, Tanya. I remember this. This subroutine is for calculating reflectivity diagnostic,
which I don't think we were doing in FIM, unless I am misinterpreting your comment.

[johnbrownsbody]

Tanya, If I am following your concern correctly, this is in the diagnostic reflectivity calculation. I suspect it is a trap to prevent too small values of mixing ratio from getting used in the subsequent reflectivity calculation, which involves exponentials.

…

On Fri, Jun 5, 2020 at 3:29 PM tanyasmirnova ***@***.***> wrote: Sorry, I did not notice that it is "ifndef" and not "ifdef". Maybe it is OK, because we set initial values before the call to calc_effectRad. do k = kts, kte re_qc1d(k) = 2.49E-6 re_qi1d(k) = 4.99E-6 re_qs1d(k) = 9.99E-6 enddo But they should be consistent with the actual values used in the computations of re_qc, re_qi and re_qs. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#34 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AFFQ33CIIVQZEM3M76O5GL3RVFPRVANCNFSM4NTAZ4PA> .

[evankalina]

Dom, does rand1 need to be added to the argument list for the calc_refl10cm subroutine and declared as a local variable there?

Regarding rand4, that was another random perturbation I added to the original version of the code for some exploratory work that I was doing. Greg disagreed and that is probably why it is commented out in the new version. I don't think you need to include rand4 here.

[climbfuji]

Dom, does rand1 need to be added to the argument list for the calc_refl10cm subroutine and declared as a local variable there?

I don't see where this is happening. Line 5383, I believe, declares it as an intent(in) variable for calc_refl10cm.

Regarding rand4, that was another random perturbation I added to the original version of the code for some exploratory work that I was doing. Greg disagreed and that is probably why it is commented out in the new version. I don't think you need to include rand4 here.

Thanks for letting me know, I'll skip it for the time being.

[evankalina]

@climbfuji Sorry, I was looking at an older version of the file without realizing it. The correct version looks good to me.

[jaymes-kenyon]

This looks good to me. If I am seeing this correctly, it looks like the effective radii initializations ("re_") will include "ifdef" conditionals, which is probably wise. This should hopefully come very close to replicating HRRRv4's microphysics (as derived from WRFV3.8.1).

[climbfuji]

Tanya, If I am following your concern correctly, this is in the diagnostic reflectivity calculation. I suspect it is a trap to prevent too small values of mixing ratio from getting used in the subsequent reflectivity calculation, which involves exponentials.
…
On Fri, Jun 5, 2020 at 3:29 PM tanyasmirnova @.***> wrote: Sorry, I did not notice that it is "ifndef" and not "ifdef". Maybe it is OK, because we set initial values before the call to calc_effectRad. do k = kts, kte re_qc1d(k) = 2.49E-6 re_qi1d(k) = 4.99E-6 re_qs1d(k) = 9.99E-6 enddo But they should be consistent with the actual values used in the computations of re_qc, re_qi and re_qs. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#34 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFFQ33CIIVQZEM3M76O5GL3RVFPRVANCNFSM4NTAZ4PA .

To be clear, we could (maybe should) implement a similar definition for WRFV381, but with the corresponding initial values,

[climbfuji]

@johnbrownsbody @joeolson42 @tanyasmirnova @evankalina @jaymes-kenyon please have a close look at the changes in commit 626ec0e. I cleaned up the initialization and bounds for the effective radii by moving those into the routine calc_effectRad. This makes much more sense, because there is only one place where these are defined.

Note that the initial values and bounds that were previously applied in the calling routines matched the internal bounds for WRF v3.8.1, hence inconsistent with those for WRFv4+. Sine the WRF v3.8.1 bounds were tighter and applied afterwards, they effectively overwrote the v4+ bounds. I mimicked and documented this behavior in calc_effectRad.

[tanyasmirnova]

Dom,
The changes in the 626ec0e commit look good to me.

[climbfuji]

Dom,
The changes in the 626ec0e commit look good to me.

Thanks, @tanyasmirnova . I haven't run the usual regression tests yet, and I still need to check the differences between the WRFv4 and WRFv3.8.1 Thompson MP code for our usual C768 November 2020 test case.