LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - pdf_utilities.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 20 | 122 | 16.4 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 18 | 5.6 %

              Line data    Source code
           1             : !-------------------------------------------------------------------------
           2             : ! $Id$
           3             : !===============================================================================
           4             : module pdf_utilities
           5             :
           6             :   implicit none
           7             :
           8             :   private ! Set default scope to private
           9             :
          10             :   public :: mean_L2N,                  &
          11             :             mean_L2N_dp,               &
          12             :             stdev_L2N,                 &
          13             :             stdev_L2N_dp,              &
          14             :             corr_NL2NN,                &
          15             :             corr_NL2NN_dp,             &
          16             :             corr_NN2NL,                &
          17             :             corr_LL2NN,                &
          18             :             corr_LL2NN_dp,             &
          19             :             corr_NN2LL,                &
          20             :             compute_mean_binormal,     &
          21             :             compute_variance_binormal, &
          22             :             calc_comp_corrs_binormal,  &
          23             :             calc_corr_chi_x,           &
          24             :             calc_corr_eta_x,           &
          25             :             calc_corr_rt_x,            &
          26             :             calc_corr_thl_x,           &
          27             :             calc_xp2
          28             :
          29             :   contains
          30             :
          31             :   !=============================================================================
          32           0 :   elemental function mean_L2N( mu_x, sigma2_on_mu2 )  &
          33             :   result( mu_x_n )
          34             :
          35             :     ! Description:
          36             :     ! For a lognormally-distributed variable x, this function finds the mean of
          37             :     ! ln x (mu_x_n) for the ith component of the PDF, given the mean of x (mu_x)
          38             :     ! and the variance of x (sigma_sqd_x) for the ith component of the PDF. The
          39             :     ! value ln x is distributed normally when x is distributed lognormally.
          40             :
          41             :     ! References:
          42             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
          43             :     !    Marcel Dekker, 401 pp.
          44             :     !  -- App. B.
          45             :     !-----------------------------------------------------------------------
          46             :
          47             :     use constants_clubb, only: &
          48             :         one  ! Constant(s)
          49             :
          50             :     use clubb_precision, only: &
          51             :         core_rknd ! Variable(s)
          52             :
          53             :     implicit none
          54             :
          55             :     ! Input Variables
          56             :     real( kind = core_rknd ), intent(in) ::  &
          57             :       mu_x,          & ! Mean of x (ith PDF component)                     [-]
          58             :       sigma2_on_mu2    ! Ratio:  sigma_x^2 / mu_x^2 (ith PDF component)    [-]
          59             :
          60             :     ! Return Variable
          61             :     real( kind = core_rknd ) ::  &
          62             :       mu_x_n  ! Mean of ln x (ith PDF component)           [-]
          63             :
          64             :
          65             :     ! Find the mean of ln x for the ith component of the PDF.
          66             :     ! The max( mu_x / sqrt( 1 + sigma_x^2 / mu_x^2 ), tiny( mu_x ) ) statement
          67             :     ! is used to prevent taking ln 0, which will produce a result of -infinity.
          68             :     ! This would happen when mu_x is 0.  However, this code should not be
          69             :     ! entered when mu_x has a value of 0.
          70           0 :     mu_x_n = log( max( mu_x / sqrt( one + sigma2_on_mu2 ), tiny( mu_x ) ) )
          71             :
          72             :
          73             :     return
          74             :
          75             :   end function mean_L2N
          76             :
          77             :   !=============================================================================
          78           0 :   elemental function mean_L2N_dp( mu_x, sigma2_on_mu2 )  &
          79             :   result( mu_x_n )
          80             :
          81             :     ! Description:
          82             :     ! For a lognormally-distributed variable x, this function finds the mean of
          83             :     ! ln x (mu_x_n) for the ith component of the PDF, given the mean of x (mu_x)
          84             :     ! and the variance of x (sigma_sqd_x) for the ith component of the PDF. The
          85             :     ! value ln x is distributed normally when x is distributed lognormally.
          86             :     ! This function uses double precision variables.
          87             :
          88             :     ! References:
          89             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
          90             :     !    Marcel Dekker, 401 pp.
          91             :     !  -- App. B.
          92             :     !-----------------------------------------------------------------------
          93             :
          94             :     use constants_clubb, only: &
          95             :         one_dp  ! Constant(s)
          96             :
          97             :     use clubb_precision, only: &
          98             :         dp ! double precision
          99             :
         100             :     implicit none
         101             :
         102             :     ! Input Variables
         103             :     real( kind = dp ), intent(in) ::  &
         104             :       mu_x,          & ! Mean of x (ith PDF component)                     [-]
         105             :       sigma2_on_mu2    ! Ratio:  sigma_x^2 / mu_x^2 (ith PDF component)    [-]
         106             :
         107             :     ! Return Variable
         108             :     real( kind = dp ) ::  &
         109             :       mu_x_n  ! Mean of ln x (ith PDF component)           [-]
         110             :
         111             :
         112             :     ! Find the mean of ln x for the ith component of the PDF.
         113           0 :     mu_x_n = log( mu_x / sqrt( one_dp + sigma2_on_mu2 ) )
         114             :
         115             :
         116             :     return
         117             :
         118             :   end function mean_L2N_dp
         119             :
         120             :   !=============================================================================
         121           0 :   elemental function stdev_L2N( sigma2_on_mu2 )  &
         122             :   result( sigma_x_n )
         123             :
         124             :     ! Description:
         125             :     ! For a lognormally-distributed variable x, this function finds the standard
         126             :     ! deviation of ln x (sigma_x_n) for the ith component of the PDF, given the
         127             :     ! mean of x (mu_x) and the variance of x (sigma_sqd_x) for the ith component
         128             :     ! of the PDF.  The value ln x is distributed normally when x is distributed
         129             :     ! lognormally.
         130             :
         131             :     ! References:
         132             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
         133             :     !    Marcel Dekker, 401 pp.
         134             :     !  -- App. B.
         135             :     !-----------------------------------------------------------------------
         136             :
         137             :     use constants_clubb, only: &
         138             :         one  ! Constant(s)
         139             :
         140             :     use clubb_precision, only: &
         141             :         core_rknd ! Variable(s)
         142             :
         143             :     implicit none
         144             :
         145             :     ! Input Variables
         146             :     real( kind = core_rknd ), intent(in) ::  &
         147             :       sigma2_on_mu2    ! Ratio:  sigma_x^2 / mu_x^2 (ith PDF component)    [-]
         148             :
         149             :     ! Return Variable
         150             :     real( kind = core_rknd ) ::  &
         151             :       sigma_x_n  ! Standard deviation of ln x (ith PDF component)   [-]
         152             :
         153             :
         154             :     ! Find the standard deviation of ln x for the ith component of the PDF.
         155           0 :     sigma_x_n = sqrt( log( one + sigma2_on_mu2 ) )
         156             :
         157             :
         158             :     return
         159             :
         160             :   end function stdev_L2N
         161             :
         162             :   !=============================================================================
         163           0 :   elemental function stdev_L2N_dp( sigma2_on_mu2 )  &
         164             :   result( sigma_x_n )
         165             :
         166             :     ! Description:
         167             :     ! For a lognormally-distributed variable x, this function finds the standard
         168             :     ! deviation of ln x (sigma_x_n) for the ith component of the PDF, given the
         169             :     ! mean of x (mu_x) and the variance of x (sigma_sqd_x) for the ith component
         170             :     ! of the PDF.  The value ln x is distributed normally when x is distributed
         171             :     ! lognormally.
         172             :     ! This function uses double precision variables.
         173             :
         174             :     ! References:
         175             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
         176             :     !    Marcel Dekker, 401 pp.
         177             :     !  -- App. B.
         178             :     !-----------------------------------------------------------------------
         179             :
         180             :     use constants_clubb, only: &
         181             :         one_dp  ! Constant(s)
         182             :
         183             :     use clubb_precision, only: &
         184             :         dp ! double precision
         185             :
         186             :     implicit none
         187             :
         188             :     ! Input Variables
         189             :     real( kind = dp ), intent(in) ::  &
         190             :       sigma2_on_mu2    ! Ratio:  sigma_x^2 / mu_x^2 (ith PDF component)    [-]
         191             :
         192             :     ! Return Variable
         193             :     real( kind = dp ) ::  &
         194             :       sigma_x_n  ! Standard deviation of ln x (ith PDF component)   [-]
         195             :
         196             :
         197             :     ! Find the standard deviation of ln x for the ith component of the PDF.
         198           0 :     sigma_x_n = sqrt( log( one_dp + sigma2_on_mu2 ) )
         199             :
         200             :
         201             :     return
         202             :
         203             :   end function stdev_L2N_dp
         204             :
         205             :   !=============================================================================
         206           0 :   elemental function corr_NL2NN( corr_x_y, sigma_y_n, y_sigma2_on_mu2 )  &
         207             :   result( corr_x_y_n )
         208             :
         209             :     ! Description:
         210             :     ! For a normally-distributed variable x and a lognormally-distributed
         211             :     ! variable y, this function finds the correlation of x and ln y (corr_x_y_n)
         212             :     ! for the ith component of the PDF, given the correlation of x and y
         213             :     ! (corr_x_y) and the standard deviation of ln y (sigma_y_n) for the ith
         214             :     ! component of the PDF.  The value ln y is distributed normally when y is
         215             :     ! distributed lognormally.
         216             :
         217             :     ! References:
         218             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
         219             :     !    Marcel Dekker, 401 pp.
         220             :     !  -- Eq. B-1.
         221             :     !-----------------------------------------------------------------------
         222             :
         223             :     use constants_clubb, only: &
         224             :         max_mag_correlation, & ! Constant(s)
         225             :         zero
         226             :
         227             :     use clubb_precision, only: &
         228             :         core_rknd ! Variable(s)
         229             :
         230             :     implicit none
         231             :
         232             :     ! Input Variables
         233             :     real( kind = core_rknd ), intent(in) :: &
         234             :       corr_x_y,        & ! Correlation of x and y (ith PDF component)       [-]
         235             :       sigma_y_n,       & ! Standard deviation of ln y (ith PDF component)   [-]
         236             :       y_sigma2_on_mu2    ! Ratio:  sigma_y^2 / mu_y^2 (ith PDF component)   [-]
         237             :
         238             :     ! Return Variable
         239             :     real( kind = core_rknd ) ::  &
         240             :       corr_x_y_n  ! Correlation of x and ln y (ith PDF component) [-]
         241             :
         242             :
         243             :     ! Find the correlation of x and ln y for the ith component of the PDF.
         244             :     ! When sigma_y = 0 and mu_y > 0, y_sigma2_on_mu2 = 0.  This results in
         245             :     ! sigma_y_n = 0.  The resulting corr_x_y_n is undefined.  However, the
         246             :     ! divide-by-zero problem needs to be addressed in the code.
         247           0 :     if ( sigma_y_n > zero ) then
         248           0 :        corr_x_y_n = corr_x_y * sqrt( y_sigma2_on_mu2 ) / sigma_y_n
         249             :     else ! sigma_y_n = 0
         250             :        ! The value of sqrt( y_sigma2_on_mu2 ) / sigma_y_n can be rewritten as:
         251             :        ! sqrt( y_sigma2_on_mu2 ) / sqrt( ln( 1 + y_sigma2_on_mu2 ) ).
         252             :        ! This can be further rewritten as:
         253             :        ! sqrt( y_sigma2_on_mu2 / ln( 1 + y_sigma2_on_mu2 ) ),
         254             :        ! which has a limit of 1 as y_sigma2_on_mu2 approaches 0 from the right.
         255             :        ! When sigma_y_n = 0, the value of corr_x_y_n is undefined, so set it
         256             :        ! to corr_x_y.
         257           0 :        corr_x_y_n = corr_x_y
         258             :     endif ! sigma_y_n > 0
         259             :
         260             :     ! Clip the magnitude of the correlation of x and ln y in the ith PDF
         261             :     ! component, just in case the correlation (ith PDF component) of x and y and
         262             :     ! the standard deviation (ith PDF component) of ln y are inconsistent,
         263             :     ! resulting in an unrealizable value for corr_x_y_n.
         264           0 :     if ( corr_x_y_n > max_mag_correlation ) then
         265             :        corr_x_y_n = max_mag_correlation
         266           0 :     elseif ( corr_x_y_n < -max_mag_correlation ) then
         267           0 :        corr_x_y_n = -max_mag_correlation
         268             :     endif
         269             :
         270             :
         271             :     return
         272             :
         273             :   end function corr_NL2NN
         274             :
         275             :   !=============================================================================
         276           0 :   elemental function corr_NL2NN_dp( corr_x_y, sigma_y_n, y_sigma2_on_mu2 )  &
         277             :   result( corr_x_y_n )
         278             :
         279             :     ! Description:
         280             :     ! For a normally-distributed variable x and a lognormally-distributed
         281             :     ! variable y, this function finds the correlation of x and ln y (corr_x_y_n)
         282             :     ! for the ith component of the PDF, given the correlation of x and y
         283             :     ! (corr_x_y) and the standard deviation of ln y (sigma_y_n) for the ith
         284             :     ! component of the PDF.  The value ln y is distributed normally when y is
         285             :     ! distributed lognormally.
         286             :     ! This function uses double precision variables.
         287             :
         288             :     ! References:
         289             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
         290             :     !    Marcel Dekker, 401 pp.
         291             :     !  -- Eq. B-1.
         292             :     !-----------------------------------------------------------------------
         293             :
         294             :     use constants_clubb, only: &
         295             :         max_mag_correlation, & ! Constant(s)
         296             :         zero_dp
         297             :
         298             :     use clubb_precision, only: &
         299             :         dp ! double precision
         300             :
         301             :     implicit none
         302             :
         303             :     ! Input Variables
         304             :     real( kind = dp ), intent(in) :: &
         305             :       corr_x_y,        & ! Correlation of x and y (ith PDF component)       [-]
         306             :       sigma_y_n,       & ! Standard deviation of ln y (ith PDF component)   [-]
         307             :       y_sigma2_on_mu2    ! Ratio:  sigma_y^2 / mu_y^2 (ith PDF component)   [-]
         308             :
         309             :     ! Return Variable
         310             :     real( kind = dp ) ::  &
         311             :       corr_x_y_n  ! Correlation of x and ln y (ith PDF component) [-]
         312             :
         313             :
         314             :     ! Find the correlation of x and ln y for the ith component of the PDF.
         315             :     ! When sigma_y = 0 and mu_y > 0, y_sigma2_on_mu2 = 0.  This results in
         316             :     ! sigma_y_n = 0.  The resulting corr_x_y_n is undefined.  However, the
         317             :     ! divide-by-zero problem needs to be addressed in the code.
         318           0 :     if ( sigma_y_n > zero_dp ) then
         319           0 :        corr_x_y_n = corr_x_y * sqrt( y_sigma2_on_mu2 ) / sigma_y_n
         320             :     else ! sigma_y_n = 0
         321             :        ! The value of sqrt( y_sigma2_on_mu2 ) / sigma_y_n can be rewritten as:
         322             :        ! sqrt( y_sigma2_on_mu2 ) / sqrt( ln( 1 + y_sigma2_on_mu2 ) ).
         323             :        ! This can be further rewritten as:
         324             :        ! sqrt( y_sigma2_on_mu2 / ln( 1 + y_sigma2_on_mu2 ) ),
         325             :        ! which has a limit of 1 as y_sigma2_on_mu2 approaches 0 from the right.
         326             :        ! When sigma_y_n = 0, the value of corr_x_y_n is undefined, so set it
         327             :        ! to corr_x_y.
         328           0 :        corr_x_y_n = corr_x_y
         329             :     endif ! sigma_y_n > 0
         330             :
         331             :     ! Clip the magnitude of the correlation of x and ln y in the ith PDF
         332             :     ! component, just in case the correlation (ith PDF component) of x and y and
         333             :     ! the standard deviation (ith PDF component) of ln y are inconsistent,
         334             :     ! resulting in an unrealizable value for corr_x_y_n.
         335           0 :     if ( corr_x_y_n > real( max_mag_correlation, kind = dp ) ) then
         336             :        corr_x_y_n = real( max_mag_correlation, kind = dp )
         337           0 :     elseif ( corr_x_y_n < -real( max_mag_correlation, kind = dp ) ) then
         338           0 :        corr_x_y_n = -real( max_mag_correlation, kind = dp )
         339             :     endif
         340             :
         341             :
         342             :     return
         343             :
         344             :   end function corr_NL2NN_dp
         345             :
         346             :   !=============================================================================
         347           0 :   subroutine corr_NN2NL( nz, ngrdcol, &
         348           0 :                          corr_x_y_n, sigma_y_n, y_sigma2_on_mu2, &
         349           0 :                          corr_x_y )
         350             :
         351             :     ! Description:
         352             :     ! For a normally-distributed variable x and a lognormally-distributed
         353             :     ! variable y, this function finds the correlation of x and y (corr_x_y) for
         354             :     ! the ith component of the PDF, given the correlation of x and ln y
         355             :     ! (corr_x_y_n) and the standard deviation of ln y (sigma_y_n) for the ith
         356             :     ! component of the PDF.  The value ln y is distributed normally when y is
         357             :     ! distributed lognormally.
         358             :
         359             :     ! References:
         360             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
         361             :     !    Marcel Dekker, 401 pp.
         362             :     !  -- Eq. B-1.
         363             :     !-----------------------------------------------------------------------
         364             :
         365             :     use constants_clubb, only: &
         366             :         max_mag_correlation, & ! Constant(s)
         367             :         zero
         368             :
         369             :     use clubb_precision, only: &
         370             :         core_rknd ! Variable(s)
         371             :
         372             :     implicit none
         373             :
         374             :     ! Input Variables
         375             :     integer, intent(in) :: &
         376             :       nz,      & ! Number of vertical levels
         377             :       ngrdcol    ! Number of grid columns
         378             :
         379             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         380             :       corr_x_y_n,      & ! Correlation of x and ln y (ith PDF component)    [-]
         381             :       sigma_y_n,       & ! Standard deviation of ln y (ith PDF component)   [-]
         382             :       y_sigma2_on_mu2    ! Ratio:  sigma_y^2 / mu_y^2 (ith PDF component)   [-]
         383             :
         384             :     ! Output Variable
         385             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
         386             :       corr_x_y  ! Correlation of x and y (ith PDF component) [-]
         387             :
         388             :     ! Local variables
         389             :     integer :: &
         390             :       j, k  ! Loop iterators
         391             :
         392             :
         393             :     ! Find the correlation of x and y for the ith component of the PDF.
         394             :     ! When sigma_y = 0 and mu_y > 0, y_sigma2_on_mu2 = 0.  This results in
         395             :     ! sigma_y_n = 0.  The resulting corr_x_y and corr_x_y_n are undefined.
         396             :     ! However, the divide-by-zero problem needs to be addressed in the code.
         397           0 :     do k = 1, nz
         398           0 :       do j = 1, ngrdcol
         399             :
         400           0 :         if ( sigma_y_n(j,k) > zero ) then
         401             :            ! Use the maximum of y_sigma2_on_mu2 and tiny( y_sigma2_on_mu2 ) instead
         402             :            ! of just y_sigma2_on_mu2.  The value of y_sigma2_on_mu2 must already be
         403             :            ! greater than 0 in order for this block of code to be entered (when
         404             :            ! y_sigma2_on_mu2 = 0, sigma_y_n = 0, and this block of code is not
         405             :            ! entered).
         406             :            corr_x_y(j,k) = corr_x_y_n(j,k) * sigma_y_n(j,k) &
         407           0 :                             / sqrt( max( y_sigma2_on_mu2(j,k), tiny( y_sigma2_on_mu2(j,k) ) ) )
         408             :         else ! sigma_y_n = 0
         409             :            ! The value of sigma_y_n / sqrt( y_sigma2_on_mu2 ) can be rewritten as:
         410             :            ! sqrt( ln( 1 + y_sigma2_on_mu2 ) ) / sqrt( y_sigma2_on_mu2 ).
         411             :            ! This can be further rewritten as:
         412             :            ! sqrt( ln( 1 + y_sigma2_on_mu2 ) / y_sigma2_on_mu2 ),
         413             :            ! which has a limit of 1 as y_sigma2_on_mu2 approaches 0 from the right.
         414             :            ! When sigma_y_n = 0, the value of corr_x_y is undefined, so set it
         415             :            ! to corr_x_y_n.
         416           0 :            corr_x_y(j,k) = corr_x_y_n(j,k)
         417             :         endif ! sigma_y_n > 0
         418             :
         419             :         ! Clip the magnitude of the correlation of x and y in the ith PDF component,
         420             :         ! just in case the correlation (ith PDF component) of x and ln y and the
         421             :         ! standard deviation (ith PDF component) of ln y are inconsistent, resulting
         422             :         ! in an unrealizable value for corr_x_y.
         423           0 :         if ( corr_x_y(j,k) > max_mag_correlation ) then
         424           0 :            corr_x_y(j,k) = max_mag_correlation
         425           0 :         elseif ( corr_x_y(j,k) < -max_mag_correlation ) then
         426           0 :            corr_x_y(j,k) = -max_mag_correlation
         427             :         endif
         428             :
         429             :       end do
         430             :     end do
         431             :
         432           0 :     return
         433             :
         434             :   end subroutine corr_NN2NL
         435             :
         436             :   !=============================================================================
         437           0 :   elemental function corr_LL2NN( corr_x_y, sigma_x_n, sigma_y_n, &
         438             :                                  x_sigma2_on_mu2, y_sigma2_on_mu2 )  &
         439             :   result( corr_x_y_n )
         440             :
         441             :     ! Description:
         442             :     ! For lognormally-distributed variables x and y, this function finds the
         443             :     ! correlation of ln x and ln y (corr_x_y_n) for the ith component of the
         444             :     ! PDF, given the correlation of x and y (corr_x_y), the standard deviation
         445             :     ! of ln x (sigma_x_n), and the standard deviation of ln y (sigma_y_n) for
         446             :     ! the ith component of the PDF.  The value of ln x (or ln y) is distributed
         447             :     ! normally when x (or y) is distributed lognormally.
         448             :
         449             :     ! References:
         450             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
         451             :     !    Marcel Dekker, 401 pp.
         452             :     !  -- Eq. C-3.
         453             :     !-----------------------------------------------------------------------
         454             :
         455             :     use constants_clubb, only: &
         456             :         one,                 & ! Constant(s)
         457             :         zero,                &
         458             :         max_mag_correlation
         459             :
         460             :     use clubb_precision, only: &
         461             :         core_rknd ! Variable(s)
         462             :
         463             :     implicit none
         464             :
         465             :     ! Input Variables
         466             :     real( kind = core_rknd ), intent(in) ::  &
         467             :       corr_x_y,        & ! Correlation of x and y (ith PDF component)       [-]
         468             :       sigma_x_n,       & ! Standard deviation of ln x (ith PDF component)   [-]
         469             :       sigma_y_n,       & ! Standard deviation of ln y (ith PDF component)   [-]
         470             :       x_sigma2_on_mu2, & ! Ratio:  sigma_x^2 / mu_x^2 (ith PDF component)   [-]
         471             :       y_sigma2_on_mu2    ! Ratio:  sigma_y^2 / mu_y^2 (ith PDF component)   [-]
         472             :
         473             :     ! Return Variable
         474             :     real( kind = core_rknd ) ::  &
         475             :       corr_x_y_n  ! Correlation of ln x and ln y (ith PDF component)  [-]
         476             :
         477             :     ! Local Variable
         478             :     real( kind = core_rknd ) ::  &
         479             :       log_arg    ! Input into the ln function    [-]
         480             :
         481             :
         482             :     ! Find the correlation of ln x and ln y for the ith component of the PDF.
         483             :     ! When sigma_x = 0 and mu_x > 0, x_sigma2_on_mu2 = 0.  This results in
         484             :     ! sigma_x_n = 0.  The resulting corr_x_y_n is undefined.  The same holds
         485             :     ! true when sigma_y = 0 and mu_y > 0.  However, the divide-by-zero problem
         486             :     ! needs to be addressed in the code.
         487           0 :     if ( sigma_x_n > zero .and. sigma_y_n > zero ) then
         488             : !      corr_x_y_n = log( one + corr_x_y * sqrt( exp( sigma_x_n**2 ) - one ) &
         489             : !                                       * sqrt( exp( sigma_y_n**2 ) - one )  ) &
         490             : !                   / ( sigma_x_n * sigma_y_n )
         491           0 :        log_arg = one + corr_x_y * sqrt( x_sigma2_on_mu2 * y_sigma2_on_mu2 )
         492           0 :        corr_x_y_n = log( log_arg ) / ( sigma_x_n * sigma_y_n )
         493             :     else ! sigma_x_n = 0 or sigma_y_n = 0
         494             :        ! The value of corr_x_y_n is undefined, so set it to corr_x_y.
         495           0 :        corr_x_y_n = corr_x_y
         496             :     endif ! sigma_x_n > 0 and sigma_y_n > 0
         497             :
         498             :     ! Clip the magnitude of the correlation of ln x and ln y in the ith PDF
         499             :     ! component, just in case the correlation (ith PDF component) of x and y,
         500             :     ! the standard deviation (ith PDF component) of ln x, and the standard
         501             :     ! deviation (ith PDF component) of ln y are inconsistent, resulting in an
         502             :     ! unrealizable value for corr_x_y_n.
         503           0 :     if ( corr_x_y_n > max_mag_correlation ) then
         504             :        corr_x_y_n = max_mag_correlation
         505           0 :     elseif ( corr_x_y_n < -max_mag_correlation ) then
         506           0 :        corr_x_y_n = -max_mag_correlation
         507             :     endif
         508             :
         509             :
         510             :     return
         511             :
         512             :   end function corr_LL2NN
         513             :
         514             :   !=============================================================================
         515           0 :   elemental function corr_LL2NN_dp( corr_x_y, sigma_x_n, sigma_y_n, &
         516             :                                     x_sigma2_on_mu2, y_sigma2_on_mu2 )  &
         517             :   result( corr_x_y_n )
         518             :
         519             :     ! Description:
         520             :     ! For lognormally-distributed variables x and y, this function finds the
         521             :     ! correlation of ln x and ln y (corr_x_y_n) for the ith component of the
         522             :     ! PDF, given the correlation of x and y (corr_x_y), the standard deviation
         523             :     ! of ln x (sigma_x_n), and the standard deviation of ln y (sigma_y_n) for
         524             :     ! the ith component of the PDF.  The value of ln x (or ln y) is distributed
         525             :     ! normally when x (or y) is distributed lognormally.
         526             :     ! This function uses double precision variables.
         527             :
         528             :     ! References:
         529             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
         530             :     !    Marcel Dekker, 401 pp.
         531             :     !  -- Eq. C-3.
         532             :     !-----------------------------------------------------------------------
         533             :
         534             :     use constants_clubb, only: &
         535             :         one_dp,              & ! Constant(s)
         536             :         zero_dp,             &
         537             :         max_mag_correlation
         538             :
         539             :     use clubb_precision, only: &
         540             :         dp ! double precision
         541             :
         542             :     implicit none
         543             :
         544             :     ! Input Variables
         545             :     real( kind = dp ), intent(in) ::  &
         546             :       corr_x_y,        & ! Correlation of x and y (ith PDF component)       [-]
         547             :       sigma_x_n,       & ! Standard deviation of ln x (ith PDF component)   [-]
         548             :       sigma_y_n,       & ! Standard deviation of ln y (ith PDF component)   [-]
         549             :       x_sigma2_on_mu2, & ! Ratio:  sigma_x^2 / mu_x^2 (ith PDF component)   [-]
         550             :       y_sigma2_on_mu2    ! Ratio:  sigma_y^2 / mu_y^2 (ith PDF component)   [-]
         551             :
         552             :
         553             :     ! Return Variable
         554             :     real( kind = dp ) ::  &
         555             :       corr_x_y_n  ! Correlation of ln x and ln y (ith PDF component)  [-]
         556             :
         557             :
         558             :     ! Find the correlation of ln x and ln y for the ith component of the PDF.
         559             :     ! When sigma_x = 0 and mu_x > 0, x_sigma2_on_mu2 = 0.  This results in
         560             :     ! sigma_x_n = 0.  The resulting corr_x_y_n is undefined.  The same holds
         561             :     ! true when sigma_y = 0 and mu_y > 0.  However, the divide-by-zero problem
         562             :     ! needs to be addressed in the code.
         563           0 :     if ( sigma_x_n > zero_dp .and. sigma_y_n > zero_dp ) then
         564             :        corr_x_y_n &
         565             :        = log( one_dp + corr_x_y * sqrt( x_sigma2_on_mu2 * y_sigma2_on_mu2 ) ) &
         566           0 :          / ( sigma_x_n * sigma_y_n )
         567             :     else ! sigma_x_n = 0 or sigma_y_n = 0
         568             :        ! The value of corr_x_y_n is undefined, so set it to corr_x_y.
         569           0 :        corr_x_y_n = corr_x_y
         570             :     endif ! sigma_x_n > 0 and sigma_y_n > 0
         571             :
         572             :     ! Clip the magnitude of the correlation of ln x and ln y in the ith PDF
         573             :     ! component, just in case the correlation (ith PDF component) of x and y,
         574             :     ! the standard deviation (ith PDF component) of ln x, and the standard
         575             :     ! deviation (ith PDF component) of ln y are inconsistent, resulting in an
         576             :     ! unrealizable value for corr_x_y_n.
         577           0 :     if ( corr_x_y_n > real( max_mag_correlation, kind = dp ) ) then
         578             :        corr_x_y_n = real( max_mag_correlation, kind = dp )
         579           0 :     elseif ( corr_x_y_n < -real( max_mag_correlation, kind = dp ) ) then
         580           0 :        corr_x_y_n = -real( max_mag_correlation, kind = dp )
         581             :     endif
         582             :
         583             :
         584             :     return
         585             :
         586             :   end function corr_LL2NN_dp
         587             :
         588             :   !=============================================================================
         589           0 :   subroutine corr_NN2LL( nz, ngrdcol, &
         590           0 :                          corr_x_y_n, sigma_x_n, sigma_y_n, &
         591           0 :                          x_sigma2_on_mu2, y_sigma2_on_mu2, &
         592           0 :                          corr_x_y )
         593             :
         594             :     ! Description:
         595             :     ! For lognormally-distributed variables x and y, this function finds the
         596             :     ! correlation of x and y (corr_x_y) for the ith component of the PDF, given
         597             :     ! the correlation of ln x and ln y (corr_x_y_n), the standard deviation of
         598             :     ! ln x (sigma_x_n), and the standard deviation of ln y (sigma_y_n) for
         599             :     ! the ith component of the PDF.  The value of ln x (or ln y) is distributed
         600             :     ! normally when x (or y) is distributed lognormally.
         601             :
         602             :     ! References:
         603             :     !  Garvey, P. R., 2000: Probability methods for cost uncertainty analysis.
         604             :     !    Marcel Dekker, 401 pp.
         605             :     !  -- Eq. C-3.
         606             :     !-----------------------------------------------------------------------
         607             :
         608             :     use constants_clubb, only: &
         609             :         one,                 & ! Constant(s)
         610             :         zero,                &
         611             :         max_mag_correlation
         612             :
         613             :     use clubb_precision, only: &
         614             :         core_rknd ! Variable(s)
         615             :
         616             :     implicit none
         617             :
         618             :     ! Input Variables
         619             :     integer, intent(in) :: &
         620             :       nz,      & ! Number of vertical levels
         621             :       ngrdcol    ! Number of grid columns
         622             :
         623             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
         624             :       corr_x_y_n,      & ! Correlation of ln x and ln y (ith PDF component) [-]
         625             :       sigma_x_n,       & ! Standard deviation of ln x (ith PDF component)   [-]
         626             :       sigma_y_n,       & ! Standard deviation of ln y (ith PDF component)   [-]
         627             :       x_sigma2_on_mu2, & ! Ratio:  sigma_x^2 / mu_x^2 (ith PDF component)   [-]
         628             :       y_sigma2_on_mu2    ! Ratio:  sigma_y^2 / mu_y^2 (ith PDF component)   [-]
         629             :
         630             :     ! Output Variable
         631             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) ::  &
         632             :       corr_x_y  ! Correlation of x and y (ith PDF component)  [-]
         633             :
         634             :     ! Local variables
         635             :     integer :: &
         636             :       j, k  ! Loop iterators
         637             :
         638             :     ! Find the correlation of x and y for the ith component of the PDF.
         639             :     ! When sigma_x = 0 and mu_x > 0, x_sigma2_on_mu2 = 0.  This results in
         640             :     ! sigma_x_n = 0.  The resulting corr_x_y and corr_x_y_n are undefined.  The
         641             :     ! same holds true when sigma_y = 0 and mu_y > 0.  However, the
         642             :     ! divide-by-zero problem needs to be addressed in the code.
         643           0 :     do k = 1, nz
         644           0 :       do j = 1, ngrdcol
         645             :
         646           0 :         if ( sigma_x_n(j,k) > zero .and. sigma_y_n(j,k) > zero ) then
         647             :     !       corr_x_y = ( exp( sigma_x_n * sigma_y_n * corr_x_y_n ) - one ) &
         648             :     !                  / ( sqrt( exp( sigma_x_n**2 ) - one ) &
         649             :     !                      * sqrt( exp( sigma_y_n**2 ) - one ) )
         650             :            corr_x_y(j,k) = ( exp( sigma_x_n(j,k) * sigma_y_n(j,k) * corr_x_y_n(j,k) ) - one ) &
         651           0 :                            / sqrt( x_sigma2_on_mu2(j,k) * y_sigma2_on_mu2(j,k) )
         652             :         else ! sigma_x_n = 0 or sigma_y_n = 0
         653             :            ! The value of corr_x_y is undefined, so set it to corr_x_y_n.
         654           0 :            corr_x_y(j,k) = corr_x_y_n(j,k)
         655             :         endif ! sigma_x_n > 0 and sigma_y_n > 0
         656             :
         657             :         ! Clip the magnitude of the correlation of x and y in the ith PDF component,
         658             :         ! just in case the correlation (ith PDF component) of ln x and ln y, the
         659             :         ! standard deviation (ith PDF component) of ln x, and the standard deviation
         660             :         ! (ith PDF component) of ln y are inconsistent, resulting in an unrealizable
         661             :         ! value for corr_x_y.
         662           0 :         if ( corr_x_y(j,k) > max_mag_correlation ) then
         663           0 :            corr_x_y(j,k) = max_mag_correlation
         664           0 :         elseif ( corr_x_y(j,k) < -max_mag_correlation ) then
         665           0 :            corr_x_y(j,k) = -max_mag_correlation
         666             :         endif
         667             :
         668             :       end do
         669             :     end do
         670             :
         671           0 :     return
         672             :
         673             :   end subroutine corr_NN2LL
         674             :
         675             :   !=============================================================================
         676           0 :   elemental function compute_mean_binormal( mu_x_1, mu_x_2, mixt_frac ) &
         677             :   result( xm )
         678             :
         679             :     ! Description:
         680             :     ! Computes the overall grid-box mean of a binormal distribution from the
         681             :     ! mean of each component
         682             :
         683             :     ! References:
         684             :     !   None
         685             :     !-----------------------------------------------------------------------
         686             :
         687             :     use clubb_precision, only: &
         688             :         core_rknd ! Constant
         689             :
         690             :     use constants_clubb, only: &
         691             :         one ! Constant
         692             :
         693             :     implicit none
         694             :
         695             :     ! Input Variables
         696             :     real( kind = core_rknd ), intent(in) :: &
         697             :       mu_x_1,    & ! First PDF component mean of 'x'                       [?]
         698             :       mu_x_2,    & ! Second PDF component mean of 'x'                      [?]
         699             :       mixt_frac    ! Weight of the first PDF component                     [-]
         700             :
         701             :     ! Output Variables
         702             :     real( kind = core_rknd ) :: &
         703             :       xm           ! Mean of 'x' (overall)                                 [?]
         704             :
         705             :     !-----------------------------------------------------------------------
         706             :
         707             :     !----- Begin Code -----
         708           0 :     xm = mixt_frac * mu_x_1 + ( one - mixt_frac ) * mu_x_2
         709             :
         710             :
         711             :     return
         712             :
         713             :   end function compute_mean_binormal
         714             :
         715             :   !=============================================================================
         716           0 :   elemental function compute_variance_binormal( xm, mu_x_1, mu_x_2, &
         717             :                                                 stdev_x_1, stdev_x_2, &
         718             :                                                 mixt_frac ) &
         719             :   result( xp2 )
         720             :
         721             :     ! Description:
         722             :     ! Computes the overall grid-box variance of a binormal distribution from the
         723             :     ! variance of each component.
         724             :
         725             :     ! References:
         726             :     !   None
         727             :     !-----------------------------------------------------------------------
         728             :
         729             :     use clubb_precision, only: &
         730             :         core_rknd ! Constant
         731             :
         732             :     use constants_clubb, only: &
         733             :         one ! Constant
         734             :
         735             :     implicit none
         736             :
         737             :     ! Input Variables
         738             :     real( kind = core_rknd ), intent(in) :: &
         739             :       xm,        & ! Overall mean of 'x'                                   [?]
         740             :       mu_x_1,    & ! First PDF component mean of 'x'                       [?]
         741             :       mu_x_2,    & ! Second PDF component mean of 'x'                      [?]
         742             :       stdev_x_1, & ! Standard deviation of 'x' in the first PDF component  [?]
         743             :       stdev_x_2, & ! Standard deviation of 'x' in the second PDF component [?]
         744             :       mixt_frac    ! Weight of the first PDF component                     [-]
         745             :
         746             :     ! Output Variables
         747             :     real( kind = core_rknd ) :: &
         748             :       xp2          ! Variance of 'x' (overall)                             [?^2]
         749             :
         750             :     !-----------------------------------------------------------------------
         751             :
         752             :     !----- Begin Code -----
         753             :     xp2 = mixt_frac * ( ( mu_x_1 - xm )**2 + stdev_x_1**2 ) &
         754           0 :           + ( one - mixt_frac ) * ( ( mu_x_2 - xm )**2 + stdev_x_2**2 )
         755             :
         756             :
         757             :     return
         758             :
         759             :   end function compute_variance_binormal
         760             :
         761             :   !=============================================================================
         762    17870112 :   subroutine calc_comp_corrs_binormal( nz, ngrdcol,    & ! In
         763    17870112 :                                        xpyp, xm, ym,   & ! In
         764    17870112 :                                        mu_x_1, mu_x_2, & ! In
         765    17870112 :                                        mu_y_1, mu_y_2, & ! In
         766    17870112 :                                        sigma_x_1_sqd,  & ! In
         767    17870112 :                                        sigma_x_2_sqd,  & ! In
         768    17870112 :                                        sigma_y_1_sqd,  & ! In
         769    17870112 :                                        sigma_y_2_sqd,  & ! In
         770    17870112 :                                        mixt_frac,      & ! In
         771    17870112 :                                        corr_x_y_1,     & ! Out
         772    17870112 :                                        corr_x_y_2      ) ! Out
         773             :
         774             :     ! Description:
         775             :     ! Calculates the PDF component correlations of variables x and y, where
         776             :     ! x and y are both distributed as two-component normals (or binormals).
         777             :     ! The PDF component correlations are set equal to each other.
         778             :     !
         779             :     ! The overall covariance of x and y, <x'y'>, can be expressed in terms of
         780             :     ! PDF parameters by integrating over the PDF:
         781             :     !
         782             :     ! <x'y'> = INT(-inf:inf) INT(-inf:inf) ( x - <x> ) ( y - <y> ) P(x,y) dy dx;
         783             :     !
         784             :     ! where <x> is the overall mean of x, <y> is the overall mean of y, and
         785             :     ! P(x,y) is the equation for the two-component normal PDF of x and y.
         786             :     !
         787             :     ! The integral is evaluated, and the equation for <x'y'> is:
         788             :     !
         789             :     ! <x'y'> = mixt_frac * ( ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
         790             :     !                        + corr_x_y_1 * sigma_x_1 * sigma_y_1 )
         791             :     !          + ( 1 - mixt_frac ) * ( ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
         792             :     !                                  + corr_x_y_2 * sigma_x_2 * sigma_y_2 );
         793             :     !
         794             :     ! where mu_x_1 is the mean of x in the 1st PDF component, mu_x_2 is the mean
         795             :     ! of x in the 2nd PDF component, mu_y_1 is the mean of y in the 1st PDF
         796             :     ! component, mu_y_2 is the mean of y in the 2nd PDF component, sigma_x_1 is
         797             :     ! the standard deviation of x in the 1st PDF component, sigma_x_2 is the
         798             :     ! standard deviation of x in the 2nd PDF component, sigma_y_1 is the
         799             :     ! standard deviation of y in the 1st PDF component, sigma_y_2 is the
         800             :     ! standard deviation of y in the 2nd PDF component, corr_x_y_1 is the
         801             :     ! correlation of x and y in the 1st PDF component, corr_x_y_2 is the
         802             :     ! correlation of x and y in the 2nd PDF component, and mixt_frac is the
         803             :     ! mixture fraction (weight of the 1st PDF component).
         804             :     !
         805             :     ! This equation can be rewritten as:
         806             :     !
         807             :     ! <x'y'> = mixt_frac * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
         808             :     !          + mixt_frac * corr_x_y_1 * sigma_x_1 * sigma_y_1
         809             :     !          + ( 1 - mixt_frac ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> )
         810             :     !          + ( 1 - mixt_frac ) * corr_x_y_2 * sigma_x_2 * sigma_y_2.
         811             :     !
         812             :     ! Setting the two PDF component correlations equal to each other
         813             :     ! (corr_x_y_1 = corr_x_y_2), the equation can be solved for the PDF
         814             :     ! component correlations:
         815             :     !
         816             :     ! corr_x_y_1 = corr_x_y_2
         817             :     ! = ( <x'y'> - mixt_frac * ( mu_x_1 - <x> ) * ( mu_y_1 - <y> )
         818             :     !            - ( 1 - mixt_frac ) * ( mu_x_2 - <x> ) * ( mu_y_2 - <y> ) )
         819             :     !   / ( mixt_frac * sigma_x_1 * sigma_y_1
         820             :     !       + ( 1 - mixt_frac ) * sigma_x_2 * sigma_y_2 );
         821             :     !
         822             :     ! where -1 <= corr_x_y_1 = corr_x_y_2 <= 1.
         823             :     !
         824             :     ! When sigma_x_1 * sigma_y_1 = 0 and sigma_x_2 * sigma_y_2 = 0, at least one
         825             :     ! of x or y are constant within each PDF component, and both PDF component
         826             :     ! correlations are undefined.
         827             :
         828             :     ! References:
         829             :     !-----------------------------------------------------------------------
         830             :
         831             :     use constants_clubb, only: &
         832             :         max_mag_correlation, & ! Variable(s)
         833             :         one, &
         834             :         zero
         835             :
         836             :     use clubb_precision, only: &
         837             :         core_rknd    ! Variable(s)
         838             :
         839             :     implicit none
         840             :
         841             :     ! ---------------- Input Variables ----------------
         842             :     integer, intent(in) :: &
         843             :       nz,      & ! Number of vertical levels
         844             :       ngrdcol    ! Number of grid columns
         845             :
         846             :     real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
         847             :       xpyp,          & ! Covariance of x and y (overall)    [(x units)(y units)]
         848             :       xm,            & ! Mean of x (overall)                [x units]
         849             :       ym,            & ! Mean of y (overall)                [y units]
         850             :       mu_x_1,        & ! Mean of x (1st PDF component)      [x units]
         851             :       mu_x_2,        & ! Mean of x (2nd PDF component)      [x units]
         852             :       mu_y_1,        & ! Mean of y (1st PDF component)      [y units]
         853             :       mu_y_2,        & ! Mean of y (2nd PDF component)      [y units]
         854             :       sigma_x_1_sqd, & ! Variance of x (1st PDF component)  [(x units)^2]
         855             :       sigma_x_2_sqd, & ! Variance of x (2nd PDF component)  [(x units)^2]
         856             :       sigma_y_1_sqd, & ! Variance of y (1st PDF component)  [(y units)^2]
         857             :       sigma_y_2_sqd, & ! Variance of y (2nd PDF component)  [(y units)^2]
         858             :       mixt_frac        ! Mixture fraction                   [-]
         859             :
         860             :     ! ---------------- Output Variables ----------------
         861             :     real ( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
         862             :       corr_x_y_1, & ! Correlation of x and y (1st PDF component)    [-]
         863             :       corr_x_y_2    ! Correlation of x and y (2nd PDF component)    [-]
         864             :
         865             :     ! ---------------- Local Variables ----------------
         866             :     integer :: i, k
         867             :
         868             :     !$acc parallel loop gang vector collapse(2) default(present)
         869  1536829632 :     do k = 1, nz
         870 25380961632 :       do i = 1, ngrdcol
         871 47688264000 :         if ( sigma_x_1_sqd(i,k) * sigma_y_1_sqd(i,k) > zero &
         872 71532396000 :              .or. sigma_x_2_sqd(i,k) * sigma_y_2_sqd(i,k) > zero ) then
         873             :
         874             :            ! Calculate corr_x_y_1 (which also equals corr_x_y_2).
         875             :            corr_x_y_1(i,k) &
         876             :            = ( xpyp(i,k) &
         877             :                - mixt_frac(i,k) * ( mu_x_1(i,k) - xm(i,k) ) * ( mu_y_1(i,k) - ym(i,k) ) &
         878             :                - ( one - mixt_frac(i,k) ) * ( mu_x_2(i,k) - xm(i,k) ) * ( mu_y_2(i,k) - ym(i,k))) &
         879             :              / ( mixt_frac(i,k) * sqrt( sigma_x_1_sqd(i,k) * sigma_y_1_sqd(i,k) ) &
         880 23699327379 :                  + ( one - mixt_frac(i,k) ) * sqrt( sigma_x_2_sqd(i,k) * sigma_y_2_sqd(i,k) ) )
         881             :
         882             :            ! The correlation must fall within the bounds of
         883             :            ! -max_mag_correlation < corr_x_y_1 (= corr_x_y_2) < max_mag_correlation
         884             :            corr_x_y_1(i,k) = max( -max_mag_correlation, &
         885 23699327379 :                              min( max_mag_correlation, corr_x_y_1(i,k) ) )
         886             :
         887             :         else ! sigma_x_1^2 * sigma_y_1^2 = 0 and sigma_x_2^2 * sigma_y_2^2 = 0.
         888             :
         889             :            ! The correlation is undefined (output as 0).
         890   144804621 :            corr_x_y_1(i,k) = zero
         891             :
         892             :         endif
         893             :         ! Set corr_x_y_2 equal to corr_x_y_1.
         894 25363091520 :         corr_x_y_2(i,k) = corr_x_y_1(i,k)
         895             :       end do
         896             :     end do
         897             :     !$acc end parallel loop
         898             :
         899    17870112 :     return
         900             :
         901             :   end subroutine calc_comp_corrs_binormal
         902             :
         903             :   !=============================================================================
         904           0 :   elemental function calc_corr_chi_x( crt_i, cthl_i,             &
         905             :                                       sigma_rt_i, sigma_thl_i,   &
         906             :                                       sigma_chi_i,               &
         907             :                                       corr_rt_x_i, corr_thl_x_i )  &
         908             :   result( corr_chi_x_i )
         909             :
         910             :     ! Description:
         911             :     ! This function calculates the correlation of extended liquid water mixing
         912             :     ! ratio, chi (old s), and a generic variable x, within the ith component of
         913             :     ! the PDF.  The variable chi can be split into mean and turbulent
         914             :     ! components, such that:
         915             :     !
         916             :     ! chi = <chi> + chi';
         917             :     !
         918             :     ! where < > denotes a mean field an ' denotes a turbulent component.
         919             :     !
         920             :     ! The linearized equation for chi' is given in Larson et al. (2001), where
         921             :     ! within the ith component of the PDF:
         922             :     !
         923             :     ! chi_(i)' = Coef_rt(i) * r_t(i)' - Coef_thl(i) * th_l(i)'.
         924             :     !
         925             :     ! The equation for chi' can be multiplied by x'.  The equation becomes:
         926             :     !
         927             :     ! chi'x'_(i) = Coef_rt(i) * r_t'x'_(i) - Coef_thl(i) * th_l'x'_(i).
         928             :     !
         929             :     ! Averaging both sides, the covariance <chi'x'> is given by the equation:
         930             :     !
         931             :     ! <chi'x'_(i)> = Coef_rt(i) * <r_t'x'_(i)> - Coef_thl(i) * <th_l'x'_(i)>.
         932             :     !
         933             :     ! This equation can be rewritten as:
         934             :     !
         935             :     ! sigma_chi(i) * sigma_x(i) * corr_chi_x(i)
         936             :     !   = Coef_rt(i) * sigma_rt(i) * sigma_x(i) * corr_rt_x(i)
         937             :     !     - Coef_thl(i) * sigma_thl(i) * sigma_x(i) * corr_thl_x(i).
         938             :     !
         939             :     ! This equation can be solved for corr_chi_x(i):
         940             :     !
         941             :     ! corr_chi_x(i)
         942             :     ! = Coef_rt(i) * ( sigma_rt(i) / sigma_chi(i) ) * corr_rt_x(i)
         943             :     !   - Coef_thl(i) * ( sigma_thl(i) / sigma_chi(i) ) * corr_thl_x(i).
         944             :     !
         945             :     ! The correlation of chi and x within the ith component of the PDF is
         946             :     ! calculated.
         947             :
         948             :     ! References:
         949             :     ! Eq. (13) and Eq. (14) of Larson, V. E., R. Wood, P. R. Field, J.-C. Golaz,
         950             :     ! T. H. Vonder Haar, W. R. Cotton, 2001:  Systematic Biases in the
         951             :     ! Microphysics and Thermodynamics of Numerical Models That Ignore
         952             :     ! Subgrid-Scale Variability. J. Atmos. Sci., 58, 1117--1128,
         953             :     ! doi:https://doi.org/10.1175/1520-0469(2001)058%3C1117:SBITMA%3E2.0.CO;2.
         954             :     !
         955             :     ! Eq. (A29) of Griffin, B. M., 2016:  Improving the Subgrid-Scale
         956             :     ! Representation of Hydrometeors and Microphysical Feedback Effects Using a
         957             :     ! Multivariate PDF.  Doctoral dissertation, University of
         958             :     ! Wisconsin -- Milwaukee, Milwaukee, WI, Paper 1144, 165 pp., URL
         959             :     ! http://dc.uwm.edu/cgi/viewcontent.cgi?article=2149&context=etd.
         960             :     !-----------------------------------------------------------------------
         961             :
         962             :     use constants_clubb, only: &
         963             :         zero,                & ! Constant(s)
         964             :         chi_tol,             &
         965             :         max_mag_correlation
         966             :
         967             :     use clubb_precision, only: &
         968             :         core_rknd ! Variable(s)
         969             :
         970             :     implicit none
         971             :
         972             :     ! Input Variables
         973             :     real( kind = core_rknd ), intent(in) :: &
         974             :       crt_i,        & ! Coefficient of r_t for chi (old s) (ith PDF comp.)   [-]
         975             :       cthl_i,       & ! Coefficient of th_l for chi (ith PDF comp.)  [(kg/kg)/K]
         976             :       sigma_rt_i,   & ! Standard deviation of r_t (ith PDF component)    [kg/kg]
         977             :       sigma_thl_i,  & ! Standard deviation of th_l (ith PDF component)       [K]
         978             :       sigma_chi_i,  & ! Standard deviation of chi (ith PDF component)    [kg/kg]
         979             :       corr_rt_x_i,  & ! Correlation of r_t and x (ith PDF component)         [-]
         980             :       corr_thl_x_i    ! Correlation of th_l and x (ith PDF component)        [-]
         981             :
         982             :     ! Return Variable
         983             :     real( kind = core_rknd ) :: &
         984             :       corr_chi_x_i  ! Correlation of chi and x (ith PDF component)   [-]
         985             :
         986             :
         987             :     ! Calculate the correlation of chi and x in the ith PDF component.
         988           0 :     if ( sigma_chi_i > chi_tol ) then
         989             :
         990             :        corr_chi_x_i = crt_i * ( sigma_rt_i / sigma_chi_i ) * corr_rt_x_i  &
         991           0 :                       - cthl_i * ( sigma_thl_i / sigma_chi_i ) * corr_thl_x_i
         992             :
         993             :     else  ! sigma_chi_i = 0
         994             :
         995             :        ! The standard deviation of chi in the ith PDF component is 0.  This
         996             :        ! means that chi is constant within the ith PDF component, and the ith
         997             :        ! PDF component covariance of chi and x is also 0.  The correlation of
         998             :        ! chi and x is undefined in the ith PDF component, so a value of 0 will
         999             :        ! be used.
        1000             :        corr_chi_x_i = zero
        1001             :
        1002             :     endif
        1003             :
        1004             :     ! Clip the magnitude of the correlation of chi and x in the ith PDF
        1005             :     ! component, just in case the correlations and standard deviations used in
        1006             :     ! calculating it are inconsistent, resulting in an unrealizable value for
        1007             :     ! corr_chi_x_i.
        1008           0 :     if ( corr_chi_x_i > max_mag_correlation ) then
        1009             :        corr_chi_x_i = max_mag_correlation
        1010           0 :     elseif ( corr_chi_x_i < -max_mag_correlation ) then
        1011           0 :        corr_chi_x_i = -max_mag_correlation
        1012             :     endif
        1013             :
        1014             :
        1015             :     return
        1016             :
        1017             :   end function calc_corr_chi_x
        1018             :
        1019             :   !=============================================================================
        1020           0 :   elemental function calc_corr_eta_x( crt_i, cthl_i,             &
        1021             :                                       sigma_rt_i, sigma_thl_i,   &
        1022             :                                       sigma_eta_i, corr_rt_x_i,  &
        1023             :                                       corr_thl_x_i )  &
        1024             :   result( corr_eta_x_i )
        1025             :
        1026             :     ! Description:
        1027             :     ! This function calculates the correlation of the variable that is
        1028             :     ! orthogonal to extended liquid water mixing ratio in a PDF transformation,
        1029             :     ! eta (old t), and a generic variable x, within the ith component of
        1030             :     ! the PDF.
        1031             :
        1032             :     ! References:
        1033             :     ! Eq. (A30) of Griffin, B. M., 2016:  Improving the Subgrid-Scale
        1034             :     ! Representation of Hydrometeors and Microphysical Feedback Effects Using a
        1035             :     ! Multivariate PDF.  Doctoral dissertation, University of
        1036             :     ! Wisconsin -- Milwaukee, Milwaukee, WI, Paper 1144, 165 pp., URL
        1037             :     ! http://dc.uwm.edu/cgi/viewcontent.cgi?article=2149&context=etd.
        1038             :     !-----------------------------------------------------------------------
        1039             :
        1040             :     use constants_clubb, only: &
        1041             :         zero,                & ! Constant(s)
        1042             :         eta_tol,             &
        1043             :         max_mag_correlation
        1044             :
        1045             :     use clubb_precision, only: &
        1046             :         core_rknd ! Variable(s)
        1047             :
        1048             :     implicit none
        1049             :
        1050             :     ! Input Variables
        1051             :     real( kind = core_rknd ), intent(in) :: &
        1052             :       crt_i,        & ! Coefficient of r_t for chi (old s) (ith PDF comp.)   [-]
        1053             :       cthl_i,       & ! Coefficient of th_l for chi (ith PDF comp.)  [(kg/kg)/K]
        1054             :       sigma_rt_i,   & ! Standard deviation of r_t (ith PDF component)    [kg/kg]
        1055             :       sigma_thl_i,  & ! Standard deviation of th_l (ith PDF component)       [K]
        1056             :       sigma_eta_i,  & ! Standard deviation of eta (ith PDF component)    [kg/kg]
        1057             :       corr_rt_x_i,  & ! Correlation of r_t and x (ith PDF component)         [-]
        1058             :       corr_thl_x_i    ! Correlation of th_l and x (ith PDF component)        [-]
        1059             :
        1060             :     ! Return Variable
        1061             :     real( kind = core_rknd ) &
        1062             :       corr_eta_x_i  ! Correlation of eta and x (ith PDF component)   [-]
        1063             :
        1064             :
        1065             :     ! Calculate the correlation of eta and x in the ith PDF component.
        1066           0 :     if ( sigma_eta_i > eta_tol ) then
        1067             :
        1068             :        corr_eta_x_i = crt_i * ( sigma_rt_i / sigma_eta_i ) * corr_rt_x_i  &
        1069           0 :                       + cthl_i * ( sigma_thl_i / sigma_eta_i ) * corr_thl_x_i
        1070             :
        1071             :     else  ! sigma_eta_i = 0
        1072             :
        1073             :        ! The standard deviation of eta in the ith PDF component is 0.  This
        1074             :        ! means that eta is constant within the ith PDF component, and the ith
        1075             :        ! PDF component covariance of eta and x is also 0.  The correlation of
        1076             :        ! eta and x is undefined in the ith PDF component, so a value of 0 will
        1077             :        ! be used.
        1078             :        corr_eta_x_i = zero
        1079             :
        1080             :     endif
        1081             :
        1082             :     ! Clip the magnitude of the correlation of eta and x in the ith PDF
        1083             :     ! component, just in case the correlations and standard deviations used in
        1084             :     ! calculating it are inconsistent, resulting in an unrealizable value for
        1085             :     ! corr_eta_x_i.
        1086           0 :     if ( corr_eta_x_i > max_mag_correlation ) then
        1087             :        corr_eta_x_i = max_mag_correlation
        1088           0 :     elseif ( corr_eta_x_i < -max_mag_correlation ) then
        1089           0 :        corr_eta_x_i = -max_mag_correlation
        1090             :     endif
        1091             :
        1092             :
        1093             :     return
        1094             :
        1095             :   end function calc_corr_eta_x
        1096             :
        1097             :   !=============================================================================
        1098           0 :   elemental function calc_corr_rt_x( crt_i, sigma_rt_i, sigma_chi_i, &
        1099             :                                      sigma_eta_i, corr_chi_x_i, corr_eta_x_i ) &
        1100             :   result( corr_rt_x_i )
        1101             :
        1102             :     ! Description:
        1103             :     ! This function calculates the correlation of rt and x based on the
        1104             :     ! correlation of chi and x and the correlation of eta and x.
        1105             :
        1106             :     ! References:
        1107             :     !-----------------------------------------------------------------------
        1108             :
        1109             :     use constants_clubb, only: &
        1110             :         two,                 & ! Constant(s)
        1111             :         zero,                &
        1112             :         rt_tol,              &
        1113             :         max_mag_correlation
        1114             :
        1115             :     use clubb_precision, only: &
        1116             :         core_rknd ! Variable(s)
        1117             :
        1118             :     implicit none
        1119             :
        1120             :     ! Input Variables
        1121             :     real( kind = core_rknd ), intent(in) :: &
        1122             :       crt_i,        & ! Coef. of r_t in chi/eta eqns. (ith PDF component)    [-]
        1123             :       sigma_rt_i,   & ! Standard deviation of r_t (ith PDF component)    [kg/kg]
        1124             :       sigma_chi_i,  & ! Standard deviation of chi (ith PDF component)    [kg/kg]
        1125             :       sigma_eta_i,  & ! Standard deviation of eta (ith PDF component)    [kg/kg]
        1126             :       corr_chi_x_i, & ! Correlation of chi and x (ith PDF component)         [-]
        1127             :       corr_eta_x_i    ! Correlation of eta and x (ith PDF component)         [-]
        1128             :
        1129             :     ! Return Variable
        1130             :     real( kind = core_rknd ) :: &
        1131             :       corr_rt_x_i   ! Correlation of rt and x (ith PDF component)            [-]
        1132             :
        1133             :
        1134             :     ! Calculate the correlation of rt and x in the ith PDF component.
        1135           0 :     if ( sigma_rt_i > rt_tol ) then
        1136             :
        1137             :        corr_rt_x_i = ( sigma_eta_i * corr_eta_x_i &
        1138             :                        + sigma_chi_i * corr_chi_x_i ) &
        1139           0 :                      / ( two * crt_i * sigma_rt_i )
        1140             :
        1141             :     else  ! sigma_rt_i = 0
        1142             :
        1143             :        ! The standard deviation of rt in the ith PDF component is 0.  This means
        1144             :        ! that rt is constant within the ith PDF component, and the ith PDF
        1145             :        ! component covariance of rt and x is also 0.  The correlation of rt and
        1146             :        ! x is undefined in the ith PDF component, so a value of 0 will be used.
        1147             :        corr_rt_x_i = zero
        1148             :
        1149             :     endif
        1150             :
        1151             :     ! Clip the magnitude of the correlation of rt and x in the ith PDF
        1152             :     ! component, just in case the correlations and standard deviations used in
        1153             :     ! calculating it are inconsistent, resulting in an unrealizable value for
        1154             :     ! corr_rt_x_i.
        1155           0 :     if ( corr_rt_x_i > max_mag_correlation ) then
        1156             :        corr_rt_x_i = max_mag_correlation
        1157           0 :     elseif ( corr_rt_x_i < -max_mag_correlation ) then
        1158           0 :        corr_rt_x_i = -max_mag_correlation
        1159             :     endif
        1160             :
        1161             :
        1162             :     return
        1163             :
        1164             :   end function calc_corr_rt_x
        1165             :
        1166             :   !=============================================================================
        1167           0 :   elemental function calc_corr_thl_x( cthl_i, sigma_thl_i, sigma_chi_i, &
        1168             :                                       sigma_eta_i, corr_chi_x_i, &
        1169             :                                       corr_eta_x_i ) &
        1170             :   result( corr_thl_x_i )
        1171             :
        1172             :     ! Description:
        1173             :     ! This function calculates the correlation of thl and x based on the
        1174             :     ! correlation of chi and x and the correlation of eta and x.
        1175             :
        1176             :     ! References:
        1177             :     !-----------------------------------------------------------------------
        1178             :
        1179             :     use constants_clubb, only: &
        1180             :         two,                 & ! Constant(s)
        1181             :         zero,                &
        1182             :         thl_tol,             &
        1183             :         max_mag_correlation
        1184             :
        1185             :     use clubb_precision, only: &
        1186             :         core_rknd ! Variable(s)
        1187             :
        1188             :     implicit none
        1189             :
        1190             :     ! Input Variables
        1191             :     real( kind = core_rknd ), intent(in) :: &
        1192             :       cthl_i,       & ! Coef. of thl:  chi/eta eqns. (ith PDF comp.) [(kg/kg)/K]
        1193             :       sigma_thl_i,  & ! Standard deviation of thl (ith PDF component)        [K]
        1194             :       sigma_chi_i,  & ! Standard deviation of chi (ith PDF component)    [kg/kg]
        1195             :       sigma_eta_i,  & ! Standard deviation of eta (ith PDF component)    [kg/kg]
        1196             :       corr_chi_x_i, & ! Correlation of chi and x (ith PDF component)         [-]
        1197             :       corr_eta_x_i    ! Correlation of eta and x (ith PDF component)         [-]
        1198             :
        1199             :     ! Return Variable
        1200             :     real( kind = core_rknd ) :: &
        1201             :       corr_thl_x_i    ! Correlation of thl and x (ith PDF component)         [-]
        1202             :
        1203             :
        1204             :     ! Calculate the correlation of thl and x in the ith PDF component.
        1205           0 :     if ( sigma_thl_i > thl_tol ) then
        1206             :
        1207             :        corr_thl_x_i = ( sigma_eta_i * corr_eta_x_i &
        1208             :                         - sigma_chi_i * corr_chi_x_i ) &
        1209           0 :                       / ( two * cthl_i * sigma_thl_i )
        1210             :
        1211             :     else  ! sigma_thl_i = 0
        1212             :
        1213             :        ! The standard deviation of thl in the ith PDF component is 0.  This
        1214             :        ! means that thl is constant within the ith PDF component, and the ith
        1215             :        ! PDF component covariance of thl and x is also 0.  The correlation of
        1216             :        ! thl and x is undefined in the ith PDF component, so a value of 0 will
        1217             :        ! be used.
        1218             :        corr_thl_x_i = zero
        1219             :
        1220             :     endif
        1221             :
        1222             :     ! Clip the magnitude of the correlation of thl and x in the ith PDF
        1223             :     ! component, just in case the correlations and standard deviations used in
        1224             :     ! calculating it are inconsistent, resulting in an unrealizable value for
        1225             :     ! corr_thl_x_i.
        1226           0 :     if ( corr_thl_x_i > max_mag_correlation ) then
        1227             :        corr_thl_x_i = max_mag_correlation
        1228           0 :     elseif ( corr_thl_x_i < -max_mag_correlation ) then
        1229           0 :        corr_thl_x_i = -max_mag_correlation
        1230             :     endif
        1231             :
        1232             :
        1233             :     return
        1234             :
        1235             :   end function calc_corr_thl_x
        1236             :
        1237             :   !=============================================================================
        1238           0 :   elemental function calc_xp2( mu_x_1, mu_x_2, &
        1239             :                                mu_x_1_n, mu_x_2_n, &
        1240             :                                sigma_x_1, sigma_x_2, &
        1241             :                                sigma_x_1_n, sigma_x_2_n, &
        1242             :                                mixt_frac, x_frac_1, x_frac_2, &
        1243             :                                x_mean )  &
        1244             :   result( xp2 )
        1245             :
        1246             :     ! Description:
        1247             :     ! Calculates the overall variance of x, <x'^2>, where the distribution of x
        1248             :     ! is a combination of a lognormal distribution and/or 0 in each PDF
        1249             :     ! component.  The fraction of each component where x is lognormally
        1250             :     ! distributed (amd greater than 0) is x_frac_i (x_frac_1 and x_frac_2 for
        1251             :     ! PDF components 1 and 2, respectively).  The fraction of each component
        1252             :     ! where x has a value of 0 is ( 1 - x_frac_i ).  This function should be
        1253             :     ! called to calculate the total variance for x when <x'^2> is not provided
        1254             :     ! by a predictive (or other) equation.
        1255             :     !
        1256             :     ! This function is used to calculate the overall variance for rain water
        1257             :     ! mixing ratio, <r_r'^2>, and the overall variance for rain drop
        1258             :     ! concentration, <N_r'^2>.  The ratio of variance to mean-value-squared is
        1259             :     ! specified for the in-precip values of r_r and N_r within each PDF
        1260             :     ! component, allowing for the calculation of sigma_rr_i and sigma_Nr_i,
        1261             :     ! as well as sigma_rr_i_n and sigma_Nr_i_n.
        1262             :
        1263             :     ! References:
        1264             :     !-----------------------------------------------------------------------
        1265             :
        1266             :     use constants_clubb, only: &
        1267             :         two,  & ! Constant(s)
        1268             :         one,  &
        1269             :         zero, &
        1270             :         eps
        1271             :
        1272             :     use clubb_precision, only: &
        1273             :         core_rknd  ! Variable(s)
        1274             :
        1275             :     implicit none
        1276             :
        1277             :     ! Input Variables
        1278             :     real( kind = core_rknd ), intent(in) :: &
        1279             :       mu_x_1,      & ! Mean of x (1st PDF comp.) in x_frac                  [-]
        1280             :       mu_x_2,      & ! Mean of x (2nd PDF comp.) in x_frac                  [-]
        1281             :       mu_x_1_n,    & ! Mean of ln x (1st PDF comp.) in x_frac               [-]
        1282             :       mu_x_2_n,    & ! Mean of ln x (2nd PDF comp.) in x_frac               [-]
        1283             :       sigma_x_1,   & ! Standard deviation of x (1st PDF comp.) in x_frac    [-]
        1284             :       sigma_x_2,   & ! Standard deviation of x (2nd PDF comp.) in x_frac    [-]
        1285             :       sigma_x_1_n, & ! Standard deviation of ln x (1st PDF comp.) in x_frac [-]
        1286             :       sigma_x_2_n, & ! Standard deviation of ln x (2nd PDF comp.) in x_frac [-]
        1287             :       mixt_frac,   & ! Mixture fraction                                     [-]
        1288             :       x_frac_1,    & ! Fraction: x distributed lognormally (1st PDF comp.)  [-]
        1289             :       x_frac_2,    & ! Fraction: x distributed lognormally (2nd PDF comp.)  [-]
        1290             :       x_mean         ! Overall mean value of x                              [-]
        1291             :
        1292             :     ! Return Variable
        1293             :     real( kind = core_rknd ) :: &
        1294             :       xp2            ! Overall variance of x                                [-]
        1295             :
        1296             :
        1297             :     ! Calculate overall variance of x, <x'^2>.
        1298           0 :     if ( abs(sigma_x_1) < eps .and. abs(sigma_x_2) < eps ) then
        1299             :
        1300             :        ! The value of x is constant within both PDF components.
        1301             :        xp2 = ( mixt_frac * x_frac_1 * mu_x_1**2 &
        1302             :              + ( one - mixt_frac ) * x_frac_2 * mu_x_2**2 &
        1303             :              ) &
        1304           0 :              - x_mean**2
        1305             :
        1306             :
        1307           0 :     elseif ( abs(sigma_x_1) < eps ) then
        1308             :
        1309             :        ! The value of x is constant within the 1st PDF component.
        1310             :        xp2 = ( mixt_frac * x_frac_1 * mu_x_1**2 &
        1311             :              + ( one - mixt_frac ) * x_frac_2 &
        1312             :                * exp( two * mu_x_2_n + two * sigma_x_2_n**2 ) &
        1313             :              ) &
        1314           0 :              - x_mean**2
        1315             :
        1316             :
        1317           0 :     elseif ( abs(sigma_x_2) < eps ) then
        1318             :
        1319             :        ! The value of x is constant within the 2nd PDF component.
        1320             :        xp2 = ( mixt_frac * x_frac_1 &
        1321             :                * exp( two * mu_x_1_n + two * sigma_x_1_n**2 ) &
        1322             :              + ( one - mixt_frac ) * x_frac_2 * mu_x_2**2 &
        1323             :              ) &
        1324           0 :              - x_mean**2
        1325             :
        1326             :
        1327             :     else  ! sigma_x_1 and sigma_x_2 > 0
        1328             :
        1329             :        ! The value of x varies within both PDF component.
        1330             :        xp2 = ( mixt_frac * x_frac_1 &
        1331             :                * exp( two * mu_x_1_n + two * sigma_x_1_n**2 ) &
        1332             :              + ( one - mixt_frac ) * x_frac_2 &
        1333             :                * exp( two * mu_x_2_n + two * sigma_x_2_n**2 ) &
        1334             :              ) &
        1335           0 :              - x_mean**2
        1336             :
        1337             :
        1338             :     endif
        1339             :
        1340             :
        1341             :     ! As a check, prevent negative values for hydrometeor variances due to
        1342             :     ! numerical loss of precision error.
        1343           0 :     if ( xp2 < zero ) then
        1344           0 :        xp2 = zero
        1345             :     endif
        1346             :
        1347             :
        1348             :     return
        1349             :
        1350             :   end function calc_xp2
        1351             :
        1352             : !===============================================================================
        1353             :
        1354             : end module pdf_utilities

Generated by: LCOV version 1.14