Refactor and cleanup of HB PBL scheme

[jimmielin]

Originator(s): @mwaxmonsky (very lightly edited by @jimmielin)

Description (include issue title and the keyword ['closes', 'fixes', 'resolves'] and issue number):

• refactor and cleanup of Holtslag-Boville PBL scheme.

List all namelist files that were added or changed: N/A

List all files eliminated and why: N/A

List all files added and what they do: N/A

List all existing files that have been modified, and describe the changes:
(Helpful git command: git diff --name-status development...<your_branch_name>)

M       phys_utils/atmos_phys_pbl_utils.F90
  - moved minimum kvf limiter to HB logic instead of in common PBL utils code.

M       schemes/holtslag_boville/holtslag_boville_diff.F90
  - various refactoring and cleanup of HB PBL scheme into pure elemental functions.

List all automated tests that failed, as well as an explanation for why they weren't fixed: None (B4B)

Is this an answer-changing PR? If so, is it a new physics package, algorithm change, tuning change, etc?
B4B as tested in CAM against 64_137

If yes to the above question, describe how this code was validated with the new/modified features:

[nusbaume]

Thanks for all of these great changes @mwaxmonsky @jimmielin! I had some clean-up and reorganization requests, but hopefully nothing that will be too big of a burden. Of course if you do run into any problems or issues with my requests then just let me know. Thanks again!

[nusbaume]

I realize this is likely deviating from the actual science text, but can we rename these to something like a_factor and b_factor? Otherwise searching for "a" is going to be a not-super-fun experience.

[nusbaume]

I believe this is the general equation for the brunt-vaisaila frequency, so can we move it to atmos_phys_pbl_utils and add the following documentation?

Suggested change

	pure elemental function calc_brunt_vaisaila_frequency(thv1, thv2, z1, z2, g) result(n2)
	real(kind_phys), intent(in) :: thv1
	real(kind_phys), intent(in) :: thv2
	real(kind_phys), intent(in) :: z1
	real(kind_phys), intent(in) :: z2
	real(kind_phys), intent(in) :: g
	real(kind_phys) :: n2
	n2 = g*2.0_kind_phys*(thv1-thv2)/((thv1+thv2)*(z1-z2))
	end function calc_brunt_vaisaila_frequency
	pure elemental function calc_brunt_vaisaila_frequency(thv1, thv2, z1, z2, g) result(n2)
	! Vallis, G. K. (2006).
	! Atmospheric and Oceanic Fluid Dynamics. Cambridge University Press.
	! https://doi.org/10.1017/CBO9780511790447
	! Page 94, equation 2.225
	real(kind_phys), intent(in) :: thv1 ! (virtual) potential temperature of top layer [K]
	real(kind_phys), intent(in) :: thv2 ! (virtual) potential temperature of bottom layer [K]
	real(kind_phys), intent(in) :: z1 ! height of top layer [m]
	real(kind_phys), intent(in) :: z2 ! height of bottom layer [m]
	real(kind_phys), intent(in) :: g ! gravitational acceleration [m s-2]
	real(kind_phys) :: n2 ! Brunt-Vaisaila frequency [s-2]
	n2 = g*2.0_kind_phys*(thv1-thv2)/((thv1+thv2)*(z1-z2))
	end function calc_brunt_vaisaila_frequency

We might also want to open an issue noting that this function needs unit tests.

[nusbaume]

I also believe this is the general equation for the bulk Richardson number, so can we move it to atmos_phys_pbl_utils and add the following documentation?

Suggested change

	pure elemental function calc_bulk_richardson_number(thv1, thv2, z1, z2, s2, g) result(ri)
	real(kind_phys), intent(in) :: thv1
	real(kind_phys), intent(in) :: thv2
	real(kind_phys), intent(in) :: z1
	real(kind_phys), intent(in) :: z2
	real(kind_phys), intent(in) :: s2
	real(kind_phys), intent(in) :: g
	real(kind_phys) :: ri
	real(kind_phys) :: n2
	n2 = calc_brunt_vaisaila_frequency(thv1, thv2, z1, z2, g)
	ri = n2/s2
	end function calc_bulk_richardson_number
	pure elemental function calc_bulk_richardson_number(thv1, thv2, z1, z2, s2, g) result(ri)
	! https://glossary.ametsoc.org/wiki/bulk-richardson-number/
	!
	! Also found in:
	! Stull, Roland B. An Introduction to Boundary Layer Meteorology. Springer Kluwer Academic Publishers, 1988. Print.
	! DOI: https://doi.org/10.1007/978-94-009-3027-8
	! Equation 5.6.3, page 177
	real(kind_phys), intent(in) :: thv1 ! (virtual) potential temperature of top layer [K]
	real(kind_phys), intent(in) :: thv2 ! (virtual) potential temperature of bottom layer [K]
	real(kind_phys), intent(in) :: z1 ! height of top layer [m]
	real(kind_phys), intent(in) :: z2 ! height of bottom layer [m]
	real(kind_phys), intent(in) :: s2 ! vertical shear squared [s-2]
	real(kind_phys), intent(in) :: g ! gravitational acceleration [m s-2]
	real(kind_phys) :: ri ! bulk Richardson number [1]
	real(kind_phys) :: n2 ! Brunt-Vaisaila frequency [s-2]
	n2 = calc_brunt_vaisaila_frequency(thv1, thv2, z1, z2, g)
	ri = n2/s2
	end function calc_bulk_richardson_number

We should also note that this function needs unit tests in a new issue as well (which could be combined with the Brunt-Vaisaila issue).

[nusbaume]

Actually, I think this equation is technically for the "gradient" Richardson number, so maybe we should rename it from bulk to gradient? There are other functions below which appear to represent a version of the "bulk" number.

[nusbaume]

I believe that this is technically the vertical shear squared, so we should update the function name:

Suggested change

	pure elemental function calc_shear_squared(u1, u2, v1, v2, z1, z2, minimum_velocity_shear_squared) result(s2)
	pure elemental function calc_vertical_shear_squared(u1, u2, v1, v2, z1, z2, minimum_velocity_shear_squared) result(s2)

I might also move this to atmos_phys_pbl_utils as well, and note that it needs unit tests in the follow-up issue.

[nusbaume]

I believe these are both technically the "bulk" Richardson number, so let's name it as such:

Suggested change

	pure elemental function calc_modified_richardson_number_at_height(uh, us, vh, vs, h, zs, thvh, thvs, tlv, ustar, g) result(rih)
	! Modified version of base richardson number equation that incorporates reference level potential temperature
	real(kind_phys), intent(in) :: uh, us, vh, vs, h, zs, thvh, thvs, tlv, ustar, g
	real(kind_phys) :: rih
	real(kind_phys), parameter :: b = 100._kind_phys ! Derivation of b comes from page 253
	real(kind_phys), parameter :: tiny = 1.e-36_kind_phys ! Prevents division by 0
	rih = g * (thvh-tlv)*(h-zs) / &
	(thvs * max((uh-us)**2 + (vh-vs)**2 + b*ustar**2, tiny))
	end function calc_modified_richardson_number_at_height
	pure elemental function calc_richardson_number_at_height(uh, us, vh, vs, h, zs, thvh, thvs, ustar, g) result(rih)
	pure elemental function calc_modified_bulk_richardson_number_at_height(uh, us, vh, vs, h, zs, thvh, thvs, tlv, ustar, g) result(rih)
	! Modified version of base richardson number equation that incorporates reference level potential temperature
	real(kind_phys), intent(in) :: uh, us, vh, vs, h, zs, thvh, thvs, tlv, ustar, g
	real(kind_phys) :: rih
	real(kind_phys), parameter :: b = 100._kind_phys ! Derivation of b comes from page 253
	real(kind_phys), parameter :: tiny = 1.e-36_kind_phys ! Prevents division by 0
	rih = g * (thvh-tlv)*(h-zs) / &
	(thvs * max((uh-us)**2 + (vh-vs)**2 + b*ustar**2, tiny))
	end function calc_modified_richardson_number_at_height
	pure elemental function calc_bulk_richardson_number_at_height(uh, us, vh, vs, h, zs, thvh, thvs, ustar, g) result(rih)

[nusbaume]

I might label this pbl_height to be specific:

Suggested change

	pure function linear_interpolate_height_wrt_richardson(z, rino) result(pblh)
	pure function linear_interpolate_pbl_height_wrt_richardson(z, rino) result(pblh)

[nusbaume]

I might update the comment here to make it clear that phih is supposed to be dimensionless:

Suggested change

	real(kind_phys) :: phih ! Vertical temperature gradient [1]
	real(kind_phys) :: phih ! dimensionless vertical temperature gradient [1]

[nusbaume]

I might update the comment here to make it clear that phiminv is supposed to be dimensionless:

Suggested change

	real(kind_phys) :: phiminv ! Wind gradient [1]
	real(kind_phys) :: phiminv ! dimensionless wind gradient [1]

This same request holds for other phi declarations below as well.