continuous_frozen_ground_index.F90

Go to the documentation of this file.

module continuous_frozen_ground_index
 
  use iso_c_binding, only          : c_short, c_int, c_float, c_double, c_long_long
  use constants_and_conversions
  use data_catalog, only           : dat
  use data_catalog_entry, only     : data_catalog_entry_t
 
  use exceptions
  implicit none
 
  private
 
  public update_continuous_frozen_ground_index, initialize_continuous_frozen_ground_index
 
contains
 
  subroutine initialize_continuous_frozen_ground_index( cfgi, cfgi_ll, cfgi_ul, active_cells )
 
    real (c_float), intent(inout)  :: cfgi(:)
    real (c_float), intent(inout)  :: cfgi_ll(:)
    real (c_float), intent(inout)  :: cfgi_ul(:)
 
    logical (c_bool), intent(in)   :: active_cells(:,:)
 
    ! [ LOCALS ]
    type (data_catalog_entry_t), pointer :: pinitial_cfgi
    type (data_catalog_entry_t), pointer :: pcfgi_lower_limit
    type (data_catalog_entry_t), pointer :: pcfgi_upper_limit
 
    ! locate the data structure associated with the gridded initial_cfgi
    pinitial_cfgi => dat%find("INITIAL_CONTINUOUS_FROZEN_GROUND_INDEX")
 
    if ( .not. associated( pinitial_cfgi ) ) then
      pinitial_cfgi => dat%find("INITIAL_FROZEN_GROUND_INDEX")
 
      if ( .not. associated( pinitial_cfgi ) ) then
 
        call warn(smessage="An INITIAL_CONTINUOUS_FROZEN_GROUND_INDEX grid (or constant) was not found.",    &
        shints="Check your control file to see that a valid INITIAL_CONTINUOUS_FROZEN_GROUND_INDEX grid or"  &
          //" constant is specified.", lfatal=false )
 
        cfgi = 0.0_c_float
 
      endif  
 
    else
 
      call pinitial_cfgi%getvalues()
 
      ! map the 2D array of INITIAL_CFGI values to the vector of active cells
      cfgi = pack( pinitial_cfgi%pGrdBase%rData, active_cells )
 
     if ( minval( cfgi ) < fzero &
        .or. maxval( cfgi ) > 300.0_c_float )  &
       call warn(smessage="One or more initial continuous frozen ground values outside of " &
         //"valid range (0 to 300)", lfatal=true )
 
    endif
 
 
 
    ! locate the data structure associated with the gridded CFGI_LOWER_LIMIT
    pcfgi_lower_limit => dat%find("CFGI_LOWER_LIMIT")
 
    ! figure 4-21 shows an alternate syntax
    if ( .not. associated( pcfgi_lower_limit ) ) then
      pcfgi_lower_limit => dat%find("LOWER_LIMIT_CFGI")
    endif  
 
    if ( .not. associated( pcfgi_lower_limit ) ) then
 
      call warn(smessage="No value supplied for CONTINUOUS_FROZEN_GROUND_INDEX_LOWER_LIMIT.",             &
      shints="Value set to the default of 56, which was appropriate for the Pacific Northwestern U.S., "  &
      //"but may be inappropriate elsewhere.", lfatal=false )
 
      cfgi_ll = 56.   ! units are degrees C-days
 
    else
 
      call pcfgi_lower_limit%getvalues()
 
      ! map the 2D array of INITIAL_CFGI values to the vector of active cells
      cfgi_ll = pack( pcfgi_lower_limit%pGrdBase%rData, active_cells )
 
     if ( minval( cfgi_ll ) < fzero                                             &
        .or. maxval( cfgi_ll ) > 300.0_c_float )                                &
       call warn(smessage="One or more CFGI lower limit values outside of "     &
         //"valid range (0 to 300)", lfatal=true )
 
    endif
 
 
    ! locate the data structure associated with the gridded CFGI_UPPER_LIMIT
    pcfgi_upper_limit => dat%find("CFGI_UPPER_LIMIT")
 
    ! figure 4-21 shows an alternate syntax; accomodate this
    if ( .not. associated( pcfgi_upper_limit ) ) then
      pcfgi_upper_limit => dat%find("UPPER_LIMIT_CFGI")
    endif  
 
    if ( .not. associated( pcfgi_upper_limit ) ) then
 
      call warn(smessage="No value supplied for CONTINUOUS_FROZEN_GROUND_INDEX_UPPER_LIMIT.",             &
      shints="Value set to the default of 83, which was appropriate for the Pacific Northwestern U.S., "  &
      //"but may be inappropriate elsewhere.", lfatal=false )
 
        cfgi_ul = 83. ! units are degree C-days
 
    else
 
      call pcfgi_upper_limit%getvalues()
 
      ! map the 2D array of INITIAL_CFGI values to the vector of active cells
      cfgi_ul = pack( pcfgi_upper_limit%pGrdBase%rData, active_cells )
 
     if ( minval( cfgi_ul ) < fzero                                            &
        .or. maxval( cfgi_ul ) > 300.0_c_float )                               &
       call warn(smessage="One or more CFGI upper limit values outside of "    &
         //"valid range (0 to 300)", lfatal=true )
 
    endif
 
    if (any(cfgi_ul - cfgi_ll <= 0.0_c_float)) then
      call warn(smessage="One or more CFGI upper limit values is less that its"     &
        //" corresponding CFGI lower limit values", lfatal=true)
    endif
 
 
  end subroutine initialize_continuous_frozen_ground_index
 
!--------------------------------------------------------------------------------------------------
 
  !> Update the continuous frozen ground index (CFGI) for a cell.
  !!
  !! Computes the continuous frozen ground index
  !! @param[inout] rCFGI   Continuous frozen ground index to be updated.
  !! @param[in]    rTAvg_F Mean daily air temperature, in \degF
  !!
  !! @note Implemented as per Molnau and Bissel (1983).
  !!
  !! @note Molnau, M. and Bissell, V.C., 1983, A continuous frozen ground index for
  !! flood forecasting: In Proceedings 51st Annual Meeting Western Snow Conference,
  !! 109-119, Canadian Water Resources Assoc. Cambridge, Ont.
  elemental subroutine update_continuous_frozen_ground_index( fCFGI, fTMax_F, fTMin_F, fSnowCover )
 
    ! [ ARGUMENTS ]
    real (c_float), intent(inout)       :: fcfgi
    real (c_float), intent(in)          :: ftmax_f
    real (c_float), intent(in)          :: ftmin_f
    real (c_float), intent(in)          :: fsnowcover
 
    ! [ LOCALS ]
    real (c_float), parameter    :: fdecay_coefficient_a                      = 0.97_c_float
    real (c_float), parameter    :: fsnow_reduction_coefficient_freezing      = 0.08_c_float
    real (c_float), parameter    :: fsnow_reduction_coefficient_thawing       = 0.5_c_float
    real (c_float), parameter    :: fcm_per_inch                              = 2.54_c_float
    real (c_float), parameter    :: freezing_point_deg_c                      = 0.0_c_float
 
    real (c_float) :: ftavg_c              ! temporary variable holding avg temp in C
    real (c_float) :: fsnowdepthcm         ! snow depth in centimeters
 
 
    ftavg_c = f_to_c( (ftmax_f + ftmin_f) / 2.0_c_float )
 
    ! assuming snow depth is 10 times the water content of the snow in inches
    fsnowdepthcm = fsnowcover * 10.0_c_float * fcm_per_inch
 
    associate( tavg => ftavg_c,                                              &
               a => fdecay_coefficient_a,                                    &
               cfgi => fcfgi,                                                &
               k_freeze => fsnow_reduction_coefficient_freezing,             &
               k_thaw => fsnow_reduction_coefficient_thawing )
 
      if( tavg > freezing_point_deg_c ) then
 
        cfgi = max( a * cfgi - tavg * exp( -0.4_c_float * k_thaw * fsnowdepthcm ), 0.0_c_float )
 
      else ! temperature is below freezing
 
        cfgi = max( a * cfgi - tavg * exp( -0.4_c_float * k_freeze * fsnowdepthcm ), 0.0_c_float )
 
      end if
 
    end associate
 
  end subroutine update_continuous_frozen_ground_index
 
end module continuous_frozen_ground_index