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 17:57:11 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 23701892776 :                  + ( 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 23701892776 :                              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   142239224 :            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