interception__gash.F90

Go to the documentation of this file.

module interception__gash
 
  use iso_c_binding, only : c_short, c_int, c_float, c_double
  use constants_and_conversions
  use data_catalog
  use data_catalog_entry
  use dictionary
  use exceptions
  use file_operations
  use parameters, only          : params
  use grid
  use netcdf4_support, only      : t_netcdf4_file
  use simulation_datetime
  use fstring
  use fstring_list
  implicit none
 
  private
 
  public :: canopy_storage_capacity_table_values,           &
            trunk_storage_capacity_table_values,            &
            stemflow_fraction_table_values,                 &
            evaporation_to_rainfall_ratio,                  &
            p_sat,                                          &
            gash_interception_storage_max_growing_season,   &
            gash_interception_storage_max_nongrowing_season
 
  public interception_gash_initialize
  public interception_gash_calculate
  public precipitation_at_saturation
 
  type (data_catalog_entry_t), pointer :: pevaporation_to_rainfall_ratio
 
  real (c_float), allocatable   :: evaporation_to_rainfall_ratio(:)
 
  real (c_float), allocatable   :: canopy_storage_capacity_table_values(:)
  real (c_float), allocatable   :: trunk_storage_capacity_table_values(:)
  real (c_float), allocatable   :: stemflow_fraction_table_values(:)
  real (c_float), allocatable   :: gash_interception_storage_max_growing_season(:)
  real (c_float), allocatable   :: gash_interception_storage_max_nongrowing_season(:)
 
  real (c_float), allocatable   :: p_sat(:)
 
  type (t_netcdf4_file), pointer     :: pncfile           ! pointer to OUTPUT NetCDF file
 
contains
 
  !> Initialize the Gash interception algorithm.
  !!
  !! Read in a canopy cover fraction and rainfall-to-evaporation ratio grids
  !!
  subroutine interception_gash_initialize( lActive, fCanopy_Cover_Fraction, iLandUseIndex )
 
    logical (c_bool), intent(in)    :: lactive(:,:)
    real (c_float), intent(in)      :: fcanopy_cover_fraction(:)
    integer (c_int), intent(in)     :: ilanduseindex(:)
 
    ! [ LOCALS ]
    integer (c_int)                 :: istat
    type (fstring_list_t)                 :: sllist
    integer (c_int)                 :: iindex
    integer (c_int)                 :: icount
    integer (c_int)                 :: inumrecs
    integer (c_int), allocatable    :: ilandusetablecodes(:)
    integer (c_int)                 :: inumberoflanduses
    logical (c_bool)                :: larelengthsequal
    real (c_float), allocatable     :: tempvals(:)
 
    icount = count( lactive )
 
    allocate( evaporation_to_rainfall_ratio( icount ), stat=istat )
    call assert( istat==0, "Problem allocating memory.", __file__, __line__ )
 
    ! locate the data structure associated with the gridded evaporation to rainfall ratio
    pevaporation_to_rainfall_ratio => dat%find("EVAPORATION_TO_RAINFALL_RATIO")
    if ( .not. associated( pevaporation_to_rainfall_ratio ) ) &
        call die("A EVAPORATION_TO_RAINFALL_RATIO grid must be supplied in order to"      &
                 //" make use of this option.", __file__, __line__)
 
 
    call pevaporation_to_rainfall_ratio%getvalues()
 
    call grid_writearcgrid("Evaporation_to_Rainfall_Ratio__echo.asc", &
      pevaporation_to_rainfall_ratio%pGrdBase )
 
    evaporation_to_rainfall_ratio = pack( pevaporation_to_rainfall_ratio%pGrdBase%rData, lactive )
 
    !! now grab the table values needed for this module
 
    ! create list of possible table headings to look for...
    call sllist%append( "LU_Code" )
    call sllist%append( "Landuse_Lookup_Code" )
 
    !> Determine how many landuse codes are present
    call params%get_parameters( slkeys=sllist, ivalues=ilandusetablecodes )
    inumberoflanduses = count( ilandusetablecodes >= 0 )
 
    call sllist%clear()
    call sllist%append("Canopy_Capacity")
    call sllist%append("Canopy_Storage_Capacity")
    call params%get_parameters( slkeys=sllist, fvalues=canopy_storage_capacity_table_values )
 
    call sllist%clear()
    call sllist%append("Trunk_Capacity")
    call sllist%append("Trunk_Storage_Capacity")
    call params%get_parameters( slkeys=sllist, fvalues=trunk_storage_capacity_table_values )
 
    call sllist%clear()
    call sllist%append("Stemflow_Fraction")
    call sllist%append("Stemflow_Frac")
    call params%get_parameters( slkeys=sllist, fvalues=stemflow_fraction_table_values )
 
    inumrecs = ubound(canopy_storage_capacity_table_values,1)
    larelengthsequal = ( inumrecs == inumberoflanduses )
 
    !> retrieve interception storage max (GROWING SEASON)
    call sllist%clear()
    call sllist%append("interception_storage_max_growing")
    call sllist%append("interception_storage_max_growing_season")
    call sllist%append("interception_storage_maximum_growing")
 
    call params%get_parameters( slkeys=sllist,                                           &
                                fvalues=gash_interception_storage_max_growing_season, &
                                lfatal=false )
 
    inumrecs = ubound(gash_interception_storage_max_growing_season,1)
    larelengthsequal = ( inumrecs == inumberoflanduses )
 
    if ( .not. larelengthsequal ) then
      call warn( smessage="The number of landuses does not match the number of interception storage "      &
                          //"maximum values for the growing season"                                        &
                          //"('interception_storage_max_growing').",                                       &
                 shints="A default value of 0.1 inches was assigned for the maximum interception storage", &
                 smodule=__file__, iline=__line__, lfatal=false )
      allocate(tempvals(inumberoflanduses), stat=istat)
      tempvals = 0.1
      call move_alloc(tempvals, gash_interception_storage_max_growing_season)
           
    endif
 
    !> retrieve interception storage max (NONGROWING SEASON)
    call sllist%clear()
    call sllist%append("interception_storage_max_nongrowing")
    call sllist%append("interception_storage_max_nongrowing_season")
    call sllist%append("interception_storage_maximum_nongrowing")
 
    call params%get_parameters( slkeys=sllist,                                           &
                                fvalues=gash_interception_storage_max_nongrowing_season, &
                                lfatal=false )
 
    inumrecs = ubound(gash_interception_storage_max_nongrowing_season,1)
    larelengthsequal = ( inumrecs == inumberoflanduses )
 
    if ( .not. larelengthsequal ) then
      call warn( smessage="The number of landuses does not match the number of interception storage "         &
                          //"maximum values for the nongrowing season"                                        &
                          //"('interception_storage_max_nongrowing').",                                       &
                 shints="A default value of 0.1 inches was assigned for the maximum interception storage",    &
                 smodule=__file__, iline=__line__, lfatal=false )
      allocate(tempvals(inumberoflanduses), stat=istat)
      tempvals = 0.1
      call move_alloc(tempvals, gash_interception_storage_max_nongrowing_season)
    endif
 
    !> @TODO add more guard code here to QA incoming data
 
    if ( .not. larelengthsequal )                                                       &
      call warn( smessage="The number of canopy storage capacity values ("              &
        //ascharacter( inumrecs )//") does not match the number of landuse values ("    &
        //ascharacter( inumberoflanduses )//").",                                       &
        smodule=__file__, iline=__line__, lfatal=.true._c_bool )
 
 
    inumrecs = ubound(trunk_storage_capacity_table_values,1)
    larelengthsequal = ( inumrecs == inumberoflanduses )
 
    if ( .not. larelengthsequal )                                                        &
      call warn( smessage="The number of trunk storage capacity values ("                &
        //ascharacter( inumrecs )//") does not match the number of landuse values ("     &
        //ascharacter( inumberoflanduses )//").",                                        &
        smodule=__file__, iline=__line__, lfatal=.true._c_bool )
 
 
    inumrecs = ubound(stemflow_fraction_table_values,1)
    larelengthsequal = ( inumrecs == inumberoflanduses )
 
    if ( .not. larelengthsequal )                                                        &
      call warn( smessage="The number of stemflow fraction values ("                &
        //ascharacter( inumrecs )//") does not match the number of landuse values ("     &
        //ascharacter( inumberoflanduses )//").",                                        &
        smodule=__file__, iline=__line__, lfatal=.true._c_bool )
 
 
    allocate( p_sat( count( lactive ) ), stat=istat )
    call assert( istat==0, "Problem allocating memory for P_SAT.", __file__, __line__ )
 
    p_sat = precipitation_at_saturation( evaporation_to_rainfall_ratio,    &
      canopy_storage_capacity_table_values( ilanduseindex ), &
      fcanopy_cover_fraction )
 
  end subroutine interception_gash_initialize
 
!--------------------------------------------------------------------------------------------------
 
  elemental function precipitation_at_saturation( E_div_P, canopy_storage_capacity, canopy_cover_fraction )  &
                                                                               result( psat )
 
    real (c_float), intent(in)    :: e_div_p
    real (c_float), intent(in)    :: canopy_storage_capacity
    real (c_float), intent(in)    :: canopy_cover_fraction
    real (c_float)                :: psat
 
    ! [ LOCALS ]
    real (c_float)                :: p_div_e
 
    if ( canopy_cover_fraction > 0.0_c_float .and. canopy_storage_capacity > 0.0_c_float ) then
 
       psat = - ( canopy_storage_capacity / ( canopy_cover_fraction * e_div_p ) )      &
                     * log( 1.0_c_float - e_div_p )
    else
 
      psat = 0.0_c_float
 
    endif
 
!                !! calc Precip needed to saturate canopy
!                 Psat=-( cancap( ilu(ip) ) / ( canfrac( ip )*cerf( ip ) ) ) &
!                       * log( ( 1 - cerf( ip ) ) / ( 1 - ( 1 - ceint2 ) * cerf( ip ) ) )
 
 
  end function precipitation_at_saturation
 
!--------------------------------------------------------------------------------------------------
 
  elemental subroutine interception_gash_calculate( fRainfall, fFog, fCanopy_Cover_Fraction,       &
                                                    fTrunk_Storage_Capacity, fStemflow_Fraction,   &
                                                    fEvaporation_to_Rainfall_Ratio,                &
                                                    fPrecipitation_at_Saturation, fInterception )
 
    real (c_float), intent(in)        :: frainfall
    real (c_float), intent(in)        :: ffog
    real (c_float), intent(in)        :: fcanopy_cover_fraction
    real (c_float), intent(in)        :: ftrunk_storage_capacity
    real (c_float), intent(in)        :: fstemflow_fraction
    real (c_float), intent(in)        :: fevaporation_to_rainfall_ratio
    real (c_float), intent(in)        :: fprecipitation_at_saturation
    real (c_float), intent(inout)     :: finterception
 
    ! [ LOCALS ]
    real (c_float)      :: frainfall_plus_fog
 
    frainfall_plus_fog = frainfall + ffog
 
    if ( fprecipitation_at_saturation .approxequal. 0.0_c_float ) then
 
      finterception = 0.0_c_float
 
     ! if(drf+dfog.lt.Psat)then
    else if ( frainfall_plus_fog < fprecipitation_at_saturation ) then
 
      !  dcanint=canfrac(ip)*(drf+dfog)
      finterception = fcanopy_cover_fraction * frainfall_plus_fog
 
!       elseif(drf+dfog.gt.tcap(ilu(ip))/tfrac(ip))then
    else if ( frainfall_plus_fog > ftrunk_storage_capacity                                 &
                                   / (fstemflow_fraction + 1.0e-6) ) then
 
      finterception =  fcanopy_cover_fraction * fprecipitation_at_saturation               &
                      + fcanopy_cover_fraction * fevaporation_to_rainfall_ratio            &
                         * ( frainfall_plus_fog - fprecipitation_at_saturation )           &
                      + ftrunk_storage_capacity
 
 
      ! dcanint = canfrac(ip) * Psat
      !           + canfrac(ip) * cerf(ip) * ( drf + dfog - Psat )
      !           + tcap( ilu(ip) )
 
    else
 
      finterception =  fcanopy_cover_fraction * fprecipitation_at_saturation               &
                      + fcanopy_cover_fraction * fevaporation_to_rainfall_ratio            &
                         * ( frainfall_plus_fog - fprecipitation_at_saturation )           &
                      + fstemflow_fraction * frainfall_plus_fog
 
      ! dcanint = canfrac(ip) * Psat
      !           + canfrac(ip) * cerf(ip) *( drf + dfog - Psat )
      !           + tfrac(ip) * ( drf + dfog )
 
    end if
 
    finterception = min( finterception, frainfall_plus_fog )
 
    !! HWB code:
!                !! calc Precip needed to saturate canopy
!                 Psat=-( cancap( ilu(ip) ) / ( canfrac( ip )*cerf( ip ) ) ) &
!                                 * log( ( 1 - cerf( ip ) ) / ( 1 - ( 1 - ceint2 ) * cerf( ip ) ) )
!
!                 if ( drf + dfog .lt. Psat )then
!                    dcanint = canfrac(ip) * ( drf + dfog )
!                 elseif ( drf + dfog .gt. tcap( ilu(ip) ) / tfrac(ip) ) then
!                    dcanint = canfrac(ip) * Psat + canfrac(ip) * cerf(ip) * ( drf + dfog - Psat ) + tcap( ilu(ip) )
!                 else
!                    dcanint = canfrac(ip) * Psat + canfrac(ip) * cerf(ip) *( drf + dfog - Psat ) + tfrac(ip) * ( drf + dfog )
!                 endif
!
!                 dpnet = drf + dfog - dcanint
 
  end subroutine interception_gash_calculate
 
end module interception__gash