LCOV - code coverage report
Current view: top level - physics/clubb/src/CLUBB_core - precipitation_fraction.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 0 226 0.0 %
Date: 2024-12-17 17:57:11 Functions: 0 3 0.0 %

          Line data    Source code
       1             : !-------------------------------------------------------------------------
       2             : ! $Id$
       3             : !===============================================================================
       4             : module precipitation_fraction
       5             : 
       6             :   ! Description:
       7             :   ! Sets overall precipitation fraction as well as the precipitation fraction
       8             :   ! in each PDF component.
       9             : 
      10             :   implicit none
      11             : 
      12             :   private
      13             : 
      14             :   public :: precip_fraction
      15             : 
      16             :   private :: component_precip_frac_weighted, &
      17             :              component_precip_frac_specify,  &
      18             :              precip_frac_assert_check
      19             : 
      20             :   integer, parameter, public :: &
      21             :     precip_frac_calc_type = 2  ! Option used to calculate component precip_frac
      22             : 
      23             :   contains
      24             : 
      25             :   !=============================================================================
      26           0 :   subroutine precip_fraction( nz, ngrdcol, hydromet_dim, &
      27           0 :                               hydromet, cloud_frac, cloud_frac_1, &
      28           0 :                               l_mix_rat_hm, l_frozen_hm, hydromet_tol, &
      29           0 :                               cloud_frac_2, ice_supersat_frac, &
      30           0 :                               ice_supersat_frac_1, ice_supersat_frac_2, &
      31           0 :                               mixt_frac, clubb_params, &
      32             :                               stats_metadata, &
      33           0 :                               stats_sfc, & 
      34           0 :                               precip_frac, &
      35           0 :                               precip_frac_1, &
      36           0 :                               precip_frac_2, &
      37           0 :                               precip_frac_tol )
      38             : 
      39             :     ! Description:
      40             :     ! Determines (overall) precipitation fraction over the horizontal domain, as
      41             :     ! well as the precipitation fraction within each PDF component, at every
      42             :     ! vertical grid level.
      43             : 
      44             :     ! References:
      45             :     !-----------------------------------------------------------------------
      46             : 
      47             :     use constants_clubb, only: &
      48             :         one,            & ! Constant(s)
      49             :         zero,           &
      50             :         cloud_frac_min, &
      51             :         fstderr
      52             : 
      53             :     use parameter_indices, only: &
      54             :         nparams, & ! Variable(s)
      55             :         iupsilon_precip_frac_rat
      56             : 
      57             :     use stats_type_utilities, only: &
      58             :         stat_update_var_pt  ! Procedure(s)
      59             : 
      60             :     use clubb_precision, only: &
      61             :         core_rknd  ! Variable(s)
      62             : 
      63             :     use error_code, only: &
      64             :         clubb_at_least_debug_level, &   ! Procedure
      65             :         err_code, &                     ! Error Indicator
      66             :         clubb_fatal_error               ! Constant
      67             : 
      68             :     use stats_type, only: &
      69             :         stats ! Type
      70             : 
      71             :     use stats_variables, only: &
      72             :         stats_metadata_type
      73             : 
      74             :     implicit none
      75             : 
      76             :     !------------------------- Input Variables -------------------------
      77             :     integer, intent(in) :: &
      78             :       nz,           & ! Number of model vertical grid levels
      79             :       ngrdcol,      & ! Number of grid columns
      80             :       hydromet_dim    ! Number of hydrometeor species
      81             : 
      82             :     real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
      83             :       hydromet    ! Mean of hydrometeor, hm (overall)           [units vary]
      84             : 
      85             :     logical, dimension(hydromet_dim), intent(in) :: &
      86             :       l_frozen_hm, & ! if true, then the hydrometeor is frozen; otherwise liquid
      87             :       l_mix_rat_hm   ! if true, then the quantity is a hydrometeor mixing ratio
      88             : 
      89             :     real( kind = core_rknd ), dimension(hydromet_dim), intent(in) :: &
      90             :       hydromet_tol    ! Tolerance values for all hydrometeors    [units vary]
      91             : 
      92             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
      93             :       cloud_frac,          & ! Cloud fraction (overall)                      [-]
      94             :       cloud_frac_1,        & ! Cloud fraction (1st PDF component)            [-]
      95             :       cloud_frac_2,        & ! Cloud fraction (2nd PDF component)            [-]
      96             :       ice_supersat_frac,   & ! Ice supersaturation fraction (overall)        [-]
      97             :       ice_supersat_frac_1, & ! Ice supersaturation fraction (1st PDF comp.)  [-]
      98             :       ice_supersat_frac_2, & ! Ice supersaturation fraction (2nd PDF comp.)  [-]
      99             :       mixt_frac              ! Mixture fraction                              [-]
     100             : 
     101             :     real( kind = core_rknd ), dimension(nparams), intent(in) :: &
     102             :       clubb_params    ! Array of CLUBB's tunable parameters    [units vary]
     103             :       
     104             :     type (stats_metadata_type), intent(in) :: &
     105             :       stats_metadata
     106             : 
     107             :     !------------------------- Inout Variables -------------------------
     108             :     type (stats), target, dimension(ngrdcol), intent(inout) :: &
     109             :       stats_sfc
     110             : 
     111             :     !------------------------- Output Variables -------------------------
     112             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
     113             :       precip_frac,   & ! Precipitation fraction (overall)               [-]
     114             :       precip_frac_1, & ! Precipitation fraction (1st PDF component)     [-]
     115             :       precip_frac_2    ! Precipitation fraction (2nd PDF component)     [-]
     116             : 
     117             :     real( kind = core_rknd ), dimension(ngrdcol), intent(out) :: &
     118             :       precip_frac_tol    ! Minimum precip. frac. when hydromet. are present  [-]
     119             : 
     120             :     !------------------------- Local Variables -------------------------
     121             : 
     122             :     ! "Maximum allowable" hydrometeor mixing ratio in-precip component mean.
     123             :     real( kind = core_rknd ), parameter :: &
     124             :       max_hm_ip_comp_mean = 0.0025_core_rknd  ! [kg/kg]
     125             : 
     126             :     real( kind = core_rknd ), parameter :: &
     127             :       precip_frac_tol_coef = 0.1_core_rknd  ! Coefficient for precip_frac_tol
     128             : 
     129             :     integer :: &
     130             :       j, k, ivar   ! Loop indices
     131             : 
     132             : 
     133             :     ! Initialize the precipitation fraction variables (precip_frac,
     134             :     ! precip_frac_1, and precip_frac_2) to 0.
     135           0 :     precip_frac(:,:)   = zero
     136           0 :     precip_frac_1(:,:) = zero
     137           0 :     precip_frac_2(:,:) = zero
     138             : 
     139             :     ! Set the minimum allowable precipitation fraction when hydrometeors are
     140             :     ! found at a grid level.
     141           0 :     if ( any( l_frozen_hm ) ) then
     142             :       ! Ice microphysics included.
     143           0 :       do j = 1, ngrdcol
     144           0 :         precip_frac_tol(j) &
     145             :         = max( precip_frac_tol_coef &
     146             :                * max( maxval( cloud_frac(j,:) ), maxval( ice_supersat_frac(j,:) ) ), &
     147           0 :                cloud_frac_min )
     148             :       end do
     149             :     else
     150             :       ! Warm microphysics.
     151           0 :       do j = 1, ngrdcol
     152           0 :         precip_frac_tol(j) = max( precip_frac_tol_coef * maxval( cloud_frac(j,:) ), &
     153           0 :                                 cloud_frac_min )
     154             :       end do
     155             :     endif
     156             :     
     157             :     !!! Find overall precipitation fraction.
     158             :     ! The precipitation fraction is the greatest cloud fraction at or above a
     159             :     ! vertical level.
     160           0 :     if ( any( l_frozen_hm ) ) then
     161             :       
     162             :       ! Ice microphysics included.
     163           0 :       precip_frac(:,nz) = max( cloud_frac(:,nz), ice_supersat_frac(:,nz) )
     164             :       
     165           0 :       do k = nz-1, 1, -1
     166           0 :         precip_frac(:,k) = max( precip_frac(:,k+1), cloud_frac(:,k), &
     167           0 :                                 ice_supersat_frac(:,k) )
     168             :       end do
     169             :       
     170             :     else
     171             :       
     172             :       ! Warm microphysics.
     173           0 :       precip_frac(:,nz) = cloud_frac(:,nz)
     174             :       
     175           0 :       do k = nz-1, 1, -1
     176           0 :         precip_frac(:,k) = max( precip_frac(:,k+1), cloud_frac(:,k) )
     177             :       end do
     178             :       
     179             :     end if
     180             : 
     181             :     !!! Special checks for overall precipitation fraction
     182           0 :     do k = 1, nz, 1
     183           0 :       do j = 1, ngrdcol
     184             : 
     185           0 :         if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) &
     186           0 :             .and. precip_frac(j,k) < precip_frac_tol(j) ) then
     187             : 
     188             :           ! In a scenario where we find any hydrometeor at this grid level, but
     189             :           ! no cloud at or above this grid level, set precipitation fraction to
     190             :           ! a minimum threshold value.
     191           0 :           precip_frac(j,k) = precip_frac_tol(j)
     192             : 
     193           0 :         elseif ( all( hydromet(j,k,:) < hydromet_tol(:) ) ) then
     194             : 
     195             :           ! The means (overall) of every precipitating hydrometeor are all less
     196             :           ! than their respective tolerance amounts.  They are all considered to
     197             :           ! have values of 0.  There are not any hydrometeor species found at
     198             :           ! this grid level.  There is also no cloud at or above this grid
     199             :           ! level, so set precipitation fraction to 0.
     200           0 :           precip_frac(j,k) = zero
     201             : 
     202             :         endif
     203             :         
     204             :       end do
     205             :     enddo ! Special checks for overall precipitation fraction loop: k = 1, nz, 1
     206             : 
     207             : 
     208             :     !!! Find precipitation fraction within each PDF component.
     209             :     !
     210             :     ! The overall precipitation fraction, f_p, is given by the equation:
     211             :     !
     212             :     ! f_p = a * f_p(1) + ( 1 - a ) * f_p(2);
     213             :     !
     214             :     ! where "a" is the mixture fraction (weight of PDF component 1), f_p(1) is
     215             :     ! the precipitation fraction within PDF component 1, and f_p(2) is the
     216             :     ! precipitation fraction within PDF component 2.  Overall precipitation
     217             :     ! fraction is found according the method above, and mixture fraction is
     218             :     ! already determined, leaving f_p(1) and f_p(2) to be solved for.  The
     219             :     ! values for f_p(1) and f_p(2) must satisfy the above equation.
     220           0 :     if ( precip_frac_calc_type == 1 ) then
     221             : 
     222             :       ! Calculatate precip_frac_1 and precip_frac_2 based on the greatest
     223             :       ! weighted cloud_frac_1 at or above a grid level.
     224             :       call component_precip_frac_weighted( nz, ngrdcol, hydromet_dim, & ! intent(in)
     225             :                                            l_frozen_hm, hydromet_tol, & ! intent(in)
     226             :                                            hydromet(:,:,:), precip_frac(:,:), & ! intent(in)
     227             :                                            cloud_frac_1(:,:), cloud_frac_2(:,:), & ! intent(in)
     228             :                                            ice_supersat_frac_1(:,:), & ! intent(in)
     229             :                                            ice_supersat_frac_2(:,:), mixt_frac(:,:), & !intent(in)
     230             :                                            precip_frac_tol(:), & ! intent(in)
     231             :                                            precip_frac_1(:,:), precip_frac_2(:,:) ) ! intent(out)
     232             :                                             
     233             :     elseif ( precip_frac_calc_type == 2 ) then
     234             : 
     235             :       ! Specified method.
     236             :       call component_precip_frac_specify( nz, ngrdcol, hydromet_dim, hydromet_tol, & ! intent(in)
     237             :                                           clubb_params(iupsilon_precip_frac_rat), & ! intent(in)
     238             :                                           hydromet(:,:,:), precip_frac(:,:), & ! intent(in)
     239             :                                           mixt_frac(:,:), precip_frac_tol(:), & ! intent(in)
     240             :                                           precip_frac_1(:,:), precip_frac_2(:,:) ) ! intent(out)
     241             : 
     242             :     else ! Invalid option selected.
     243             : 
     244             :        write(fstderr,*) "Invalid option to calculate precip_frac_1 " &
     245             :                         // "and precip_frac_2."
     246             :        err_code = clubb_fatal_error
     247             :        return
     248             : 
     249             :     endif ! precip_frac_calc_type
     250             : 
     251             : 
     252             :     ! Increase Precipiation Fraction under special conditions.
     253             :     !
     254             :     ! There are scenarios that sometimes occur that require precipitation
     255             :     ! fraction to be boosted.  Precipitation fraction is calculated from cloud
     256             :     ! fraction and ice supersaturation fraction.  For numerical reasons, CLUBB's
     257             :     ! PDF may become entirely subsaturated with respect to liquid and ice,
     258             :     ! resulting in both a cloud fraction of 0 and an ice supersaturation
     259             :     ! fraction of 0.  When this happens, precipitation fraction drops to 0 when
     260             :     ! there aren't any hydrometeors present at that grid level, or to
     261             :     ! precip_frac_tol when there is at least one hydrometeor present at that
     262             :     ! grid level.  However, sometimes there are large values of hydrometeors
     263             :     ! found at that grid level.  When this occurs, the PDF component in-precip
     264             :     ! mean of a hydrometeor can become ridiculously large.  This is because the
     265             :     ! ith PDF component in-precip mean of a hydrometeor, mu_hm_i,  is given by
     266             :     ! the equation:
     267             :     !
     268             :     ! mu_hm_i = hm_i / precip_frac_i;
     269             :     !
     270             :     ! where hm_i is the overall ith PDF component mean of the hydrometeor, and
     271             :     ! precip_frac_i is the ith PDF component precipitation fraction.  When
     272             :     ! precip_frac_i has a value of precip_frac_tol and hm_i is large, mu_hm_i
     273             :     ! can be huge.  This can cause enormous microphysical process rates and
     274             :     ! result in numerical instability.  It is also very inaccurate.
     275             :     !
     276             :     ! In order to limit this problem, the ith PDF component precipitation
     277             :     ! fraction is increased in order to decrease mu_hm_i.  First, an "upper
     278             :     ! limit" is set for mu_hm_i when the hydrometeor is a mixing ratio.  This is
     279             :     ! called max_hm_ip_comp_mean.  At every vertical level and for every
     280             :     ! hydrometeor mixing ratio, a check is made to try to prevent mu_hm_i from
     281             :     ! exceeding the "upper limit".  The check is:
     282             :     !
     283             :     ! hm_i / precip_frac_i ( which = mu_hm_i )  >  max_hm_ip_comp_mean,
     284             :     !
     285             :     ! which can be rewritten:
     286             :     !
     287             :     ! hm_i > precip_frac_i * max_hm_ip_comp_mean.
     288             :     !
     289             :     ! Since hm_i has not been calculated yet, the assumption for this check is
     290             :     ! that all of the hydrometeor is found in one PDF component, which is the
     291             :     ! worst-case scenario in violating this limit.  The check becomes:
     292             :     !
     293             :     ! <hm> / ( mixt_frac * precip_frac_1 )  >  max_hm_ip_comp_mean;
     294             :     !    in PDF comp. 1; and
     295             :     ! <hm> / ( ( 1 - mixt_frac ) * precip_frac_2 )  >  max_hm_ip_comp_mean;
     296             :     !    in PDF comp. 2.
     297             :     !
     298             :     ! These limits can be rewritten as:
     299             :     !
     300             :     ! <hm>  >  mixt_frac * precip_frac_1 * max_hm_ip_comp_mean;
     301             :     !    in PDF comp. 1; and
     302             :     ! <hm>  >  ( 1 - mixt_frac ) * precip_frac_2 * max_hm_ip_comp_mean;
     303             :     !    in PDF comp. 2.
     304             :     !
     305             :     ! When component precipitation fraction is found to be in excess of the
     306             :     ! limit, precip_frac_i is increased to:
     307             :     !
     308             :     ! <hm> / ( mixt_frac * max_hm_ip_comp_mean );
     309             :     !    when the limit is exceeded in PDF comp. 1; and
     310             :     ! <hm> / ( ( 1 - mixt_frac ) * max_hm_ip_comp_mean );
     311             :     !    when the limit is exceeded in PDF comp. 2.
     312             :     !
     313             :     ! Of course, precip_frac_i is not allowed to exceed 1, so when
     314             :     ! <hm> / mixt_frac (or <hm> / ( 1 - mixt_frac )) is already greater than
     315             :     ! max_hm_ip_comp_mean, mu_hm_i will also have to be greater than
     316             :     ! max_hm_ip_comp_mean.  However, the value of mu_hm_i is still reduced when
     317             :     ! compared to what it would have been using precip_frac_tol.  In the event
     318             :     ! that multiple hydrometeor mixing ratios violate the check, the code is set
     319             :     ! up so that precip_frac_i is increased based on the highest hm_i.
     320             :     
     321             : 
     322           0 :    do ivar = 1, hydromet_dim
     323           0 :       if ( l_mix_rat_hm(ivar) ) then
     324           0 :         do k = 1, nz
     325           0 :           do j = 1, ngrdcol
     326             : 
     327             :             ! The hydrometeor is a mixing ratio.
     328             : 
     329           0 :             if ( hydromet(j,k,ivar) >= hydromet_tol(ivar) .and. &
     330             :                  hydromet(j,k,ivar) > mixt_frac(j,k) * precip_frac_1(j,k) &
     331             :                                      * max_hm_ip_comp_mean ) then
     332             : 
     333             :                 ! Increase precipitation fraction in the 1st PDF component.
     334             :                 precip_frac_1(j,k) &
     335             :                 = min( hydromet(j,k,ivar) &
     336           0 :                        / ( mixt_frac(j,k) * max_hm_ip_comp_mean ), one )
     337             : 
     338             :                 ! The value of precip_frac_1 must be at least precip_frac_tol
     339             :                 ! when precipitation is found in the 1st PDF component.
     340           0 :                 precip_frac_1(j,k) = max( precip_frac_1(j,k), precip_frac_tol(j) )
     341             : 
     342             :             endif ! <hm>/(mixt_frac*precip_frac_1) > max_hm_ip_comp_mean
     343             : 
     344           0 :             if ( hydromet(j,k,ivar) >= hydromet_tol(ivar) .and. &
     345             :                   hydromet(j,k,ivar) > ( one - mixt_frac(j,k) ) * precip_frac_2(j,k) &
     346           0 :                                      * max_hm_ip_comp_mean ) then
     347             : 
     348             :                 ! Increase precipitation fraction in the 2nd PDF component.
     349             :                 precip_frac_2(j,k) &
     350             :                 = min( hydromet(j,k,ivar) &
     351           0 :                        / ( ( one - mixt_frac(j,k) ) * max_hm_ip_comp_mean ), one )
     352             : 
     353             :                 ! The value of precip_frac_2 must be at least precip_frac_tol
     354             :                 ! when precipitation is found in the 2nd PDF component.
     355           0 :                 precip_frac_2(j,k) = max( precip_frac_2(j,k), precip_frac_tol(j) )
     356             : 
     357             :             endif ! <hm>/((1-mixt_frac)*precip_frac_2) > max_hm_ip_comp_mean
     358             : 
     359             :           end do
     360             :         end do ! k = 1, nz
     361             :       end if ! l_mix_rat_hm(ivar)
     362             :     end do ! ivar = 1, hydromet_dim
     363             : 
     364             :     ! Recalculate overall precipitation fraction for consistency.
     365             :     precip_frac(:,:) = mixt_frac(:,:) * precip_frac_1(:,:) &
     366           0 :                        + ( one - mixt_frac(:,:) ) * precip_frac_2(:,:)
     367             : 
     368             :     ! Double check that precip_frac_tol <= precip_frac <= 1 when hydrometeors
     369             :     ! are found at a grid level.
     370             :     ! PLEASE DO NOT ALTER precip_frac, precip_frac_1, or precip_frac_2 anymore
     371             :     ! after this point in the code.
     372           0 :     do k = 1, nz, 1
     373           0 :       do j = 1, ngrdcol
     374           0 :         if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
     375           0 :           precip_frac(j,k) = min( max( precip_frac(j,k), precip_frac_tol(j) ), one )
     376             :         end if ! any( hydromet(k,:) >= hydromet_tol(:) )
     377             :       end do
     378             :     end do ! k = 1, nz, 1
     379             : 
     380             : 
     381             :     ! Statistics
     382           0 :     if ( stats_metadata%l_stats_samp ) then
     383           0 :       if ( stats_metadata%iprecip_frac_tol > 0 ) then
     384           0 :         do j = 1, ngrdcol
     385           0 :           call stat_update_var_pt( stats_metadata%iprecip_frac_tol, 1, precip_frac_tol(j), & ! intent(in)
     386           0 :                                    stats_sfc(j) ) ! intent(inout)
     387             :         end do
     388             :       end if ! stats_metadata%iprecip_frac_tol
     389             :     end if ! stats_metadata%l_stats_samp
     390             : 
     391             : 
     392             :     ! Assertion check for precip_frac, precip_frac_1, and precip_frac_2.
     393           0 :     if ( clubb_at_least_debug_level( 2 ) ) then
     394           0 :       do j = 1, ngrdcol
     395             :         call precip_frac_assert_check( nz, hydromet_dim, hydromet_tol,                    & ! intent(in)
     396           0 :                                        hydromet(j,:,:), mixt_frac(j,:), precip_frac(j,:), & ! in
     397           0 :                                        precip_frac_1(j,:), precip_frac_2(j,:),            & ! intent(in)
     398           0 :                                        precip_frac_tol(j) )                                 ! intent(in)
     399             :       end do
     400             :     endif
     401             : 
     402           0 :     return
     403             : 
     404             :   end subroutine precip_fraction
     405             : 
     406             :   !=============================================================================
     407             :   subroutine component_precip_frac_weighted( nz, ngrdcol, hydromet_dim, &
     408             :                                              l_frozen_hm, hydromet_tol, &
     409             :                                              hydromet, precip_frac, &
     410             :                                              cloud_frac_1, cloud_frac_2, &
     411             :                                              ice_supersat_frac_1, &
     412             :                                              ice_supersat_frac_2, mixt_frac, &
     413             :                                              precip_frac_tol, &
     414             :                                              precip_frac_1, precip_frac_2 )
     415             : 
     416             :     ! Description:
     417             :     ! Set precipitation fraction in each component of the PDF.  The weighted 1st
     418             :     ! PDF component precipitation fraction (weighted_pfrac_1) at a grid level is
     419             :     ! calculated by the greatest value of mixt_frac * cloud_frac_1 at or above
     420             :     ! the relevant grid level.  Likewise, the weighted 2nd PDF component
     421             :     ! precipitation fraction (weighted_pfrac_2) at a grid level is calculated by
     422             :     ! the greatest value of ( 1 - mixt_frac ) * cloud_frac_2 at or above the
     423             :     ! relevant grid level.
     424             :     !
     425             :     ! The fraction weighted_pfrac_1 / ( weighted_pfrac_1 + weighted_pfrac_2 ) is
     426             :     ! the weighted_pfrac_1 fraction.  Multiplying this fraction by overall
     427             :     ! precipitation fraction and then dividing by mixt_frac produces the 1st PDF
     428             :     ! component precipitation fraction (precip_frac_1).  Then, calculate the 2nd
     429             :     ! PDF component precipitation fraction (precip_frac_2) accordingly.
     430             : 
     431             :     ! References:
     432             :     !-----------------------------------------------------------------------
     433             : 
     434             :     use constants_clubb, only: &
     435             :         one,  & ! Constant(s)
     436             :         zero
     437             : 
     438             :     use clubb_precision, only: &
     439             :         core_rknd  ! Variable(s)
     440             : 
     441             :     implicit none
     442             : 
     443             :     ! Input Variables
     444             :     integer, intent(in) :: &
     445             :       nz,           & ! Number of model vertical grid levels
     446             :       ngrdcol,      & ! Number of grid columns
     447             :       hydromet_dim    ! Number of hydrometeor species
     448             : 
     449             :     logical, dimension(hydromet_dim), intent(in) :: &
     450             :       l_frozen_hm  ! if true, then the hydrometeor is frozen; otherwise liquid
     451             : 
     452             :     real( kind = core_rknd ), dimension(hydromet_dim), intent(in) :: &
     453             :       hydromet_tol    ! Tolerance values for all hydrometeors    [units vary]
     454             : 
     455             :     real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
     456             :       hydromet    ! Mean of hydrometeor, hm (overall)           [units vary]
     457             : 
     458             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
     459             :       precip_frac,         & ! Precipitation fraction (overall)              [-]
     460             :       cloud_frac_1,        & ! Cloud fraction (1st PDF component)            [-]
     461             :       cloud_frac_2,        & ! Cloud fraction (2nd PDF component)            [-]
     462             :       ice_supersat_frac_1, & ! Ice supersaturation fraction (1st PDF comp.)  [-]
     463             :       ice_supersat_frac_2, & ! Ice supersaturation fraction (2nd PDF comp.)  [-]
     464             :       mixt_frac              ! Mixture fraction                              [-]
     465             : 
     466             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
     467             :       precip_frac_tol    ! Minimum precip. frac. when hydromet. are present  [-]
     468             : 
     469             :     ! Output Variables
     470             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
     471             :       precip_frac_1, & ! Precipitation fraction (1st PDF component)     [-]
     472             :       precip_frac_2    ! Precipitation fraction (2nd PDF component)     [-]
     473             : 
     474             :     ! Local Variables
     475             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
     476             :       weighted_pfrac_1, & ! Product of mixt_frac and cloud_frac_1           [-]
     477             :       weighted_pfrac_2    ! Product of ( 1 - mixt_frac ) and cloud_frac_2   [-]
     478             : 
     479             :     integer :: k, j  ! Loop index
     480             : 
     481             : 
     482             :     !!! Find precipitation fraction within PDF component 1.
     483             :     ! The method used to find overall precipitation fraction will also be to
     484             :     ! find precipitation fraction within PDF component 1 and PDF component 2.
     485             :     ! In order to do so, it is assumed (poorly) that PDF component 1 overlaps
     486             :     ! PDF component 1 and that PDF component 2 overlaps PDF component 2 at every
     487             :     ! vertical level in the vertical profile.  
     488             :     
     489             :     ! The weighted precipitation fraction in PDF component 1 is the greatest
     490             :     ! value of the product of mixture fraction (mixt_frac) and 1st PDF
     491             :     ! component cloud fraction at or above a vertical level.  Likewise, the
     492             :     ! weighted precipitation fraction in PDF component 2 is the greatest
     493             :     ! value of the product of ( 1 - mixt_frac ) and 2nd PDF component cloud
     494             :     ! fraction at or above a vertical level.
     495             :     if ( any( l_frozen_hm ) ) then
     496             :       
     497             :       ! Ice microphysics included.
     498             : 
     499             :       ! Weighted precipitation fraction in PDF component 1.
     500             :       weighted_pfrac_1(:,nz) &
     501             :       = max( mixt_frac(:,nz) * cloud_frac_1(:,nz), &
     502             :              mixt_frac(:,nz) * ice_supersat_frac_1(:,nz) )
     503             : 
     504             :       ! Weighted precipitation fraction in PDF component 2.
     505             :       weighted_pfrac_2(:,nz) &
     506             :       = max( ( one - mixt_frac(:,nz) ) * cloud_frac_2(:,nz), &
     507             :              ( one - mixt_frac(:,nz) ) * ice_supersat_frac_2(:,nz) )
     508             :              
     509             :       do k = nz, 1, -1
     510             :         do j = 1, ngrdcol
     511             :           ! Weighted precipitation fraction in PDF component 1.
     512             :           weighted_pfrac_1(j,k) &
     513             :           = max( weighted_pfrac_1(j,k+1), &
     514             :                  mixt_frac(j,k) * cloud_frac_1(j,k), &
     515             :                  mixt_frac(j,k) * ice_supersat_frac_1(j,k) )
     516             : 
     517             :           ! Weighted precipitation fraction in PDF component 2.
     518             :           weighted_pfrac_2(j,k) &
     519             :           = max( weighted_pfrac_2(j,k+1), &
     520             :                  ( one - mixt_frac(j,k) ) * cloud_frac_2(j,k), &
     521             :                  ( one - mixt_frac(j,k) ) * ice_supersat_frac_2(j,k) )
     522             :         end do
     523             :       end do
     524             :       
     525             :       
     526             :     else
     527             :       
     528             :        ! Warm microphysics.
     529             : 
     530             :        ! Weighted precipitation fraction in PDF component 1.
     531             :        weighted_pfrac_1(:,nz) = mixt_frac(:,nz) * cloud_frac_1(:,nz)
     532             : 
     533             :        ! Weighted precipitation fraction in PDF component 2.
     534             :        weighted_pfrac_2(:,nz) = ( one - mixt_frac(:,nz) ) * cloud_frac_2(:,nz)
     535             :        
     536             :       do k = nz, 1, -1
     537             :         do j = 1, ngrdcol
     538             :           ! Weighted precipitation fraction in PDF component 1.
     539             :           weighted_pfrac_1(j,k) &
     540             :           = max( weighted_pfrac_1(j,k+1), &
     541             :                  mixt_frac(j,k) * cloud_frac_1(j,k) )
     542             : 
     543             :           ! Weighted precipitation fraction in PDF component 2.
     544             :           weighted_pfrac_2(j,k) &
     545             :           = max( weighted_pfrac_2(j,k+1), &
     546             :                  ( one - mixt_frac(j,k) ) * cloud_frac_2(j,k) )
     547             :         end do
     548             :       end do
     549             :       
     550             :     end if
     551             :     
     552             : 
     553             :     ! Calculate precip_frac_1 and special cases for precip_frac_1.
     554             :     do k = 1, nz, 1
     555             :       do j = 1, ngrdcol
     556             :         ! Calculate precipitation fraction in the 1st PDF component.
     557             :         if ( weighted_pfrac_1(j,k) + weighted_pfrac_2(j,k) > zero ) then
     558             : 
     559             :           ! Adjust weighted 1st PDF component precipitation fraction by
     560             :           ! multiplying it by a factor.  That factor is overall precipitation
     561             :           ! fraction divided by the sum of the weighted 1st PDF component
     562             :           ! precipitation fraction and the weighted 2nd PDF component
     563             :           ! precipitation fraction.  The 1st PDF component precipitation
     564             :           ! fraction is then found by dividing the adjusted weighted 1st PDF
     565             :           ! component precipitation fraction by mixture fraction.
     566             :           precip_frac_1(j,k) &
     567             :           = weighted_pfrac_1(j,k) &
     568             :             * ( precip_frac(j,k) &
     569             :                 / ( weighted_pfrac_1(j,k) + weighted_pfrac_2(j,k) ) ) &
     570             :             / mixt_frac(j,k)
     571             :         else
     572             : 
     573             :           ! Usually, the sum of the weighted 1st PDF component precipitation
     574             :           ! fraction and the 2nd PDF component precipitation fraction go to 0
     575             :           ! when overall precipitation fraction goes to 0.  Since 1st PDF
     576             :           ! component weighted precipitation fraction is 0, 1st PDF component
     577             :           ! precipitation fraction also 0.
     578             :           precip_frac_1(j,k) = zero
     579             : 
     580             :         end if
     581             :       end do
     582             :     end do
     583             :       
     584             :     do k = 1, nz, 1
     585             :       do j = 1, ngrdcol
     586             :         ! Special cases for precip_frac_1.
     587             :         if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) &
     588             :             .and. precip_frac_1(j,k) &
     589             :                   > min( one, precip_frac(j,k) / mixt_frac(j,k) ) ) then
     590             : 
     591             :           ! Using the above method, it is possible for precip_frac_1 to be
     592             :           ! greater than 1.  For example, the mixture fraction at level k+1 is
     593             :           ! 0.10 and the cloud_frac_1 at level k+1 is 1, resulting in a
     594             :           ! weighted_pfrac_1 of 0.10.  This product is greater than the product
     595             :           ! of mixt_frac and cloud_frac_1 at level k.  The mixture fraction at
     596             :           ! level k is 0.05, resulting in a precip_frac_1 of 2.  The value of
     597             :           ! precip_frac_1 is limited at 1.  The leftover precipitation fraction
     598             :           ! (a result of the decreasing weight of PDF component 1 between the
     599             :           ! levels) is applied to PDF component 2.
     600             :           ! Additionally, when weighted_pfrac_1 at level k is greater than
     601             :           ! overall precipitation fraction at level k, the resulting calculation
     602             :           ! of precip_frac_2 at level k will be negative.
     603             :           precip_frac_1(j,k) = min( one, precip_frac(j,k) / mixt_frac(j,k) )
     604             : 
     605             :         elseif ( any( hydromet(j,k,:) >= hydromet_tol(:) ) &
     606             :                 .and. precip_frac_1(j,k) > zero &
     607             :                 .and. precip_frac_1(j,k) < precip_frac_tol(j) ) then
     608             : 
     609             :           ! In a scenario where we find precipitation in the 1st PDF component
     610             :           ! (it is allowed to have a value of 0 when all precipitation is found
     611             :           ! in the 2nd PDF component) but it is tiny (less than tolerance
     612             :           ! level), boost 1st PDF component precipitation fraction to tolerance
     613             :           ! level.
     614             :           precip_frac_1(j,k) = precip_frac_tol(j)
     615             : 
     616             :         elseif ( all( hydromet(j,k,:) < hydromet_tol(:) ) ) then
     617             : 
     618             :           ! The means (overall) of every precipitating hydrometeor are all less
     619             :           ! than their respective tolerance amounts.  They are all considered to
     620             :           ! have values of 0.  There are not any hydrometeor species found at
     621             :           ! this grid level.  There is also no cloud at or above this grid
     622             :           ! level, so set 1st component precipitation fraction to 0.
     623             :           precip_frac_1(j,k) = zero
     624             : 
     625             :         end if
     626             :       end do
     627             :     end do
     628             : 
     629             : 
     630             :     !!! Find precipitation fraction within PDF component 2.
     631             :     ! The equation for precipitation fraction within PDF component 2 is:
     632             :     !
     633             :     ! f_p(2) = ( f_p - a * f_p(1) ) / ( 1 - a );
     634             :     !
     635             :     ! given the overall precipitation fraction, f_p (calculated above), the
     636             :     ! precipitation fraction within PDF component 1, f_p(1) (calculated above),
     637             :     ! and mixture fraction, a.  Any leftover precipitation fraction from
     638             :     ! precip_frac_1 will be included in this calculation of precip_frac_2.
     639             :     do k = 1, nz, 1
     640             :       do j = 1, ngrdcol
     641             : 
     642             :         if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
     643             : 
     644             :           ! Calculate precipitation fraction in the 2nd PDF component.
     645             :           precip_frac_2(j,k) &
     646             :           = max( ( precip_frac(j,k) - mixt_frac(j,k) * precip_frac_1(j,k) ) &
     647             :                  / ( one - mixt_frac(j,k) ), &
     648             :                  zero )
     649             : 
     650             :           ! Special cases for precip_frac_2.
     651             :           if ( precip_frac_2(j,k) > one ) then
     652             : 
     653             :              ! Again, it is possible for precip_frac_2 to be greater than 1.
     654             :              ! For example, the mixture fraction at level k+1 is 0.10 and the
     655             :              ! cloud_frac_1 at level k+1 is 1, resulting in a weighted_pfrac_1
     656             :              ! of 0.10.  This product is greater than the product of mixt_frac
     657             :              ! and cloud_frac_1 at level k.  Additionally, precip_frac (overall)
     658             :              ! is 1 for level k.  The mixture fraction at level k is 0.5,
     659             :              ! resulting in a precip_frac_1 of 0.2.  Using the above equation,
     660             :              ! precip_frac_2 is calculated to be 1.8.  The value of
     661             :              ! precip_frac_2 is limited at 1.  The leftover precipitation
     662             :              ! fraction (as a result of the increasing weight of component 1
     663             :              ! between the levels) is applied to PDF component 1.
     664             :              precip_frac_2(j,k) = one
     665             : 
     666             :              ! Recalculate precipitation fraction in the 1st PDF component.
     667             :              precip_frac_1(j,k) &
     668             :              = ( precip_frac(j,k) - ( one - mixt_frac(j,k) ) ) / mixt_frac(j,k)
     669             : 
     670             :              ! Double check precip_frac_1
     671             :              if ( precip_frac_1(j,k) > one ) then
     672             : 
     673             :                 precip_frac_1(j,k) = one
     674             : 
     675             :                 precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
     676             :                                    / ( one - mixt_frac(j,k) )
     677             : 
     678             :              elseif ( precip_frac_1(j,k) > zero .and. precip_frac_1(j,k) < precip_frac_tol(j) ) then
     679             : 
     680             :                 precip_frac_1(j,k) = precip_frac_tol(j)
     681             : 
     682             :                 ! fp = a*fp1+(1-a)*fp2 solving for fp2
     683             :                 precip_frac_2(j,k) = precip_frac_1(j,k) &
     684             :                                      * ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
     685             :                                          - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
     686             :                 
     687             :              end if
     688             : 
     689             :           elseif ( precip_frac_2(j,k) > zero &
     690             :                    .and. precip_frac_2(j,k) < precip_frac_tol(j) ) then
     691             : 
     692             :             ! In a scenario where we find precipitation in the 2nd PDF
     693             :             ! component (it is allowed to have a value of 0 when all
     694             :             ! precipitation is found in the 1st PDF component) but it is tiny
     695             :             ! (less than tolerance level), boost 2nd PDF component
     696             :             ! precipitation fraction to tolerance level.
     697             :             precip_frac_2(j,k) = precip_frac_tol(j)
     698             : 
     699             :             ! Recalculate precipitation fraction in the 1st PDF component.
     700             :             precip_frac_1(j,k) &
     701             :              = ( precip_frac(j,k) - ( one - mixt_frac(j,k) ) * precip_frac_2(j,k) ) &
     702             :                / mixt_frac(j,k)
     703             : 
     704             :             ! Double check precip_frac_1
     705             :             if ( precip_frac_1(j,k) > one ) then
     706             : 
     707             :               precip_frac_1(j,k) = one
     708             : 
     709             :               precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
     710             :                                    / ( one - mixt_frac(j,k) )
     711             : 
     712             :             elseif ( precip_frac_1(j,k) > zero .and. precip_frac_1(j,k) < precip_frac_tol(j) ) then
     713             : 
     714             :               precip_frac_1(j,k) = precip_frac_tol(j) 
     715             : 
     716             :               ! fp = a*fp1+(1-a)*fp2 solving for fp2
     717             :               precip_frac_2(j,k) = precip_frac_1(j,k) &
     718             :                                    * ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
     719             :                                        - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
     720             : 
     721             :             end if
     722             : 
     723             :           end if ! Special cases for precip_frac_2
     724             : 
     725             :         else ! all( hydromet(k,:) < hydromet_tol(:) )
     726             : 
     727             :           ! The means (overall) of every precipitating hydrometeor are all less
     728             :           ! than their respective tolerance amounts.  They are all considered to
     729             :           ! have values of 0.  There are not any hydrometeor species found at
     730             :           ! this grid level.  There is also no cloud at or above this grid
     731             :           ! level, so set 2nd component precipitation fraction to 0.
     732             :           precip_frac_2(j,k) = zero
     733             : 
     734             :         end if ! any( hydromet(k,:) > hydromet_tol(:) )
     735             :       end do
     736             :     end do
     737             : 
     738             :     return
     739             : 
     740             :   end subroutine component_precip_frac_weighted
     741             : 
     742             :   !=============================================================================
     743           0 :   subroutine component_precip_frac_specify( nz, ngrdcol, hydromet_dim, hydromet_tol, &
     744             :                                             upsilon_precip_frac_rat, &
     745           0 :                                             hydromet, precip_frac, &
     746           0 :                                             mixt_frac, precip_frac_tol, &
     747           0 :                                             precip_frac_1, precip_frac_2 )
     748             : 
     749             :     ! Description:
     750             :     ! Calculates the precipitation fraction in each PDF component.
     751             :     !
     752             :     ! The equation for precipitation fraction is:
     753             :     !
     754             :     ! f_p = mixt_frac * f_p(1) + ( 1 - mixt_frac ) * f_p(2);
     755             :     !
     756             :     ! where f_p is overall precipitation fraction, f_p(1) is precipitation
     757             :     ! fraction in the 1st PDF component, f_p(2) is precipitation fraction in the
     758             :     ! 2nd PDF component, and mixt_frac is the mixture fraction.  Using this
     759             :     ! method, a new specified parameter is introduced, upsilon, where:
     760             :     !
     761             :     ! upsilon = mixt_frac * f_p(1) / f_p; and where 0 <= upsilon <= 1.
     762             :     !
     763             :     ! In other words, upsilon is the ratio of mixt_frac * f_p(1) to f_p.  Since
     764             :     ! f_p and mixt_frac are calculated previously, and upsilon is specified,
     765             :     ! f_p(1) can be calculated by:
     766             :     !
     767             :     ! f_p(1) = upsilon * f_p / mixt_frac;
     768             :     !
     769             :     ! and has an upper limit of 1.  The value of f_p(2) can then be calculated
     770             :     ! by:
     771             :     !
     772             :     ! f_p(2) = ( f_p - mixt_frac * f_p(1) ) / ( 1 - mixt_frac );
     773             :     !
     774             :     ! and also has an upper limit of 1.  When upsilon = 1, all of the
     775             :     ! precipitation is found in the 1st PDF component (as long as
     776             :     ! f_p <= mixt_frac, otherwise it would cause f_p(1) to be greater than 1).
     777             :     ! When upsilon = 0, all of the precipitation is found in the 2nd PDF
     778             :     ! component (as long as f_p <= 1 - mixt_frac, otherwise it would cause
     779             :     ! f_p(2) to be greater than 1).  When upsilon is between 0 and 1,
     780             :     ! precipitation is split between the two PDF components accordingly.
     781             : 
     782             :     ! References:
     783             :     !-----------------------------------------------------------------------
     784             : 
     785             :     use constants_clubb, only: &
     786             :         one,  & ! Constant(s)
     787             :         zero, &
     788             :         eps
     789             : 
     790             :     use clubb_precision, only: &
     791             :         core_rknd  ! Variable(s)
     792             : 
     793             :     implicit none
     794             : 
     795             :     ! Input Variables
     796             :     integer, intent(in) :: &
     797             :       nz,           & ! Number of model vertical grid levels
     798             :       ngrdcol,      & ! Number of grid columns
     799             :       hydromet_dim    ! Number of hydrometeor species
     800             : 
     801             :     real( kind = core_rknd ), dimension(hydromet_dim), intent(in) :: &
     802             :       hydromet_tol    ! Tolerance values for all hydrometeors    [units vary]
     803             : 
     804             :     real( kind = core_rknd ), intent(in) :: &
     805             :       upsilon_precip_frac_rat    ! ratio mixt_frac*precip_frac_1/precip_frac [-]
     806             : 
     807             :     real( kind = core_rknd ), dimension(ngrdcol,nz,hydromet_dim), intent(in) :: &
     808             :       hydromet    ! Mean of hydrometeor, hm (overall)           [units vary]
     809             : 
     810             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
     811             :       precip_frac, & ! Precipitation fraction (overall)                      [-]
     812             :       mixt_frac      ! Mixture fraction                                      [-]
     813             : 
     814             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) :: &
     815             :       precip_frac_tol    ! Minimum precip. frac. when hydromet. are present  [-]
     816             : 
     817             :     ! Output Variables
     818             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
     819             :       precip_frac_1, & ! Precipitation fraction (1st PDF component)     [-]
     820             :       precip_frac_2    ! Precipitation fraction (2nd PDF component)     [-]
     821             : 
     822             :     integer :: k, j  ! Loop index.
     823             : 
     824             : 
     825             :     ! There are hydrometeors found at this grid level.
     826           0 :     if ( abs(upsilon_precip_frac_rat - one) < &
     827             :           abs(upsilon_precip_frac_rat + one) / 2 * eps ) then
     828             : 
     829             : 
     830             :       ! Loop over all vertical levels.
     831           0 :       do k = 1, nz, 1
     832           0 :         do j = 1, ngrdcol
     833             : 
     834           0 :           if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
     835             :              
     836           0 :             if ( precip_frac(j,k) <= mixt_frac(j,k) ) then
     837             : 
     838             :               ! All the precipitation is found in the 1st PDF component.
     839           0 :               precip_frac_1(j,k) = precip_frac(j,k) / mixt_frac(j,k)
     840           0 :               precip_frac_2(j,k) = zero
     841             : 
     842             :             else ! precip_frac(k) > mixt_frac(k)
     843             : 
     844             :               ! Some precipitation is found in the 2nd PDF component.
     845           0 :               precip_frac_1(j,k) = one
     846             :               precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
     847           0 :                                  / ( one - mixt_frac(j,k) )
     848             : 
     849             :               if ( precip_frac_2(j,k) > one &
     850           0 :                    .and. abs(precip_frac(j,k) - one) < abs(precip_frac(j,k) + one) / 2 * eps ) then
     851             : 
     852             :                  ! Set precip_frac_2 = 1.
     853           0 :                  precip_frac_2(j,k) = one
     854             : 
     855           0 :               elseif ( precip_frac_2(j,k) < precip_frac_tol(j) ) then
     856             : 
     857             :                 ! Since precipitation is found in the 2nd PDF component, it
     858             :                 ! must have a value of at least precip_frac_tol.
     859           0 :                 precip_frac_2(j,k) = precip_frac_tol(j)
     860             : 
     861             :                 ! Recalculate precip_frac_1
     862             :                 precip_frac_1(j,k) &
     863             :                  = ( precip_frac(j,k) &
     864             :                      - ( one - mixt_frac(j,k) ) * precip_frac_2(j,k) ) &
     865           0 :                    / mixt_frac(j,k)
     866             : 
     867             :                 ! Double check precip_frac_1
     868           0 :                 if ( precip_frac_1(j,k) > one ) then
     869             : 
     870           0 :                   precip_frac_1(j,k) = one
     871             : 
     872             :                   precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
     873           0 :                                        / ( one - mixt_frac(j,k) )
     874             : 
     875           0 :                 elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
     876             : 
     877           0 :                   precip_frac_1(j,k) = precip_frac_tol(j)
     878             : 
     879             :                   ! fp = a*fp1+(1-a)*fp2 solving for fp2
     880             :                   precip_frac_2(j,k) = precip_frac_1(j,k) * &
     881             :                                        ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
     882           0 :                                        - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
     883             : 
     884             :                 end if
     885             : 
     886             :               end if ! precip_frac_2(k) < precip_frac_tol
     887             : 
     888             :             endif ! precip_frac(k) <= mixt_frac(k)
     889             :             
     890             :           else ! all( hydromet(k,:) < hydromet_tol(:) )
     891             :             ! There aren't any hydrometeors found at the grid level.
     892           0 :             precip_frac_1(j,k) = zero
     893           0 :             precip_frac_2(j,k) = zero
     894             :           end if
     895             :           
     896             :         end do
     897             :       end do
     898             : 
     899           0 :     elseif ( abs(upsilon_precip_frac_rat - zero) < &
     900             :           abs(upsilon_precip_frac_rat + zero) / 2 * eps ) then
     901             : 
     902           0 :       do k = 1, nz, 1
     903           0 :         do j = 1, ngrdcol
     904             : 
     905           0 :           if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
     906             :             
     907           0 :             if ( precip_frac(j,k) <= ( one - mixt_frac(j,k) ) ) then
     908             : 
     909             :               ! All the precipitation is found in the 2nd PDF component.
     910           0 :               precip_frac_1(j,k) = zero
     911           0 :               precip_frac_2(j,k) = precip_frac(j,k) / ( one - mixt_frac(j,k) )
     912             : 
     913             :             else ! precip_frac(k) > ( 1 - mixt_frac(k) )
     914             : 
     915             :               ! Some precipitation is found in the 1st PDF component.
     916             :               precip_frac_1(j,k) = ( precip_frac(j,k) - ( one - mixt_frac(j,k) ) ) &
     917           0 :                                     / mixt_frac(j,k)
     918           0 :               precip_frac_2(j,k) = one
     919             : 
     920             :               if ( precip_frac_1(j,k) > one &
     921           0 :                    .and. abs(precip_frac(j,k) - one) < abs(precip_frac(j,k) + one) / 2 * eps ) then
     922             : 
     923             :                 ! Set precip_frac_1 = 1.
     924           0 :                 precip_frac_1(j,k) = one
     925             : 
     926           0 :               elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
     927             : 
     928             :                 ! Since precipitation is found in the 1st PDF component, it
     929             :                 ! must have a value of at least precip_frac_tol.
     930           0 :                 precip_frac_1(j,k) = precip_frac_tol(j)
     931             : 
     932             :                 ! Recalculate precip_frac_2
     933             :                 precip_frac_2(j,k) = ( precip_frac(j,k) &
     934             :                                        - mixt_frac(j,k) * precip_frac_1(j,k) ) &
     935           0 :                                       / ( one - mixt_frac(j,k) )
     936             : 
     937             :                 ! Double check precip_frac_2
     938           0 :                 if ( precip_frac_2(j,k) > one ) then
     939             : 
     940           0 :                   precip_frac_2(j,k) = one
     941             : 
     942             :                   precip_frac_1(j,k) = ( ( precip_frac(j,k) - one ) + mixt_frac(j,k) ) &
     943           0 :                                        / mixt_frac(j,k)
     944             : 
     945           0 :                 elseif ( precip_frac_2(j,k) < precip_frac_tol(j) ) then
     946             : 
     947           0 :                   precip_frac_2(j,k) = precip_frac_tol(j)
     948             : 
     949             :                   ! fp = a*fp1+(1-a)*fp2 solving for fp1
     950             :                   precip_frac_1(j,k) = ( precip_frac(j,k) - precip_frac_2(j,k) ) / mixt_frac(j,k) &
     951           0 :                                        + precip_frac_2(j,k)
     952             : 
     953             :                 endif
     954             : 
     955             :               end if ! precip_frac_1(k) < precip_frac_tol
     956             : 
     957             :             endif ! precip_frac(k) <= ( 1 - mixt_frac(k) )
     958             :             
     959             :           else ! all( hydromet(k,:) < hydromet_tol(:) )
     960             :             ! There aren't any hydrometeors found at the grid level.
     961           0 :             precip_frac_1(j,k) = zero
     962           0 :             precip_frac_2(j,k) = zero
     963             :           end if
     964             :           
     965             :         end do
     966             :       end do
     967             : 
     968             :     else  ! 0 < upsilon_precip_frac_rat < 1
     969             : 
     970           0 :       do k = 1, nz, 1
     971           0 :         do j = 1, ngrdcol
     972             : 
     973           0 :           if ( any( hydromet(j,k,:) >= hydromet_tol(:) ) ) then
     974             :             ! Precipitation is found in both PDF components.  Each component
     975             :             ! must have a precipitation fraction that is at least
     976             :             ! precip_frac_tol and that does not exceed 1.
     977             : 
     978             :             ! Calculate precipitation fraction in the 1st PDF component.
     979             :             precip_frac_1(j,k) &
     980           0 :             = upsilon_precip_frac_rat * precip_frac(j,k) / mixt_frac(j,k)
     981             : 
     982             :             ! Special cases for precip_frac_1
     983           0 :             if ( precip_frac_1(j,k) > one ) then
     984           0 :               precip_frac_1(j,k) = one
     985           0 :             elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
     986           0 :               precip_frac_1(j,k) = precip_frac_tol(j)
     987             :             endif
     988             : 
     989             :             ! Calculate precipitation fraction in the 2nd PDF component.
     990             :             precip_frac_2(j,k) = ( precip_frac(j,k) &
     991             :                                    - mixt_frac(j,k) * precip_frac_1(j,k) ) &
     992           0 :                                  / ( one - mixt_frac(j,k) )
     993             : 
     994             :             ! Special case for precip_frac_2
     995           0 :             if ( precip_frac_2(j,k) > one ) then
     996             : 
     997             :               ! Set precip_frac_2 to 1.
     998           0 :               precip_frac_2(j,k) = one
     999             : 
    1000             :               ! Recalculate precipitation fraction in the 1st PDF component.
    1001             :               precip_frac_1(j,k) &
    1002           0 :               = ( precip_frac(j,k) - ( one - mixt_frac(j,k) ) ) / mixt_frac(j,k)
    1003             : 
    1004             :               ! Double check precip_frac_1
    1005           0 :               if ( precip_frac_1(j,k) > one ) then
    1006             : 
    1007           0 :                 precip_frac_1(j,k) = one
    1008             : 
    1009             :                 precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
    1010           0 :                                      / ( one - mixt_frac(j,k) )
    1011             : 
    1012           0 :               elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
    1013             : 
    1014           0 :                 precip_frac_1(j,k) = precip_frac_tol(j)
    1015             : 
    1016             :                 ! fp = a*fp1+(1-a)*fp2 solving for fp2
    1017             :                 precip_frac_2(j,k) = precip_frac_1(j,k) &
    1018             :                                      * ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
    1019           0 :                                      - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
    1020             : 
    1021             :               endif
    1022             : 
    1023           0 :             elseif ( precip_frac_2(j,k) < precip_frac_tol(j) ) then
    1024             : 
    1025             :               ! Set precip_frac_2 to precip_frac_tol.
    1026           0 :               precip_frac_2(j,k) = precip_frac_tol(j)
    1027             : 
    1028             :               ! Recalculate precipitation fraction in the 1st PDF component.
    1029             :               precip_frac_1(j,k) = ( precip_frac(j,k) &
    1030             :                                      - ( one - mixt_frac(j,k) ) * precip_frac_2(j,k) ) &
    1031           0 :                                    / mixt_frac(j,k)
    1032             : 
    1033             :               ! Double check precip_frac_1
    1034           0 :               if ( precip_frac_1(j,k) > one ) then
    1035             : 
    1036           0 :                 precip_frac_1(j,k) = one
    1037             : 
    1038             :                 precip_frac_2(j,k) = ( precip_frac(j,k) - mixt_frac(j,k) ) &
    1039           0 :                                      / ( one - mixt_frac(j,k) )
    1040             : 
    1041           0 :               elseif ( precip_frac_1(j,k) < precip_frac_tol(j) ) then
    1042             : 
    1043           0 :                 precip_frac_1(j,k) = precip_frac_tol(j)
    1044             : 
    1045             :                 ! fp = a*fp1+(1-a)*fp2 solving for fp2
    1046             :                 precip_frac_2(j,k) = precip_frac_1(j,k) &
    1047             :                                      * ( ( ( precip_frac(j,k) / precip_frac_1(j,k)) &
    1048           0 :                                      - mixt_frac(j,k) ) / ( one - mixt_frac(j,k) ) )
    1049             : 
    1050             :               end if
    1051             : 
    1052             :             end if ! Special cases for precip_frac_2
    1053             :          
    1054             :           else ! all( hydromet(k,:) < hydromet_tol(:) )
    1055             :             ! There aren't any hydrometeors found at the grid level.
    1056           0 :             precip_frac_1(j,k) = zero
    1057           0 :             precip_frac_2(j,k) = zero
    1058             :           end if
    1059             :           
    1060             :         end do
    1061             :       end do
    1062             : 
    1063             :     end if  ! upsilon_precip_frac_rat
    1064             : 
    1065           0 :     return
    1066             : 
    1067             :   end subroutine component_precip_frac_specify
    1068             : 
    1069             :   !=============================================================================
    1070           0 :   subroutine precip_frac_assert_check( nz, hydromet_dim, hydromet_tol, &
    1071           0 :                                        hydromet, mixt_frac, precip_frac, &
    1072           0 :                                        precip_frac_1, precip_frac_2, &
    1073             :                                        precip_frac_tol )
    1074             : 
    1075             :     ! Description:
    1076             :     ! Assertion check for the precipitation fraction code.
    1077             : 
    1078             :     ! References:
    1079             :     !-----------------------------------------------------------------------
    1080             : 
    1081             :     use constants_clubb, only: &
    1082             :         one,     & ! Constant(s)
    1083             :         zero,    &
    1084             :         fstderr, &
    1085             :         eps
    1086             : 
    1087             :     use clubb_precision, only: &
    1088             :         core_rknd  ! Variable(s)
    1089             : 
    1090             :     use error_code, only: &
    1091             :         err_code, &                     ! Error Indicator
    1092             :         clubb_fatal_error               ! Constant
    1093             : 
    1094             :     implicit none
    1095             : 
    1096             :     ! Input Variables
    1097             :     integer, intent(in) :: &
    1098             :       nz,           & ! Number of model vertical grid levels
    1099             :       hydromet_dim    ! Number of hydrometeor species
    1100             : 
    1101             :     real( kind = core_rknd ), dimension(hydromet_dim), intent(in) :: &
    1102             :       hydromet_tol    ! Tolerance values for all hydrometeors    [units vary]
    1103             : 
    1104             :     real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(in) :: &
    1105             :       hydromet    ! Mean of hydrometeor, hm (overall)           [units vary]
    1106             : 
    1107             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
    1108             :       mixt_frac,     & ! Mixture fraction                               [-]
    1109             :       precip_frac,   & ! Precipitation fraction (overall)               [-]
    1110             :       precip_frac_1, & ! Precipitation fraction (1st PDF component)     [-]
    1111             :       precip_frac_2    ! Precipitation fraction (2nd PDF component)     [-]
    1112             : 
    1113             :     real( kind = core_rknd ), intent(in) :: &
    1114             :       precip_frac_tol    ! Minimum precip. frac. when hydromet. are present  [-]
    1115             : 
    1116             :     ! Local Variables
    1117             :     integer :: k  ! Loop index
    1118             : 
    1119             : 
    1120             :     ! Loop over all vertical levels.
    1121           0 :     do k = 1, nz, 1
    1122             : 
    1123           0 :        if ( any( hydromet(k,:) >= hydromet_tol(:) ) ) then
    1124             : 
    1125             :           ! Overall precipitation fraction cannot be less than precip_frac_tol
    1126             :           ! when a hydrometeor is present at a grid level.
    1127           0 :           if ( precip_frac(k) < precip_frac_tol ) then
    1128             :              write(fstderr,*) "precip_frac < precip_frac_tol when " &
    1129           0 :                               // "a hydrometeor is present"
    1130           0 :              write(fstderr,*) "level = ", k
    1131           0 :              write(fstderr,*) "precip_frac = ", precip_frac(k), &
    1132           0 :                               "precip_frac_tol = ", precip_frac_tol
    1133           0 :              err_code = clubb_fatal_error
    1134           0 :              return
    1135             :           endif
    1136             : 
    1137             :           ! Overall precipitation fraction cannot exceed 1.
    1138           0 :           if ( precip_frac(k) > one ) then
    1139           0 :              write(fstderr,*) "precip_frac > 1"
    1140           0 :              write(fstderr,*) "level = ", k
    1141           0 :              write(fstderr,*) "precip_frac = ", precip_frac(k)
    1142           0 :              err_code = clubb_fatal_error
    1143           0 :              return
    1144             :           endif
    1145             : 
    1146             :           ! Precipitation fraction in the 1st PDF component is allowed to be 0
    1147             :           ! when all the precipitation is found in the 2nd PDF component.
    1148             :           ! Otherwise, it cannot be less than precip_frac_tol when a hydrometeor
    1149             :           ! is present at a grid level.  In other words, it cannot have a value
    1150             :           ! that is greater than 0 but less than precip_frac_tol
    1151             :           if ( precip_frac_1(k) > zero &
    1152           0 :                .and. precip_frac_1(k) < precip_frac_tol ) then
    1153           0 :              write(fstderr,*) "0 < precip_frac_1 < precip_frac_tol"
    1154           0 :              write(fstderr,*) "level = ", k
    1155           0 :              write(fstderr,*) "precip_frac_1 = ", precip_frac_1(k), &
    1156           0 :                               "precip_frac_tol = ", precip_frac_tol
    1157           0 :              err_code = clubb_fatal_error
    1158           0 :              return
    1159             :           endif
    1160             : 
    1161             :           ! Precipitation fraction in the 1st PDF component cannot exceed 1.
    1162           0 :           if ( precip_frac_1(k) > one ) then
    1163           0 :              write(fstderr,*) "precip_frac_1 > 1"
    1164           0 :              write(fstderr,*) "level = ", k
    1165           0 :              write(fstderr,*) "precip_frac_1 = ", precip_frac_1(k)
    1166           0 :              err_code = clubb_fatal_error
    1167           0 :              return
    1168             :           endif
    1169             : 
    1170             :           ! Precipiation fraction in the 1st PDF component cannot be negative.
    1171           0 :           if ( precip_frac_1(k) < zero ) then
    1172           0 :              write(fstderr,*) "precip_frac_1 < 0"
    1173           0 :              write(fstderr,*) "level = ", k
    1174           0 :              write(fstderr,*) "precip_frac_1 = ", precip_frac_1(k)
    1175           0 :              err_code = clubb_fatal_error
    1176           0 :              return
    1177             :           endif
    1178             : 
    1179             :           ! Precipitation fraction in the 2nd PDF component is allowed to be 0
    1180             :           ! when all the precipitation is found in the 1st PDF component.
    1181             :           ! Otherwise, it cannot be less than precip_frac_tol when a hydrometeor
    1182             :           ! is present at a grid level.  In other words, it cannot have a value
    1183             :           ! that is greater than 0 but less than precip_frac_tol
    1184             :           if ( precip_frac_2(k) > zero &
    1185           0 :                .and. precip_frac_2(k) < precip_frac_tol ) then
    1186           0 :              write(fstderr,*) "0 < precip_frac_2 < precip_frac_tol"
    1187           0 :              write(fstderr,*) "level = ", k
    1188           0 :              write(fstderr,*) "precip_frac_2 = ", precip_frac_2(k), &
    1189           0 :                               "precip_frac_tol = ", precip_frac_tol
    1190           0 :              err_code = clubb_fatal_error
    1191           0 :              return
    1192             :           endif
    1193             : 
    1194             :           ! Precipitation fraction in the 2nd PDF component cannot exceed 1.
    1195           0 :           if ( precip_frac_2(k) > one ) then
    1196           0 :              write(fstderr,*) "precip_frac_2 > 1"
    1197           0 :              write(fstderr,*) "level = ", k
    1198           0 :              write(fstderr,*) "precip_frac_2 = ", precip_frac_2(k)
    1199           0 :              err_code = clubb_fatal_error
    1200           0 :              return
    1201             :           endif
    1202             : 
    1203             :           ! Precipiation fraction in the 2nd PDF component cannot be negative.
    1204           0 :           if ( precip_frac_2(k) < zero ) then
    1205           0 :              write(fstderr,*) "precip_frac_2 < 0"
    1206           0 :              write(fstderr,*) "level = ", k
    1207           0 :              write(fstderr,*) "precip_frac_2 = ", precip_frac_2(k)
    1208           0 :              err_code = clubb_fatal_error
    1209           0 :              return
    1210             :           endif
    1211             : 
    1212             :        else  ! all( hydromet(k,:) < hydromet_tol(:) )
    1213             : 
    1214             :           ! Overall precipitation fraction must be 0 when no hydrometeors are
    1215             :           ! found at a grid level.
    1216           0 :           if ( abs(precip_frac(k)) > eps) then
    1217           0 :              write(fstderr,*) "precip_frac /= 0 when no hydrometeors are found"
    1218           0 :              write(fstderr,*) "level = ", k
    1219           0 :              write(fstderr,*) "precip_frac = ", precip_frac(k)
    1220           0 :              err_code = clubb_fatal_error
    1221           0 :              return
    1222             :           endif
    1223             : 
    1224             :           ! Precipitation fraction in the 1st PDF component must be 0 when no
    1225             :           ! hydrometeors are found at a grid level.
    1226           0 :           if ( abs(precip_frac_1(k)) > eps) then
    1227             :              write(fstderr,*) "precip_frac_1 /= 0 when no hydrometeors " &
    1228           0 :                               // "are found"
    1229           0 :              write(fstderr,*) "level = ", k
    1230           0 :              write(fstderr,*) "precip_frac_1 = ", precip_frac_1(k)
    1231           0 :              err_code = clubb_fatal_error
    1232           0 :              return
    1233             :           endif
    1234             : 
    1235             :           ! Precipitation fraction in the 2nd PDF component must be 0 when no
    1236             :           ! hydrometeors are found at a grid level.
    1237           0 :           if ( abs(precip_frac_2(k)) > eps) then
    1238             :              write(fstderr,*) "precip_frac_2 /= 0 when no hydrometeors " &
    1239           0 :                               // "are found"
    1240           0 :              write(fstderr,*) "level = ", k
    1241           0 :              write(fstderr,*) "precip_frac_2 = ", precip_frac_2(k)
    1242           0 :              err_code = clubb_fatal_error
    1243           0 :              return
    1244             :           endif
    1245             : 
    1246             :        endif  ! any( hydromet(k,:) >= hydromet_tol(:) )
    1247             : 
    1248             :        ! The precipitation fraction equation is:
    1249             :        !
    1250             :        ! precip_frac
    1251             :        !    = mixt_frac * precip_frac_1 + ( 1 - mixt_frac ) * precip_frac_2;
    1252             :        !
    1253             :        ! which means that:
    1254             :        !
    1255             :        ! precip_frac
    1256             :        ! - ( mixt_frac * precip_frac_1 + ( 1 - mixt_frac ) * precip_frac_2 )
    1257             :        ! = 0.
    1258             :        !
    1259             :        ! Check that this is true with numerical round off.
    1260           0 :        if ( ( precip_frac(k) &
    1261             :               - ( mixt_frac(k) * precip_frac_1(k) &
    1262             :                   + ( one - mixt_frac(k) ) * precip_frac_2(k) ) ) &
    1263           0 :             > ( epsilon( precip_frac(k) ) * precip_frac(k) ) ) then
    1264             :           write(fstderr,*) "mixt_frac * precip_frac_1 " &
    1265             :                            // "+ ( 1 - mixt_frac ) * precip_frac_2 " &
    1266           0 :                            // "/= precip_frac within numerical roundoff"
    1267           0 :           write(fstderr,*) "level = ", k
    1268             :           write(fstderr,*) "mixt_frac * precip_frac_1 " &
    1269           0 :                            // "+ ( 1 - mixt_frac ) * precip_frac_2 = ", &
    1270           0 :                            mixt_frac(k) * precip_frac_1(k) &
    1271           0 :                            + ( one - mixt_frac(k) ) * precip_frac_2(k)
    1272           0 :           write(fstderr,*) "precip_frac = ", precip_frac(k)
    1273           0 :           err_code = clubb_fatal_error
    1274           0 :           return
    1275             :        endif
    1276             : 
    1277             :     enddo  ! k = 1, nz, 1
    1278             : 
    1279             : 
    1280             :     return
    1281             : 
    1282             :   end subroutine precip_frac_assert_check
    1283             : 
    1284             : !===============================================================================
    1285             : 
    1286             : end module precipitation_fraction

Generated by: LCOV version 1.14