LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - fill_holes.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 65 | 272 | 23.9 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 7 | 14.3 %

              Line data    Source code
           1             : !-----------------------------------------------------------------------
           2             : ! $Id: fill_holes.F90 8738 2018-07-19 19:58:53Z bmg2@uwm.edu $
           3             : !===============================================================================
           4             : module fill_holes
           5             :
           6             :   implicit none
           7             :
           8             :   public :: fill_holes_driver, &
           9             :             fill_holes_vertical, &
          10             :             hole_filling_hm_one_lev, &
          11             :             fill_holes_hydromet, &
          12             :             fill_holes_wv, &
          13             :             clip_hydromet_conc_mvr, &
          14             :             setup_stats_indices
          15             :
          16             :   private ! Set Default Scope
          17             :
          18             :   contains
          19             :
          20             :   !=============================================================================
          21   268051680 :   subroutine fill_holes_vertical( nz, ngrdcol, threshold, upper_hf_level, &
          22   268051680 :                                   dz, rho_ds, &
          23   268051680 :                                   field )
          24             :
          25             :     ! Description:
          26             :     !   This subroutine clips values of 'field' that are below 'threshold' as much
          27             :     !   as possible (i.e. "fills holes"), but conserves the total integrated mass
          28             :     !   of 'field'.  This prevents clipping from acting as a spurious source.
          29             :     !
          30             :     !   Mass is conserved by reducing the clipped field everywhere by a constant
          31             :     !   multiplicative coefficient.
          32             :     !
          33             :     !   This subroutine does not guarantee that the clipped field will exceed
          34             :     !   threshold everywhere; blunt clipping is needed for that.
          35             :     !
          36             :     !   The lowest level (k=1) should not be included, as the hole-filling scheme
          37             :     !   should not alter the set value of 'field' at the surface (for momentum
          38             :     !   level variables), or consider the value of 'field' at a level below the
          39             :     !   surface (for thermodynamic level variables).
          40             :     !
          41             :     !   For momentum level variables only, the hole-filling scheme should not
          42             :     !   alter the set value of 'field' at the upper boundary level (k=nz).
          43             :     !   So for momemtum level variables, call with upper_hf_level=nz-1, and
          44             :     !   for thermodynamic level variables, call with upper_hf_level=nz.
          45             :     !
          46             :     ! References:
          47             :     !   ``Numerical Methods for Wave Equations in Geophysical Fluid
          48             :     !     Dynamics'', Durran (1999), p. 292.
          49             :     !-----------------------------------------------------------------------
          50             :
          51             :     use grid_class, only: &
          52             :         grid ! Type
          53             :
          54             :     use clubb_precision, only: &
          55             :         core_rknd ! Variable(s)
          56             :
          57             :     use constants_clubb, only: &
          58             :         eps, &
          59             :         one, &
          60             :         num_hf_draw_points ! The number of points on either side of the hole;
          61             :                            ! Mass is drawn from these points to fill the hole
          62             :
          63             :     implicit none
          64             :
          65             :     ! --------------------- Input variables ---------------------
          66             :     integer, intent(in) :: &
          67             :       nz, &
          68             :       ngrdcol
          69             :
          70             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
          71             :       dz      ! Spacing between thermodynamic grid levels; centered over
          72             :               ! momentum grid levels
          73             :               ! OR
          74             :               ! Spcaing between momentum grid levels; centered over
          75             :               ! thermodynamic grid levels
          76             :
          77             :     integer, intent(in) :: &
          78             :       upper_hf_level  ! Upper grid level of global hole-filling range      []
          79             :
          80             :     real( kind = core_rknd ), intent(in) :: &
          81             :       threshold  ! A threshold (e.g. w_tol*w_tol) below which field must not
          82             :                  ! fall                           [Units vary; same as field]
          83             :
          84             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) ::  &
          85             :       rho_ds       ! Dry, static density on thermodynamic or momentum levels    [kg/m^3]
          86             :
          87             :     ! --------------------- Input/Output variable ---------------------
          88             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(inout) :: &
          89             :       field  ! The field (e.g. wp2) that contains holes [Units same as threshold]
          90             :
          91             :     ! --------------------- Local Variables ---------------------
          92             :     integer :: &
          93             :       i,             & ! Loop index for column                           []
          94             :       k,             & ! Loop index for absolute grid level              []
          95             :       j, &
          96             :       k_start,     & ! Lower grid level of local hole-filling range    []
          97             :       k_end          ! Upper grid level of local hole-filling range    []
          98             :
          99             :     real( kind = core_rknd ), dimension(ngrdcol,nz)  ::  &
         100   536103360 :       rho_ds_dz,            & ! rho_ds * dz
         101   536103360 :       invrs_denom_integral, & ! Inverse of the integral in the denominator (see description)
         102   536103360 :       field_clipped           ! The raw field (e.g. wp2) that contains no holes
         103             :                               !                          [Units same as threshold]
         104             :
         105             :     real( kind = core_rknd ) ::  &
         106             :       field_avg,          & ! Vertical average of field [Units of field]
         107             :       field_clipped_avg,  & ! Vertical average of clipped field [Units of field]
         108             :       mass_fraction         ! Coefficient that multiplies clipped field
         109             :                             ! in order to conserve mass.                      []
         110             :
         111             :     real( kind = core_rknd ), dimension(ngrdcol) ::  &
         112   536103360 :       denom_integral_global,  & ! Integral in the denominator for global filling
         113   536103360 :       numer_integral_global,  & ! Integral in the numerator for global filling
         114   536103360 :       field_avg_global,       & ! Vertical average of field [Units of field]
         115   536103360 :       mass_fraction_global      ! Coefficient that multiplies clipped field
         116             :                                 ! in order to conserve mass.                      []
         117             :
         118             :     logical :: &
         119             :       l_field_below_threshold
         120             :
         121             :     real( kind = core_rknd ) ::  &
         122             :       rho_k_sum
         123             :
         124             :     ! --------------------- Begin Code ---------------------
         125             :
         126   268051680 :     l_field_below_threshold = .false.
         127             :
         128             :     !$acc parallel loop gang vector collapse(2) default(present) &
         129             :     !$acc          reduction(.or.:l_field_below_threshold)
         130 23052444480 :     do k = 1, nz
         131 >38071*10^7 :       do i = 1, ngrdcol
         132 >38044*10^7 :         if ( field(i,k) < threshold ) then
         133  8912262547 :           l_field_below_threshold = .true.
         134             :         end if
         135             :       end do
         136             :     end do
         137             :     !$acc end parallel loop
         138             :
         139             :     ! If all field values are above the specified threshold, no hole filling is required
         140   268051680 :     if ( .not. l_field_below_threshold ) then
         141             :       return
         142             :     end if
         143             :
         144             :     !$acc enter data create( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
         145             :     !$acc                    numer_integral_global, field_avg_global, mass_fraction_global )
         146             :
         147             :     !$acc parallel loop gang vector collapse(2) default(present)
         148  2411213304 :     do k = 1, nz
         149 39834224419 :       do i = 1, ngrdcol
         150 39806187055 :         rho_ds_dz(i,k) = rho_ds(i,k) * dz(i,k)
         151             :       end do
         152             :     end do
         153             :     !$acc end parallel loop
         154             :
         155             :
         156             :     ! denom_integral does not change throughout the hole filling algorithm
         157             :     ! so we can calculate it before hand. This results in unneccesary computations,
         158             :     ! but is parallelizable and reduces the cost of the serial k loop
         159             :     !$acc parallel loop gang vector collapse(2) default(present)
         160   468308083 :     do i = 1, ngrdcol
         161 35255156657 :       do k = 2+num_hf_draw_points, upper_hf_level-num_hf_draw_points
         162             :
         163             :         ! This loop and division could be written more compactly as
         164             :         !   invrs_denom_integral(i,k) = one / sum(rho_ds_dz(i,k-num_hf_draw_points:k+num_hf_draw_points))
         165             :         ! but has been manually written in loop form to improve performance
         166             :         ! when using OpenMP target offloading
         167             :         ! See: https://github.com/larson-group/clubb/issues/1138#issuecomment-1974918151
         168 34786848574 :         rho_k_sum = 0.0_core_rknd
         169             :
         170 >20872*10^7 :         do j = k - num_hf_draw_points, k + num_hf_draw_points
         171 >20872*10^7 :           rho_k_sum = rho_k_sum + rho_ds_dz(i,j)
         172             :         end do
         173             :
         174 35227119293 :         invrs_denom_integral(i,k) = one / rho_k_sum
         175             :       end do
         176             :     end do
         177             :     !$acc end parallel loop
         178             :
         179             :     !$acc parallel loop gang vector default(present)
         180   468308083 :     do i = 1, ngrdcol
         181             :
         182             :       ! Make one pass up the profile, filling holes as much as we can using
         183             :       ! nearby mass, ignoring the first level.
         184             :       !
         185             :       ! This loop can only be done in serial due to the field values required for the next
         186             :       ! iteration potentially changing. We could in theory expose more parallelism in cases
         187             :       ! where there are large enough gaps between vertical levels which need hole-filling,
         188             :       ! but levels which require hole-filling are often close or consecutive.
         189 35255156657 :       do k = 2+num_hf_draw_points, upper_hf_level-num_hf_draw_points
         190             :
         191 34786848574 :         k_start = k - num_hf_draw_points
         192 34786848574 :         k_end   = k + num_hf_draw_points
         193             :
         194 >17012*10^7 :         if ( any( field(i,k_start:k_end) < threshold ) ) then
         195             :
         196             :           ! Compute the field's vertical average cenetered at k, which we must conserve,
         197             :           ! see description of the vertical_avg function in advance_helper_module
         198             :           field_avg = sum( rho_ds_dz(i,k_start:k_end) * field(i,k_start:k_end) ) &
         199 55432201566 :                       * invrs_denom_integral(i,k)
         200             :
         201             :           ! Clip small or negative values from field.
         202  9238700261 :           if ( field_avg >= threshold ) then
         203             :             ! We know we can fill in holes completely
         204  7504203996 :             field_clipped(i,k_start:k_end) = max( threshold, field(i,k_start:k_end) )
         205             :           else
         206             :             ! We can only fill in holes partly;
         207             :             ! to do so, we remove all mass above threshold.
         208 47927997570 :             field_clipped(i,k_start:k_end) = min( threshold, field(i,k_start:k_end) )
         209             :           endif
         210             :
         211             :           ! Compute the clipped field's vertical integral.
         212             :           ! clipped_total_mass >= original_total_mass,
         213             :           ! see description of the vertical_avg function in advance_helper_module
         214             :           field_clipped_avg = sum( rho_ds_dz(i,k_start:k_end) * field_clipped(i,k_start:k_end) ) &
         215 55432201566 :                               * invrs_denom_integral(i,k)
         216             :
         217             :           ! Avoid divide by zero issues by doing nothing if field_clipped_avg ~= threshold
         218  9238700261 :           if ( abs(field_clipped_avg-threshold) > abs(field_clipped_avg+threshold)*eps/2) then
         219             :             ! Compute coefficient that makes the clipped field have the same mass as the
         220             :             ! original field.  We should always have mass_fraction > 0.
         221             :             mass_fraction = ( field_avg - threshold ) &
         222  9238700143 :                             / ( field_clipped_avg - threshold )
         223             :
         224             :             ! Calculate normalized, filled field
         225             :             field(i,k_start:k_end) = threshold &
         226 55432200858 :                         + mass_fraction * ( field_clipped(i,k_start:k_end) - threshold )
         227             :           endif
         228             :
         229             :         endif
         230             :
         231             :       end do
         232             :
         233             :     end do
         234             :     !$acc end parallel loop
         235             :
         236             :     l_field_below_threshold = .false.
         237             :
         238             :     !$acc parallel loop gang vector collapse(2) default(present) reduction(.or.:l_field_below_threshold)
         239  2411213304 :     do k = 1, nz
         240 39834224419 :       do i = 1, ngrdcol
         241 39806187055 :         if ( field(i,k) < threshold ) then
         242  7527583684 :           l_field_below_threshold = .true.
         243             :         end if
         244             :       end do
         245             :     end do
         246             :     !$acc end parallel loop
         247             :
         248             :     ! If all field values are above the threshold, no further hole filling is required
         249    28037364 :     if ( .not. l_field_below_threshold ) then
         250             :       !$acc exit data delete( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
         251             :       !$acc                   numer_integral_global, field_avg_global, mass_fraction_global )
         252             :       return
         253             :     end if
         254             :
         255             :
         256             :     ! Now we fill holes globally to maximize the chance that all holes are filled.
         257             :     ! To improve parallelism we assume that global hole filling needs to be done
         258             :     ! for each grid column, perform all calculations required, and only check
         259             :     ! if any holes need filling before the final step of updating the field.
         260             :
         261             :     ! Compute the numerator and denominator integrals
         262             :     !$acc parallel loop gang vector default(present)
         263   188450527 :     do i = 1, ngrdcol
         264   177169619 :       numer_integral_global(i) = 0.0_core_rknd
         265   188450527 :       denom_integral_global(i) = 0.0_core_rknd
         266             :     end do
         267             :     !$acc end parallel loop
         268             :
         269             :     !$acc parallel loop gang vector default(present)
         270   188450527 :     do i = 1, ngrdcol
         271 14893607254 :       do k = 2, upper_hf_level
         272 14705156727 :         numer_integral_global(i) = numer_integral_global(i) + rho_ds_dz(i,k) * field(i,k)
         273             :
         274 14882326346 :         denom_integral_global(i) = denom_integral_global(i) + rho_ds_dz(i,k)
         275             :       end do
         276             :     end do
         277             :     !$acc end parallel loop
         278             :
         279             :
         280             :     !$acc parallel loop gang vector default(present)
         281   188450527 :     do i = 1, ngrdcol
         282             :
         283             :       ! Find the vertical average of field, using the precomputed numerator and denominator,
         284             :       ! see description of the vertical_avg function in advance_helper_module
         285   177169619 :       field_avg_global(i) = numer_integral_global(i) / denom_integral_global(i)
         286             :
         287             :       ! Clip small or negative values from field.
         288   188450527 :       if ( field_avg_global(i) >= threshold ) then
         289             :         ! We know we can fill in holes completely
         290 14841272858 :         field_clipped(i,2:upper_hf_level) = max( threshold, field(i,2:upper_hf_level) )
         291             :       else
         292             :         ! We can only fill in holes partly;
         293             :         ! to do so, we remove all mass above threshold.
         294    41053488 :         field_clipped(i,2:upper_hf_level) = min( threshold, field(i,2:upper_hf_level) )
         295             :       endif
         296             :
         297             :     end do
         298             :     !$acc end parallel loop
         299             :
         300             :     ! To compute the clipped field's vertical integral we only need to recompute the numerator
         301             :     !$acc parallel loop gang vector default(present)
         302   188450527 :     do i = 1, ngrdcol
         303   177169619 :       numer_integral_global(i) = 0.0_core_rknd
         304 14893607254 :       do k = 2, upper_hf_level
         305 14882326346 :         numer_integral_global(i) = numer_integral_global(i) + rho_ds_dz(i,k) * field_clipped(i,k)
         306             :       end do
         307             :     end do
         308             :     !$acc end parallel loop
         309             :
         310             :     !$acc parallel loop gang vector default(present)
         311   188450527 :     do i = 1, ngrdcol
         312             :
         313             :       ! Do not complete calculations or update field values for this
         314             :       ! column if there are no holes that need filling
         315  6072990059 :       if ( any( field(i,2:upper_hf_level) < threshold ) ) then
         316             :
         317             :         ! Compute the clipped field's vertical integral,
         318             :         ! see description of the vertical_avg function in advance_helper_module
         319   145255352 :         field_clipped_avg = numer_integral_global(i) / denom_integral_global(i)
         320             :
         321             :         ! Avoid divide by zero issues by doing nothing if field_clipped_avg ~= threshold
         322   145255352 :         if ( abs(field_clipped_avg-threshold) > abs(field_clipped_avg+threshold)*eps/2) then
         323             :
         324             :           ! Compute coefficient that makes the clipped field have the same mass as the
         325             :           ! original field.  We should always have mass_fraction > 0.
         326             :           mass_fraction_global(i) = ( field_avg_global(i) - threshold ) &
         327   145255352 :                                     / ( field_clipped_avg - threshold )
         328             :
         329             :           ! Calculate normalized, filled field
         330             :           field(i,2:upper_hf_level) = threshold + mass_fraction_global(i) &
         331             :                                                   * ( field_clipped(i,2:upper_hf_level) &
         332 12201449568 :                                                       - threshold )
         333             :         end if
         334             :       end if
         335             :
         336             :     end do
         337             :     !$acc end parallel loop
         338             :
         339             :     !$acc exit data delete( invrs_denom_integral, field_clipped, denom_integral_global, rho_ds_dz, &
         340             :     !$acc                   numer_integral_global, field_avg_global, mass_fraction_global )
         341             :
         342             :     return
         343             :
         344             :   end subroutine fill_holes_vertical
         345             :
         346             :   !===============================================================================
         347           0 :   subroutine hole_filling_hm_one_lev( num_hm_fill, hm_one_lev, & ! Intent(in)
         348           0 :                                    hm_one_lev_filled ) ! Intent(out)
         349             :
         350             :   ! Description:
         351             :   ! Fills holes between same-phase (i.e. either liquid or frozen) hydrometeors for
         352             :   ! one height level.
         353             :   !
         354             :   ! Warning: Do not input hydrometeors of different phases, e.g. liquid and frozen.
         355             :   ! Otherwise heat will not be conserved.
         356             :   !
         357             :   ! References:
         358             :   !
         359             :   ! None
         360             :   !-----------------------------------------------------------------------
         361             :
         362             :     use constants_clubb, only: &
         363             :         one, & ! Variable(s)
         364             :         zero, &
         365             :         eps
         366             :
         367             :     use clubb_precision, only: &
         368             :         core_rknd ! Variable(s)
         369             :
         370             :     use error_code, only: &
         371             :         clubb_at_least_debug_level  ! Procedure
         372             :
         373             :     implicit none
         374             :
         375             :     ! Input Variables
         376             :     integer, intent(in) :: num_hm_fill ! number of hydrometeors involved
         377             :
         378             :     real(kind = core_rknd), dimension(num_hm_fill), intent(in) :: hm_one_lev
         379             :
         380             :     ! Output Variables
         381             :     real(kind = core_rknd), dimension(num_hm_fill), intent(out) :: hm_one_lev_filled
         382             :
         383             :     ! Local Variables
         384             :     integer :: num_neg_hm ! number of holes
         385             :
         386             :     real(kind = core_rknd) :: &
         387             :       total_hole, & ! Size of the hole ( missing mass, less than 0 )
         388             :       total_mass    ! Total mass to fill the hole
         389             :       ! total mass of water substance = total_mass + total_hole
         390             :
         391             :     integer :: i ! loop iterator
         392             :
         393             :   !-----------------------------------------------------------------------
         394             :
         395             :     !----- Begin Code -----
         396             :
         397             :     ! Initialization
         398           0 :     hm_one_lev_filled = 0._core_rknd
         399             :     total_hole = 0._core_rknd
         400             :     total_mass = 0._core_rknd
         401             :     num_neg_hm = 0
         402             :
         403             :     ! Determine the total size of the hole and the number of neg. hydrometeors
         404             :     ! and the total mass of hole filling material
         405           0 :     do i=1, num_hm_fill
         406             : !       print *, "hm_one_lev(",i,") = ", hm_one_lev(i)
         407           0 :        if ( hm_one_lev(i) < zero ) then
         408           0 :           total_hole = total_hole + hm_one_lev(i) ! less than zero
         409             :           num_neg_hm = num_neg_hm + 1
         410             :        else
         411           0 :           total_mass = total_mass + hm_one_lev(i)
         412             :        endif
         413             :
         414             :     enddo
         415             :
         416             : !    print *, "total_hole = ", total_hole
         417             : !    print *, "total_mass = ", total_mass
         418             : !    print *, "num_neg_hm = ", num_neg_hm
         419             :
         420             :     ! There is no water substance at all to fill the hole
         421           0 :     if ( abs(total_mass) < eps ) then
         422             :
         423           0 :        if ( clubb_at_least_debug_level( 2 ) ) then
         424           0 :           print *, "Warning: One-level hole filling was not successful! total_mass ~= 0"
         425             :        endif
         426             :
         427           0 :        hm_one_lev_filled = hm_one_lev
         428             :
         429             :        return
         430             :     endif
         431             :
         432             :     ! Fill the holes and adjust the remaining quantities:
         433             :     ! hm_filled(i) = 0, if hm(i) < 0
         434             :     ! or
         435             :     ! hm_filled(i) = (1 + total_hole/total_mass)*hm(i), if hm(i) > 0
         436           0 :     do i=1, num_hm_fill
         437             :
         438             :        ! if there is not enough material, fill the holes partially with all the material available
         439           0 :        if ( abs(total_hole) > total_mass ) then
         440             :
         441           0 :           if ( clubb_at_least_debug_level( 2 ) ) then
         442             :              print *, "Warning: One-level hole filling was not able to fill holes completely!" // &
         443           0 :                       " The holes were filled partially. |total_hole| > total_mass"
         444             :           endif
         445             :
         446           0 :           hm_one_lev_filled(i) = min(hm_one_lev(i), zero) * ( one + total_mass / total_hole )
         447             :
         448             :        else ! fill holes completely
         449           0 :           hm_one_lev_filled(i) = max(hm_one_lev(i), zero) * ( one + total_hole / total_mass )
         450             :
         451             :        endif
         452             :
         453             :     enddo
         454             :
         455             :     ! Assertion checks (water substance conservation, non-negativity)
         456           0 :     if ( clubb_at_least_debug_level( 2 ) ) then
         457             :
         458           0 :        if ( abs(sum( hm_one_lev ) - sum(hm_one_lev_filled)) > &
         459             :             abs(sum( hm_one_lev ) + sum(hm_one_lev_filled)) * eps/2 ) then
         460           0 :           print *, "Warning: Hole filling was not conservative!"
         461             :        endif
         462             :
         463           0 :        if ( any( hm_one_lev_filled < zero ) ) then
         464           0 :           print *, "Warning: Hole filling failed! A hole could not be filled."
         465             :        endif
         466             :
         467             :     endif
         468             :
         469             :     return
         470             :
         471             :   end subroutine hole_filling_hm_one_lev
         472             :   !-----------------------------------------------------------------------
         473             :
         474             :   !-----------------------------------------------------------------------
         475           0 :   subroutine fill_holes_hydromet( nz, hydromet_dim, hydromet, & ! Intent(in)
         476           0 :                                   l_frozen_hm, l_mix_rat_hm, & ! Intent(in)
         477           0 :                                   hydromet_filled ) ! Intent(out)
         478             :
         479             :   ! Description:
         480             :   ! Fills holes between same-phase hydrometeors(i.e. for frozen hydrometeors).
         481             :   ! The hole filling conserves water substance between all same-phase (frozen or liquid)
         482             :   ! hydrometeors at each height level.
         483             :   !
         484             :   ! Attention: The hole filling for the liquid phase hydrometeors is not yet implemented
         485             :   !
         486             :   ! Attention: l_frozen_hm and l_mix_rat_hm need to be set up before this subroutine is called!
         487             :   !
         488             :   ! References:
         489             :   !
         490             :   ! None
         491             :   !-----------------------------------------------------------------------
         492             :
         493             :     use clubb_precision, only: &
         494             :         core_rknd
         495             :
         496             :     use constants_clubb, only: &
         497             :         zero
         498             :
         499             :     implicit none
         500             :
         501             :     ! Input Variables
         502             :     integer, intent(in) :: hydromet_dim, nz
         503             :
         504             :     real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(in) :: &
         505             :       hydromet
         506             :
         507             :     logical, dimension(hydromet_dim), intent(in) :: &
         508             :       l_frozen_hm, & ! if true, then the hydrometeor is frozen; otherwise liquid
         509             :       l_mix_rat_hm   ! if true, then the quantity is a hydrometeor mixing ratio
         510             :
         511             :     ! Output Variables
         512             :     real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(out) :: &
         513             :       hydromet_filled
         514             :
         515             :     ! Local Variables
         516             :     integer :: i,j ! Loop iterators
         517             :
         518             :     integer :: num_frozen_hm ! Number of frozen hydrometeor mixing ratios
         519             :
         520             :     real( kind = core_rknd ), dimension(:,:), allocatable :: &
         521           0 :       hydromet_frozen,       & ! Frozen hydrometeor mixing ratios
         522           0 :       hydromet_frozen_filled   ! Frozen hydrometeor mixing ratios after hole filling
         523             :
         524             :   !-----------------------------------------------------------------------
         525             :
         526             :     !----- Begin Code -----
         527             :
         528             :     ! Determine the number of frozen hydrometeor mixing ratios
         529           0 :     num_frozen_hm = 0
         530           0 :     do i=1,hydromet_dim
         531           0 :        if ( l_frozen_hm(i) .and. l_mix_rat_hm(i) ) then
         532           0 :           num_frozen_hm = num_frozen_hm + 1
         533             :        endif
         534             :     enddo
         535             :
         536             :     ! Allocation
         537           0 :     allocate( hydromet_frozen(nz,num_frozen_hm) )
         538           0 :     allocate( hydromet_frozen_filled(nz,num_frozen_hm) )
         539             :
         540             :     ! Determine frozen hydrometeor mixing ratios
         541           0 :     j = 1
         542           0 :     do i = 1,hydromet_dim
         543           0 :        if ( l_frozen_hm(i) .and. l_mix_rat_hm(i) ) then
         544           0 :           hydromet_frozen(:,j) = hydromet(:,i)
         545           0 :           j = j+1
         546             :        endif
         547             :     enddo
         548             :
         549             :     ! Fill holes for the frozen hydrometeors
         550           0 :     do i=1,nz
         551           0 :        if ( any( hydromet_frozen(i,:) < zero ) ) then
         552             :           call hole_filling_hm_one_lev( num_frozen_hm, hydromet_frozen(i,:), & ! Intent(in)
         553           0 :                                      hydromet_frozen_filled(i,:) ) ! Intent(out)
         554             :        else
         555           0 :           hydromet_frozen_filled(i,:) = hydromet_frozen(i,:)
         556             :        endif
         557             :     enddo
         558             :
         559             :     ! Setup the filled hydromet array
         560             :     j = 1
         561           0 :     do i=1, hydromet_dim
         562           0 :        if ( l_frozen_hm(i) .and. l_mix_rat_hm(i) ) then
         563           0 :           hydromet_filled(:,i) = hydromet_frozen_filled(:,j)
         564           0 :           j = j+1
         565             :        else
         566           0 :           hydromet_filled(:,i) = hydromet(:,i)
         567             :        endif
         568             :     enddo
         569             :
         570             :     !!! Here we could do the same hole filling for all the liquid phase hydrometeors
         571             :
         572           0 :     return
         573           0 :   end subroutine fill_holes_hydromet
         574             :   !-----------------------------------------------------------------------
         575             :
         576             :     !-----------------------------------------------------------------------
         577           0 :   subroutine fill_holes_wv( nz, dt, exner, hydromet_name, & ! Intent(in)
         578           0 :                             rvm_mc, thlm_mc, hydromet )! Intent(inout)
         579             :
         580             :   ! Description:
         581             :   ! Fills holes using the cloud water mixing ratio from the current height level.
         582             :   !
         583             :   ! References:
         584             :   !
         585             :   ! None
         586             :   !-----------------------------------------------------------------------
         587             :
         588             :     use clubb_precision, only: &
         589             :         core_rknd
         590             :
         591             :     use constants_clubb, only: &
         592             :         zero_threshold, &
         593             :         Lv, &
         594             :         Ls, &
         595             :         Cp
         596             :
         597             :     implicit none
         598             :
         599             :     ! Input Variables
         600             :     integer, intent(in) :: nz
         601             :
         602             :     real( kind = core_rknd ), intent(in) ::  &
         603             :       dt           ! Timestep         [s]
         604             :
         605             :     character(len=10), intent(in) :: hydromet_name
         606             :
         607             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         608             :       exner        ! Exner function                            [-]
         609             :
         610             :     ! Input/Output Variables
         611             :     real( kind = core_rknd ), dimension(nz), intent(inout) :: &
         612             :       hydromet, &  ! Hydrometeor array                         [units vary]
         613             :       rvm_mc, &
         614             :       thlm_mc
         615             :
         616             :     ! Local Variables
         617             :     integer :: k ! Loop iterator
         618             :
         619             :     real( kind = core_rknd ) :: rvm_clip_tndcy
         620             :   !-----------------------------------------------------------------------
         621             :
         622             :     !----- Begin Code -----
         623             :
         624           0 :     do k = 2, nz, 1
         625             :
         626           0 :        if ( hydromet(k) < zero_threshold ) then
         627             :
         628             :           ! Set rvm_clip_tndcy to the time tendency applied to vapor and removed
         629             :           ! from the hydrometeor.
         630           0 :           rvm_clip_tndcy = hydromet(k) / dt
         631             :
         632             :           ! Adjust the tendency rvm_mc accordingly
         633           0 :           rvm_mc(k) = rvm_mc(k) + rvm_clip_tndcy
         634             :
         635             :           ! Adjust the tendency of thlm_mc according to whether the
         636             :           ! effect is an evaporation or sublimation tendency.
         637           0 :           select case ( trim( hydromet_name ) )
         638             :           case( "rrm" )
         639           0 :              thlm_mc(k) = thlm_mc(k) - rvm_clip_tndcy * ( Lv / ( Cp*exner(k) ) )
         640             :           case( "rim", "rsm", "rgm" )
         641           0 :              thlm_mc(k) = thlm_mc(k) - rvm_clip_tndcy * ( Ls / ( Cp*exner(k) ) )
         642             :           case default
         643           0 :              error stop "Fatal error in microphys_driver"
         644             :           end select
         645             :
         646             :           ! Set the mixing ratio to 0
         647           0 :           hydromet(k) = zero_threshold
         648             :
         649             :        endif ! hydromet(k,i) < 0
         650             :
         651             :     enddo ! k = 2..gr%nz
         652             :
         653           0 :     return
         654             :   end subroutine fill_holes_wv
         655             :   !-----------------------------------------------------------------------
         656             :
         657             :   !-----------------------------------------------------------------------
         658           0 :   subroutine fill_holes_driver( gr, nz, dt, hydromet_dim, hm_metadata,   & ! Intent(in)
         659             :                                 l_fill_holes_hm,                            & ! Intent(in)
         660           0 :                                 rho_ds_zm, rho_ds_zt, exner,                & ! Intent(in)
         661             :                                 stats_metadata,                             & ! Intent(in)
         662             :                                 stats_zt,                                   & ! intent(inout)
         663           0 :                                 thlm_mc, rvm_mc, hydromet )                   ! Intent(inout)
         664             :
         665             :   ! Description:
         666             :   ! Fills holes between same-phase hydrometeors(i.e. for frozen hydrometeors).
         667             :   ! The hole filling conserves water substance between all same-phase (frozen or liquid)
         668             :   ! hydrometeors at each height level.
         669             :   !
         670             :   ! Attention: The hole filling for the liquid phase hydrometeors is not yet implemented
         671             :   !
         672             :   ! Attention: l_frozen_hm and l_mix_rat_hm need to be set up before this subroutine is called!
         673             :   !
         674             :   ! References:
         675             :   !
         676             :   ! None
         677             :   !-----------------------------------------------------------------------
         678             :
         679             :     use grid_class, only: &
         680             :         grid ! Type
         681             :
         682             :     use clubb_precision, only: &
         683             :         core_rknd   ! Variable(s)
         684             :
         685             :     use constants_clubb, only: &
         686             :         zero,            &
         687             :         zero_threshold,  &
         688             :         Lv,              &
         689             :         Ls,              &
         690             :         Cp,              &
         691             :         fstderr,         &
         692             :         num_hf_draw_points
         693             :
         694             :     use corr_varnce_module, only: &
         695             :         hm_metadata_type
         696             :
         697             :     use stats_type_utilities, only: &
         698             :         stat_begin_update, & ! Subroutines
         699             :         stat_end_update
         700             :
         701             :     use stats_variables, only: &
         702             :         stats_metadata_type
         703             :
         704             :     use error_code, only: &
         705             :         clubb_at_least_debug_level  ! Procedure
         706             :
         707             :     use stats_type, only: &
         708             :       stats ! Type
         709             :
         710             :     implicit none
         711             :
         712             :     !----------------------- Input Variables -----------------------
         713             :     type (grid), target, intent(in) :: &
         714             :       gr
         715             :
         716             :     integer, intent(in) :: &
         717             :       hydromet_dim, &
         718             :       nz
         719             :
         720             :     type (hm_metadata_type), intent(in) :: &
         721             :       hm_metadata
         722             :
         723             :     logical, intent(in) :: &
         724             :       l_fill_holes_hm
         725             :
         726             :     real( kind = core_rknd ), intent(in) ::  &
         727             :       dt           ! Timestep         [s]
         728             :
         729             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         730             :       rho_ds_zm, & ! Dry, static density on momentum levels   [kg/m^3]
         731             :       rho_ds_zt    ! Dry, static density on thermo. levels    [kg/m^3]
         732             :
         733             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         734             :       exner  ! Exner function                                       [-]
         735             :
         736             :     type (stats_metadata_type), intent(in) :: &
         737             :       stats_metadata
         738             :
         739             :     !----------------------- Input/Output Variables -----------------------
         740             :     type (stats), target, intent(inout) :: &
         741             :       stats_zt
         742             :
         743             :     real( kind = core_rknd ), dimension(nz, hydromet_dim), intent(inout) :: &
         744             :       hydromet    ! Mean of hydrometeor fields    [units vary]
         745             :
         746             :     real( kind = core_rknd ), dimension(nz), intent(inout) :: &
         747             :       rvm_mc,  & ! Microphysics contributions to vapor water           [kg/kg/s]
         748             :       thlm_mc    ! Microphysics contributions to liquid potential temp [K/s]
         749             :
         750             :     !----------------------- Local Variables -----------------------
         751             :     integer :: i, k ! Loop iterators
         752             :
         753             :     real( kind = core_rknd ), dimension(nz, hydromet_dim) :: &
         754           0 :       hydromet_filled,  & ! Frozen hydrometeor mixing ratios after hole filling
         755           0 :       hydromet_clipped    ! Clipped mean of hydrometeor fields    [units vary]
         756             :
         757             :     character( len = 10 ) :: hydromet_name
         758             :
         759             :     real( kind = core_rknd ) :: &
         760             :       max_velocity    ! Maximum sedimentation velocity                     [m/s]
         761             :
         762             :     integer :: ixrm_hf, ixrm_wvhf, ixrm_cl, &
         763             :                ixrm_bt, ixrm_mc
         764             :
         765             :     logical :: l_hole_fill = .true.
         766             :
         767             :     !----------------------- Begin Code -----------------------
         768             :
         769             :     ! Start stats output for the _hf variables (changes in the hydromet array
         770             :     ! due to fill_holes_hydromet and fill_holes_vertical)
         771           0 :     if ( stats_metadata%l_stats_samp ) then
         772             :
         773           0 :        do i = 1, hydromet_dim
         774             :
         775             :           ! Set up the stats indices for hydrometeor at index i
         776             :           call setup_stats_indices( i, stats_metadata, hydromet_dim,  & ! Intent(in)
         777             :                                     hm_metadata%hydromet_list,     & ! Intent(in)
         778             :                                     ixrm_bt, ixrm_hf, ixrm_wvhf,      & ! Intent(inout)
         779             :                                     ixrm_cl, ixrm_mc,                 & ! Intent(inout)
         780           0 :                                     max_velocity )                      ! Intent(inout)
         781             :
         782           0 :           call stat_begin_update( gr%nz, ixrm_hf, hydromet(:,i) / dt, & ! intent(in)
         783           0 :                                   stats_zt ) ! intent(inout)
         784             :
         785             :        enddo ! i = 1, hydromet_dim
         786             :
         787             :     endif ! stats_metadata%l_stats_samp
         788             :
         789             :     ! If we're dealing with negative hydrometeors, we first try to fill the
         790             :     ! holes proportionally from other same-phase hydrometeors at each height
         791             :     ! level.
         792           0 :     if ( any( hydromet < zero_threshold ) .and. l_fill_holes_hm ) then
         793             :
         794             :        call fill_holes_hydromet( nz, hydromet_dim, hydromet, & ! Intent(in)
         795             :                                  hm_metadata%l_frozen_hm, hm_metadata%l_mix_rat_hm, & ! Intent(in)
         796           0 :                                  hydromet_filled ) ! Intent(out)
         797             :
         798           0 :        hydromet = hydromet_filled
         799             :
         800             :     endif ! any( hydromet < zero ) .and. l_fill_holes_hm
         801             :
         802           0 :     hydromet_filled = zero
         803             :
         804           0 :     do i = 1, hydromet_dim
         805             :
         806             :       ! Set up the stats indices for hydrometeor at index i
         807             :       call setup_stats_indices( i, stats_metadata, hydromet_dim,  & ! Intent(in)
         808             :                                 hm_metadata%hydromet_list,     & ! Intent(in)
         809             :                                 ixrm_bt, ixrm_hf, ixrm_wvhf,      & ! Intent(inout)
         810             :                                 ixrm_cl, ixrm_mc,                 & ! Intent(inout)
         811           0 :                                 max_velocity )                      ! Intent(inout)
         812             :
         813           0 :       hydromet_name = hm_metadata%hydromet_list(i)
         814             :
         815             :       ! Print warning message if any hydrometeor species has a value < 0.
         816           0 :       if ( clubb_at_least_debug_level( 1 ) ) then
         817           0 :          if ( any( hydromet(:,i) < zero_threshold ) ) then
         818             :
         819           0 :             do k = 1, nz
         820           0 :                if ( hydromet(k,i) < zero_threshold ) then
         821           0 :                   write(fstderr,*) trim( hydromet_name ) //" < ", &
         822           0 :                                    zero_threshold, &
         823           0 :                                    " in fill_holes_driver at k= ", k
         824             :                endif ! hydromet(k,i) < 0
         825             :             enddo ! k = 1, nz
         826             :
         827             :          endif ! hydromet(:,i) < 0
         828             :       endif ! clubb_at_least_debug_level( 1 )
         829             :
         830             :
         831             :       ! Store the previous value of the hydrometeor for the effect of the
         832             :       ! hole-filling scheme.
         833             : !      if ( stats_metadata%l_stats_samp ) then
         834             : !         call stat_begin_update( ixrm_hf, hydromet(:,i) &
         835             : !                                          / dt, stats_zt )
         836             : !      endif
         837             :
         838             :       ! If we're dealing with a mixing ratio and hole filling is enabled,
         839             :       ! then we apply the hole filling algorithm
         840           0 :       if ( any( hydromet(:,i) < zero_threshold ) ) then
         841             :
         842           0 :          if ( hydromet_name(1:1) == "r" .and. l_hole_fill ) then
         843             :
         844             :             !$acc data copyin( gr, gr%dzt, rho_ds_zt ) &
         845             :             !$acc        copy( hydromet(:,i) )
         846             :
         847             :             ! Apply the hole filling algorithm
         848             :             ! upper_hf_level = nz since we are filling the zt levels
         849             :             call fill_holes_vertical( gr%nz, 1, zero_threshold, gr%nz, & ! In
         850             :                                       gr%dzt, rho_ds_zt,               & ! In
         851           0 :                                       hydromet(:,i) )                    ! InOut
         852             :
         853             :             !$acc end data
         854             :
         855             :          endif ! Variable is a mixing ratio and l_hole_fill is true
         856             :
         857             :       endif ! hydromet(:,i) < 0
         858             :
         859             :       ! Enter the new value of the hydrometeor for the effect of the
         860             :       ! hole-filling scheme.
         861           0 :       if ( stats_metadata%l_stats_samp ) then
         862           0 :          call stat_end_update( gr%nz, ixrm_hf, hydromet(:,i) / dt, & ! intent(in)
         863           0 :                                stats_zt ) ! intent(inout)
         864             :       endif
         865             :
         866             :       ! Store the previous value of the hydrometeor for the effect of the water
         867             :       ! vapor hole-filling scheme.
         868           0 :       if ( stats_metadata%l_stats_samp ) then
         869           0 :          call stat_begin_update( gr%nz, ixrm_wvhf, hydromet(:,i) / dt, & ! intent(in)
         870           0 :                                  stats_zt ) ! intent(inout)
         871             :       endif
         872             :
         873           0 :       if ( any( hydromet(:,i) < zero_threshold ) ) then
         874             :
         875           0 :          if ( hydromet_name(1:1) == "r" .and. l_hole_fill ) then
         876             :
         877             :             ! If the hole filling algorithm failed, then we attempt to fill
         878             :             ! the missing mass with water vapor mixing ratio.
         879             :             ! We noticed this is needed for ASEX A209, particularly if Latin
         880             :             ! hypercube sampling is enabled.  -dschanen 11 Nov 2010
         881             :             call fill_holes_wv( nz, dt, exner, hydromet_name, & ! Intent(in)
         882           0 :                                 rvm_mc, thlm_mc, hydromet(:,i) )   ! Intent(out)
         883             :
         884             :          endif ! Variable is a mixing ratio and l_hole_fill is true
         885             :
         886             :       endif ! hydromet(:,i) < 0
         887             :
         888             :       ! Enter the new value of the hydrometeor for the effect of the water vapor
         889             :       ! hole-filling scheme.
         890           0 :       if ( stats_metadata%l_stats_samp ) then
         891           0 :          call stat_end_update( gr%nz, ixrm_wvhf, hydromet(:,i) / dt, & ! intent(in)
         892           0 :                                stats_zt ) ! intent(inout)
         893             :       endif
         894             :
         895             :       ! Clipping for hydrometeor mixing ratios.
         896           0 :       if ( hm_metadata%l_mix_rat_hm(i) ) then
         897             :
         898             :          ! Store the previous value of the hydrometeor for the effect of
         899             :          ! clipping.
         900           0 :          if ( stats_metadata%l_stats_samp ) then
         901           0 :             call stat_begin_update( gr%nz, ixrm_cl, hydromet(:,i) / dt, & ! intent(in)
         902           0 :                                     stats_zt ) ! intent(inout)
         903             :          endif
         904             :
         905           0 :          if ( any( hydromet(:,i) < zero_threshold ) ) then
         906             :
         907             :             ! Clip any remaining negative values of precipitating hydrometeor
         908             :             ! mixing ratios to 0.
         909           0 :             where ( hydromet(:,i) < zero_threshold )
         910             :                hydromet(:,i) = zero_threshold
         911             :             end where
         912             :
         913             :          endif ! hydromet(:,i) < 0
         914             :
         915             :          ! Eliminate very small values of mean precipitating hydrometeor mixing
         916             :          ! ratios by setting them to 0.
         917           0 :          do k = 2, gr%nz, 1
         918             :
         919           0 :             if ( hydromet(k,i) <= hm_metadata%hydromet_tol(i) ) then
         920             :
         921             :                rvm_mc(k) &
         922             :                = rvm_mc(k) &
         923           0 :                  + ( hydromet(k,i) / dt )
         924             :
         925           0 :                if ( .not. hm_metadata%l_frozen_hm(i) ) then
         926             :
         927             :                   ! Rain water mixing ratio
         928             :
         929             :                   thlm_mc(k) &
         930             :                   = thlm_mc(k) &
         931             :                     - ( Lv / ( Cp * exner(k) ) ) &
         932           0 :                       * ( hydromet(k,i) / dt )
         933             :
         934             :                else ! Frozen hydrometeor mixing ratio
         935             :
         936             :                   thlm_mc(k) &
         937             :                   = thlm_mc(k) &
         938             :                     - ( Ls / ( Cp * exner(k) ) ) &
         939           0 :                       * ( hydromet(k,i) / dt )
         940             :
         941             :                endif ! l_frozen_hm(i)
         942             :
         943           0 :                hydromet(k,i) = zero
         944             :
         945             :             endif ! hydromet(k,i) <= hydromet_tol(i)
         946             :
         947             :          enddo ! k = 2, gr%nz, 1
         948             :
         949             :
         950             :          ! Enter the new value of the hydrometeor for the effect of clipping.
         951           0 :          if ( stats_metadata%l_stats_samp ) then
         952             :             call stat_end_update( gr%nz, ixrm_cl, hydromet(:,i) / dt, & ! intent(in)
         953           0 :                                   stats_zt ) ! intent(inout)
         954             :          endif
         955             :
         956             :       endif ! l_mix_rat_hm(i)
         957             :
         958             :     enddo ! i = 1, hydromet_dim, 1
         959             :
         960             :     ! Calculate clipping for hydrometeor concentrations.
         961             :     call clip_hydromet_conc_mvr( gr%nz, hydromet_dim, hm_metadata, & ! Intent(in)
         962             :                                  hydromet,                            & ! Intent(in)
         963           0 :                                  hydromet_clipped )                     ! Intent(out)
         964             :
         965             :     ! Clip hydrometeor concentrations and output stats.
         966           0 :     do i = 1, hydromet_dim
         967             :
         968           0 :        if ( .not. hm_metadata%l_mix_rat_hm(i) ) then
         969             :
         970           0 :           if ( stats_metadata%l_stats_samp ) then
         971             :
         972             :              ! Set up the stats indices for hydrometeor at index i
         973             :              call setup_stats_indices( i, stats_metadata, hydromet_dim,   & ! Intent(in)
         974             :                                        hm_metadata%hydromet_list,      & ! Intent(in)
         975             :                                        ixrm_bt, ixrm_hf, ixrm_wvhf,       & ! In/Out
         976             :                                        ixrm_cl, ixrm_mc,                  & ! In/Out
         977           0 :                                        max_velocity )                       ! In/Out
         978             :
         979             :              ! Store the previous value of the hydrometeor for the effect of
         980             :              ! clipping.
         981           0 :              call stat_begin_update( gr%nz, ixrm_cl, hydromet(:,i) / dt, & ! intent(in)
         982           0 :                                      stats_zt ) ! intent(inout)
         983             :
         984             :           endif ! stats_metadata%l_stats_samp
         985             :
         986             :           ! Apply clipping of hydrometeor concentrations.
         987           0 :           hydromet(:,i) = hydromet_clipped(:,i)
         988             :
         989             :           ! Enter the new value of the hydrometeor for the effect of clipping.
         990           0 :           if ( stats_metadata%l_stats_samp ) then
         991             :              call stat_end_update( gr%nz, ixrm_cl, hydromet(:,i) / dt, & ! intent(in)
         992           0 :                                    stats_zt ) ! intent(inout)
         993             :           endif
         994             :
         995             :        endif ! .not. l_mix_rat_hm(i)
         996             :
         997             :     enddo ! i = 1, hydromet_dim, 1
         998             :
         999             :
        1000           0 :     return
        1001             :
        1002             :   end subroutine fill_holes_driver
        1003             :
        1004             :   !=============================================================================
        1005           0 :   subroutine clip_hydromet_conc_mvr( nz, hydromet_dim, hm_metadata,  & ! Intent(in)
        1006           0 :                                      hydromet,                          & ! Intent(in)
        1007           0 :                                      hydromet_clipped )                   ! Intent(out)
        1008             :
        1009             :     ! Description:
        1010             :     ! Increases the value of a hydrometeor concentration when it is too small,
        1011             :     ! according to the hydrometeor mixing ratio (which remains unchanged) and
        1012             :     ! the maximum acceptable drop or particle mean volume radius for that
        1013             :     ! hydrometeor species.
        1014             :
        1015             :     ! References:
        1016             :     !-----------------------------------------------------------------------
        1017             :
        1018             :     use grid_class, only: &
        1019             :         grid ! Type
        1020             :
        1021             :     use corr_varnce_module, only: &
        1022             :         hm_metadata_type
        1023             :
        1024             :     use constants_clubb, only: &
        1025             :         pi,          & ! Variable(s)
        1026             :         four_thirds, &
        1027             :         one,         &
        1028             :         zero,        &
        1029             :         rho_lw,      &
        1030             :         rho_ice
        1031             :
        1032             :     use index_mapping, only: &
        1033             :         Nx2rx_hm_idx, & ! Procedure(s)
        1034             :         mvr_hm_max
        1035             :
        1036             :     use clubb_precision, only: &
        1037             :         core_rknd    ! Variable(s)
        1038             :
        1039             :     implicit none
        1040             :
        1041             :     ! Input Variables
        1042             :     integer, intent(in) :: &
        1043             :       nz
        1044             :
        1045             :     integer, intent(in) :: &
        1046             :       hydromet_dim    ! Number of hydrometeor fields
        1047             :
        1048             :     type (hm_metadata_type), intent(in) :: &
        1049             :           hm_metadata
        1050             :
        1051             :     real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(in) :: &
        1052             :       hydromet    ! Mean of hydrometeor fields    [units vary]
        1053             :
        1054             :     ! Output Variable
        1055             :     real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(out) :: &
        1056             :       hydromet_clipped    ! Clipped mean of hydrometeor fields    [units vary]
        1057             :
        1058             :     ! Local Variables
        1059             :     real( kind = core_rknd ) :: &
        1060             :       Nxm_min_coef    ! Coefficient for min mean value of a concentration [1/kg]
        1061             :
        1062             :     integer :: &
        1063             :       k,   & ! Vertical grid index
        1064             :       idx    ! Hydrometeor species index
        1065             :
        1066             :
        1067             :     ! Clipping for hydrometeor concentrations.
        1068           0 :     do idx = 1, hydromet_dim
        1069             :
        1070           0 :        if ( .not. hm_metadata%l_mix_rat_hm(idx) ) then
        1071             :
        1072             :           ! The field is a hydrometeor concentration.
        1073             :
        1074           0 :           if ( .not. hm_metadata%l_frozen_hm(idx) ) then
        1075             :
        1076             :              ! Clipping for mean rain drop concentration, <Nr>.
        1077             :              !
        1078             :              ! When mean rain water mixing ratio, <rr>, is found at a grid
        1079             :              ! level, mean rain drop concentration must be at least a minimum
        1080             :              ! value so that average rain drop mean volume radius stays within
        1081             :              ! an upper bound.  Otherwise, mean rain drop concentration is 0.
        1082             :              ! This can also be applied to values at a sample point, rather than
        1083             :              ! a grid-box mean.
        1084             :
        1085             :              ! The minimum mean rain drop concentration is given by:
        1086             :              !
        1087             :              ! <Nr> = <rr> / ( (4/3) * pi * rho_lw * mvr_rain_max^3 );
        1088             :              !
        1089             :              ! where mvr_rain_max is the maximum acceptable average rain drop
        1090             :              ! mean volume radius.
        1091             :
        1092             :              Nxm_min_coef &
        1093           0 :              = one / ( four_thirds * pi * rho_lw * mvr_hm_max(idx,hm_metadata)**3 )
        1094             :
        1095             :           else ! l_frozen_hm(idx)
        1096             :
        1097             :              ! Clipping for mean frozen hydrometeor concentration, <Nx>, where x
        1098             :              ! stands for any frozen hydrometeor species.
        1099             :              !
        1100             :              ! When mean frozen hydrometeor mixing ratio, <rx>, is found at a
        1101             :              ! grid level, mean frozen hydrometeor concentration must be at
        1102             :              ! least a minimum value so that average frozen hydrometeor mean
        1103             :              ! volume radius stays within an upper bound.  Otherwise, mean
        1104             :              ! frozen hydrometeor concentration is 0.  This can also be applied
        1105             :              ! to values at a sample point, rather than a grid-box mean.
        1106             :
        1107             :              ! The minimum mean frozen hydrometeor concentration is given by:
        1108             :              !
        1109             :              ! <Nx> = <rx> / ( (4/3) * pi * rho_ice * mvr_x_max^3 );
        1110             :              !
        1111             :              ! where mvr_x_max is the maximum acceptable average frozen
        1112             :              ! hydrometeor mean volume radius for frozen hydrometeor species, x.
        1113             :
        1114             :              Nxm_min_coef &
        1115           0 :              = one / ( four_thirds * pi * rho_ice * mvr_hm_max(idx,hm_metadata)**3 )
        1116             :
        1117             :           endif ! .not. l_frozen_hm(idx)
        1118             :
        1119             :           ! Loop over vertical levels and increase hydrometeor concentrations
        1120             :           ! when necessary.
        1121           0 :           do k = 2, nz, 1
        1122             :
        1123           0 :              if ( hydromet(k,Nx2rx_hm_idx(idx,hm_metadata)) > zero ) then
        1124             :
        1125             :                 ! Hydrometeor mixing ratio, <rx>, is found at the grid level.
        1126           0 :                 hydromet_clipped(k,idx) &
        1127             :                 = max( hydromet(k,idx), &
        1128           0 :                        Nxm_min_coef * hydromet(k,Nx2rx_hm_idx(idx,hm_metadata)) )
        1129             :
        1130             :              else ! <rx> = 0
        1131             :
        1132           0 :                 hydromet_clipped(k,idx) = zero
        1133             :
        1134             :              endif ! hydromet(k,Nx2rx_hm_idx(idx)) > 0
        1135             :
        1136             :           enddo ! k = 2, nz, 1
        1137             :
        1138             :           ! The lowest thermodynamic level is below the model's lower boundary.
        1139           0 :           hydromet_clipped(1,idx) = hydromet(1,idx)
        1140             :
        1141             :        else ! l_mix_rat_hm(idx)
        1142             :
        1143             :           ! The field is a hydrometeor mixing ratio.
        1144           0 :           hydromet_clipped(:,idx) = hydromet(:,idx)
        1145             :
        1146             :        endif ! .not. l_mix_rat_hm(idx)
        1147             :
        1148             :     enddo ! idx = 1, hydromet_dim, 1
        1149             :
        1150             :
        1151           0 :     return
        1152             :
        1153             :   end subroutine clip_hydromet_conc_mvr
        1154             :
        1155             :   !-----------------------------------------------------------------------
        1156           0 :   subroutine setup_stats_indices( ihm, stats_metadata, hydromet_dim,  & ! Intent(in)
        1157           0 :                                   hydromet_list,                      & ! Intent(in)
        1158             :                                   ixrm_bt, ixrm_hf, ixrm_wvhf,        & ! Intent(inout)
        1159             :                                   ixrm_cl, ixrm_mc,                   & ! Intent(inout)
        1160             :                                   max_velocity )                        ! Intent(inout)
        1161             :
        1162             :   ! Description:
        1163             :   !
        1164             :   ! Determines the stats output indices depending on the hydrometeor.
        1165             :
        1166             :   ! Attention: hydromet_list needs to be set up before this routine is called.
        1167             :   !
        1168             :   ! Bogus example
        1169             :   ! References:
        1170             :   !
        1171             :   ! None
        1172             :   !-----------------------------------------------------------------------
        1173             :
        1174             :     use clubb_precision, only: &
        1175             :         core_rknd
        1176             :
        1177             :     use constants_clubb, only: &
        1178             :         zero
        1179             :
        1180             :     use stats_variables, only: &
        1181             :         stats_metadata_type
        1182             :
        1183             :     implicit none
        1184             :
        1185             :     ! Input Variables
        1186             :     integer, intent(in) :: ihm, hydromet_dim
        1187             :
        1188             :     type (stats_metadata_type), intent(in) :: &
        1189             :       stats_metadata
        1190             :
        1191             :     character(len=10), dimension(hydromet_dim), intent(in) :: &
        1192             :       hydromet_list
        1193             :
        1194             :     ! Input/Output Variables
        1195             :     real( kind = core_rknd ), intent(inout) :: &
        1196             :       max_velocity ! Maximum sedimentation velocity [m/s]
        1197             :
        1198             :     integer, intent(inout) :: ixrm_hf, ixrm_wvhf, ixrm_cl, &
        1199             :                               ixrm_bt, ixrm_mc
        1200             :
        1201             :   !-----------------------------------------------------------------------
        1202             :
        1203             :     !----- Begin Code -----
        1204             :
        1205             :     ! Initializing max_velocity in order to avoid a compiler warning.
        1206             :     ! Regardless of the case, it will be reset in the 'select case'
        1207             :     ! statement immediately below.
        1208           0 :     max_velocity = zero
        1209             :
        1210           0 :     select case ( trim( hydromet_list(ihm) ) )
        1211             :       case ( "rrm" )
        1212           0 :         ixrm_bt   = stats_metadata%irrm_bt
        1213           0 :         ixrm_hf   = stats_metadata%irrm_hf
        1214           0 :         ixrm_wvhf = stats_metadata%irrm_wvhf
        1215           0 :         ixrm_cl   = stats_metadata%irrm_cl
        1216           0 :         ixrm_mc   = stats_metadata%irrm_mc
        1217             :
        1218           0 :         max_velocity = -9.1_core_rknd ! m/s
        1219             :
        1220             :       case ( "rim" )
        1221           0 :         ixrm_bt   = stats_metadata%irim_bt
        1222           0 :         ixrm_hf   = stats_metadata%irim_hf
        1223           0 :         ixrm_wvhf = stats_metadata%irim_wvhf
        1224           0 :         ixrm_cl   = stats_metadata%irim_cl
        1225           0 :         ixrm_mc   = stats_metadata%irim_mc
        1226             :
        1227           0 :         max_velocity = -1.2_core_rknd ! m/s
        1228             :
        1229             :       case ( "rsm" )
        1230           0 :         ixrm_bt   = stats_metadata%irsm_bt
        1231           0 :         ixrm_hf   = stats_metadata%irsm_hf
        1232           0 :         ixrm_wvhf = stats_metadata%irsm_wvhf
        1233           0 :         ixrm_cl   = stats_metadata%irsm_cl
        1234           0 :         ixrm_mc   = stats_metadata%irsm_mc
        1235             :
        1236             :         ! Morrison limit
        1237             : !         max_velocity = -1.2_core_rknd ! m/s
        1238             :         ! Made up limit.  The literature suggests that it is quite possible
        1239             :         ! that snow flake might achieve a terminal velocity of 2 m/s, and this
        1240             :         ! happens in the COAMPS microphysics -dschanen 29 Sept 2009
        1241           0 :         max_velocity = -2.0_core_rknd ! m/s
        1242             :
        1243             :       case ( "rgm" )
        1244           0 :         ixrm_bt   = stats_metadata%irgm_bt
        1245           0 :         ixrm_hf   = stats_metadata%irgm_hf
        1246           0 :         ixrm_wvhf = stats_metadata%irgm_wvhf
        1247           0 :         ixrm_cl   = stats_metadata%irgm_cl
        1248           0 :         ixrm_mc   = stats_metadata%irgm_mc
        1249             :
        1250           0 :         max_velocity = -20._core_rknd ! m/s
        1251             :
        1252             :       case ( "Nrm" )
        1253           0 :         ixrm_bt   = stats_metadata%iNrm_bt
        1254           0 :         ixrm_hf   = 0
        1255           0 :         ixrm_wvhf = 0
        1256           0 :         ixrm_cl   = stats_metadata%iNrm_cl
        1257           0 :         ixrm_mc   = stats_metadata%iNrm_mc
        1258             :
        1259           0 :         max_velocity = -9.1_core_rknd ! m/s
        1260             :
        1261             :       case ( "Nim" )
        1262           0 :         ixrm_bt   = stats_metadata%iNim_bt
        1263           0 :         ixrm_hf   = 0
        1264           0 :         ixrm_wvhf = 0
        1265           0 :         ixrm_cl   = stats_metadata%iNim_cl
        1266           0 :         ixrm_mc   = stats_metadata%iNim_mc
        1267             :
        1268           0 :         max_velocity = -1.2_core_rknd ! m/s
        1269             :
        1270             :       case ( "Nsm" )
        1271           0 :         ixrm_bt   = stats_metadata%iNsm_bt
        1272           0 :         ixrm_hf   = 0
        1273           0 :         ixrm_wvhf = 0
        1274           0 :         ixrm_cl   = stats_metadata%iNsm_cl
        1275           0 :         ixrm_mc   = stats_metadata%iNsm_mc
        1276             :
        1277             :         ! Morrison limit
        1278             : !         max_velocity = -1.2_core_rknd ! m/s
        1279             :         ! Made up limit.  The literature suggests that it is quite possible
        1280             :         ! that snow flake might achieve a terminal velocity of 2 m/s, and this
        1281             :         ! happens in the COAMPS microphysics -dschanen 29 Sept 2009
        1282           0 :         max_velocity = -2.0_core_rknd ! m/s
        1283             :
        1284             :       case ( "Ngm" )
        1285           0 :         ixrm_bt   = stats_metadata%iNgm_bt
        1286           0 :         ixrm_hf   = 0
        1287           0 :         ixrm_wvhf = 0
        1288           0 :         ixrm_cl   = stats_metadata%iNgm_cl
        1289           0 :         ixrm_mc   = stats_metadata%iNgm_mc
        1290             :
        1291           0 :         max_velocity = -20._core_rknd ! m/s
        1292             :
        1293             :       case ( "Ncm" )
        1294           0 :         ixrm_bt   = stats_metadata%iNcm_bt
        1295           0 :         ixrm_hf   = 0
        1296           0 :         ixrm_wvhf = 0
        1297           0 :         ixrm_cl   = stats_metadata%iNcm_cl
        1298           0 :         ixrm_mc   = stats_metadata%iNcm_mc
        1299             :
        1300             :         ! Use the rain water limit, since Morrison has no explicit limit on
        1301             :         ! cloud water.  Presumably these numbers are never large.
        1302             :         ! -dschanen 28 Sept 2009
        1303           0 :         max_velocity = -9.1_core_rknd ! m/s
        1304             :
        1305             :       case default
        1306           0 :         ixrm_bt   = 0
        1307           0 :         ixrm_hf   = 0
        1308           0 :         ixrm_wvhf = 0
        1309           0 :         ixrm_cl   = 0
        1310           0 :         ixrm_mc   = 0
        1311             :
        1312           0 :         max_velocity = -9.1_core_rknd ! m/s
        1313             :
        1314             :     end select
        1315             :
        1316             :
        1317           0 :     return
        1318             :
        1319             :   end subroutine setup_stats_indices
        1320             :   !-----------------------------------------------------------------------
        1321             :
        1322             : end module fill_holes

Generated by: LCOV version 1.14