pro ZTREND, $
        IERROR,N,M,NUM,T,TMIN,TMAX, $
        M1,M2,M3,M4,M5,M6, $
        QBO,SOLAR,EXTRA1,EXTRA2, $
        ALPHA,BETA,GAMMA,DELTA,EPS1,EPS2, $
        VAR_ALPHA,VAR_BETA,VAR_GAMMA,VAR_DELTA,VAR_EPS1,VAR_EPS2, $
        ALPHA_FIT,BETA_FIT,GAMMA_FIT,DELTA_FIT,EPS1_FIT,EPS2_FIT, $
        RESIDUAL_FIT,VAR_FACTOR, $
        F,WK,A,AINV,S,C,COV,CS
;+
; PROCEDURE: Performs multiple linear regression trend analysis of
; an arbitrary time series.  OPTIONAL: error analysis for regression
; coefficients (uses standard multivariate noise model).
;
; _/_/_/_/_/_/_/_/_/_/ S H O R T   T U T O R I A L _/_/_/_/_/_/_/_/_/_/
; Form of general regression trend model used in this procedure
; (t = time index = 0,1,2,3,...,N-1):
;
;            T(t)=ALPHA(t)
;                    +BETA(t)*t
;                        +GAMMA(t)*QBO(t)
;                            +DELTA(t)*SOLAR(t)
;                                +EPS1(t)*EXTRA1(t)
;                                    +EPS2(t)*EXTRA2(t)
;                                        +RESIDUAL_FIT(t),
;
; where ALPHA represents the 12-month seasonal fit, BETA is the 
; 12-month seasonal trend coefficient, RESIDUAL_FIT(t) represents the 
; error time series, and GAMMA, DELTA, EPS1, and EPS2 are 12-month 
; coefficients corresponding to the ozone driving quantities QBO 
; (quasi-biennial oscillation), SOLAR (solar-UV proxy), and proxies
; EXTRA1 and EXTRA2 (for example, these latter two might be ENSO, 
; vorticity, geopotential heights, or temperature), respectively.
; The general model above assumes simple linear relationships 
; between T(t) and surrogates which is hopefully valid as a
; first approximation.  Note that for total ozone trends based on 
; chemical species such as involving Chlorine, the trend term 
; BETA(t)*t could be replaced (ignored by setting m2=0 in the 
; procedure call), with EPS1(t)*EXTRA1(t) where EXTRA1(t) is 
; the chemical proxy time series.
;
; This procedure assumes the following form for the coefficients
; ALPHA, BETA, GAMMA,...) in effort to approximate realistic
; seasonal dependence of sensitivity between T(t) and surrogate
; [The expansion shown below is for ALPHA(t) - similar expansions
; for BETA(t), GAMMA(t), DELTA(t), EPS1(t), and EPS2(t)]:
;
; ALPHA(t) = A0                                           <== Constant
;  +A1*cos(1*2*PI*t*dt/365.25)+A2*sin(1*2*PI*t*dt/365.25  <== 12-month
;  +A3*cos(2*2*PI*t*dt/365.25)+A4*sin(2*2*PI*t*dt/365.25) <== 6-month
;  +A5*cos(3*2*PI*t*dt/365.25)+A6*sin(3*2*pi*t*dt/365.25) <== 4-month
;  +          .                      .                  .
;  +          .                      .                  .
;  +          .                      .                  .
;
; where A0, A1, A2,... are constants, and t (t=1,2,...,n) is the
; time index (NO UNITS).
;
; Trend models often use a particular harmonic expansion to
; represent the seasonalility of the action between T(t) and a
; particular surrogate.  Shown below are two EXAMPLES of such
; chosen harmonic expansions for total ozone trend analyses [T(t)
; is a total ozone time series]:
;
;   1) Stolarski, Bloomfield, McPeters [1991] model:
;
;     ALPHA(t): A0, A1, A2, A3, A4, A5, A6, A7, A8 (const+12,6,4,3mo)
;     BETA(t) : A0, A1, A2, A3, A4 (const+12mo+6mo)
;     GAMMA(t): A0, A1, A2 (const+12mo)
;     DELTA(t): A0, A1, A2 (  "    "  )
;              (no proxies other than QBO and SOLAR)
;
;   2) Randel and Cobb [1994] zonally asymmetric model:
;
;     ALPHA(t): A0, A1, A2, A3, A4, A5, A6 (const+12,6,4mo)
;     BETA(t) : A0, A1, A2, A3, A4, A5, A6 (       "      )
;     GAMMA(t): A0, A1, A2, A3, A4, A5, A6 (       "      )
;     DELTA(t): A0, A1, A2, A3, A4, A5, A6 (       "      )
;     EPS1(t):  A0, A1, A2, A3, A4, A5, A6 (       "      )
;      (no proxies other than QBO, SOLAR, and an "EXTRA" proxy ENSO)
;
; _/_/_/_/_/_/ P A R A M E T E R   I N F O R M A T I O N _/_/_/_/_/_/_/
;
; INPUT: (1st three lines of parameter list - all MUST be declared
;         prior to calling ZTREND.PRO)
;
;       IERROR <== If preset to 1, standard error analysis is
;                  performed for regression coefficients.  Any
;                  other value SKIPS it
;       N <== Length (dimension) of array "T" to be modelled
;       M <== Total number of regression constants in the model
;             (this is EQUIVALENT to M1+M2+M3+...+M6 (see below)
;       NUM <== Number of time series measurements PER YEAR. NUM is
;               used to introduce seasonal dependence of the 
;               coefficients ALPHA(t), BETA(t), ..., EPS2(t)
;       T <== The arbitrary time series for which regression trend
;             analysis is performed
;       TMIN,TMAX <== Values of input series T that are less than TMAX
;                     or greater than TMAX will NOT be included in
;                     the regression (these are the flagging values)
;       M1, M2, M3, M4, M5, M6 <== Number (1,3,5,7,...) of constants 
;                                  to include in seasonal, trend, qbo
;                                  solar, extra1, and extra2 proxy
;                                  coefficients, respectively. These
;                                  must all be ODD numbers if not
;                                  ZERO. When you set any of M1,
;                                  M2,...,M6 to ZERO, this will
;                                  IGNORE the respective surrogate
;                                  during regression
;       QBO <== Quasi-biennial oscillation time series proxy
;       SOLAR <== Solar time series proxy
;       EXTRA1 <== An extra time series proxy (e.g., ENSO)
;       EXTRA2 <== Another extra time series proxy (e.g., vorticity)
;       
; OUTPUT (lines 4-8 in parameter list)
;
;       ALPHA <== Seasonal coefficient (see tutorial above)
;       BETA  <== Seasonal trend coefficient in units
;                 reciprocal to sampling interval.  Example: If T(t)
;                 is total ozone in Dobson units and the number of
;                 samples per year "NUM" is 52, then trends (i.e.,
;                 BETA(t)) will be in Dobson units per week
;       GAMMA <== Seasonal qbo coefficient
;       DELTA <== Seasonal solar coefficient
;       EPS1 <== Seasonal coefficient for 1st extra proxy
;       EPS2 <== Seasonal coefficient for 2n2 extra proxy
;
;       VAR_ALPHA, VAR_BETA, VAR_GAMMA, <== Seasonal variances
;       VAR_DELTA, VAR_EPS1, VAR_EPS2   <== of the coefficients
;
;       ALPHA_FIT    <== Seasonal fit (same length "N" as array "T(t)")
;       BETA_FIT     <== Trend fit                     "
;       GAMMA_FIT    <== QBO fit                       "
;       DELTA_FIT    <== Solar fit                     "
;       EPS1_FIT     <== 1st extra proxy fit           "
;       EPS2_FIT     <== 2nd extra proxy fit           "
;       RESIDUAL_FIT <== Residual fit                  "
;       VAR_FACTOR   <== Seasonal variability factor for coefficient
;                        variance series.  This seasonal variability
;                        IS INVOKED AUTOMATICALLY in this program
;                        as follows:
;
;                        VAR_ALPHA(*)=VAR_ALPHA(*)*VAR_FACTOR(*)
;                        VAR_BETA(*)=VAR_BETA(*)*VAR_FACTOR(*)
;                        VAR_GAMMA(*)=VAR_GAMMA(*)*VAR_FACTOR(*)
;                        VAR_DELTA(*)=VAR_DELTA(*)*VAR_FACTOR(*)
;                        VAR_EPS1(*)=VAR_EPS1(*)*VAR_FACTOR(*)
;                        VAR_EPS2(*)=VAR_EPS2(*)*VAR_FACTOR(*)
;
;       F,WK,A,AINV,S,C,COV,CS <== Internal arrays to this procedure
;                                  (see "zregr.pro" for meaning)
;
; IMPORTANT NOTES:
;
;  1) The time series length N must be greater than M, otherwise
;     you are solving for more constants in the regression than
;     the number of time series values.  Also, the multivariate
;     statistical analysis in this procedure is valid only if N>M.
;
;  2) EACH of M1, M2, ..., M6 must be LESS THAN the number of data
;     samples per year NUM.  Otherwise the linear system becomes
;     singular.
;
;  3) The above linear model can be written as Y = XB + e, where
;     Y is the time series vector (T(t) above) being modeled, 
;     X is a matrix of data values, and B is the coefficient
;     vector to be derived.  Solving for B via least squares
;     yields (prime denotes transpose):
;
;                                  -1
;                          B = (X'X) X'Y
;
;  4) This procedure uses a standard multivariate statistical
;     model that yields
;                                                2      -1
;        covariance matrix for B = cov[B] = sigma * (X'X)
;
;                              2
;     where data variance sigma is estimated from the vector
;     dot product e*e divided by N-M.
;
;  This routine also internally invokes a more realistic
;  seasonality of variances for the time-dependent regression
;  coefficients (i.e., VAR_ALPHA, VAR_BETA, ...)  by modulating
;  their values by the annual climatology of e*e/(N-M) (see plot
;  below):
;
;    Annual    |
;  Climatology |
;      of      |         *
;     e*e      |
;    -----     |
;    N - M     |             *
;              |     *
;              |
;              | *                                           *
;              |                 *                       *
;              |
;              |                     *               *
;              |                         *       *
;              |                             *
;             -+---+---+---+---+---+---+---+---+---+---+---+---+
;                J   F   M   A   M   J   J   A   S   O   N   D
;  For example,
;
;            VAR_ALPHA(*)=VAR_ALPHA(*)*VAR_FACTOR(*)
;
;  where VAR_FACTOR is equivalent to the above annual climatology
;  function for e*e/(N-M) divided by it's annual mean.
;
;  Physical reason for using this modulation approach:
;
;  Values for e*e/(N-M) are generally highly seasonally-dependent.
;  For example, taking annual samplings such as sampling only
;  January monthly means will produce different estimates of
;  e*e/(N-M) as opposed to annual sampling using July monthly means.
;  Several papers in the past on ozone trends have used annual
;  sampling.  A caveat with this approach is that annual sampling
;  produces aliasing effects where all variabilities with periods
;  shorter than two years (Nyquist folding period) will not be
;  detected properly.
;
;  The modulation factor VAR_FACTOR provides an extra step that
;  goes beyond standard multivariate statistics.  After calling
;  "ztrend.pro" if you want you can once again regain the original
;  multivariate coefficient variances by dividing by VAR_FACTOR:
;
;             VAR_ALPHA(*)=VAR_ALPHA(*)/VAR_FACTOR(*)
;             VAR_BETA(*)=VAR_BETA(*)/VAR_FACTOR(*)
;             VAR_GAMMA(*)=VAR_GAMMA(*)/VAR_FACTOR(*)
;             VAR_DELTA(*)=VAR_DELTA(*)/VAR_FACTOR(*)
;             VAR_EPS1(*)=VAR_EPS1(*)/VAR_FACTOR(*)
;             VAR_EPS2(*)=VAR_EPS2(*)/VAR_FACTOR(*)
;
; Address questions/comments to:
;   Dr. Jerry Ziemke
;   NASA/Goddard Space Flight Center            PH: (301) 614-6034
;   Code 916                                   FAX: (301) 614-5903
;   Greenbelt, MD  20771
;   (Affil: UMBC GEST, Baltimore, MD)
;   e-mail: ziemke@jwocky.gsfc.nasa.gov

      ICHECK=0 & IPRINT=0
      F=FLTARR(N,M)
      WK=FLTARR(N,2*N) & A=FLTARR(M,M)
      AINV=FLTARR(M,M) & S=FLTARR(M) & C=FLTARR(M)
      COV=FLTARR(M,M) & CS=FLTARR(52,N)
      ALPHA=FLTARR(NUM) & BETA=FLTARR(NUM) & GAMMA=FLTARR(NUM)
      DELTA=FLTARR(NUM) & EPS1=FLTARR(NUM) & EPS2=FLTARR(NUM)
      VAR_ALPHA=FLTARR(NUM) & VAR_BETA=FLTARR(NUM)
      VAR_GAMMA=FLTARR(NUM) & VAR_DELTA=FLTARR(NUM)
      VAR_EPS1=FLTARR(NUM) & VAR_EPS2=FLTARR(NUM)
      ALPHA_FIT=FLTARR(N) & BETA_FIT=FLTARR(N) & GAMMA_FIT=FLTARR(N)
      DELTA_FIT=FLTARR(N) & EPS1_FIT=FLTARR(N) & EPS2_FIT=FLTARR(N)
      RESIDUAL_FIT=FLTARR(N)
      VAR_FACTOR=FLTARR(NUM)


      CT=0
      FOR I=1,N DO BEGIN
        IF (T(I-1) GT TMIN) AND (T(I-1) LT TMAX) THEN BEGIN
          CT=CT+1
        ENDIF
      ENDFOR
      IF CT EQ 0 THEN BEGIN
        PRINT,'** All data either missing or bad - SKIPPING routine...'
        GOTO,JMP2
      ENDIF
      DT=365.25/NUM
      W1=2.0*3.14159265*DT/365.25
      FOR I=1,N DO BEGIN
        CS(0,I-1)=1.0
        FOR J=1,25 DO BEGIN
          CS(2*J-1,I-1)=COS(J*I*W1)
          CS(2*J,I-1)=SIN(J*I*W1)
        ENDFOR
      ENDFOR
      IC=0
      IF M1 GT 0 THEN BEGIN
        FOR J=1,M1 DO BEGIN
          IC=IC+1
          FOR I=1,N DO BEGIN
            F(I-1,IC-1)=CS(J-1,I-1)
          ENDFOR
        ENDFOR
      ENDIF
      IF M2 GT 0 THEN BEGIN
        FOR J=1,M2 DO BEGIN
          IC=IC+1
          FOR I=1,N DO BEGIN
            F(I-1,IC-1)=I*CS(J-1,I-1)
          ENDFOR
        ENDFOR
      ENDIF
      IF M3 GT 0 THEN BEGIN
        FOR J=1,M3 DO BEGIN
          IC=IC+1
          FOR I=1,N DO BEGIN
            F(I-1,IC-1)=QBO(I-1)*CS(J-1,I-1)
          ENDFOR
        ENDFOR
      ENDIF
      IF M4 GT 0 THEN BEGIN
        FOR J=1,M4 DO BEGIN
          IC=IC+1
          FOR I=1,N DO BEGIN
            F(I-1,IC-1)=SOLAR(I-1)*CS(J-1,I-1)
          ENDFOR
        ENDFOR
      ENDIF
      IF M5 GT 0 THEN BEGIN
        FOR J=1,M5 DO BEGIN
          IC=IC+1
          FOR I=1,N DO BEGIN
            F(I-1,IC-1)=EXTRA1(I-1)*CS(J-1,I-1)
          ENDFOR
        ENDFOR
      ENDIF
      IF M6 GT 0 THEN BEGIN
        FOR J=1,M6 DO BEGIN
          IC=IC+1
          FOR I=1,N DO BEGIN
            F(I-1,IC-1)=EXTRA2(I-1)*CS(J-1,I-1)
          ENDFOR
        ENDFOR
      ENDIF
      ;_/_/_/_ BEGIN _/_/_/_/_/zregr.pro ROUTINE_/_/_/_ BEGIN _/_/_/_/_/
      NGOOD=0
      FOR I=1,N DO BEGIN
        IF (T(I-1) GT TMIN) AND (T(I-1) LT TMAX) THEN NGOOD=NGOOD+1
      ENDFOR
      FOR J=1,M DO BEGIN
        S(J-1)=0.
        FOR I=1,N DO BEGIN
          IF (T(I-1) GT TMIN) AND (T(I-1) LT TMAX) THEN BEGIN
            S(J-1)=S(J-1)+T(I-1)*F(I-1,J-1)         ; <== S=F'T
          ENDIF
        ENDFOR
      ENDFOR
      ;Construct matrix A:
      FOR K=1,M DO BEGIN
        FOR J=1,M DO BEGIN
          A(J-1,K-1)=0.
          FOR I=1,N DO BEGIN
            IF (T(I-1) GT TMIN) AND (T(I-1) LT TMAX) THEN BEGIN
              A(J-1,K-1)=A(J-1,K-1)+F(I-1,J-1)*F(I-1,K-1) ; <== A=F'F
            ENDIF
          ENDFOR
        ENDFOR
      ENDFOR
;      AINV=INVERT(A)
      FOR I=1,M DO BEGIN
        FOR J=1,M DO BEGIN
          WK(I-1,J-1)=A(I-1,J-1)
        ENDFOR
      ENDFOR
      FOR I=1,M DO BEGIN
        FOR J=M+1,2*M DO BEGIN
          IF M+I EQ J THEN BEGIN
            WK(I-1,J-1)=1.
          ENDIF ELSE BEGIN
            WK(I-1,J-1)=0.
          ENDELSE
        ENDFOR
      ENDFOR
      FOR K=1,M DO BEGIN
        FOR I=1,M DO BEGIN
          IF I NE K THEN BEGIN
            IF WK(K-1,K-1) EQ 0 THEN BEGIN
              L=1
JMP:          IF WK(L-1,K-1) EQ 0 THEN BEGIN
                L=L+1
                GOTO, JMP
              ENDIF
              FOR J=K,2*M DO BEGIN
                WK(K-1,J-1)=WK(K-1,J-1)+WK(L-1,J-1)
              ENDFOR
            ENDIF
            U=-WK(I-1,K-1)/WK(K-1,K-1)
            FOR J=K+1,2*M DO BEGIN
              WK(I-1,J-1)=WK(I-1,J-1)+U*WK(K-1,J-1)
            ENDFOR
          ENDIF
        ENDFOR
      ENDFOR
      FOR J=1,M DO BEGIN
        FOR I=1,M DO BEGIN
          AINV(I-1,J-1)=WK(I-1,J+M-1)/WK(I-1,I-1)
        ENDFOR
      ENDFOR
      IF  IPRINT EQ 1 THEN BEGIN
        FOR J=1,M DO BEGIN
          PRINT,'Column ',J,' of inverse matrix:',FORMAT='(1X,A7,I3,A19)'
          FOR I=1,M DO BEGIN
            PRINT, A(I-1,J-1)
          ENDFOR
        ENDFOR
      ENDIF
      IF ICHECK EQ 1 THEN begin
        PRINT,'-----------------------------------------------'
        PRINT,'Simple error check by multiplying A by AINV:'
        PRINT,'(this should ideally yield the identity matrix)'
        FOR J=1,M DO BEGIN
          PRINT,'Column ',J,' of product matrix A*AINV:', $
                FORMAT='(1X,A7,I3,A27)'
          FOR I=1,M DO BEGIN
            SUM=0.
            FOR K=1,M DO BEGIN
              SUM=SUM+A(I-1,K-1)*AINV(K-1,J-1)
            ENDFOR
            PRINT, SUM
          ENDFOR
        ENDFOR
      ENDIF
      ;Coefficient vector solution:
      FOR I=1,M DO BEGIN
        C(I-1)=0.
        FOR J=1,M DO BEGIN
          C(I-1)=C(I-1)+AINV(I-1,J-1)*S(J-1)
        ENDFOR
      ENDFOR
      IF  IERROR EQ 1 THEN BEGIN
        FOR I=1,N DO BEGIN
          IF (T(I-1) GT TMIN) AND (T(I-1) LT TMAX) THEN BEGIN
            SUM=0
            FOR J=1,M DO BEGIN
              SUM=SUM+C(J-1)*F(I-1,J-1)
            ENDFOR
            RESIDUAL_FIT(I-1)=T(I-1)-SUM
          ENDIF
        ENDFOR
        SUM=0
        FOR I=1,N DO BEGIN
          IF (T(I-1) GT TMIN) AND (T(I-1) LT TMAX) THEN BEGIN
            SUM=SUM+RESIDUAL_FIT(I-1)^2
          ENDIF
        ENDFOR
        SIGSQR=SUM/(NGOOD-M)
        FOR K=1,M DO BEGIN
          FOR J=1,M DO BEGIN
            COV(J-1,K-1)=SIGSQR*AINV(J-1,K-1)
          ENDFOR
        ENDFOR
      ENDIF
      ;_/_/_/_/_ END _/_/_/_/_/zregr.pro ROUTINE_/_/_/_/_ END _/_/_/_/_/
      ;Individual seasonal coeffs (ALPHA, BETA,...)
      FOR I=1,NUM DO BEGIN
        IF M1 GT 0 THEN BEGIN
          SUM=0.
          FOR J=1,M1 DO BEGIN
            SUM=SUM+C(J-1)*CS(J-1,I-1)
          ENDFOR
          ALPHA(I-1)=SUM
        ENDIF
        IF M2 GT 0 THEN BEGIN
          SUM=0.
          FOR J=1,M2 DO BEGIN
            SUM=SUM+C(J+M1-1)*CS(J-1,I-1)
          ENDFOR
          BETA(I-1)=SUM
        ENDIF
        IF M3 GT 0 THEN BEGIN
          SUM=0.
          FOR J=1,M3 DO BEGIN
            SUM=SUM+C(J+M1+M2-1)*CS(J-1,I-1)
          ENDFOR
          GAMMA(I-1)=SUM
        ENDIF
        IF M4 GT 0 THEN BEGIN
          SUM=0.
          FOR J=1,M4 DO BEGIN
            SUM=SUM+C(J+M1+M2+M3-1)*CS(J-1,I-1)
          ENDFOR
          DELTA(I-1)=SUM
        ENDIF
        IF M5 GT 0 THEN BEGIN
          SUM=0.
          FOR J=1,M5 DO BEGIN
            SUM=SUM+C(J+M1+M2+M3+M4-1)*CS(J-1,I-1)
          ENDFOR
          EPS1(I-1)=SUM
        ENDIF
        IF M6 GT 0 THEN BEGIN
          SUM=0.
          FOR J=1,M6 DO BEGIN
            SUM=SUM+C(J+M1+M2+M3+M4+M5-1)*CS(J-1,I-1)
          ENDFOR
          EPS2(I-1)=SUM
        ENDIF
      ENDFOR
      IF IERROR EQ 1 THEN BEGIN
        ;Seasonal factor (VAR_FACTOR) for coeff variances
        FOR I=1,NUM DO BEGIN
          SUM=0.
          SUMCT=0.
          FOR IYEAR=1,N/NUM DO BEGIN
            MM=I+NUM*(IYEAR-1)
            IF (T(MM-1) GT TMIN) AND (T(MM-1) LT TMAX) THEN BEGIN
              SUM=SUM+RESIDUAL_FIT(MM-1)^2
              SUMCT=SUMCT+1.
            ENDIF
          ENDFOR
          IF (SUMCT GT 0.) THEN VAR_FACTOR(I-1)=SUM/SUMCT
          IF (SUMCT EQ 0.) THEN VAR_FACTOR(I-1)=0.
        ENDFOR
        SUM=0.
        FOR I=1,NUM DO SUM=SUM+VAR_FACTOR(I-1)
        SUM=SUM/NUM
        IF SUM GT 0 THEN BEGIN
          FOR I=1,NUM DO VAR_FACTOR(I-1)=VAR_FACTOR(I-1)/SUM
        ENDIF
        IF M1 GT 0 THEN BEGIN
          FOR I=1,NUM DO BEGIN
            SUM=0.
            FOR K=1,M1 DO BEGIN
              FOR J=1,M1 DO BEGIN
                SUM=SUM+COV(J-1,K-1)*CS(J-1,I-1)*CS(K-1,I-1)
              ENDFOR
            ENDFOR
            VAR_ALPHA(I-1)=SUM*VAR_FACTOR(I-1)
          ENDFOR
        ENDIF
        IF M2 GT 0 THEN BEGIN
          FOR I=1,NUM DO BEGIN
            SUM=0.
            FOR K=1,M2 DO BEGIN
              FOR J=1,M2 DO BEGIN
                SUM=SUM+COV(J+M1-1,K+M1-1)*CS(J-1,I-1)*CS(K-1,I-1)
              ENDFOR
            ENDFOR
            VAR_BETA(I-1)=SUM*VAR_FACTOR(I-1)
          ENDFOR
        ENDIF
        IF M3 GT 0 THEN BEGIN
          FOR I=1,NUM DO BEGIN
            SUM=0.
            FOR K=1,M3 DO BEGIN
              FOR J=1,M3 DO BEGIN
                SUM=SUM+COV(J+M1+M2-1,K+M1+M2-1)*CS(J-1,I-1)*CS(K-1,I-1)
              ENDFOR
            ENDFOR
            VAR_GAMMA(I-1)=SUM*VAR_FACTOR(I-1)
          ENDFOR
        ENDIF
        IF M4 GT 0 THEN BEGIN
          FOR I=1,NUM DO BEGIN
            SUM=0.
            FOR K=1,M4 DO BEGIN
              FOR J=1,M4 DO BEGIN
                SUM=SUM+COV(J+M1+M2+M3-1,K+M1+M2+M3-1) $
                   *CS(J-1,I-1)*CS(K-1,I-1)
              ENDFOR
            ENDFOR
            VAR_DELTA(I-1)=SUM*VAR_FACTOR(I-1)
          ENDFOR
        ENDIF
        IF M5 GT 0 THEN BEGIN
          FOR I=1,NUM DO BEGIN
            SUM=0.
            FOR K=1,M5 DO BEGIN
              FOR J=1,M5 DO BEGIN
                SUM=SUM+COV(J+M1+M2+M3+M4-1,K+M1+M2+M3+M4-1) $
                   *CS(J-1,I-1)*CS(K-1,I-1)
              ENDFOR
            ENDFOR
            VAR_EPS1(I-1)=SUM*VAR_FACTOR(I-1)
          ENDFOR
        ENDIF
        IF M6 GT 0 THEN BEGIN
          FOR I=1,NUM DO BEGIN
            SUM=0.
            FOR K=1,M6 DO BEGIN
              FOR J=1,M6 DO BEGIN
                SUM=SUM+COV(J+M1+M2+M3+M4+M5-1,K+M1+M2+M3+M4+M5-1) $
                   *CS(J-1,I-1)*CS(K-1,I-1)
              ENDFOR
            ENDFOR
            VAR_EPS2(I-1)=SUM*VAR_FACTOR(I-1)
          ENDFOR
        ENDIF
      ENDIF
      ;Individual time series fits:
      FOR I=1,N DO BEGIN
        II=(I-1) MOD NUM + 1
        ALPHA_FIT(I-1)=ALPHA(II-1)
        BETA_FIT(I-1)=BETA(II-1)*I
        GAMMA_FIT(I-1)=GAMMA(II-1)*QBO(I-1)
        DELTA_FIT(I-1)=DELTA(II-1)*SOLAR(I-1)
        EPS1_FIT(I-1)=EPS1(II-1)*EXTRA1(I-1)
        EPS2_FIT(I-1)=EPS2(II-1)*EXTRA2(I-1)
      ENDFOR

JMP2: RETURN
      END