C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
      PROGRAM EXAMPLE1

C Program: Gives an example of linear regression trend analysis 
C          of a time series using the Fortran procedure "ZTREND.F".
C          The trend model includes seasonal cycle, linear trend,
C          QBO, and SOLAR proxy terms.  (please read documentation
C          in subroutine "ZTREND.F")
C The data:
C          Nimbus 7 TOMS version 7 1979-1992 monthly zonal mean
C          total column ozone time series averaged between 40N and 45N.
C          (Data Source:  see http://toms.gsfc.nasa.gov/)
C To run this example program:
C          1) Make sure the following files are in your current
C             directory:
C                         ztrend.f
C                           data.txt
C                       example1.f    (this program)
C          2) Compile:
C                  f77 -o example1 example1.f ztrend.f
C          3) Run program:
C                  example1

      !Constant parameters necessary to run subroutine ZTREND.F:
      PARAMETER(IERROR=1,N=168,NUM=12,TMIN=50.,TMAX=600.,
     &          M1=9,M2=7,M3=7,M4=7,M5=0,M6=0,M=M1+M2+M3+M4+M5+M6)
      !Note meaning of the constants in the parameter statement above:
      !IERROR=1 ; <== Do multivariate statistical analysis
      !N=168    ; <== Total number of months in time series (1979-92)
      !NUM=12   ; <== Number of data samples per year
      !TMIN=50  ; <== Ozone values < 50 DU are flagged as bad data
      !TMAX=600 ; <== Ozone values > 600 DU are flagged as bad data
      !M1=9     ; <== 9 constants for seasonal coeff
      !M2=7     ; <== 7     "      "  linear trend coefficient
      !M3=7     ; <== 7     "      "  QBO coefficient
      !M4=7     ; <== 7     "      "  SOLAR coefficient
      !M5=0     ; <== 0     "      "  EXTRA1 coeff (EXTRA1 not used)
      !M6=0     ; <== 0     "      "  EXTRA2 coeff (EXTRA2 not used)
      !M=M1+M2+M3+M4+M5+M6  ; <== Total number of constants in model

      !Define arrays for procedure ztrend.pro:
      REAL T(N),QBO(N),SOLAR(N),EXTRA1(N),EXTRA2(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)
      REAL RESIDUAL_FIT(N),VAR_FACTOR(NUM)
      REAL F(N,M),WK(N,2*N),A(M,M),AINV(M,M),S(M),C(M),COV(M,M),
     &     CS(52,N)

      !Define a 12-element array for printing out derived coefficients:
      REAL DUMMY(12)

      !Read in Ozone time series, and QBO and SOLAR proxies:
      CHARACTER*80 STUFF
      OPEN(12,FORM='FORMATTED',STATUS='OLD',FILE='data.txt')
      READ(12,*) STUFF
      DO I=1,N
        READ(12,*) T(I)
      ENDDO
      READ(12,*) STUFF
      DO I=1,N
        READ(12,*) QBO(I)
      ENDDO
      READ(12,*) STUFF
      DO I=1,N
        READ(12,*) SOLAR(I)
      ENDDO
      CLOSE(12)

      !Now do the trend analysis:
      CALL 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)

      !Print out the coefficients and their 2-sigma uncertainties:
      WRITE(*,*)'Seasonal Cycle Coefficient (Dobson Units):'
      WRITE(*,*)
     & ' Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec'
      WRITE(*,'(12F5.0)') ALPHA
      DO I=1,NUM
        DUMMY(I)=2*SQRT(VAR_ALPHA(I))
      ENDDO
      WRITE(*,*)'Seasonal Cycle Coefficient 2-SIGMA (Dobson Units):'
      WRITE(*,*)
     & ' Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec'
      WRITE(*,'(12F5.1)') DUMMY
      WRITE(*,*)' '

      WRITE(*,*)'Trend Coefficient (% per decade):'
      DO I=1,NUM
        DUMMY(I)=BETA(I)*12000./ALPHA(I)
      ENDDO
      WRITE(*,*)
     & ' Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec'
      WRITE(*,'(12F5.1)') DUMMY
      WRITE(*,*)'Trend Coefficient 2-SIGMA (% per decade):'
      DO I=1,NUM
        DUMMY(I)=2*SQRT(VAR_BETA(I))*12000./ALPHA(I)
      ENDDO
      WRITE(*,*)
     & ' Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec'
      WRITE(*,'(12F5.1)') DUMMY
      WRITE(*,*)' '

      WRITE(*,*)'QBO Coefficient (DU per 10m/s):'
      WRITE(*,*)
     & ' Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec'
      DO I=1,NUM
        DUMMY(I)=GAMMA(I)*10.
      ENDDO
      WRITE(*,'(12F5.1)') DUMMY
      WRITE(*,*)'QBO Coefficient 2-SIGMA (DU per 10m/s):'
      DO I=1,NUM
        DUMMY(I)=2*SQRT(VAR_GAMMA(I))*10.
      ENDDO
      WRITE(*,*)
     & ' Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec'
      WRITE(*,'(12F5.1)') DUMMY
      WRITE(*,*)' '

      WRITE(*,*)'SOLAR Coefficient (DU per 100F10.7):'
      WRITE(*,*)
     & ' Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec'
      DO I=1,NUM
        DUMMY(I)=DELTA(I)*100.
      ENDDO
      WRITE(*,'(12F5.1)') DUMMY
      WRITE(*,*)'SOLAR Coefficient 2-SIGMA (DU per 100F10.7):'
      DO I=1,NUM
        DUMMY(I)=2*SQRT(VAR_DELTA(I))*100.
      ENDDO
      WRITE(*,*)
     & ' Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec'
      WRITE(*,'(12F5.1)') DUMMY


      STOP
      END