et__hargreaves_samani.F90

Go to the documentation of this file.

!> @file
!>  Contains a single module, et_hargreaves, which
!>  calculates potential evapotranspiration by means of the Hargreaves-Samani (1985) method.
 
 
!>  Calculates potential evapotranspiration by means of the
!>  Hargreaves-Samani (1985) method.
module et__hargreaves_samani
!!****h* SWB/et_hargreaves
! NAME
!   et_hargreaves.f95 - Evapotranspiration calculation using the
!   Hargreaves method.
!
! SYNOPSIS
!   This module calculates evapotranspiration using the Hargreaves
!   method.
!
! NOTES
!
!  Reference:
!
!   Allen, R.G., and others, 2006, FAO Irrigation and Drainage Paper No. 56,
!    "Crop Evapotranspiration (Guidelines for computing crop water
!    requirements)", Food and Agriculture Organization, Rome, Italy.
!
!!***
 
  use iso_c_binding, only : c_short, c_int, c_float, c_double
  use constants_and_conversions
  use fstring, only                   : ascharacter
  use logfiles, only                  : logs, log_all
  use meteorological_calculations
  use parameters, only                : params
  use solar_calculations
 
  implicit none
 
 
  ! ET parameters -- default values are from Hargreaves and Samani (1985)
  real (c_double) :: et_slope         ! = 0.0023     
  real (c_double) :: et_exponent      ! = 0.5
  real (c_double) :: et_constant      ! = 17.8
 
contains
 
subroutine et_hargreaves_initialize( ) !pConfig, sRecord )
  !! Configures the module, using the command line in 'sRecord'
  ! [ ARGUMENTS ]
!  type (T_MODEL_CONFIGURATION), pointer :: pConfig ! pointer to data structure that contains
                                                   ! model options, flags, and other settings
!  character (len=*),intent(inout) :: sRecord
 
  ! [ LOCALS ]
 
  real (kind=c_float), allocatable :: fet_slope(:)
  real (kind=c_float), allocatable :: fet_exponent(:)
  real (kind=c_float), allocatable :: fet_constant(:)
 
  call params%get_parameters( skey="Hargreaves_ET_slope", fvalues=fet_slope)
  call params%get_parameters( skey="Hargreaves_ET_exponent", fvalues=fet_exponent)
  call params%get_parameters( skey="Hargreaves_ET_constant", fvalues=fet_constant)
 
  if ( fet_slope(1) > ftinyval) then
    et_slope = real(fet_slope(1), c_double)
  else
    et_slope = 0.0023_c_double
  endif
 
  if ( fet_exponent(1) > ftinyval) then
    et_exponent = real(fet_exponent(1), c_double)
  else
    et_exponent = 0.5_c_double
  endif
 
  if ( fet_constant(1) > ftinyval) then
    et_constant = real(fet_constant(1), c_double)
  else
    et_constant = 17.8_c_double
  endif
 
  call logs%write("****************************************************************************")
  call logs%write("initializing Hargreaves ET parameters")
  call logs%write("****************************************************************************")
  call logs%write("  ET_SLOPE     = "//ascharacter(real(et_slope, c_float)))
  call logs%write("  ET_EXPONENT  = "//ascharacter(real(et_exponent, c_float)))
  call logs%write("  ET_CONSTANT  = "//ascharacter(real(et_constant, c_float)))
 
  !  write(UNIT=LU_LOG,FMT=*) "Configuring Hargreaves PET model"
 
!   if (pConfig%rSouthernLatitude <= rNO_DATA_NCDC &
!     .or. pConfig%rNorthernLatitude <= rNO_DATA_NCDC) then
 
!     call Chomp( sRecord,sOption )
!     read ( unit=sOption, fmt=*, iostat=iStat ) rValue
!     call Assert( iStat == 0, "Could not read the southerly latitude" )
!     pConfig%rSouthernLatitude = dpTWOPI * rValue / 360.0_c_float
 
!     call Chomp( sRecord,sOption )
!     read ( unit=sOption, fmt=*, iostat=iStat ) rValue
!     call Assert( iStat == 0, "Could not read the northerly latitude" )
!     pConfig%rNorthernLatitude = dpTWOPI * rValue / 360.0_c_float
 
!   else
 
!     call echolog("Southern and northern latitude values have been determined" &
!       //"~from project grid bounds and projection parameters. The values supplied" &
!       //"~with the Hargreaves PET option will be ignored...")
 
!   endif
 
end subroutine et_hargreaves_initialize
 
!------------------------------------------------------------------------------
 
impure elemental function et_hargreaves_calculate( iDayOfYear, iNumDaysInYear, fLatitude, fTMin, fTMax )  result(fReferenceET0)
  !! Computes the potential ET for each cell, based on TMIN and TMAX.
  !! Stores cell-by-cell PET values in the model grid.
 
  integer (c_int),intent(in) :: idayofyear
  integer (c_int),intent(in) :: inumdaysinyear
  real (c_float), intent(in) :: flatitude
  real (c_float), intent(in) :: ftmin
  real (c_float), intent(in) :: ftmax
  real (c_double)             :: freferenceet0
 
  ! [ LOCALS ]
  real (c_double) :: fdelta, fomega_s, fd_r, fra
  real (c_double) :: dlatitude_radians
 
  dlatitude_radians = flatitude * degrees_to_radians
 
  fd_r =relative_earth_sun_distance__d_r(idayofyear,inumdaysinyear)
  fdelta = solar_declination__delta(idayofyear, inumdaysinyear)
 
  fomega_s = sunrise_sunset_angle__omega_s(dlatitude_radians, fdelta)
 
  ! NOTE that the following equation returns extraterrestrial radiation in
  ! MJ / m**2 / day.  The Hargreaves equation requires extraterrestrial
  ! radiation to be expressed in units of mm / day.
  fra = extraterrestrial_radiation__ra(dlatitude_radians, fdelta, fomega_s, fd_r)
 
  freferenceet0 = et0_hargreaves( equivalent_evaporation(fra), ftmin, ftmax)
 
end function et_hargreaves_calculate
 
 
elemental function et0_hargreaves( rRa, rTMinF, rTMaxF )   result(rET_0)
 
  ! [ ARGUMENTS ]
  real (c_double),intent(in) :: rra
  real (c_float),intent(in)  :: rtminf
  real (c_float),intent(in)  :: rtmaxf
 
  ! [ RETURN VALUE ]
  real (c_double) :: ret_0
 
  ! [ LOCALS ]
  real (c_double) :: rtdelta
  real (c_double) :: rtavg
 
  rtavg = (rtminf + rtmaxf) / 2.0_c_double
 
  rtdelta = abs(f_to_k(real(rtmaxf, c_double)) - f_to_k(real(rtminf, c_double)))
 
  ret_0 = max(dzero, &
                mm_to_in( et_slope * rra * (f_to_c(rtavg) + et_constant) * (rtdelta**et_exponent) ) )
!                mm_to_in( 0.0023_c_float * rRa * (F_to_C(rTavg) + 17.8_c_float) * sqrt(rTDelta)) )
 
end function et0_hargreaves
 
 
end module et__hargreaves_samani