my silly penman-monteith code

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