Skip to content

Commit

Permalink
Adding clubb files for use with this tag.
Browse files Browse the repository at this point in the history
  • Loading branch information
@bstephens82
bstephens82 committed Nov 6, 2023
1 parent 10d4adb commit e5ee1a2
Show file tree
Hide file tree
Showing 92 changed files with 112,602 additions and 0 deletions.
279 changes: 279 additions & 0 deletions code_for_ben_clubb_silhs_files/clubb/LY93_pdf.F90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,279 @@
! $Id$
!===============================================================================
module LY93_pdf

! Description:
! The multivariate, two-component PDF of Lewellen and Yoh (1993).

! References:
! Lewellen, W. S. and Yoh, S., 1993. Binormal Model of Ensemble Partial
! Cloudiness. J. Atmos. Sci., 50, 9, 1228--1237.
!-------------------------------------------------------------------------

implicit none

public :: LY93_driver, & ! Procedure(s)
calc_mixt_frac_LY93, &
calc_params_LY93

private ! default scope

contains

!=============================================================================
subroutine LY93_driver( nz, wm, rtm, thlm, wp2, rtp2, & ! In
thlp2, Skw, Skrt, Skthl, & ! In
mu_w_1, mu_w_2, & ! Out
mu_rt_1, mu_rt_2, & ! Out
mu_thl_1, mu_thl_2, & ! Out
sigma_w_1_sqd, sigma_w_2_sqd, & ! Out
sigma_rt_1_sqd, sigma_rt_2_sqd, & ! Out
sigma_thl_1_sqd, sigma_thl_2_sqd, & ! Out
mixt_frac ) ! Out

! Description:
! Calculates the mixture fraction and the PDF component means and PDF
! component variances of w, rt, and theta-l following Lewellen and Yoh.

! References:
! Lewellen, W. S. and Yoh, S., 1993. Binormal Model of Ensemble Partial
! Cloudiness. J. Atmos. Sci., 50, 9, 1228--1237.
!-----------------------------------------------------------------------

use clubb_precision, only: &
core_rknd ! Variable(s)

implicit none

integer, intent(in) :: &
nz

! Input Variables
real( kind = core_rknd), dimension(nz), intent(in) :: &
wm, & ! Mean of w (overall) [m/s]
wp2, & ! Variance of w (overall) [m^2/s^2]
Skw, & ! Skewness of w (overall) [-]
rtm, & ! Mean of rt (overall) [kg/kg]
rtp2, & ! Variance of rt (overall) [kg^2/kg^2]
Skrt, & ! Skewness of rt (overall) [-]
thlm, & ! Mean of thl (overall) [K]
thlp2, & ! Variance of thl (overall) [K^2]
Skthl ! Skewness of thl (overall) [-]

! Output Variables
real( kind = core_rknd), dimension(nz), intent(out) :: &
mu_w_1, & ! Mean of w (1st PDF component) [m/s]
mu_w_2, & ! Mean of w (2nd PDF component) [m/s]
mu_rt_1, & ! Mean of rt (1st PDF component) [kg/kg]
mu_rt_2, & ! Mean of rt (2nd PDF component) [kg/kg]
mu_thl_1, & ! Mean of thl (1st PDF component) [K]
mu_thl_2, & ! Mean of thl (2nd PDF component) [K]
sigma_w_1_sqd, & ! Variance of w (1st PDF component) [m^2/s^2]
sigma_w_2_sqd, & ! Variance of w (2nd PDF component) [m^2/s^2]
sigma_rt_1_sqd, & ! Variance of rt (1st PDF component) [m^2/s^2]
sigma_rt_2_sqd, & ! Variance of rt (2nd PDF component) [m^2/s^2]
sigma_thl_1_sqd, & ! Variance of thl (1st PDF component) [m^2/s^2]
sigma_thl_2_sqd, & ! Variance of thl (2nd PDF component) [m^2/s^2]
mixt_frac ! Mixture fraction [-]

! Local Variables
real( kind = core_rknd), dimension(nz) :: &
Sk_max ! Maximum of magnitudes of skewness [-]


! Find the maximum of the magnitudes of skewness.
Sk_max = max( abs( Skw ), abs( Skrt ), abs( Skthl ) )

! Calculate mixture fraction.
mixt_frac = calc_mixt_frac_LY93( nz, Sk_max )

! Calculate the PDF parameters for w.
call calc_params_LY93( nz, wm, wp2, Skw, mixt_frac, & ! In
mu_w_1, mu_w_2, & ! Out
sigma_w_1_sqd, sigma_w_2_sqd ) ! Out

! Calculate the PDF parameters for rt.
call calc_params_LY93( nz, rtm, rtp2, Skrt, mixt_frac, & ! In
mu_rt_1, mu_rt_2, & ! Out
sigma_rt_1_sqd, sigma_rt_2_sqd ) ! Out

! Calculate the PDF parameters for thl.
call calc_params_LY93( nz, thlm, thlp2, Skthl, mixt_frac, & ! In
mu_thl_1, mu_thl_2, & ! Out
sigma_thl_1_sqd, sigma_thl_2_sqd ) ! Out


return

end subroutine LY93_driver

!=============================================================================
function calc_mixt_frac_LY93( nz, Sk_max ) &
result( mixt_frac )

! Description:
! Calculates mixture fraction iteratively according to Lewellen and Yoh.

! References:
! Eq. (21) of Lewellen, W. S. and Yoh, S., 1993. Binormal Model of Ensemble
! Partial Cloudiness. J. Atmos. Sci., 50, 9, 1228--1237.
!-----------------------------------------------------------------------

use constants_clubb, only: &
one, & ! Constant(s)
three_fourths, &
one_half, &
zero

use clubb_precision, only: &
core_rknd ! Variable(s)

implicit none

integer, intent(in) :: &
nz

! Input Variable
real( kind = core_rknd), dimension(nz), intent(in) :: &
Sk_max ! Maximum of magnitudes of skewness [-]

! Return Variable
real( kind = core_rknd), dimension(nz) :: &
mixt_frac ! Mixture fraction [-]

! Local Variables
real( kind = core_rknd) :: &
mixt_frac_low, & ! Low value of mixture frac. in iterative solver [-]
mixt_frac_high, & ! High value of mixture frac.in iterative solver [-]
expr_equal_zero ! Expr. mixt_frac^6 - Sk_max * ( 1 - mixt_frac ) [-]

! Tolerance for mixture fraction in solver [-]
real( kind = core_rknd) :: &
LY_mixt_frac_tol = 1.0e-4_core_rknd

integer :: k ! Vertical level index


do k = 1, nz, 1

if ( Sk_max(k) > 0.84_core_rknd ) then

mixt_frac_low = one_half
mixt_frac_high = one

do ! solve iteratively for mixture fraction

mixt_frac(k) = one_half * ( mixt_frac_low + mixt_frac_high )

expr_equal_zero &
= mixt_frac(k)**6 - Sk_max(k)**2 * ( one - mixt_frac(k) )

if ( abs( expr_equal_zero ) < LY_mixt_frac_tol ) then
! Mixture fraction has been solved for within the specificed
! tolerance.
exit
else
if ( expr_equal_zero > zero ) then
mixt_frac_high = mixt_frac(k)
else ! expr_equal_zero < 0
mixt_frac_low = mixt_frac(k)
endif
endif

enddo ! solve iteratively for mixture fraction

else ! Sk_max <= 0.84

mixt_frac(k) = three_fourths

endif

enddo ! k = 1, nz, 1


return

end function calc_mixt_frac_LY93

!=============================================================================
subroutine calc_params_LY93( nz, xm, xp2, Skx, mixt_frac, & ! In
mu_x_1, mu_x_2, & ! Out
sigma_x_1_sqd, sigma_x_2_sqd ) ! Out

! Description:
! Calculates the PDF component means and PDF component variances for
! variable x according to Lewellen and Yoh.

! References:
! Eq. (14), Eq. (15), Eq. (16), Eq. (17), and Eq. (18) of
! Lewellen, W. S. and Yoh, S., 1993. Binormal Model of Ensemble Partial
! Cloudiness. J. Atmos. Sci., 50, 9, 1228--1237.
!-----------------------------------------------------------------------

use constants_clubb, only: &
three, & ! Constant(s)
one, &
one_third, &
zero

use clubb_precision, only: &
core_rknd ! Variable(s)

implicit none

integer, intent(in) :: &
nz

! Input Variables
real( kind = core_rknd), dimension(nz), intent(in) :: &
xm, & ! Mean of x (overall) [units vary]
xp2, & ! Variance of x (overall) [(units vary)^2]
Skx, & ! Skewness of x (overall) [-]
mixt_frac ! Mixture fraction [-]

! Output Variables
real( kind = core_rknd), dimension(nz), intent(out) :: &
mu_x_1, & ! Mean of x (1st PDF component) [units vary]
mu_x_2, & ! Mean of x (2nd PDF component) [units vary]
sigma_x_1_sqd, & ! Variance of x (1st PDF component) [(units vary)^2]
sigma_x_2_sqd ! Variance of x (2nd PDF component) [(units vary)^2]

! Local Variables
real( kind = core_rknd), dimension(nz) :: &
sgn_Skx, & ! Sign of Skx [-]
B_x ! Spread of the PDF component means function [units vary]


! Find the sign of Skx
where ( Skx >= zero )
sgn_Skx = one
elsewhere ! Skx < 0
sgn_Skx = -one
endwhere

! Calculate B_x, the LY function for the spread of the PDF component means.
B_x = sgn_Skx * sqrt( xp2 ) &
* ( abs( Skx ) / ( one - mixt_frac ) )**one_third

! Calculate the mean of x in the 1st PDF component.
mu_x_1 = xm - B_x * ( one - mixt_frac )

! Calculate the mean of x in the 2nd PDF component.
mu_x_2 = xm + B_x * mixt_frac

! Calculate the variance of x in the 1st PDF component.
sigma_x_1_sqd = xp2 - B_x**2 * ( one - mixt_frac ) &
* ( one + mixt_frac + mixt_frac**2 ) &
/ ( three * mixt_frac )

! Calculate the variance of x in the 2nd PDF component.
sigma_x_2_sqd = xp2 + B_x**2 * ( one - mixt_frac )**2 / three


return

end subroutine calc_params_LY93

!=============================================================================

end module LY93_pdf
Loading

0 comments on commit e5ee1a2

Please sign in to comment.