Index exceeds matrix dimensions : Error

I have a long for loop with many IF statements. I keep arriving at the error:

Index exceeds matrix dimensions.
  Error in ESCIMOloop (line )
      ShortRadBal(t)=(1 - Albedo(t))*GlobalRad(t);
  Attempted to access Albedo(2); index out of bounds because numel(Albedo)=1.

Albedo is calculated earlier in the for-loop, so I don't understand why following equations are saying it is out of bounds. The Albedo at T(1)=0, but all following albedos should be calculated starting with line 79. I'm new to matlab and very confused.

        %ESCIMO energy balance method for snow melting.
        %Data from the ... Switzerland
        clear all

        x=load('C:\Users\wesser\Desktop\P2\ESCIMO\ESCIMOdata.dat'); %Load the file with data

        %Field Observed Parameters
        td=x(:,4);              %hour of the day
        T=x(:,5);                %Air Temperature (K)
        Relhum=x(:,6);       %Relative Humidity (%)
        Wind=x(:,7);          %Wind Speed (m/s)
        P=x(:,8);                %Precipitation (mm/h)
        GlobalRad=x(:,9);          %Global Radiation (W/m2)
        Inlong=x(:,10);         %Incoming Longwave Radiation (W/m2)
        Snow_obs=x(:,11);           %SWE Observed

        %List parameters that start with a known value at t=1
        VapPresS(1)=0;      %Surface Vapor Pressure is 0 a the start
        Sublim(1)=0;     %(Re-)Sublimation is 0 at the start
        Melt(1)=0;      %Melt is 0 at the start
        ModeledSWE(1)=0;    %Modeled SWE is 0 at the start
        SnowAge(1)=0;       %Snow age is 0 at the start
        Albedo(1)=0;        %Albedo is 0 at the start

        %Parameters (to be CALIBRATED)
        PM=importdata('C:\Users\wesser\Desktop\P2\ESCIMO\parameter.txt'); % model parameters input

        Amin=PM(1,:);           %Minimum Albedo
        Aadd=PM(2,:);           %Additive Albedo
        DPp=PM(3,:);            %Decline parameter (Positive temperatures)
        DPn=PM(4,:);            %Decline parameter (negative temperatures)
        SS=PM(5,:);             %Significant snowfall (mm/h)
        Tp=PM(6,:);             %Phase transition Temperature (K)
        Se=PM(7,:);             %Snow Emissivity 
        SHF=PM(8,:);            %Soil Heat Flux (W/m2)

        %Calibrated parameter ranges ------TBD
        %Amin - # to #         %Minimum Albedo
        %Aadd - # to #          %Additive Albedo
        %DPp - # to #            %Decline parameter (Positive temperatures)
        %DPn - # to #            %Decline parameter (negative temperatures)
        %SS - # to #             %Significant snowfall (mm/h)
        %Tp - # to #             %Phase transition Temperature (K)
        %Se - # to #             %Snow Emissivity 
        %SHF - # to #            %Soil Heat Flux (W/m2)

        %PARAMETERS (used default)
        SHCW=4180.0;            %Specific heat capacity of water (J/(kgK))
        SHCS=2100.0;            %Specific heat capacity of snow (J/(kgK))
        MH=337500.0;             %Melting heat (J/kg)
        SBC=0.0000000567;       %Stefan-Boltzmann-Constant (W/(m2K^4))
        CSH=2835500.0;          %Condensation/Sublimation heat (J/kg)

        %Precipitation Loop
        for t=2 : length(P);             %time series loop
            if T(t)>=Tp;       %Precip that is rain       
               Pr(t)=P(t);
            else
               Pr(t)=0;
            end

            if T(t)<Tp;      %Precip that is snow
               Ps(t)=P(t);
            else
               Ps(t)=0;
            end

            if Ps(t)>SS;         %Snow age
                SnowAge(t)=0;
                if ModeledSWE(t-1)>0;
                    SnowAge(t)=SnowAge(t-1)+0.04167;
                else
                    SnowAge(t)=0;
                end
            end

            if T(t)>273.16;
                Z(t)=DPp;
            else
                Z(t)=DPn;
            end
            if Ps(t)>SS;    %Albedo
                Albedo(t)=Amin+Aadd;
                if SnowAge(t)>0;
                    Albedo(t)=Amin+(Albedo(t-1)-Amin)*exp(Z(t)*0.04167);
                else   
                   Albedo(t)= 0;
                end
            end

            SurfTemp(t)=min(T(t),273.16);

            VapPresA(t)=6.1078*exp((17.08085*(T(t)-273.16))/(234.175+(T(t)-273.16)))*(Relhum(t)/100);

            if ModeledSWE(t-1)>0;
                VapPresS(t)=6.1078*exp((17.088085*(SurfTemp(t)-273.16))/(234.175+(SurfTemp(t)-273.16)));
            else
                VapPresS(t)=0;
            end

            ShortRadBal(t)=(1- Albedo(t))*GlobalRad(t);
            LongRadBal(t)=Inlong(t)-(Se*SBC*SurfTemp(t)^4);
            SenHeatFlux(t)=18.85*(0.18+0.098*Wind(t))*(T(t)-SurfTemp(t));
            LatHeatFlux(t)=32.82*(0.18+0.098*Wind(t))*(VapPresA(t)-VapPresS(t));

            if Pr(t)>0;
                AdvFluxR(t)=SHCW*(T(t)-273.16)*Pr(t)/3600;
            else
                AdvFluxR(t)=0;
            end

            if Ps(t)>0;
                AdvFluxS(t)=SHCS*(T(t)-273.16)*Ps(t)/3600;
            else
                AdvFluxS(t)=0;
            end

            EB(t)=ShortRadBal(t)+LongRadBal(t)+SenHeatFlux(t)+LatHeatFlux(t)+AdvFluxR(t)+AdvFluxS(t)+SHF;

            if ModeledSWE(t-1)>0;
               Sublim(t)=3600*LatHeatFlux(t)/CSH;
            else
               Sublim(t)=0;
            end

            if T(t)>273.16;
                if ModeledSWE(t-1)>0 & EB(t)>0;
                   Melt(t)= min(EB(t)*3600/MH ,ModeledSWE(t-1));
                else
                    Melt(t)=0;
                end
            else
                Melt(t)=0;
            end

            if ModeledSWE(t-1)>0;
                max((ModeledSWE(t-1)+Pr(t)+Ps(t)+Sublim(t)-Melt(t)),0);
            else
                max((Ps(t)+Sublim(t)-Melt(t)),0);
            end
        end

Answer 1

[Roger Wohlwend]

Albedo(t) only exist if Ps(t) > SS. It seems that for k = 2 this condition is not met, the code in the if-clause is not execute and Albedo(2) is never calculated. That is why Albedo(2) does not exist and the error occurs. You can solve the Problem by introducting an else-clause

if Ps(t)>SS
   Albedo(k) = something
else
   Albedo(k) = something else
end

With every Iteration you increase the Albedo vector and several other vectors. That is not a good programming technique. I recommend preallocating Memory for all those vectors. That makes the code run faster. For more Information search for "Preallocation" in the Matlab documentary. In your case insert the following line before the Loop

Albedo = Zeros(length(P),1)

and repeat that for all other vectors (VapPresS, Sublim, Melt, ..). Then you can delete Albedo(1) = 0 (and Melt(1) = 0, Sublim(1) = 0, ...).

[Image Analyst]

You're lucky Roger is trying to help you. Many people (like me) are not going to examine lengthy code like that unless it can be copied and pasted and run in the debugger, which it can't because you forgot to attach ESCIMOdata.dat.

[hydrogall]

You're right. I am lucky to have had Roger's help as now my code is running perfectly. I am very appreciative of his help and grateful of those who take the time to teach a newbie. We all go through the learning curve.