function ETmodel
%LAST MODIFIED: 22AUG07
clear all

imax=360;
jmax=720;

daynum=[15 46 74 105 135 166 196 227 258 288 319 349];                                      %based on the middle of each month
kPAR=0.30;                                                                                  %kPAR
kRn=0.60;                                                                                   %kRn

%Clear the Matrices
for y=1986:1995,%1995,%for 1994, change all k=1:x to 1:9
    year=num2str(y)
    year2dig=y-1900;
    
    %Read in the data and store them into 3D arrays    
    mmax=12;
    for m=1:mmax, 
        if m<10,
            Taname=['raw/cru5_temp_mean_hd_',num2str(year),'0',num2str(m),'.asc'];
            Tdiname=['raw/cru5_diurnal_tmp_hd_',num2str(year),'0',num2str(m),'.asc'];            
            NDVIname=['raw/fasir_ndvi413_hd_',num2str(year),'0',num2str(m),'.asc'];
            Rnname=['raw/srb_ntf_av_1d_',num2str(year),'0',num2str(m),'00.asc'];
            eaname=['raw/cru5_vapor_pres_hd_',num2str(year),'0',num2str(m),'.asc'];
            ch1aname=['raw/avhrrpf.ch1.1ntegl.',num2str(year2dig),'0',num2str(m),'01.txt'];
            ch1bname=['raw/avhrrpf.ch1.1ntegl.',num2str(year2dig),'0',num2str(m),'11.txt'];
            ch1cname=['raw/avhrrpf.ch1.1ntegl.',num2str(year2dig),'0',num2str(m),'21.txt'];
            ch2aname=['raw/avhrrpf.ch2.1ntegl.',num2str(year2dig),'0',num2str(m),'01.txt'];
            ch2bname=['raw/avhrrpf.ch2.1ntegl.',num2str(year2dig),'0',num2str(m),'11.txt'];
            ch2cname=['raw/avhrrpf.ch2.1ntegl.',num2str(year2dig),'0',num2str(m),'21.txt'];
            Ta(:,:,m)=load(Taname);
            Tdi(:,:,m)=load(Tdiname);
            Tmax(:,:,m)=Ta(:,:,m)+Tdi(:,:,m)/2;
            NDVI(:,:,m)=load(NDVIname);
            Rn0(:,:,m)=load(Rnname);
            ea(:,:,m)=load(eaname);
            ch1a(:,:,m)=load(ch1aname);
            ch1b(:,:,m)=load(ch1bname);
            ch1c(:,:,m)=load(ch1cname);
            ch2a(:,:,m)=load(ch2aname);
            ch2b(:,:,m)=load(ch2bname);
            ch2c(:,:,m)=load(ch2cname);
            ch1avg(:,:,m)=(((ch1a(:,:,m)+ch1b(:,:,m)+ch1c(:,:,m))/3) - 10) *.002*.01;                       % averages channel and scales
            ch2avg(:,:,m)=(((ch2a(:,:,m)+ch2b(:,:,m)+ch2c(:,:,m))/3) - 10) *.002*.01;
            
        else
            Taname=['raw/cru5_temp_mean_hd_',num2str(year),num2str(m),'.asc'];
            Tdiname=['raw/cru5_diurnal_tmp_hd_',num2str(year),num2str(m),'.asc'];
            NDVIname=['raw/fasir_ndvi413_hd_',num2str(year),num2str(m),'.asc'];   
            Rnname=['raw/srb_ntf_av_1d_',num2str(year),num2str(m),'00.asc'];
            eaname=['raw/cru5_vapor_pres_hd_',num2str(year),num2str(m),'.asc'];
            ch1aname=['raw/avhrrpf.ch1.1ntegl.',num2str(year2dig),num2str(m),'01.txt'];
            ch1bname=['raw/avhrrpf.ch1.1ntegl.',num2str(year2dig),num2str(m),'11.txt'];
            ch1cname=['raw/avhrrpf.ch1.1ntegl.',num2str(year2dig),num2str(m),'21.txt'];
            ch2aname=['raw/avhrrpf.ch2.1ntegl.',num2str(year2dig),num2str(m),'01.txt'];
            ch2bname=['raw/avhrrpf.ch2.1ntegl.',num2str(year2dig),num2str(m),'11.txt'];
            ch2cname=['raw/avhrrpf.ch2.1ntegl.',num2str(year2dig),num2str(m),'21.txt'];
            Ta(:,:,m)=load(Taname);
            Tdi(:,:,m)=load(Tdiname);
            Tmax(:,:,m)=Ta(:,:,m)+Tdi(:,:,m)/2;
            NDVI(:,:,m)=load(NDVIname);
            Rn0(:,:,m)=load(Rnname);
            ea(:,:,m)=load(eaname);
            ch1a(:,:,m)=load(ch1aname);
            ch1b(:,:,m)=load(ch1bname);
            ch1c(:,:,m)=load(ch1cname);
            ch2a(:,:,m)=load(ch2aname);
            ch2b(:,:,m)=load(ch2bname);
            ch2c(:,:,m)=load(ch2cname);
            ch1avg(:,:,m)=(((ch1a(:,:,m)+ch1b(:,:,m)+ch1c(:,:,m))/3) - 10) *.002*.01;
            ch2avg(:,:,m)=(((ch2a(:,:,m)+ch2b(:,:,m)+ch2c(:,:,m))/3) - 10) *.002*.01;
        end
    end
    
    for i=1:180,
        for j=1:360,
            Rn(2*i-1,2*j-1,:)=Rn0(i,j,:);
            Rn(2*i-1,2*j,:)=Rn0(i,j,:);
            Rn(2*i,2*j-1,:)=Rn0(i,j,:);
            Rn(2*i,2*j,:)=Rn0(i,j,:);
        end
    end
    
    for i=1:imax,
        lat(i)=-1.0056*i+91.006;
        
        for j=1:jmax,
            for m=1:mmax,
                SAVI(i,j,m)=0.45*NDVI(i,j,m)+0.132;  unadjusted SAVI
                SAVI_AVHRR(i,j,m)=1.5*(ch2avg(i,j,m)-ch1avg(i,j,m))/(ch2avg(i,j,m)+ch1avg(i,j,m)+0.5);  %unadjusted SAVI
                NDVI_FASIR0(i,j,m)=NDVI(i,j,m);                                                         %NDVI FASIR
                NDVI_AVHRR(i,j,m)=NDVI(i,j,m);                                                          %NDVI AVHRR (used to correct SAVI AVHRR)
                NDVI_AVHRR(i,j,m)=(ch2avg(i,j,m)-ch1avg(i,j,m))/(ch2avg(i,j,m)+ch1avg(i,j,m));          %NDVI AVHRR (used to correct SAVI AVHRR)
            end
        end
    end
    SAVImax(i,j)=max(SAVI0(i,j,:));
    SAVImin(i,j)=min(SAVI0(i,j,:));
    reshapeNDVI_FASIR=reshape(NDVI_FASIR0,imax*jmax*mmax,1); 
    NDVI_FASIRmax=max(reshapeNDVI_FASIR);
    NDVImin=min(reshapeNDVI_AVHRR);
    if NDVImin<0, NDVImin=0; end
    for i=1:imax,
        for j=1:jmax,
            for m=1:mmax,
                if NDVI_FASIRmax==0, NDVI_FASIR(i,j,m)=0;
                else NDVI_FASIR(i,j,m)=NDVI(i,j,m)/NDVI_FASIRmax; end                                   %NDVI FASIR stretched to 0-1
                if NDVI(i,j,m)<=0, SAVI(i,j,m)=0;
                else SAVI(i,j,m)=SAVI_AVHRR(i,j,m)*NDVI_FASIR(i,j,m)/NDVI_AVHRR(i,j,m); end             %adjusted SAVI AVHRR to FASIR
                
                %Biophysical Parameters
                daylength(i,j,m)=(exp(7.42+0.045*lat(i))/3600)*(sin((daynum(m)-79)*0.01721))+12; 
                if daylength(i,j,m) > 24, daylength(i,j,m)=24; end
                if daylength(i,j,m) <= 0, daylength(i,j,m)=0; end                                       %Daylength bound to 0 and 24
                
                es_mid(i,j,m)=.61121*exp(17.502*Tmax(i,j,m)/(Tmax(i,j,m)+240.97));                      %es: Saturation Pressure of Water Vapor
                es_avg(i,j,m)=.61121*exp(17.502*Ta(i,j,m)/(Ta(i,j,m)+240.97));              
                s_mid(i,j,m)=17.502*240.97*es_mid(i,j,m)/((Tmax(i,j,m)+240.97)^2);                      %s: Rate of Change of Vapor Pressure with Temperature
                s_avg(i,j,m)=17.502*240.97*es_avg(i,j,m)/((Ta(i,j,m)+240.97)^2);                
                EFeq_mid(i,j,m)=s_mid(i,j,m)/(s_mid(i,j,m)+.066);
                EFeq_avg(i,j,m)=s_avg(i,j,m)/(s_avg(i,j,m)+.066);
                VPD_mid(i,j,m)=es_mid(i,j,m)-.1*ea(i,j,m);                                              %VPD: Vapor Pressure Deficit
                VPD_avg(i,j,m)=es_avg(i,j,m)-.1*ea(i,j,m);                                   
                RH_mid(i,j,m)=min(1,max(0,.1*ea(i,j,m)/es_mid(i,j,m)));                                 %RH: Relative Humidity
                RH_avg(i,j,m)=min(1,max(0,.1*ea(i,j,m)/es_avg(i,j,m)));                         
                
                if daylength(i,j,m)==0, RnDay(i,j,m)=0; RnDaymm(i,j,m)=0; fIPAR(i,j,m)=0; fAPAR(i,j,m)=0; LAI(i,j,m)=0; Rns(i,j,m)=0; Rnc(i,j,m)=0;
                else
                    RnDay(i,j,m)=Rn(i,j,m);                                                             %RnDay: daytime Net Radiation
                    RnDaymm(i,j,m)=RnDay(i,j,m)*1.038;                                                  %RnDaymm: converts Rn W/m2 to mm
                    fIPAR(i,j,m)=NDVI(i,j,m)-0.05;                                                      %fIPAR: fraction of Intercepted Photosynthetically Active Radiation
                    fAPAR(i,j,m)=1.2*(1.136*SAVI(i,j,m)-0.04);                                          %fAPAR: fraction of Absorbed Photosynthetically Active Radiation
                    LAI(i,j,m)=-(1/kPAR)*log(1-fIPAR(i,j,m));                                           %LAI: Leaf Area Index 
                    Rns(i,j,m)=RnDaymm(i,j,m)*exp(-kRn*LAI(i,j,m));                                     %Rns: Soil Net Radiation
                    Rnc(i,j,m)=RnDaymm(i,j,m)-Rns(i,j,m);                                               %Rnc: Canopy Net Radiation
                end
                growthopt(i,j,m)=RnDaymm(i,j,m)*SAVI(i,j,m)*Tmax(i,j,m)/VPD_avg(i,j,m);                 %Growthopt or Phen
            end
            growthoptmax(i,j)=max(growthopt(i,j,:));                                                    %Max Growthopt
            fAPARmax(i,j)=max(fAPAR(i,j,:));                                                            %Max fAPAR
            for m=1:mmax,
                if growthopt(i,j,m)==growthoptmax(i,j), Topt(i,j)=Tmax(i,j,m); end                      %Topt: Ta at Max Growthopt
                if growthopt(i,j,m)==growthoptmax(i,j), JDopt(i,j)=daynum(m); end                       %JDopt: JD at Max Growthopt
            end
            
            %Constraint Functions       
            for m=1:mmax,
                fwet(i,j,m)=RH_avg(i,j,m)^4;                                                            %fwet: Relative Surface Wetness
                fc(i,j,m)=1*NDVI(i,j,m)-0.05;                                                           %fc: Fractional Total Vegetation Cover
                if fc(i,j,m)==0, fg(i,j,m)=0;
                else fg(i,j,m)=min(1,max(0,fAPAR(i,j,m)/fc(i,j,m))); end                                %fg: Green Canopy Fraction
                if Topt(i,j)==0, fT(i,j,m)=0;
                else fT(i,j,m)=min(1,1.14/(1+exp(0.2*(Topt(i,j)-10-Tmax(i,j,m))))); end                 %fT: Plant Temperature Constraint
                if fAPARmax(i,j)==0, fM(i,j,m)=0; 
                elseif daynum(m)<JDopt(i,j), fM(i,j,m)=1; 
                else fM(i,j,m)=fAPAR(i,j,m)/fAPARmax(i,j,:); end                                        %fM: Plant Moisture Constraint
                if RH_mid(i,j,m)==0, fSM(i,j,m)=0;
                else fSM(i,j,m)=RH_mid(i,j,m)^(VPD_mid(i,j,m)); end%/RH_mid(i,j,m)); end                %fSM: Soil Moisture Constraint
                if fSM(i,j,m)>1, fSM(i,j,m)=1; end
                if fg(i,j,m)>1, fg(i,j,m)=1; end
                
                %Final LE Equations
                %PT(i,j,m)=1.26*EFeq_mid(i,j,m)*RnDaymm(i,j,m);%*daylength(i,j,m)*0.0000004*3600*30;
                LEsdry(i,j,m)=(1-fwet(i,j,m))*fSM(i,j,m)*1.26*EFeq_mid(i,j,m)*Rns(i,j,m);%*daylength(i,j,m)*0.0000004*3600*30; %LEsdry
                LEswet(i,j,m)=fwet(i,j,m)*1.26*EFeq_mid(i,j,m)*Rns(i,j,m);%*daylength(i,j,m)*0.0000004*3600*30; %LEswet
                LEcdry(i,j,m)=(1-fwet(i,j,m))*fT(i,j,m)*fg(i,j,m)*fM(i,j,m)*1.26*EFeq_mid(i,j,m)*Rnc(i,j,m);%*daylength(i,j,m)*0.0000004*3600*30; %LEcdry
                LEcwet(i,j,m)=fwet(i,j,m)*1.26*EFeq_mid(i,j,m)*Rnc(i,j,m);%*daylength(i,j,m)*0.0000004*3600*30; %LEcwet aka Interception Evaporation or LEi
                LE(i,j,m)=LEsdry(i,j,m)+LEswet(i,j,m)+LEcdry(i,j,m)+LEcwet(i,j,m);
                %Uniform negative values: -999
                if land_mask(i,j)==0, LE(i,j,m)=-999; end
                if LE(i,j,m)<0, LE(i,j,m)=0; end
                
                %Write out variables
                for m=1:mmax,
                    if m<10,
                        %             fullPath=['output/','LEsdry_',num2str(year),'_0',num2str(m),'.txt'];
                        %             varName=LEsdry(:,:,m);
                        %             str1 = sprintf ('save %s varName -ascii', fullPath);
                        %             eval (str1);
                        %             fullPath=['output/','LEswet_',num2str(year),'_0',num2str(m),'.txt'];
                        %             varName=LEswet(:,:,m);
                        %             str1 = sprintf ('save %s varName -ascii', fullPath);
                        %             eval (str1);
                        %             fullPath=['output/','LEcwet_',num2str(year),'_0',num2str(m),'.txt'];
                        %             varName=LEcwet(:,:,m);
                        %             str1 = sprintf ('save %s varName -ascii', fullPath);
                        %             eval (str1);
                        %             fullPath=['output/','LEcdry_',num2str(year),'_0',num2str(m),'.txt'];
                        %             varName=LEcdry(:,:,m);
                        %             str1 = sprintf ('save %s varName -ascii', fullPath);
                        %             eval (str1);
                        %                                     fullPath=['output/','LE_',num2str(year),'_0',num2str(m),'.txt'];
                        %                                     varName=LE(:,:,m);
                        %                                     str1 = sprintf ('save %s varName -ascii', fullPath);
                        %                                     eval (str1);
                        %             fullPath=['output/','PT_',num2str(year),'_0',num2str(m),'.txt'];
                        %             varName=PT(:,:,m);
                        %             str1 = sprintf ('save %s varName -ascii', fullPath);
                        %             eval (str1);
                    else
                        %                         fullPath=['output/','LEsdry_',num2str(year),'_',num2str(m),'.txt'];
                        %                         varName=LEsdry(:,:,m);
                        %                         str1 = sprintf ('save %s varName -ascii', fullPath);
                        %                         eval (str1);
                        %                         fullPath=['output/','LEswet_',num2str(year),'_',num2str(m),'.txt'];
                        %                         varName=LEswet(:,:,m);
                        %                         str1 = sprintf ('save %s varName -ascii', fullPath);
                        %                         eval (str1);
                        %                         fullPath=['output/','LEcwet_',num2str(year),'_',num2str(m),'.txt'];
                        %                         varName=LEcwet(:,:,m);
                        %                         str1 = sprintf ('save %s varName -ascii', fullPath);
                        %                         eval (str1);
                        %                         fullPath=['output/','LEcdry_',num2str(year),'_',num2str(m),'.txt'];
                        %                         varName=LEcdry(:,:,m);
                        %                         str1 = sprintf ('save %s varName -ascii', fullPath);
                        %                         eval (str1);
                        %                                     fullPath=['output/','LE_',num2str(year),'_',num2str(m),'.txt'];
                        %                                     varName=LE(:,:,m);
                        %                                     str1 = sprintf ('save %s varName -ascii', fullPath);
                        %                                     eval (str1);
                        %             fullPath=['output/','PT_',num2str(year),'_',num2str(m),'.txt'];
                        %             varName=PT(:,:,m);
                        %             str1 = sprintf ('save %s varName -ascii', fullPath);
                        %             eval (str1);
                    end
                end
            end