runoff__curve_number.F90

Go to the documentation of this file.

module runoff__curve_number
 
  use iso_c_binding, only                    : c_int, c_float, c_double, c_bool
  use constants_and_conversions, only        : false, true
  use datetime
  use exceptions
  use simulation_datetime
  use fstring
  use fstring_list
  use parameters, only                    : params, params_dict
  implicit none
 
  private
 
  real (c_float), allocatable    :: cn_arciii(:,:)
  real (c_float), allocatable    :: cn_arcii(:,:)
  real (c_float), allocatable    :: cn_arci(:,:)
  real (c_float), allocatable    :: prev_5_days_rain(:,:)
  integer (c_int), allocatable   :: ilandusecodes(:)
  integer (c_int)                :: daycount
  integer (c_int), parameter     :: five_day_sum = 6
 
 
  public :: runoff_curve_number_initialize
  public :: runoff_curve_number_calculate
  public :: update_previous_5_day_rainfall
  public :: update_curve_number_fn
  public :: cn_arci, cn_arcii, cn_arciii
  public :: five_day_sum, prev_5_days_rain
 
  real (c_float), parameter :: amc_dry_growing = 1.40_c_float
  real (c_float), parameter :: amc_dry_dormant = 0.50_c_float
  real (c_float), parameter :: amc_wet_growing = 2.10_c_float
  real (c_float), parameter :: amc_wet_dormant = 1.10_c_float
 
  type (datetime_t) :: date_last_updated
 
contains
 
  subroutine runoff_curve_number_initialize( cell_is_active )
 
   use ieee_arithmetic, only : ieee_is_nan, ieee_is_finite
 
    logical (c_bool), intent(in)   :: cell_is_active(:,:)
 
    ! [ LOCALS ]
    integer (c_int)              :: inumberoflanduses
    integer (c_int)              :: inumberofsoilgroups
    integer (c_int), allocatable :: icurvenumberseqnums(:)
    type (fstring_list_t)              :: sllist
    type (fstring_list_t)              :: slcurvenumber
    integer (c_int)              :: istat
    integer (c_int)              :: isoilsindex
    character (len=:), allocatable    :: stext
    real (c_float), allocatable  :: cn(:)
 
    call sllist%append("LU_Code")
    call sllist%append("Landuse_Lookup_Code")
 
    !> Determine how many soil groups are present
 
    ! retrieve a string list of all keys associated with curve number (i.e. "CN_1", "CN_2", "CN_3", etc)
    slcurvenumber = params%grep_name( "CN", lfatal=true )
 
    ! Convert the string list to an vector of integers; this call strips off the "CN_" part of label
    icurvenumberseqnums = slcurvenumber%get_integer()
 
    ! count how many items are present in the vector; this should equal the number of soils groups
    inumberofsoilgroups = count( icurvenumberseqnums > 0 )
 
    !> Determine how many landuse codes are present
    call params%get_parameters( slkeys=sllist, ivalues=ilandusecodes )
    inumberoflanduses = count( ilandusecodes >= 0 )
 
    allocate( cn_arciii(inumberoflanduses, inumberofsoilgroups), stat=istat )
    call assert( istat == 0, "Failed to allocate memory for curve number table - AMC III", &
      __file__, __line__)
 
    allocate( cn_arcii(inumberoflanduses, inumberofsoilgroups), stat=istat )
    call assert( istat == 0, "Failed to allocate memory for curve number table - AMC II", &
      __file__, __line__)
 
    allocate( cn_arci(inumberoflanduses, inumberofsoilgroups), stat=istat )
    call assert( istat == 0, "Failed to allocate memory for curve number table - AMC I", &
      __file__, __line__)
 
    ! we should have the curve number table fully filled out following this block
    do isoilsindex = 1, inumberofsoilgroups
      stext = "CN_"//ascharacter(isoilsindex)
      call params%get_parameters( skey=stext, fvalues=cn )
      cn_arcii(:, isoilsindex) = cn
    enddo
 
    allocate( prev_5_days_rain( count(cell_is_active), 6 ), stat=istat )
    call assert( istat == 0, "Failed to allocate memory for curve number PREV_5_DAYS_RAIN table", &
      __file__, __line__)
 
    prev_5_days_rain = 0.0_c_float
 
    ! DAYCOUNT will vary from 1 to 5 and serve as the second index value for
    ! PREV_5_DAYS_RAIN
    daycount = 0
 
    ! calculate curve numbers for antecedent runof conditions I and III
    cn_arci = cn_ii_to_cn_i( cn_arcii )
    cn_arciii = cn_ii_to_cn_iii( cn_arcii )
 
    call date_last_updated%parsedate("01/01/1000")
 
  end subroutine runoff_curve_number_initialize
 
!--------------------------------------------------------------------------------------------------
 
  elemental function return_landuse_index_fn(iLanduseCode)      result(iLanduseIndex)
 
    integer (c_int), intent(in)         :: ilandusecode
    integer (c_int)                     :: ilanduseindex
 
    ilanduseindex = 0
 
    do while (ilanduseindex < ubound( ilandusecodes, 1) )
 
      ilanduseindex = ilanduseindex + 1
 
      if (ilandusecode == ilandusecodes(ilanduseindex) ) exit
 
    end do
 
  end function return_landuse_index_fn
 
!--------------------------------------------------------------------------------------------------
 
  elemental function prob_runoff_enhancement( fCFGI, fCFGI_UL, fCFGI_LL )    result(fPf)
 
    real (c_float), intent(in) :: fcfgi
    real (c_float), intent(in) :: fcfgi_ul
    real (c_float), intent(in) :: fcfgi_ll
 
    real (c_float)             :: fpf
 
    if(fcfgi <= fcfgi_ll ) then
      fpf = 0_c_float
    elseif(fcfgi >= fcfgi_ul) then
      fpf = 1.0_c_float
    else
      fpf = ( fcfgi - fcfgi_ll ) / ( fcfgi_ul - fcfgi_ll )
    end if
 
  end function prob_runoff_enhancement
 
!--------------------------------------------------------------------------------------------------
 
  subroutine update_previous_5_day_rainfall( infil, indx )
 
    real (c_float), intent(in)            :: infil
    integer (c_int), intent(in)           :: indx
 
    ! we only want to increment DAYCOUNT 1x per day!
    if ( .not. date_last_updated == sim_dt%curr ) then
 
      if ( daycount < 5 ) then
        daycount = daycount + 1
      else
        daycount = 1
      endif
 
      date_last_updated = sim_dt%curr
 
    endif
 
    prev_5_days_rain( indx, daycount ) = infil
    prev_5_days_rain( indx, five_day_sum ) = sum( prev_5_days_rain( indx, 1:5 ) )
 
  end subroutine
!--------------------------------------------------------------------------------------------------
 
!   elemental function update_curve_number_fn( iLanduseIndex, iSoilsIndex, fSoilStorage, &
!                           fSoilStorage_Max, fCFGI )  result( CN_adj )
 
!     integer (c_int), intent(in)  :: iLanduseIndex
!     integer (c_int), intent(in)  :: iSoilsIndex
!     real (c_float), intent(in)   :: fSoilStorage
!     real (c_float), intent(in)   :: fSoilStorage_Max
!     real (c_float), intent(in)   :: fCFGI
!     real (c_float)               :: CN_adj
 
!     ! [ LOCALS ]
!     real (c_float) :: Pf
!     real (c_float) :: CN_base
!     real (c_float) :: fFraction_FC   ! fraction of field capacity
!     real (c_float) :: frac
 
!     fFraction_FC = 0.0
!     if (fSoilStorage_Max > 0.0) fFraction_FC = min( 1.0_c_float, fSoilStorage / fSoilStorage_Max )
 
!     ! Correct the curve number...
!     !
!     ! current thinking on this: if the soil moisture is at average levels or *greater*
!     ! *and* the CFGI indicates frozen ground conditions, assume that runoff enhancement
!     ! is probable. If the soils froze under relatively dry conditions, assume that some pore
!     ! space is open in the frozen ground.
!     if( fCFGI > CFGI_LL .and. fFraction_FC > 0.5 ) then
 
!       Pf = prob_runoff_enhancement( fCFGI )
 
!       ! use probability of runoff enhancement to calculate a weighted
!       ! average of curve number under Type II vs Type III antecedent
!       ! runoff conditions
!       CN_adj = CN_ARCII(iLanduseIndex, iSoilsIndex) * ( 1.0_c_float - Pf ) &
!                 + CN_ARCIII(iLanduseIndex, iSoilsIndex) * Pf
 
!     elseif ( fFraction_FC > 0.5_c_float ) then   ! adjust upward toward AMC III
 
!       ! we want a number ranging from 0 to 1 indicating the relative split
!       ! between CN AMCII and CN AMCIII
!       frac = ( fFraction_FC - 0.5_c_float ) / 0.5_c_float
 
!       CN_adj = CN_ARCII(iLanduseIndex, iSoilsIndex) * ( 1.0_c_float - frac ) &
!                 + CN_ARCIII(iLanduseIndex, iSoilsIndex) * frac
 
!     elseif ( fFraction_FC < 0.5_c_float ) then   ! adjust downward toward AMC I
 
!       ! we want a number ranging from 0 to 1 indicating the relative split
!       ! between CN AMCII and CN AMCI
!       !
!       ! ex. fFraction_FC = 0.48; frac = 2.0 * 0.48 = 0.96
!       !     CN_adj = 0.96 ( CN_ARCII ) + ( 1.0 - 0.96 ) * CN_ARCIII
!       !
!       frac = fFraction_FC / 0.5_c_float
 
!       CN_adj = CN_ARCI(iLanduseIndex, iSoilsIndex) * ( 1.0_c_float - frac ) &
!                 + CN_ARCII(iLanduseIndex, iSoilsIndex) * frac
 
!     else
 
!       CN_adj = CN_ARCII(iLanduseIndex, iSoilsIndex)
 
!     endif
 
 
!     ! ensure that whatever modification have been made to the curve number
!     ! remain within reasonable bounds
!     CN_adj = MIN( CN_adj, 100.0_c_float )
!     CN_adj = MAX( CN_adj, 0.0_c_float )
 
!   end function update_curve_number_fn
 
!--------------------------------------------------------------------------------------------------
 
  elemental function update_curve_number_fn( iLanduseIndex, iSoilsIndex, cell_index,   &
                                             it_is_growing_season, fSoilStorage_Max,   &
                                             fCFGI, fCFGI_LL, fCFGI_UL )                 result( CN_adj )
 
  integer (c_int), intent(in)  :: ilanduseindex
  integer (c_int), intent(in)  :: isoilsindex
  integer (c_int), intent(in)  :: cell_index
  logical (c_bool), intent(in) :: it_is_growing_season
  real (c_float), intent(in)   :: fsoilstorage_max
  real (c_float), intent(in)   :: fcfgi
  real (c_float), intent(in)   :: fcfgi_ll
  real (c_float), intent(in)   :: fcfgi_ul
  real (c_float)               :: cn_adj
 
  ! [ LOCALS ]
  real (c_float) :: pf
  real (c_float) :: frac
  real (c_float) :: finflow
 
  finflow = prev_5_days_rain( cell_index, five_day_sum )
 
  associate( cn_i   => cn_arci( ilanduseindex, isoilsindex ),        &
              cn_ii  => cn_arcii( ilanduseindex, isoilsindex ),       &
              cn_iii => cn_arciii( ilanduseindex, isoilsindex ) )
 
    ! Correct the curve number...
 
    if( ( fcfgi > fcfgi_ll ) .and. ( fsoilstorage_max > 0.0_c_float ) ) then
 
       pf = prob_runoff_enhancement( fcfgi, fcfgi_ll, fcfgi_ul )
 
       ! use probability of runoff enhancement to calculate a weighted
       ! average of curve number under Type II vs Type III antecedent
       ! runoff conditions
       cn_adj = cn_ii * (1-pf) +  cn_iii * pf
 
    else if ( it_is_growing_season ) then
 
      if ( finflow < amc_dry_growing ) then
 
        cn_adj = cn_i
 
      else if ( finflow >= amc_dry_growing                                  &
          .and. finflow < amc_wet_growing ) then
 
         cn_adj = cn_ii
 
      else
 
        cn_adj = cn_iii
 
      end if
 
    else ! dormant (non-growing) season
 
      if ( finflow < amc_dry_dormant ) then
 
        cn_adj = cn_i
 
      else if ( finflow >= amc_dry_dormant                                  &
          .and. finflow < amc_wet_dormant ) then
 
        cn_adj = cn_ii
 
      else
 
        cn_adj = cn_iii
 
      end if
 
    end if
 
  end associate
 
  ! ensure that whatever modification have been made to the curve number
  ! remain within reasonable bounds
  cn_adj = min(cn_adj,100.0_c_float)
  cn_adj = max(cn_adj,30.0_c_float)
 
 
  end function update_curve_number_fn
 
!--------------------------------------------------------------------------------------------------
 
  elemental function cn_ii_to_cn_iii(CN_II)   result(CN_III)
 
    real (c_float), intent(in)  :: cn_ii
    real (c_float)              :: cn_iii
 
    cn_iii = cn_ii / ( 0.427_c_float + 0.00573_c_float * cn_ii )
 
    cn_iii = max( cn_iii, 30.0_c_float )
    cn_iii = min( cn_iii, 100.0_c_float )
 
  end function cn_ii_to_cn_iii
 
!--------------------------------------------------------------------------------------------------
 
  elemental function cn_ii_to_cn_i(CN_II)   result(CN_I)
 
    real (c_float), intent(in)  :: cn_ii
    real (c_float)              :: cn_i
 
    ! The following comes from page 192, eq. 3.145 of "SCS Curve Number Methodology"
    cn_i = cn_ii / (2.281_c_float - 0.01281_c_float * cn_ii )
 
    cn_i = max( cn_i, 30.0_c_float )
    cn_i = min( cn_i, 100.0_c_float )
 
  end function cn_ii_to_cn_i
 
!--------------------------------------------------------------------------------------------------
  !> Calculate the runoff by means of the curve number method.
  !!
  !! Runoff is calculated using a modification of the curve number method.
  !! Specifically, the initial abstraction is defined as 0.05 * SMax rather than the
  !! standard 0.20 * SMax of the original method. This redefinition of the initial abstraction
  !! term was found to be more appropriate for long-term continuous simulations than the
  !! standard 0.2 * SMax of the original.
  !! @param[in]  iLanduseIndex  Pre-configured index number corresponding to a line in the landuse lookup table.
  !! @param[in]  iSoilsGroup   The numerical index associated with the soils group of interest.
  !! @param[in]  fInflow   The sum of net rainfall, snowmelt, and runon from upslope cells (and possibly irrigation).
  !! @param[in]  fCFGI   The current value of the continuous frozen-ground index.
  !! @param[in]  lIsGrowingSeason  Logical value indicating whether dormant season or growing season
  !!             criteria values are to be used when calculating runoff.
  !!
  !! @retval   fRunoff  Daily runoff calculated by means of the SCS Curve Number method.
  !!
  !! @note Reference: Woodward, D. E., R. H. Hawkins, R. Jiang, A. Hjelmfeldt Jr, J. Van Mullem,
  !!       and Q. D. Quan. "Runoff Curve Number Method: Examination of the Initial Abstraction Ratio."
  !!       In Conference Proceeding Paper, World Water and Environmental Resources Congress, 2003.
  elemental subroutine runoff_curve_number_calculate(runoff,                              &
                                                     curve_num_adj,                       &
                                                     cell_index,                          &
                                                     landuse_index,                       &
                                                     soil_group,                          &
                                                     it_is_growing_season,                &
                                                     inflow,                              &
                                                     soil_storage_max,                    &
                                                     continuous_frozen_ground_index,      &
                                                     cfgi_lower_limit,                    &
                                                     cfgi_upper_limit )
 
    real (c_float), intent(inout)  :: runoff
    real (c_float), intent(inout)  :: curve_num_adj
    integer (c_int), intent(in)    :: cell_index
    integer (c_int), intent(in)    :: landuse_index
    integer (c_int), intent(in)    :: soil_group
    logical (c_bool), intent(in)   :: it_is_growing_season
    real (c_float), intent(in)     :: inflow
    real (c_float), intent(in)     :: soil_storage_max
    real (c_float), intent(in)     :: continuous_frozen_ground_index
    real (c_float), intent(in)     :: cfgi_lower_limit
    real (c_float), intent(in)     :: cfgi_upper_limit
 
    ! [ LOCALS ]
!    real (c_float) :: CN_05
    real (c_float) :: smax
    real (c_float) :: cn_adj
 
    curve_num_adj = update_curve_number_fn( landuse_index, soil_group,                        &
                                            cell_index,                                       &
                                            it_is_growing_season,                             &
                                            soil_storage_max,                                 &
                                            continuous_frozen_ground_index,                   &
                                            cfgi_lower_limit,                                 &
                                            cfgi_upper_limit )
 
    smax = ( 1000.0_c_float / curve_num_adj ) - 10.0_c_float
 
!     ! Equation 9, Hawkins and others, 2002
!     CN_05 = 100_c_float / &
!             ((1.879_c_float * ((100.0_c_float / CN_adj ) - 1.0_c_float )**1.15_c_float) + 1.0_c_float)
 
    ! Equation 8, Hawkins and others, 2002
    ! adjust Smax for alternate initial abstraction amount
    smax = 1.33_c_float * ( smax**1.15_c_float )
 
    ! now consider runoff if Ia ~ 0.05S
    if ( inflow > 0.05_c_float * smax ) then
      runoff = ( inflow - 0.05_c_float * smax )**2  / ( inflow + 0.95_c_float * smax )
    else
      runoff = 0.0_c_float
    end if
 
  end subroutine runoff_curve_number_calculate
 
end module runoff__curve_number