function [ET_ro Ao] = CIMIS_ET(time,T,U,RH,Rsd,Lat,Lon,Z)
% INPUT
% time: matlab time vector (midpoint) (HAS TO BE A COLUMN VECTOR!)
% T: measured Temperature [oC]
% U: measured wind speed [m/s]
% RH: measured relative humidity [%]
% Rsd: measured shortwave downwelling solar radiation [W/m2]
% Z: station altitude above sea level [m]
% Lat,Lon: Latitude and Longitude of station [degrees, West is negative]
% OUTPUT
% ET_ro: ASCE-ET (correct, confirmed by worksheet of Rick Snyder)
% Ao: aerodynamic contribution to ET
% ET: CIMIS ET (not correct, algorithm sketchy)

% REFERENCES
% ASCE Method:
% http://www.fao.org/docrep/X0490E/x0490e07.htm#calculation%20procedures
% http://biomet.ucdavis.edu/evapotranspiration.html
% CIMIS Method:
% http://wwwcimis.water.ca.gov/cimis/infoEtoCimisEquation.jsp

dt = median(diff(time));
time_vec = datevec(time);
DoY = ceil( time - datenum(time_vec(1),1,1) );
ToD = time_vec(:,4)+time_vec(:,5)/60;
Lz = 120; %  longitude of the local time meridian (degrees West), 120° for Pacific time zone 
Lon = - Lon;

Tk = T + 273.16; % not needed for CIMIS
Rsd = Rsd * 3600 / 1e6 * (dt*24); % convert from W/m2 to MJ/m2/h; equs are written for hourly 
% Saturation vapor pressure 
es = 0.6108 * exp(T * 17.27./ (T + 237.3));
ea = RH/100 .* es;

% find net longwave radiation
GSC = 0.082; % MJ m -2 min-1 (solar constant)
% Ra: extraterrestrial radiation
% dr = inverse relative distance Earth-Sun (correction for eccentricity of Earth's orbit around the sun) 
% omega1 = Solar time angle 1/2 hour before omega, that is, at the beginning of period (radians) 
% omega2 = Solar time angle 1/2 hour after  omega, at end of period (radians) 
% fe = Station latitude (radians) 
% delta = Declination of the sun above the celestial equator (radians)
% theta = solar altitude 
dr = 1+0.033*cos(2*pi/365 * DoY); 
        b = 2*pi/364*(DoY-81);
    Sc = 0.1645 * sin(2*b) - 0.1255 * cos(b) - 0.025*sin(b);
omega = pi/12 * ( ToD + 4/60*(Lz-Lon)+Sc -12  ) ; % CIMIS: t must be center of time interval; 
omega = pi/12 * ( ToD + 1/15*(Lz-Lon)+Sc -12  ) ; % ASCE

omega1 = omega - (dt/2*24)*pi/12;
omega2 = omega + (dt/2*24)*pi/12;
    phi = pi*Lat/180;
delta = 0.409 * sin(2*pi/365*DoY-1.39);

%     omegas = acos(-tan(phi).*tan(delta));
% Ra = 24*60/pi*GSC*dr .* ( omegas *sin(phi).*sin(delta) +cos(phi)*cos(delta).*sin(omegas) ); % FAO daily
theta = asin(sin(phi)*sin(delta) + cos(phi)*cos(delta).*cos(omega) )*180/pi; % CIMIS & ASCE
sin_theta = (omega2-omega1).*sin(phi).*sin(delta) + cos(phi)*cos(delta).*(sin(omega2)-sin(omega1)); % ASCE

Ra = 6*60/pi*GSC*dr .* ( (omega2-omega1) *sin(phi).*sin(delta) + cos(phi)*cos(delta).*(sin(omega2)-sin(omega1)) ); % CIMIS
Ra =12*60/pi*GSC*dr .* sin_theta; % ASCE

% Rso  = clear sky total global radiation at the earth surface (MJ m-2 h-1) 
% Rspm = solar radiation (MJm-2 h-1) during the last hour for which ?>10° 
% Rsopm= clear sky solar radiation (MJ m-2h-1) during the last hour for which ?>10 ° 
Rso = Ra * (0.75+2e-5*Z);

emissivity = 0.34 - 0.14 * sqrt(ea);
Rsdpm = repmat(NaN,length(time),1); Rsopm = repmat(NaN,length(time),1); % need to determine these later (nigthtime onlY)

f(theta>10) = 1.35 * Rsd(theta>10)./Rso(theta>10) - 0.35;
f(theta<=10)= 1; % assume it was sunny during last daylight episode
for i = 2:length(f)
if theta(i)<10, f(i) = f(i-1); end  
end

f = f';

f(f<0)=0.595; % CIMIS
f(f>1)=1;

% Long wave Radiation Equation (MJ m-2 h-1) 
Rnl = emissivity .* f * 2.04e-10 .* Tk.^4 *(dt*24);

% f = cloudiness factor 
% epsilon = net emissivity of the surface 
% ? = Steffan-Boltzman constant (2.04*10-10) MJ m -2 h-1 K-4 

% net radiation for albedo = 0.23
Rn = (0.77 * Rsd - Rnl)*1e6/3600;
 
% VPD - Vapor pressure deficit 
VPD = es - ea; % (kPa) 
% DEL - Slope of the saturation vapor pressure vs. air temperature curve at the average hourly air temperature 
DEL = (4099 * es)./(T + 237.3).^2;
% DEL_= es./Tk.*(6790.5./Tk-5.028); % given in Ventura, D. Spano, P. Duce, R.L. Snyder (Irrigaiton Science 1999) --> SAME RESULT

% Barometric pressure 
P = 101.3 - 0.0115 * Z + 5.44 * 10^-7 * Z^2; % CIMIS
% GAM - Psychrometer constant (kPa C-1) 
GAM = 0.000646 * (1 + 0.000946*T) .* P; % CIMIS
% W - Weighting function 
W = DEL./(DEL + GAM);
% FU2 - Wind function 
FU2 = repmat(NaN,length(time),1);   
FU2(Rn<=0)= 0.125 + 0.0439 * U(Rn<=0);
FU2(Rn>0) = 0.030 + 0.0576 * U(Rn>0) ;
% NR - Convert Rn from Wm-2 to mm 
NR = Rn./(694.5*(1-0.000946*T));
% Hourly ETo is approximately equal to RET 
ET = W.*NR + (1-W).*VPD .* FU2; %[mm/day]

% ASCE Method % http://biomet.ucdavis.edu/evapotranspiration.html 
Rn = 0.77 * Rsd - Rnl; % MJ/m2/h
P = 101.3 * ((293-0.0065*Z)/293)^5.26; % ASCE
lambda = 2.501-2.361*0.001*T;
GAM = 0.00163 ./ lambda .* P; % ASCE
% day:  ra = 208/u2; rs = 50; 
% night:ra = 208/u2; rs =200; 
% for u<0.5: u = 0.5 m/s;
U(U<0.5) = 0.5;
G = repmat(NaN,length(time),1); GAM_star = repmat(NaN,length(time),1);
GAM_star(Rn>0) = GAM(Rn>0).*(1+0.24*U(Rn>0));  
GAM_star(Rn<=0)= GAM(Rn<=0).*(1+0.96*U(Rn<=0));
G(Rn>0) = 0.1 * Rn(Rn>0); G(Rn<=0) = 0.5 * Rn(Rn<=0);

Ro = DEL./lambda.*(Rn-G)./(DEL+GAM_star);
Ao = (37*GAM./Tk).*U.*(es-ea) ./ (DEL+GAM_star) * (dt*24);
ET_ro = Ro+Ao;

% TEST ALGORITHM USING RIMC FIELD DATA
% [datastr, trash,data(:,2),data(:,3),data(:,4),data(:,5),data(:,6),data(:,7),data(:,8),data(:,9),data(:,10),...
%     data(:,11),data(:,12),data(:,13),data(:,14),data(:,15),data(:,16),data(:,17),data(:,18),data(:,19),data(:,20),data(:,21)] = ...
%     textread('C:\Documents and Settings\Jan Kleissl\My Documents\1ucsd\web\demroes\data\RIMCSoil\RIMC_20090305.dat',...
%     '%21c%u%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f','delimiter',',','headerlines',4); clear trash
% time_RIMC = datenum( str2num(datastr(:,2:5)) , str2num(datastr(:,7:8)),str2num(datastr(:,10:11)) ,...
%                 str2num(datastr(:,13:14)) , str2num(datastr(:,16:17)) ,str2num(datastr(:,19:20)));
% Rsd = interp1(time_RIMC,data(:,5),time_RIMF);
% load temp.mat
% time = time_RIMF; Lat = 32.885248; Lon = -117.240184; Z = 424*.3048;
% U([1:117 404:574]) = []; RH([1:117 404:574]) = []; Rsd([1:117 404:574]) = []; T([1:117 404:574]) = []; time_RIMF([1:117 404:574]) = [];

% TEST ALGORITHM BY COMPARING AGAINST CIMIS FILES
% from file: CIMIS135Blythe_Aug25-Oct052008_daily.csv
% CIMIS_ET_ = [5.334 6.096 6.604 5.08 4.572 4.572 5.842 6.35 5.842 6.096 6.096 6.35 5.588 6.096 6.096 ...
%      5.08 3.556 5.08 5.588 5.334 5.588 6.35 4.318 5.08 5.334 5.588 5.842 5.334 5.08 5.334 5.842 5.334 ...
%      5.588 6.096 6.858 5.588 5.08 5.334 5.842 3.81 4.318]';
% CIMIS_Rsd= [587 602 557 473 445 441 585 606 585 614 604 606 486 523 595 549 399 591 587 589 576 560 441 ...
%     543 538 554 572 575 563 543 539 525 518 498 516 500 510 489 478 321 513]'* .484 * 24 / 694.5 *.77; % Ly/day --> mm, Rs,net
% daystarts = find(diff(time_vec(:,3))~=0);
% figure(10), hold on,figure(11), hold on
% for i = 1:length(daystarts)
%     trash = ET(daystarts(i):daystarts(i)+round(1/dt));
%     trash1 = ET_ro(daystarts(i):daystarts(i)+round(1/dt));
%     trash2= W(daystarts(i):daystarts(i)+round(1/dt)).*NR(daystarts(i):daystarts(i)+round(1/dt)); 
%     trash3 = (1-W(daystarts(i):daystarts(i)+round(1/dt))).*VPD(daystarts(i):daystarts(i)+round(1/dt)) .* FU2(daystarts(i):daystarts(i)+round(1/dt));
%     trash4= Rnl(daystarts(i):daystarts(i)+round(1/dt))*1e6/3600/694.5; 
%     fprintf('ET on %s =%5.2f, ASCE: %5.2f, W*NR=%5.2f, (1-W)*VPD*FU2=%5.2f, Rnl =%5.2f, CIMIS:ET =%5.2f, Rsd =%5.2f, Rsd-ET =%5.2f\n',datestr(time(daystarts(i)+1)),sum(trash(~isnan(trash)==1)),sum(trash1(~isnan(trash1)==1)),...
%         sum(trash2(~isnan(trash2)==1)), sum(trash3(~isnan(trash3)==1)), sum(trash4(~isnan(trash4)==1)), CIMIS_ET_(i), CIMIS_Rsd(i), -CIMIS_ET_(i)+CIMIS_Rsd(i) ); % sum(CIMIS2(daystarts(i):daystarts(i)+round(1/dt),8))
%     figure(10),plot(sum(trash2(~isnan(trash2)==1)),CIMIS_Rsd(i),'k.')
%     figure(11),plot(sum(trash(~isnan(trash)==1)),CIMIS_ET_(i),'k.')
% end
% figure(10),xlabel('W*Rn'),ylabel('Rsd'),axis([2.5 8 2.5 8])
% figure(11),xlabel('ET'),ylabel('ET [CIMIS]'),axis([2.5 8 2.5 8])
% 
% figure,subplot(2,1,1),plot(Ra,'y-'),hold on,plot(Rso,'k-'),plot(Rsd,'r-'),plot(Rnl,'b-')
% axis([0 48 0 5])
% subplot(2,1,2),plot(W,'k:'),hold on,plot(NR,'b-'),plot(VPD .* FU2,'r-');axis([0 48 -inf inf])
% legend('W','NR','VPD*FU2')
% 
% XTicks=[datenum(2008,7,1):datenum(2009,1,1)]; XTickLabels=datestr(XTicks,6);
% figure,plot(time,Rsd,'k-')
% set(gca,'XTick',XTicks,'XTickLabel',XTickLabels,'XGrid','on','Box','on')
% axis([datenum(2008,9,16) datenum(2008,9,19) 0 inf])