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

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

      DT=365.25/NUM
      W1=2.0*3.14159265*DT/365.25
      DO I=1,N
        CS(1,I)=1.0
        DO J=1,25
          CS(2*J,I)=COS(J*W1*I)
          CS(2*J+1,I)=SIN(J*W1*I)
        ENDDO
      ENDDO

      IC=0
      IF (M1.GT.0) THEN
        DO J=1,M1
          IC=IC+1
          DO I=1,N
            F(I,IC)=CS(J,I)
          ENDDO
        ENDDO
      ENDIF
      IF (M2.GT.0) THEN
        DO J=1,M2
          IC=IC+1
          DO I=1,N
            F(I,IC)=I*CS(J,I)
          ENDDO
        ENDDO
      ENDIF
      IF (M3.GT.0) THEN
        DO J=1,M3
          IC=IC+1
          DO I=1,N
            F(I,IC)=QBO(I)*CS(J,I)
          ENDDO
        ENDDO
      ENDIF
      IF (M4.GT.0) THEN
        DO J=1,M4
          IC=IC+1
          DO I=1,N
            F(I,IC)=SOLAR(I)*CS(J,I)
          ENDDO
        ENDDO
      ENDIF
      IF (M5.GT.0) THEN
        DO J=1,M5
          IC=IC+1
          DO I=1,N
            F(I,IC)=EXTRA1(I)*CS(J,I)
          ENDDO
        ENDDO
      ENDIF
      IF (M6.GT.0) THEN
        DO J=1,M6
          IC=IC+1
          DO I=1,N
            F(I,IC)=EXTRA2(I)*CS(J,I)
          ENDDO
        ENDDO
      ENDIF
      !_/_/_/_ BEGIN _/_/_/_/_/zregr.f ROUTINE_/_/_/_ BEGIN _/_/_/_/_/
      NGOOD=0
      DO I=1,N
        IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX)) NGOOD=NGOOD+1
      ENDDO
      DO J=1,M
        S(J)=0.
        DO I=1,N
          IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX)) THEN
            S(J)=S(J)+T(I)*F(I,J)         ! <== S=F'T
          ENDIF
        ENDDO
      ENDDO
      !Construct matrix A:
      DO K=1,M
        DO J=1,M
          A(J,K)=0.
          DO I=1,N
            IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX)) THEN
              A(J,K)=A(J,K)+F(I,J)*F(I,K) ! <== A=F'F
            ENDIF
          ENDDO
        ENDDO
      ENDDO
      DO I=1,M
        DO J=1,M
          WK(I,J)=A(I,J)
        ENDDO
      ENDDO
      DO I=1,M
        DO J=M+1,2*M
          IF (M+I.EQ.J) THEN
            WK(I,J)=1.0
          ELSE
            WK(I,J)=0.0
          ENDIF
        ENDDO
      ENDDO
      DO K=1,M
        DO I=1,M
          IF (I.NE.K) THEN
            IF (WK(K,K).EQ.0.0) THEN
              L=1
 30           IF (WK(L,K).EQ.0.0) THEN
                L=L+1
                GOTO 30
              ENDIF
              DO J=K,2*M
                WK(K,J)=WK(K,J)+WK(L,J)
              ENDDO
            ENDIF
            U=-WK(I,K)/WK(K,K)
            DO J=K+1,2*M
              WK(I,J)=WK(I,J)+U*WK(K,J)
            ENDDO
          ENDIF
        ENDDO
      ENDDO
      DO J=1,M
        DO I=1,M
          AINV(I,J)=WK(I,J+M)/WK(I,I)
        ENDDO
      ENDDO
      IF (IPRINT.EQ.1) THEN
        DO J=1,M
          WRITE(*,'(1X,A7,I3,A19)')'Column ',J,' of inverse matrix:'
          DO I=1,M
            WRITE(*,*) A(I,J)
          ENDDO
        ENDDO
      ENDIF
      IF (ICHECK.EQ.1) THEN
        WRITE(*,*)'-----------------------------------------------'
        WRITE(*,*)'Simple error check by multiplying A by AINV:'
        WRITE(*,*)'(this should ideally yield the identity matrix)'
        DO J=1,M
          WRITE(*,'(1X,A7,I3,A27)')
     &      'Column ',J,' of product matrix A*AINV:'
          DO I=1,M
            SUM=0.
            DO K=1,M
              SUM=SUM+A(I,K)*AINV(K,J)
            ENDDO
            WRITE(*,*) SUM
          ENDDO
        ENDDO
      ENDIF
      !Coefficient vector solution:
      DO I=1,M
        C(I)=0.
        DO J=1,M
          C(I)=C(I)+AINV(I,J)*S(J)
        ENDDO
      ENDDO
      IF (IERROR.EQ.1) THEN
        DO I=1,N
          IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX)) THEN
            SUM=0
            DO J=1,M
              SUM=SUM+C(J)*F(I,J)
            ENDDO
            RESIDUAL_FIT(I)=T(I)-SUM
          ENDIF
        ENDDO
        SUM=0
        DO I=1,N
          IF ((T(I).GT.TMIN).AND.(T(I).LT.TMAX))
     &      SUM=SUM+RESIDUAL_FIT(I)**2
        ENDDO
        SIGSQR=SUM/(NGOOD-M)
        DO K=1,M
          DO J=1,M
            COV(J,K)=SIGSQR*AINV(J,K)
          ENDDO
        ENDDO
      ENDIF
      !_/_/_/_/_ END _/_/_/_/_/zregr.f ROUTINE_/_/_/_/_ END _/_/_/_/_/
      !Individual seasonal coeffs (ALPHA, BETA,...)
      DO I=1,NUM
        IF (M1.GT.0) THEN
          SUM=0.
          DO J=1,M1
            SUM=SUM+C(J)*CS(J,I)
          ENDDO
          ALPHA(I)=SUM
        ENDIF
        IF (M2.GT.0) THEN
          SUM=0.
          DO J=1,M2
            SUM=SUM+C(J+M1)*CS(J,I)
          ENDDO
          BETA(I)=SUM
        ENDIF
        IF (M3.GT.0) THEN
          SUM=0.
          DO J=1,M3
            SUM=SUM+C(J+M1+M2)*CS(J,I)
          ENDDO
          GAMMA(I)=SUM
        ENDIF
        IF (M4.GT.0) THEN
          SUM=0.
          DO J=1,M4
            SUM=SUM+C(J+M1+M2+M3)*CS(J,I)
          ENDDO
          DELTA(I)=SUM
        ENDIF
        IF (M5.GT.0) THEN
          SUM=0.
          DO J=1,M5
            SUM=SUM+C(J+M1+M2+M3+M4)*CS(J,I)
          ENDDO
          EPS1(I)=SUM
        ENDIF
        IF (M6.GT.0) THEN
          SUM=0.
          DO J=1,M6
            SUM=SUM+C(J+M1+M2+M3+M4+M5)*CS(J,I)
          ENDDO
          EPS2(I)=SUM
        ENDIF
      ENDDO
      IF (IERROR.EQ.1) THEN
        !Seasonal factor (VAR_FACTOR) for coeff variances
        DO I=1,NUM
          SUM=0.
          SUMCT=0.
          DO IYEAR=1,N/NUM
            MM=I+NUM*(IYEAR-1)
            IF ((T(MM).GT.TMIN).AND.(T(MM).LT.TMAX)) THEN
              SUM=SUM+RESIDUAL_FIT(MM)**2
              SUMCT=SUMCT+1.
            ENDIF
          ENDDO
          IF (SUMCT.GT.0.) VAR_FACTOR(I)=SUM/SUMCT
          IF (SUMCT.EQ.0.) VAR_FACTOR(I)=0.
        ENDDO
        SUM=0.
        DO I=1,NUM
          SUM=SUM+VAR_FACTOR(I)
        ENDDO
        SUM=SUM/NUM
        IF (SUM.GT.0.) THEN
          DO I=1,NUM
            VAR_FACTOR(I)=VAR_FACTOR(I)/SUM
          ENDDO
        ENDIF
        IF (M1.GT.0) THEN
          DO I=1,NUM
            SUM=0.
            DO K=1,M1
              DO J=1,M1
                SUM=SUM+COV(J,K)*CS(J,I)*CS(K,I)
              ENDDO
            ENDDO
            VAR_ALPHA(I)=SUM*VAR_FACTOR(I)
          ENDDO
        ENDIF
        IF (M2.GT.0) THEN
          DO I=1,NUM
            SUM=0.
            DO K=1,M2
              DO J=1,M2
                SUM=SUM+COV(J+M1,K+M1)*CS(J,I)*CS(K,I)
              ENDDO
            ENDDO
            VAR_BETA(I)=SUM*VAR_FACTOR(I)
          ENDDO
        ENDIF
        IF (M3.GT.0) THEN
          DO I=1,NUM
            SUM=0.
            DO K=1,M3
              DO J=1,M3
                SUM=SUM+COV(J+M1+M2,K+M1+M2)*CS(J,I)*CS(K,I)
              ENDDO
            ENDDO
            VAR_GAMMA(I)=SUM*VAR_FACTOR(I)
          ENDDO
        ENDIF
        IF (M4.GT.0) THEN
          DO I=1,NUM
            SUM=0.
            DO K=1,M4
              DO J=1,M4
                SUM=SUM+COV(J+M1+M2+M3,K+M1+M2+M3)
     &              *CS(J,I)*CS(K,I)
              ENDDO
            ENDDO
            VAR_DELTA(I)=SUM*VAR_FACTOR(I)
          ENDDO
        ENDIF
        IF (M5.GT.0) THEN
          DO I=1,NUM
            SUM=0.
            DO K=1,M5
              DO J=1,M5
                SUM=SUM+COV(J+M1+M2+M3+M4,K+M1+M2+M3+M4)
     &              *CS(J,I)*CS(K,I)
              ENDDO
            ENDDO
            VAR_EPS1(I)=SUM*VAR_FACTOR(I)
          ENDDO
        ENDIF
        IF (M6.GT.0) THEN
          DO I=1,NUM
            SUM=0.
            DO K=1,M6
              DO J=1,M6
                SUM=SUM+COV(J+M1+M2+M3+M4+M5,K+M1+M2+M3+M4+M5)
     &              *CS(J,I)*CS(K,I)
              ENDDO
            ENDDO
            VAR_EPS2(I)=SUM*VAR_FACTOR(I)
          ENDDO
        ENDIF
      ENDIF
      !Individual time series fits:
      DO I=1,N
        II=MOD(I-1,NUM)+1
        ALPHA_FIT(I)=ALPHA(II)
        BETA_FIT(I)=BETA(II)*I
        GAMMA_FIT(I)=GAMMA(II)*QBO(I)
        DELTA_FIT(I)=DELTA(II)*SOLAR(I)
        EPS1_FIT(I)=EPS1(II)*EXTRA1(I)
        EPS2_FIT(I)=EPS2(II)*EXTRA2(I)
      ENDDO

      RETURN
      END