Warning! NONE OF THESE PARAMETERS ARE EVEN CLOSE TO REALISTIC! YOU NEED TO GO FIND THEM!
I need to learn Penman-Monteith if I want to be a wind person.
About what I knew prior to this silly exercise was "it calculates net ET!" how? where? to what extent?
There was only one solution... write a MatLab code to try it out. So, without understanding exactly what half of the parameters meant, I made a little PM toy. I want to put it here not because it's a great example of code, but maybe because I missed something, or maybe because someone else might want to not type all of PM into a code.
Source of the PM equation I used was the FAO website (
http://www.fao.org/docrep/x0490e/x0490e06.htm#penman monteith equation
I used their examples,too:
%------------------------------------------BEGIN!
%penman.m
% calculates net ET using PM equation
%parameters and surrogate functions
%----------------grass set--------------------
%d = 2/3; % zero plane displacement height
%zm = 2; % height of wind measurements in m
%zom = 0.123; % roughness length govern momentum transfer
%zh = 2; % height of humidity measures
%zoh = 0.1*zom; % roughness length hv if not given
%uz = 1; % wind speed at height
%h =0.12; % crop height
%k = 0.41;
%rl= 70; % single grass leaf stomatal resistance
% note that the FAO grass reference has h = 0.12, rl=70,and albedo =0.23.
%-----------------------------------------
%----------------- test set---------------------------
Rn = 2000; % net radiation (MJ/M2/d)
G = 500; % soil heat flux (MJ/M2/day)
rhoA = 2; % mean air density
cP = 4; % specific heat air
es = 10;
ea = 7; % es-ea vpd of air
%ra = 5; % aerdyn resistance
d = 2/3; % zero plane displacement height
zm = 2; % height of wind measurements in m
zom = 0.123; % roughness length govern momentum transfer
zh = 2; % height of humidity measures
%zoh = 3; % roughness length governing trans of heat/vap
zoh = 0.1*zom; % roughness length hv if not given
k = 0.41; % vonkarmans constant
uz = 1; % wind speed at height
h =0.12; % crop height
ra = (log((zm-d)*h/zom*h)*log((zh-d)*h/zoh*h))/k^2*uz;
gradient = 2; % slope of the saturation vp/t relate
gamma = 7; % psychrometric constant
% note: the bulk surface resistance describes the resistence of the vapour
% flow through the transpiring crop and soil evaporation. when veg does not
% cover the soil completely, the resistence factor should include the
% effects of evap from the soil surface.if the crop is not transpiring at
% its potential rate then teh resistence depends on the water status of the
% veg. we can do an approximation
%rs = 2; % bulk surface resistance, if applied in full
rl= 100; % bulk stomatal resist of illuminated leaf
%LAIa = 3; % active (lit) LAI
%----------------------LAI curve---------------------------------
t = 1:180; % length growing season for crop (age for tree)
LAImax = 7; % max feasible LAI
r = 0.05; % dampening coefficient for LAI curve
Lai(t) = LAImax./(1+exp(-r*t)); % generate LAI curve
% plot Lai if you want.
% plot(t, Lai(:));
% xlabel('time')
% ylabel('LAI');
% title ('Leaf Area Index with respect to time')
z = 100; % use z to specify your LAI from curve
LAIa = 0.5*(Lai(z));
% alternate route for grasses
% LAIa =12*h;
% rs approximation
rs = rl/LAIa;
%-------------------------end test set--------------------------
%penman equation for lambdaET
lambdaET= (gradient*(Rn-G)+rhoA*cP*(es-ea)/ra)/(gradient+gamma*(rs/ra));
% ------------------further useful substitutions--------------
% from the grass example I am working this from
% ra = 208/u; % aerdyn resistance
% rs = 70; % bulk surface resistance
% u = rs/(0.34*ra);
R = 0.287; % specific gas constant Kj/kg/K
% Tkelv = 1.01(T+273); % temp in kelv
% rhoA = P/(Tkelv*R); % mean air density
epsi =0.622; % molecular weight of watervap
P = 1000; % atmospheric pressure
lambda = 2.45; % latent heat of vaporization
cp = (gamma*epsi*lambda)/P;
% we can approximate Cp this way
%---------------------penman ET0------------------------------
delta = 5; % slope vapor press curve kPa/C
Rn = 1000; % net rads at surve MJ/m2/d
% note that if mean is used daily estimate of ET will be underest.
T = 5; % mean daily temp at 2m h
u = 0.5; % wind speed at 2 m height
ET0= (0.408*delta*(Rn-G)+gamma*(900/T+273)*u*(es-ea))/(delta+gamma*(1+0.34*u));
%---------------------------------------------END!!!!
I hope this is helpful for someone. FEEL FREE TO TELL ME IF I MISSED SOMETHING IN THE EQUATION. Ultimately, MatLab "toys" like to grow into bigger things, I feel like, so lest I mess up... come help!
In case you want to know more about PM ET function:
Wiki
UF Extension-- go extension!
UCD more complicated
