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 17:57:11 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   178701120 :   subroutine fill_holes_vertical( nz, ngrdcol, threshold, upper_hf_level, &
      22   178701120 :                                   dz, rho_ds, &
      23   178701120 :                                   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   357402240 :       rho_ds_dz,            & ! rho_ds * dz
     101   357402240 :       invrs_denom_integral, & ! Inverse of the integral in the denominator (see description)
     102   357402240 :       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   357402240 :       denom_integral_global,  & ! Integral in the denominator for global filling
     113   357402240 :       numer_integral_global,  & ! Integral in the numerator for global filling
     114   357402240 :       field_avg_global,       & ! Vertical average of field [Units of field]
     115   357402240 :       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   178701120 :     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 15368296320 :     do k = 1, nz
     131 >25380*10^7 :       do i = 1, ngrdcol
     132 >25363*10^7 :         if ( field(i,k) < threshold ) then
     133  8966831621 :           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   178701120 :     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  2380439580 :     do k = 1, nz
     149 39324789090 :       do i = 1, ngrdcol
     150 39297109560 :         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   462318936 :     do i = 1, ngrdcol
     161 34805633527 :       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 34343314591 :         rho_k_sum = 0.0_core_rknd
     169             : 
     170 >20605*10^7 :         do j = k - num_hf_draw_points, k + num_hf_draw_points
     171 >20605*10^7 :           rho_k_sum = rho_k_sum + rho_ds_dz(i,j)
     172             :         end do
     173             : 
     174 34777953997 :         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   462318936 :     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 34805633527 :       do k = 2+num_hf_draw_points, upper_hf_level-num_hf_draw_points
     190             : 
     191 34343314591 :         k_start = k - num_hf_draw_points
     192 34343314591 :         k_end   = k + num_hf_draw_points
     193             : 
     194 >16718*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 55769172702 :                       * invrs_denom_integral(i,k)
     200             : 
     201             :           ! Clip small or negative values from field.
     202  9294862117 :           if ( field_avg >= threshold ) then
     203             :             ! We know we can fill in holes completely
     204  7520717778 :             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 48248454924 :             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 55769172702 :                               * invrs_denom_integral(i,k)
     216             : 
     217             :           ! Avoid divide by zero issues by doing nothing if field_clipped_avg ~= threshold
     218  9294862117 :           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  9294862046 :                             / ( field_clipped_avg - threshold )
     223             : 
     224             :             ! Calculate normalized, filled field
     225             :             field(i,k_start:k_end) = threshold &
     226 55769172276 :                         + 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  2380439580 :     do k = 1, nz
     240 39324789090 :       do i = 1, ngrdcol
     241 39297109560 :         if ( field(i,k) < threshold ) then
     242  7580956695 :           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    27679530 :     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   187199796 :     do i = 1, ngrdcol
     264   175993909 :       numer_integral_global(i) = 0.0_core_rknd
     265   187199796 :       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   187199796 :     do i = 1, ngrdcol
     271 14794806837 :       do k = 2, upper_hf_level
     272 14607607041 :         numer_integral_global(i) = numer_integral_global(i) + rho_ds_dz(i,k) * field(i,k)
     273             : 
     274 14783600950 :         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   187199796 :     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   175993909 :       field_avg_global(i) = numer_integral_global(i) / denom_integral_global(i)
     286             : 
     287             :       ! Clip small or negative values from field.
     288   187199796 :       if ( field_avg_global(i) >= threshold ) then
     289             :         ! We know we can fill in holes completely
     290 14741288806 :         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    42312144 :         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   187199796 :     do i = 1, ngrdcol
     303   175993909 :       numer_integral_global(i) = 0.0_core_rknd
     304 14794806837 :       do k = 2, upper_hf_level
     305 14783600950 :         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   187199796 :     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  5908083612 :       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   145463245 :         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   145463245 :         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   145463245 :                                     / ( 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 12218912580 :                                                       - 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