changing time scale in cross correlation

Hi all, I'm running a 1D cross correlation on two timeseries that are given in years, however I want to show a lag in months not years. I'm not sure where to change the script to do this. The following is my script to create the correlation and lag:

    t1m=time(10:341);           %Excludes ends of time length to
    s1m=ssha_area(10:341);      %better fit data line
    t2m=ctime(10:160);
    s2m=chlo_area(10:160);
    N1=length(t1m);            % sample size
    N2=length(t2m);            %
    N21=ceil(N1/2);            % half the number of obs(rounded)
    N22=ceil(N2/2);            %
    ts1 = t1m(1);              %
    ts2 = t2m(1);              %
    tlast1 = t1m(end);         % last data point
    tlast2 = t2m(end);         %
    ti1 = (tlast1-ts1)/(N1-1); % mean sampling interval
    ti2 = (tlast2-ts2)/(N2-1); % 
    te1 = tlast1+ti1;          % end plus one increment
    te2 = tlast2+ti2;          %
    %
    % frequency vector
    %
    fmin1=0;                   % the min freq is always 0
    fmin2=0;                   %
    fi1=1/(te1-ts1);           % the lowest non-zero freq in data
                               % is always 1 divided by the length 
                               % of the data set
    fi2=1/(te2-ts2);           %
    fs1=1/ti1;                 % the sampling frequency
    fs2=1/ti2;                 %
    fmax1=fs1-fi1;             % before applying the fftshift, the 
                               % maximum frequency is always the
                               % sampling frequency minus the
    fmax2=fs2-fi2;             % frequency increment
    fpre1=(fmin1:fi1:fmax1)';  % frequency grid before fftshift
    fpre2=(fmin2:fi2:fmax2)';  %
    f1=[fpre1(N21+1:N1)-fs1;fpre1(1:N21)];

    f2=[fpre2(N22+1:N2)-fs2;fpre2(1:N22)]; 
    % Now apply the Fourier transform
    %
    S1=fft(s1m,N1);       
    S2=fft(s2m,N2);
    %
    % Then apply fftshift to S1 and S2:
    %
    S1shift=fftshift(S1);
    S2shift=fftshift(S2);
    %
    % and compute the magnitude
    %
    S1mag=fftshift(abs(S1));
    S2mag=fftshift(abs(S2));

    % Now plot the original data.
    %
      figure(1), clf
      subplot(2,1,1), grid on, box on, hold on
      plot(t1m,s1m,'b--.')
      legend('signal 1 unfiltered')
      xlabel('time (years)'), ylabel('signal (m)')
      subplot(2,1,2), grid on, box on, hold on
      plot(f1,S1mag,'b--.')
      legend('signal 1 unfiltered signal')
      xlabel('frequency (years^{-1})'), ylabel('magnitude (m)')
      figure(2), clf
      subplot(2,1,1), grid on, box on, hold on
      plot(t2m,s2m,'b--.')
      legend('signal 2 unfiltered')
      xlabel('time (years)'), ylabel('signal (m)')
      subplot(2,1,2), grid on, box on, hold on
      plot(f2,S2mag,'b--.')
      legend('signal 2 unfiltered')
      xlabel('frequency (hrs^{-1})'), ylabel('amplitude (m)')
    end

This is where the cross correlation begins really

    % Cross correlate the two time-series, and find the dominant frequency at which they are correlated, and the phase difference at that frequency. 

    s1m = s1m-mean(s1m);
    s2m = s2m-mean(s2m);
    s1m=s1m(1:N2);

    [x12pre,lag12pre]=xcorr(s1m,s2m)

Ok I think this is where I would change the script but that would seem to be dependent on my sampling interval(ti1) as well so I'm not sure if there are multiple parts I have to change.

    % Convert lag12pre into values of unit years.
    %
    lag12=((lag12pre*ti1));
    %
    % Normalise the cross-correlation to fall between -1 and 1
    %
    [a11,lag11]=xcorr(s1m);
    [a22,lag22]=xcorr(s2m);
    lagind11=find(lag11==0);
    lagind22=find(lag22==0);
    x12=x12pre/sqrt(a11(lagind11).*a22(lagind22));
    %
    % Calculate the cross-spectrum by Fourier transforming the
    % cross-correlation. 
    %
    % (ifftshift).
    %
    x12shift=ifftshift(x12);
    lag12shift=lag12-min(lag12);
    %
    % Create a new frequency vector, corresponding to the
    % time-lag axis of the cross-correlation of the signals.
    %
    N12=length(lag12);         % number of observations
    NT12=lag12(end)-lag12(1); % length of record (years)
    dT12=NT12/(N12-1);         % mean sample space (years).
    fs12=1/dT12;               % sampling frequency (max frequency)
    df12=(1/NT12)*(N12-1)/N12; % minimum non-zero frequency.
    %
    f12pre=0:df12:fs12-df12;   % the frequency vector
    f12=fftshift(f12pre);      % and the vector that corresponds
                               % to the fftshifted signal
    f12(f12>(fs12-df12/2)/2)=f12(f12>(fs12-df12/2)/2)-fs12;
    %
    % Calculate the Fourier transform, and fftshift it:
    %
    X12=fft(x12shift);
    X12shift=fftshift(X12);
    %
    % Normalise the magnitude by the length of the record.
    %
    X12mag=abs(X12shift);
    %
    % Calculate the angle
    X12phs=angle(X12shift);
    %
    if 1
       figure(3), clf
       subplot(3,1,1), grid on, box on, hold on
       title('Lag Between SSHA and Chlorophyll')
       plot(lag12,x12,'b--.')
       xlabel('lag (years)'), ylabel('correlation (m^2)')
       subplot(3,1,2), grid on, box on, hold on
       plot(f12,X12mag,'b--.')
       xlabel('frequency (years^{-1})'), ylabel('magnitude of    xcorr (m^2)')
       subplot(3,1,3), grid on, box on, hold on
       plot(f12,X12phs,'b--.')
       xlabel('frequency (years^{-1})'), ylabel('lag (rads)')
    end

[Daniel Shub]

You need to work through your problem enough that you can give a much smaller MWE. Don't give us anything that we don't need ... For example, the question doesn't see to deal with plotting, so get rid of all of that.