function [PET] = penmont_PET(tas,Rnet,vp,ps,u)

% [PET] = penmont_PET(tas,Rnet,vp,ps)
%
% This function will use the Penman-Monteith method to calculate reference
% (potential) evapotranspiration from monthly model data.

% PET  = potential evapotranspiration (mm day-1)
%
% tas       = temperature at 2m, degrees C
% vp        = Vapor pressure, Pa
% Rnet      = surface net radiation, W/m2
% ps        = surface pressure, Pa
% u         = surface wind at 2m, m/2
%
% Written by Benjamin I. Cook, modified by Nathan Steiger
%
% Based on: Xu and Singh (2002), from Water Resources Management
%           Lawrence (2005), BAMS
%
%           FAO Document:
%           http://www.fao.org/docrep/X0490E/x0490e00.htm#Contents
%
% For Tetens, both above and below zero:
%       http://cires.colorado.edu/~voemel/vp.html

if isempty(u)
    u=1; % assumed uniform wind speed at 2m (m/s)
else
    %u(u<=0.5)=0.5;%Limit very low wind speeds for the calculations (m/s)
end

% Calculate the latent heat of vaporization (MJ kg-1)
lambda_lv=2.501-(2.361e-3).*tas;

% Convert Pressure to kPa
ps=ps./1000;

% Saturation vapor pressure (kPa)
es=0.611.*exp((17.27.*tas)./(tas+237.3));

% Actual vapor pressure (kPa)
%vp=0.611.*exp((17.27.*Td)./(Td+237.3));
vp=vp./1000;

% Dewpoint Temperature (degrees C)
%Td=((tas+273.15).*(1-(((tas+273.15).*log(rh./100))./(2.501e6./461.5))).^(-1))-273.15;

% Relative Humidity
%rh=vp./es;

% Slope of the vapor pressure curve (kPa C-1)
%delta_vpc=(4098.*es.*tas)./((tas+237.3).^2);
delta_vpc=(4098.*es)./((tas+237.3).^2);

% Psychometric constant (kPa C-1)
psych_const=0.00163.*(ps./lambda_lv);

% Net radiation
Rnet=Rnet*0.0864;  % convert W/m2 to MJ/m2/d

% Soil Heat Flux
% soil heat capacity*daily temp change*effective soil depth
% gflux=0.7.*tdiff.*0.15;
% at daily time steps, we can ignore ground heat flux
gflux=0;

%VPD=es-vp;

% Potential Evapotranspiration (mm)
PET=(0.408.*delta_vpc.*(Rnet-gflux)+psych_const.*(900./(tas+273.15)).*u.*(es-vp))./(delta_vpc+psych_const.*(1+0.34.*u));
% For negative (mean) temperatures in Celsius: PET=0.
PET(tas<0)=0;









