Create a vector to identify storm time steps

I am trying to index storms from a matlab script I received from a colleague which generates storm duration statistics by identifying "holes". Essentially I want to create a vector which identifies which storm each time step belongs to. The function would number storms based on the arrival date and end date. Forgive me I am not very familiar with Matlab functions.

Data is daily precipitation like this: 
39240  0
39241  0
39242  0
39243  0
39244  0
39245  16.2000000000000
39246  0.200000000000000
39247  0
39248  0
39249  0
39250  0
39251  0.200000000000000
39252  0
39253  12.4000000000000
39254  0.200000000000000
39255  0
39256  0
39257  0
39258  28.6000000000000
39259  0.200000000000000
39260  0
39261  0
39262  34
39263  0.800000000000000
39264  0

The function itself I am trying to alter:

% Conducts analysis of a rainfall time-series; calculating storm frequency,
% average storm intensity, and average depth of rainfall per storm

% created by xxxx, last updated 10/16/12

% input data should be in the format of date/time (:,1) and total rainfall
% per observation (:,2)
% output units will be equal to input units

% lines 34 and 43 index and identify storm start and end dates, if needed,
% currently commented-out

function [out, arrival_dates, end_date] = storm_stat(input_rain)

% Function outputs:
% out(:,1) = st_num = number of storms during record
% out(:,2) = st_arrival = arrival time of storms (start time to start time)
% out(:,3) = st_inter = length between storms (end time to start time)
% out(:,4) = st_dur = storm duration
% out(:,5) = st_depth = depth of rainfall during each storm
% out(:,6) = st_inten = average depth of rain per unit time

rain_date = input_rain(:,1);
rain_obs = input_rain(:,2);

% Index location of storms within the input data
[~,L,N] = bwboundaries(rain_obs>1,'noholes');

storm_num = N;  % number of storm events in record
out(:,1) = storm_num;

% index storm arrival rows and calculate storm arrival interval
aa = find(diff(L)>0);
aa = aa+1;
% if the first time step has precipitation it does not get counted with the 2 lines above;
% the following if loop adjusts for that
if rain_obs(1) > 0 % check if first time step is > 0
  aaa = nan(N,1); % preallocate new temp variable with the length of N
  aaa(2:end) = aa; % copy aa to new temp variable, starting at 2nd row
  aaa(1) = 1; % set first time step to 1 (indicating rain on the first time step)
  aa = aaa; % convert back to aa
end
% arrival_dates = rain_date(aa);

bb = diff(aa);
bb = [NaN; bb];
storm_arrival = bb; % arrival times are for the start of one storm relative to the start of the previous storm
out(:,2) = nanmean(storm_arrival);

% index storm end rows 
cc = find(diff(L)<0);
cc = cc+1;
% end_date = rain_date(bb);

% calculate inter-storm period (storm_interval)
% time from end of one storm to beginning of the next storm
% first row = NaN because there is no inter-storm interval for the first
% storm
for j=2:N
  storm_interval(j,1) = aa(j,:) - cc(j-1,:);
end
out(:,3) = nanmean(storm_interval);

% determine length of each storm
for j=1:N
  storm_duration(j,1) = length(rain_obs(L==j));
end     
out(:,4) = nanmean(storm_duration);

% calculate depth of each storm
for j=1:N
  storm_depth(j,1) = nansum(rain_obs(L==j));
end
out(:,5) = nanmean(storm_depth);

% calculate intensity of each storm
storm_intensity = storm_depth./storm_duration;
out(:,6) = nanmean(storm_intensity);

% identify each storm period
storm_period = 
out(