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