program multiproxy
c
c     Copyright (C) 1996, revised 1997 (some additional comments added 2005)
c
c     Michael E. Mann
c
c     multivariate proxy based instrumental surface temperature pattern
c     reconstruction
c
c     uses the following approach:
c
c     (1) uses monthly instrumental surface temperature "training" data
c         consisting of all quasi-continueous 5 deg x 5 deg gridpoint
c         records between 1902 and 1993
c
c     (2) normalize the above by variance of the detrended
c         series--this insures that non-stationarity of recent
c         instrumental period will not influence the estimate of
c         variance upon renormalization of reconstructed data.
c         variance is calculated for annual/season average selected
c         for reconstruction.
c
c     (3) perform EOF decomposition of the above into M (typically 
c         about 40) dominant modes that will potentially be used.
c
c     (4) form the annual or seasonal (warm/cold) averages of 
c         the associated principal  components over calibration or 'training'
c         interval (1902-19XX where XX=1980 in most analyses--both ends of the 
c         interval can be varied, e.g. XX=1971.
c
c         similarly, standardize each of the p proxy indicators
c         by their detrended variance over the calibration interval
c
c     (5) find transfer function between each of the p proxy indicators and
c         the first n < p PCs
c
c     (6) invert the transfer function to find the optimal weights on the
c         n PCs given the values of the p proxy indicators during the
c         preinstrumental interval
c
c         Use these weights to calculate the associated surface temperature 
c         patterns associated distribution of proxy indicator values. Restore 
c         the variance removed in step (2) to the proxy-reconstructed data.
c
c     (7) perform 2 types of verification:
c
c         b) verification using subset of quasi-continuous gridpoints available 
c            back to 1854 
c         c) verification against isolated long instrumental records, and Quinn 
c            El Nino chronology
c
c     output provided:
c
c     instrumental data EOF decomposition:
c
c     * preference/option selection information
c     * instrumental surface temperature eigenvalue spectrum
c     * first 10 teofs of instrumental surface temp
c     * first 10 pcs of seasonally/annually averaged pcs
c
c     proxy/instrumental pc training:
c
c     * first 10 ceofs/pcs of proxy/intr pc training
c     * spatial patterns corresponding to the combination of
c       teofs represented by each ceof
c      
c     proxy reconstruction
c
c     * proxy reconstructed data
c     * loadings on each of the pcs for the reconstructions
c
c     indices (raw, training, reconstructed, verification)
c
c     * global ave
c     * nhem ave
c     * nino3
c
c     calibration/verfication statistics:
c
c     * global ave
c     * nhem ave
c     * nino3
c     * multivariate
c     
c
c     ipc represents an upper limit on how many EOFs we
c     compute for potential use in the calibration procedure
c     typically we set to 40
c
      parameter (ipcmax=40)
      parameter (itermax=1000)
c
      parameter (maxmonth=1680)
      parameter (maxstat=1082)
      parameter (maxlong=240)
      parameter (maxproxy=450)
      parameter (mmax = maxstat)
      parameter (max0 = mmax+20)
      parameter (nlarge=8000)
      parameter (nmax1 = 1104)
      parameter (nmax0=1500)
      parameter (nmax2=150)
      parameter (ntrnmax = 80)
      parameter (pi=3.14159265359d0)
      parameter (half=0.5d0)
      parameter (outofbounds=-999.d0)
c
      real alat1,alat2,alon1,alon2
      real lon1,lon2,lat1,lat2
      real *8 lon(max0),lat(max0)
      integer imask(max0),imatch(max0)
      integer iwhich(ipcmax)
      real *8 newlon(max0),newlat(max0)
      real *8 quinn(nmax0),qrankraw(nmax0),
     $         qrankcal(nmax0),
     $         quinnsort(nmax0),
     $            pquinnraw(nmax0),
     $            pquinncal(nmax0),
     $            quinnthresh(10)
      real *8 pquinnrand(nmax0)
      real *8 corrrand(itermax)
      integer iseq(nmax0)
      real *8 series1(nmax0),series2(nmax0),series3(nmax0)
      integer nkeep(max0)
      real *8 cs(max0),cslat(max0),area,areanh
      real *8 sigma(max0),sigma0(max0)
      real *8 aprox(nmax0,maxproxy),ahist(nmax0,maxproxy)
      real *8 alath(maxproxy),alonh(maxproxy)
      real *8 anom(maxmonth,max0),standard(maxmonth,max0)
      real *8 nhemr(nmax2),globr(nmax2),ninor(nmax2),detrr(nmax2)
      real *8 ninosort(nmax2),ninothresh(10)
      real *8 nhemt(nmax2),globt(nmax2),ninot(nmax2),detrt(nmax2)
      real *8 nhemc(nmax0),globc(nmax0),ninoc(nmax0),detrc(nmax0)
      real *8 soir(nmax0)
      real *8 ninod(nmax0)
      real *8 nhemv(nmax0),globv(nmax0),soi(nmax0),detrv(nmax0)
      real *8 enso1(nmax0),enso2(nmax0)
      real *8 rpckeep(max0,ipcmax)
      real *8 pccalib(ipcmax)
      real *8 varrawgp(mmax),varrawgp0(mmax),amsetrngp(mmax),
     $          amsecalgp(mmax),amsevergp(mmax)
      real *8 varcalgp(mmax),varcalpc(ipcmax)
      real *8 corrgp(mmax),corrpc(ipcmax)
      real *8 varrawpc(ipcmax),amsecalpc(ipcmax)
      real *8 annual(nmax2,max0),anew(nmax2,max0)
      real *8 anom0(nmax0,maxlong)
      real *8 a(nmax2,mmax),weight(mmax),weight0(max0)
      real *8 weightprx(maxproxy)
      real *8 aint(nmax0,max0)
      real *8 aver0(nmax0,max0)
      real *8 sd(max0),datave(max0)
      real *8 proxave(maxproxy),sdprox(maxproxy)
      real *8 mean(max0),annmean(max0)
      real *8 S(mmax),S0(mmax),partsum(mmax+1),sum(mmax+1)
      real *8 SS(ipcmax)
      integer index0(max0)
      integer istart(maxproxy),ikeep(maxproxy)
      integer igood(max0)
      real *8 eof(mmax,ipcmax),tpc(nmax0,ipcmax)
      real *8 rpc(ipcmax)
      complex *16 AA(nmax1,mmax),UU(nmax1,nmax1),VV(mmax,mmax)
      real *8 B0(ntrnmax),x0(ipcmax),
     $        work0(ipcmax)
c
      real *8 beta(maxproxy,ipcmax)
c
      real *8 zero,dt,one,small,adum
      real *8 sumtot
      character *83 name
      character *90 name1,name2
      character *90 flnm(maxproxy)
      character *90 nome0
      character *90 nome(maxproxy)
      character *94 nome2(maxproxy)
      character *90 flnm2
c
      zero = 0.d0
      one = 1.d0
      small = 1e-8
c
c     dt    = time interval spacing in units of years
c     nscan = length of time series in units of "dt"
c     iabv  = number of time series in analysis (ie, "spatial" array)
c
c     define limits for entire joint proxy/instrumental dataset
c
c     "nlow" will define reference year #1
c
      IP=0
      IP2 = 0
      nlow = 500
      nhigh = 1980
      ntot = nhigh-nlow+1
c
c     define limits for raw data and training interval
c
      itrainmin0 = 1902
      itrainmax0 = nhigh
      iyearmax = 1993
      irawmin = 1902
      iearlyinst = 1820
      irawmax = nhigh
      nraw = irawmax-irawmin+1
      nskip0 = 0
      nkeep1 = nmax1-nskip0
c      
      do i=1,maxproxy
         istart(i)=99999
      end do
      do i=1,nmax0
         do j=1,max0
            aint(i,j)=zero
            aver0(i,j)=zero
         end do
      end do
      do i=1,nmax2
         do j=1,max0
            annual(i,j)=zero 
         end do
      end do
      do i=1,maxmonth
         do j=1,max0
            anom(i,j)=zero
         end do
      end do
c
      dt = 1.d0
c
c     set defaults
c
      iproxy = 0
      itrainmin = itrainmin0
      itrainmax = itrainmax0
      iverif = 0
      imonte = 1
      ivermin = 1854
      ivermax = itrainmin-1
      nverif = ivermax-ivermin+1
      ireconmin = 1854
      ireconmax = itrainmax0
      iset = 1
      iunif = 1
      iproxmin = 1820
      iproxmax = 1980
      iysoi1 = 1865
      iysoi2 = ireconmax
c
      ntrain = itrainmax-itrainmin+1
      nstart1 = itrainmin-nlow+1
      nend1 = itrainmax-nlow+1
      nscan = nend1-nstart1+1
      ipc = 40
      neofs = 10
      ifilter = 0
      nfilter = neofs
      do i=1,nfilter
         iwhich(i)=i
      end do
      iseason = 0
      itrain=1
      itype = 0
      alat1=-90.0
      alat2=90.0
      alon1=0.0
      alon2=360.0
      iteofs = 0
c
c
c
8888  write (6,*) 
     $  '1) train against ',itrainmin,
     $      '-',itrainmax,' temperature gridpoint data'
      if (iseason.eq.0) then
         write (6,*) '    seasonality: all seasons'
      else
      if (iseason.eq.1) then
         write (6,*) '    seasonality: cold season'
      else
         write (6,*) '    seasonality: warm season'
      endif
      endif
c
c
      if (iproxy.eq.0) then
        write (6,*) '2) train with proxy+long instr data:'
      else
       if (iproxy.eq.1) then
         write (6,*) '2) train with proxy data:'
       else
        if (iproxy.eq.2) then
         write (6,*) '2) train with long instr data:'
        else
        endif
       endif
      endif
      if (ifilter.ne.2) then
      write (6,*) '   (train with first',neofs,' pcs)'
      endif
      if (ifilter.eq.1) then
         write (6,*) '*   use specified subset of the first ',neofs,
     $      'eofs: ',(iwhich(j),j=1,nfilter)
      else
        if (ifilter.eq.2) then
           write (6,*) ' use "rule N" criterion applied to proxies'
        endif
      endif
      if (iset.eq.1) then
       write (6,*) '  (use all records available for any given year)'
      else
       write (6,*) '  (use all records dating back to ',iproxmin,' )'
      endif
      if (iunif.eq.1) then
       write (6,*) '   * use uniform weights on all proxies'
      else
       write (6,*) '   * use pre-specified weights for proxies'
      endif
c
      write (6,*) '3) training interval: ',itrainmin,itrainmax
      write (6,*) '4) reconstruction interval: ',ireconmin,ireconmax
      if (iverif.eq.0) then
        write (6,*) '5) verification: '
        write (6,*) '    (a)',ivermin,'-',ivermax,
     $     ' independent sub-period'
        write (6,*) '    (b)',iysoi1,'-',ivermax,
     $        ' soi/nino3 correlation'
        write (6,*) '    (c) long early historical records'
        write (6,*) '    (d) Quinn index'
        if (imonte.eq.1) then
           write (6,*) '         * compute bootstrap confidences'
        endif
      else
         if (iverif.eq.1) then
          write (6,*) '5) verification: ' 
          write (6,*) '   (a)',ivermin,'-',
     $          ivermax,' independent sub-period'
          write (6,*) '   (b)',iysoi1,'-',ivermax,
     $         ' soi/nino3 correlation'
         else
          write (6,*) '5) verification: ' 
          write (6,*) '   (a) long early historical records'
          write (6,*) '   (b) Quinn index'
          if (imonte.eq.1) then
           write (6,*) '         * compute bootstrap confidences'
          endif
         endif
      endif
c
c      
      write (6,*) 'option to change (0=continue)'
      read (5,*) ioption
c
      if (ioption.eq.1) goto 1111
      if (ioption.eq.2) goto 2222
      if (ioption.eq.3) goto 3333
      if (ioption.eq.4) goto 4444
      if (ioption.eq.5) goto 5555
c
      goto 6666
c
c
c
1111  write (6,*) 'seasonality'
      write (6,*) '(0) all seasons'
      write (6,*) '(1) cold season'
      write (6,*) '(2) warm season'
      read (5,*) iseason
      goto 8888
c
2222  write (6,*) 'train with: '
      write (6,*) '(0) all proxies + long histor data'
      write (6,*) '(1) proxy data only'
      write (6,*) '(2) long histor data only' 
      read (5,*) iproxy
      write (6,*) 'records to be used'
      write (6,*) 'earliest proxy sampling year is ',nlow
      write (6,*) '(0) fixed set dating back to specified year'
      write (6,*) '(1) all available for any given year'
      read (5,*) iset
      if (iset.eq.0) then
         write (6,*) 'year back to which records must date'
        read (5,*) iproxmin
      endif
      write (6,*) '  weights on proxies'
      write (6,*) ' (0) pre-specified'
      write (6,*) ' (1) uniform'
      read (5,*) iunif
      write (6,*) 'keep first how many pcs for training (1<pcs< ',
     $      min(ipc,maxproxy),')'
      write (6,*) '(enter 0 for specified subset,'
      write (6,*) ' enter -1 for automatic rule N application)'
      read (5,*) neofs
      if (neofs.eq.-1) then
         ifilter = 2
         neofs = ipc
      endif
      if (neofs.eq.0) then
        ifilter = 1
        write (6,*) 
     $     'keep first how many pcs for training (1<pcs< ',
     $      min(ipc,maxproxy),')'
        read (5,*) neofs
        nfilter = 0
        write (6,*) 'eof #s to use (0=done)'
        do j=1,ipcmax
           read (5,*) ii
           if (ii.eq.0) goto 8888
           nfilter = nfilter+1
           iwhich(nfilter)=ii
        end do
      else
       nfilter = neofs
       do i=1,nfilter
          iwhich(i)=i
       end do
      endif
      goto 8888
c
c
3333  write (6,*) 'training data available for the interval'
      write (6,*) itrainmin0,' to ',itrainmax0
      write (6,*)
      write (6,*) 'select training interval:'
      write (6,*) 'beginning year: '
      read (5,*) imin
      write (6,*) 'ending year: '
      read (5,*) imax
      if  ((imin.lt.itrainmin0).or.(imax.gt.itrainmax0)) then
         write (6,*) 'invalid subset'
         goto 3333
      endif
      itrainmax = imax
      itrainmin = imin 
      nskip0 = 12*(itrainmin-irawmin)
      nkeep1 = nmax1-nskip0
      ivermax = itrainmin-1
      ntrain = imax-imin+1
      write (6,*) 'training inverval: ',itrainmin,' to ',
     $             itrainmax
      goto 8888
c
c
4444  write (6,*) 'select reconstruction interval:'
      write (6,*) 'beginning year: '
      read (5,*) ireconmin
      write (6,*) 'ending year: '
      read (5,*) ireconmax
      goto 8888
c
c
5555  write (6,*) 'verification type: '
      write (6,*) '(0) ',ivermin,'-',ivermax,' sub-period + long histor'
      write (6,*) '(1) ',ivermin,'-',ivermax,' sub-period only'
      write (6,*) '(2) long histor only'
      read (5,*) iverif
      if (iverif.ne.1) then
        write (6,*) 'minimum verification year:'
        read (5,*) ivermin
        write (6,*) 'maximum verification year:'
        read (5,*) ivermax
        write (6,*) 'calc bootstrap stats for g^2?'
        write (6,*) '(0) no'
        write (6,*) '(1) yes'
        read (5,*) imonte
      endif
      goto 8888
c
6666  continue
      ntrain = itrainmax-itrainmin+1
      nmax = ntrain
      nstart1 = itrainmin-nlow+1
      nend1 = itrainmax-nlow+1
      nscan = nend1-nstart1+1
      iearly = ireconmin 
      ilate = ireconmax
c
c     read in instrumental surface temperature data (these correspond to the 1082 quasi-continuous
c          instrumental gridpoint records available from 1902-1993 as described in Mann et al 1998
c          supplementary information)
c
      write (6,*) 'reading in monthly instrumental data...'
c
      do i=7,maxstat/20+7
        ix1 = i-6
        id1 = ix1/10
        id2 = ix1-id1*10
        name = '/p0/tape/big2/DATA/JONESBRIFFA/MONTHLY/glb-train-month'
     $     //char(48+id1)//char(48+id2)//'.int'
        open (unit=1,file=name,status='old')
        do imonth=1,maxmonth
          read (1,*,end=44) ayr,(anom(imonth,k+20*(ix1-1)),k=1,20)
        end do
44      close (unit=1)
      end do
c
      nmonlong = imonth-1
      write (6,*) nmonlong,' months read in...'
c
      im1 = (itrainmin-1854)*12+1
      im3 = (itrainmax-1854+1)*12
      im2 = (iyearmax-1854+1)*12
      ncount = im2-im1+1
      nmonths = im3-im1+1
      write (6,*) 'number of months for TPCA decomposition: ',
     $     ncount
c
c     read in locations of instrumental gridpoint data and assign
c     areal weightings
c
      open (unit=1,file=
     $   '/p0/tape/big2/DATA/JONESBRIFFA/MONTHLY/globe-1902.dat',
     $    status='old')
      do i=1,maxstat
         igood(i)=0
         read (1,*) idum1,idum2,lon(i),lat(i)
         cs(i) = dcos(lat(i)*pi/180.d0)
         if ( ((lon(i).ge.alon1).and.(lon(i).le.alon2)).and.
     $        ((lat(i).ge.alat1).and.(lat(i).le.alat2)) ) 
     $          igood(i)=1       
      end do
      close (unit=1)
c
c
c     standardize instrumental data by detrended training period
c     variance for indicated seasonality. 
c
c     note: we remove calibration mean from monthly data
c
c
      jc = 0
      area = 0.0
      areanh = 0.0
      do j=1,maxstat
c
c      first determine if this gridpoint is in selected
c      subdomain
c
       if (igood(j).eq.1) then
c
         jc = jc+1
         cslat(jc)=cs(j)
         area = area + cs(j)
         if (lat(j).gt.zero) areanh = areanh + cs(j)
         newlon(jc)=lon(j)
         newlat(jc)=lat(j)
         do i=1,im2
            standard(i,jc)=anom(i,j)
         end do
c
       endif
      end do
      ispat = jc
c
      open (unit=15,file='TPCA/tpca-standard.out',status='unknown')
c
c     note: im1 and im3 bound the desired training interval 
c     of "nmonths" duration
c     
c     however: the interval im1 and im2  of "ncount" duration
c     may have to be longer in order to insure an overconstrained PCA 
c     (ncount > ngridpts)
c
c     base reference period mean on calibration period+update
c     through 1993 (im1 to im3)
c
c
      do  j=1,ispat
c
       asum = 0.0
       amean = 0.0
       do i=im1,im3
          amean = amean + standard(i,j)
       end do
       mean(j) = amean/float(nmonths)
c
       sum1 = zero
       sum2 = zero
       sum5 = zero
       sum6 = zero
       do i=im1,im2
          anom(i,j)=standard(i,j)-mean(j)
          sum1 = sum1 + float(i)
          sum2 = sum2 + float(i**2)
          sum5 = sum5 + anom(i,j)
          sum6 = sum6 + anom(i,j)*float(i)
       end do
       slope = (ncount*sum6-sum1*sum5)
     $       /(ncount*sum2-sum1**2)
       ainter = (sum2*sum5-sum1*sum6)/
     $      (ncount*sum2-sum1**2)
c
c
c      standardize monthly data by detrended variance
c
       amean = 0.0
       sigma(j)=0.0
       do i=im1,im2
         standard(i,j)=anom(i,j)-slope*float(i)
     $      - ainter
         sigma(j)=sigma(j)+standard(i,j)**2
       end do
       sigma(j)=sqrt(sigma(j)/float(ncount))
       do i=im1,im2
          anom(i,j)=anom(i,j)/sigma(j)
       end do
       write (15,*) j,mean(j),sigma(j)
      end do
c
      close (unit=15)
      iabv = ispat
c
c
      write (6,*) 'total of ',iabv,' gridpoints of training data'
      write (6,*)
c
c
c     FIRST, WE PERFORM THE EOF DECOMPOSITION ON THE
c     STANDARDIZED MONTHLY DATA
c
c
      M = iabv
      NU = M
      NV = M
      do i=1,iabv+1
        partsum(i)=zero
        sum(i)=zero
      end do
      sumtot = 0.d0
      niter = 0
c
c
c     set gridpoint weight on instrumental data
c
c     default weight = cos(latitude)
c
      do j=1,iabv
         weight(j)=cslat(j)
      end do
c
c
c     calculate EOFs of surface temperature record
c 
      i = 0
      do i2=im1,im2
         i = i+1
         do j=1,iabv
            adum = weight(j)*anom(i2,j)
            AA(i,j)=dcmplx(adum,zero)
         end do
      end do
      N = i
c
c
c     calculate EOFs using (complex) SVD
c
c     note: this is an expensive step--if the same instrumental
c           eigenvectors are to be used in multiple analyses, these
c           should be saved and read in, rather than re-calculated each
c           time
c
c
      call CSVD(AA,nmax1,mmax,N,M,IP,NU,NV,S,UU,VV)
c
c     collect eigenvalues, fractional variances explained..
c
      do i=1,M
           sum(i) = sum(i) + S(i)**2
           sumtot = sumtot + S(i)**2
           do j=i,M
             partsum(i)=partsum(i)+S(j)**2
           end do
      end do
c
c
      open (unit=16,file='TPCA/tpca-eigenvals.out',status='unknown')
c
      do i=1,M
          write (16,*) i,S(i),
     $      real(sum(i)/sumtot),
     $      one-real(partsum(i+1)/sumtot)
      end do
c
      do i=1,ipc
        id1 = i/10
        id2 = i-id1*10
        name1 = 'TPCA/eof'
     $     //char(48+id1)//char(48+id2)//'.out'
        open (unit=30+i,file=name1,status='unknown')
        do j=1,iabv
           eof(j,i)=real(VV(j,i))
           write (30+i,*) j,real(VV(j,i))
        end do
        close (unit=30+i)
      end do
c
      do i=1,ipc
        id1 = i/10
        id2 = i-id1*10
        name2 = 'TPCA/pc'
     $     //char(48+id1)//char(48+id2)//'.out'
        open (unit=30+i,file=name2,status='unknown')
        do j=1,nmax1
           write (30+i,*) j,real(UU(j,i))
        end do
        close (unit=30+i)
      end do
c
      sumtot0 = sumtot
      do j=1,ipc
      S0(j)=S(j)
      end do
c
      write (6,*) 'done calculating TPCA ...'
c
c 
c     now we will calculate the appropriately seasonally/annually
c     averaged pcs
c
c     only use the specified training interval bounded by 
c     im1 and im3 for calibration, but calculate annual mean
c     pcs between im1 and im2  
c
c
      do j=1,ipc
c
      S0(j)=S(j)
      j0 = j
c
      amean = 0.0
      icount = 0
      nyear = 0
      asum = 0.0
      if (iseason.eq.0) then
         mlim = 12
      else
        if (iseason.eq.1) then
           mlim = 6
           icount = 3
           asum = 3.0*real(UU(1,j))
        else
           mlim = 6
           icount = 0
           asum = 0.0
        endif
      endif
      do i=im1,im2
         iwhere = i-im1+1
         ii = mod(i-1,12)+1
         iflag = 0
         if (iseason.eq.0) then
            iflag=1
         else
          if (iseason.eq.1) then
            if ((ii.ge.10).or.(ii.le.3)) iflag=1
          else
            if ((ii.le.9).and.(ii.ge.4)) iflag=1
          endif
         endif
         if (iflag.eq.1) then
            icount = icount+1
            asum = asum + real(UU(iwhere,j)) 
            if (icount.eq.mlim) then
               nyear = nyear+1
               anew(nyear,j0)=asum/float(mlim)
               icount = 0
               asum = 0.0
            endif
         endif
      end do
c     
      end do
c
c     calculate appropriate seasonal/annual average for
c     the raw gridpoint data
c
      amean = 0.0
      icount = 0
      nyear = 0
      asum = 0.0
c
      open (unit=15,file='sigmas.out',status='unknown')
c
      do j=1,iabv
c
      sigma0(j)=0.0
      amean = 0.0
      icount = 0
      nyear = 0
      asum = 0.0
      if (iseason.eq.0) then
         mlim = 12
      else
        if (iseason.eq.1) then
c          use first JFM to form first cold-season half-year mean,
c          assuming persistence for previous OND
           mlim = 6
           icount = 3
           asum = 3.0*anom(1,j)
        else
           mlim = 6
        endif
      endif
      nyear0 = 0
      do i=im1,im2
         iwhere = i-im1+1
         ii = mod(i-1,12)+1
         iflag = 0
         if (iseason.eq.0) then
            iflag=1
         else
          if (iseason.eq.1) then
            if ((ii.ge.10).or.(ii.le.3)) iflag=1
          else
            if ((ii.le.9).and.(ii.ge.4)) iflag=1
          endif
         endif
         if (iflag.eq.1) then
            icount = icount+1
            asum = asum + anom(i,j)*sigma(j)
            if (icount.eq.mlim) then
               if (i.le.im3) nyear0=nyear0+1
               nyear = nyear+1
               annual(nyear,j)=asum/float(mlim)
               if (i.le.im3) amean = amean+annual(nyear,j)
               icount = 0
               asum = 0.0
            endif
         endif
      end do
      annmean(j)=amean/float(nyear0)
c
c     calculated detrended seasonal/annual variance
c
      sum1 = 0.0
      sum2 = 0.0
      sum5 = 0.0
      sum6 = 0.0
      do i=1,nyear0
         sum1 = sum1 + float(i)
         sum2 = sum2 + float(i**2)
         sum5 = sum5 + annual(i,j)
         sum6 = sum6 + annual(i,j)*float(i)
      end do
c
      slope = (nyear0*sum6-sum1*sum5)
     $       /(nyear0*sum2-sum1**2)
      ainter = (sum2*sum5-sum1*sum6)/
     $      (nyear0*sum2-sum1**2)
c
      do i=1,nyear
         detr=annual(i,j)-slope*float(i)
     $      - ainter
         sigma0(j)=sigma0(j)+detr**2
      end do
      sigma0(j)=sqrt(sigma0(j)/float(nyear))
c
      write (15,*) j,sigma(j),sigma0(j),annmean(j)
      end do
      close (unit=15)
c     
      nyears = nyear0
c
      do j=1,ipc
c 
c
c      remove calibration mean from pcs
c
       datave(j)=0.0
       do i=1,nyear0
          datave(j)=datave(j)+anew(i,j)
       end do
       datave(j)=datave(j)/float(nyear0)


       sum1 = 0.0
       sum2 = 0.0
       sum5 = 0.0
       sum6 = 0.0
       do i=1,nyear
          anew(i,j)=anew(i,j)-datave(j)
       end do
       do i=1,nyear0
          sum1 = sum1 + float(i)
          sum2 = sum2 + float(i**2)
          sum5 = sum5 + anew(i,j)
          sum6 = sum6 + anew(i,j)*float(i)
       end do
c
       slope = (nyear0*sum6-sum1*sum5)
     $       /(nyear0*sum2-sum1**2)
       ainter = (sum2*sum5-sum1*sum6)/
     $      (nyear0*sum2-sum1**2)
c
c      now standardize the seasonal pcs (based on 
c      calibration period variance)
c
       amean = 0.0
       sd(j)=0.0
       do i=1,nyear
         standard(i,j)=anew(i,j)-slope*float(i)
     $      - ainter
       end do
       do i=1,nyear0
         sd(j)=sd(j)+standard(i,j)**2
       end do
       sd(j)=sqrt(sd(j)/float(nyear0))
       do i=1,nyear
          anew(i,j)=anew(i,j)/sd(j)
       end do
      end do
c
      nyearf = nyear
      nyearf2 = nyear+irawmax-itrainmax
      nyear = nyear0
      write (6,*) 'nyear nearf nyears'
      write (6,*) nyear, nyearf, nyears
c
      write (6,*) 'reading in training data...'
c
c     read in proxy data
c
555   format (f5.2,1a40)
c
c     read in list of proxy data to be used
c
      if (iproxy.eq.0) then
        open (unit=1,
     $  file='/p0/tape/big2/MULTIPROXY/DATA/multiproxy.dat',
     $  status='old')
      else
        if (iproxy.eq.1) then
         open (unit=1,
     $   file='/p0/tape/big2/MULTIPROXY/DATA/multiproxy-proxy.dat',
     $   status='old')
        else
          if (iproxy.eq.2) then
           open (unit=1,
     $     file='/p0/tape/big2/MULTIPROXY/DATA/multiproxy-instr.dat',
     $     status='old')
          endif
        endif
      endif
c
c
      read (1,*)
      i=0
1010  i=i+1
         if (i.gt.maxproxy) goto 415
         read (1,555,end=415) weightprx(i),nome(i)
         if (weightprx(i).lt.zero) then
            nome0=nome(i)
            i=i-1
            goto 1010
         else
           if (weightprx(i).lt.small) then
               weightprx(i)=small
           else
              if (iunif.eq.1) weightprx(i)=one
           endif
         endif
         if (weightprx(i).ge.zero) then
c
         flnm(i) = '/p0/tape/big2/MULTIPROXY/DATA/'
     $       //nome0(1:index(nome0,' ')-1)//'/'//nome(i)
         write (6,*) flnm(i)
         endif
      goto 1010
415   continue
      nproxy = i-1
      close (unit=1)
c     
c     now read in actual proxy records
c
      do j=1,nproxy
         ii = j
         open (unit=9,file=flnm(j),status='old')
         iflag = 0
         do i=1,nmax0
            aprox(i,ii)=-99999999.0
         end do
         do i=1,nlarge
            read (9,*,end=88) ayear,adat
            iyr = int(ayear+0.5)
            if (iflag.eq.0) then
                istart(j)=iyr
                iflag = 1
            endif
            if ((iyr.ge.nlow).and.(iyr.le.nhigh)) then
               iyear = iyr-nlow+1
               aprox(iyear,ii)=adat
            endif
         end do
88       close (unit=9)
c
c        extend series based on assumed persistence of final value if series ends 
c        before end of training interval
c
         iyear = nhigh-nlow+2
737      iyear = iyear-1
         if (aprox(iyear,ii).lt.-99999990.0) goto 737
         do iyr = iyear+1,nhigh-nlow+1
            aprox(iyr,ii)=aprox(iyear,ii)
         end do
c       
      end do
c
c
      write (6,*) 'nyear nearf nyears'
      write (6,*) nyear, nyearf, nyears
      write (6,*) 'read in: '
      write (6,*) nproxy, ' time series'
c
c
c
9999  continue
c
c     standardize training data by detrended reference
c     period variance
c
c
      write (6,*) 'ntrain: ',ntrain
c
      do j=1,nproxy
c 
       proxave(j)=zero
       do i=1,ntrain
          iyear = itrainmin+i-nlow
          proxave(j)=proxave(j)+aprox(iyear,j)
       end do
       proxave(j)=proxave(j)/float(ntrain)
c
c      remove 1902-19XX mean from training data
c
       do i=nlow,nhigh
          iyear = i-nlow+1
          aprox(iyear,j)=aprox(iyear,j)-proxave(j)
       end do
      end do
c
c     standardize proxy data using detrended variance
c
      do j=1,nproxy
c 
       sum1 = zero
       sum2 = zero
       sum5 = zero
       sum6 = zero
       do i=1,ntrain
          iyear = itrainmin+i-nlow
          sum1 = sum1 + float(i)
          sum2 = sum2 + float(i**2)
          sum5 = sum5 + aprox(iyear,j)
          sum6 = sum6 + aprox(iyear,j)*float(i)
       end do
c
       slope = (ntrain*sum6-sum1*sum5)
     $       /(ntrain*sum2-sum1**2)
       ainter = (sum2*sum5-sum1*sum6)/
     $      (ntrain*sum2-sum1**2)
c
c
       amean = 0.0
       sdprox(j)=0.0
       do i=1,ntrain
         iyear = itrainmin+i-nlow
         standprx=aprox(iyear,j)-slope*float(i)
     $      - ainter
         sdprox(j)=sdprox(j)+standprx**2
       end do
       sdprox(j)=sqrt(sdprox(j)/float(ntrain))
       do i=nlow,nhigh
          iyear = i-nlow+1
          aprox(iyear,j)=aprox(iyear,j)/sdprox(j)
       end do
      end do
c
c     select proxy records to be used
c
      if (iset.eq.0) then
        mkeep = 0
        do i=1,nproxy
         if (istart(i).le.iproxmin) then
            mkeep = mkeep + 1
            ikeep(mkeep)=iabv+i
         endif
        end do
        write (6,*) 'total of ',mkeep,
     $      ' training records back to ',iproxmin
      endif
c 
c     now determine suggested number of EOFs in training
c     based on rule N applied to the proxy data alone
c     during the interval t > iproxmin (the minimum
c     year by which each proxy is required to have started,
c     note that default is iproxmin = 1820 if variable
c     proxy network is allowed (latest begin date
c     in network)
c
c     we seek the n first eigenvectors whose eigenvalues
c     exceed 1/nproxy'
c
c     nproxy' is the effective climatic spatial degrees of freedom
c     spanned by the proxy network (typically an appropriate
c     estimate is 20-40)
c
      write (6,*) 'determining optimal # EOFs for training'
      im1 = iproxmin
      im2 = iproxmax
      N0 = iproxmax-iproxmin+1
      jj = 0
      do j=1,nproxy
         if (istart(j).le.iproxmin) then
           jj = jj+1
           i = 0
           do i2=im1,im2
              i = i+1
              iyear = i2-nlow+1
              adum = aprox(iyear,j)
              AA(i,jj)=dcmplx(adum,zero)
           end do
         end if
      end do
      M0 = jj
      NU0 = M0
      NV0 = M0
c
c     taking into account spatial dependence between
c     different proxy records, estimate number
c     of climatic spatial degrees of freedom sampled
c        (actually, estimate reasonable range)
      np1 = 60
      np2 = 40
      np3 = M0/2
      write (6,*)  np1,np2,np3
      ruleN1 = one/float(np1)
      ruleN2 = one/float(np2)
      ruleN3 = one/float(np3)
c
      call CSVD(AA,nmax1,mmax,N0,M0,IP0,NU0,NV0,S,UU,VV)
c
c     collect eigenvalues, fractional variances explained..
c
      sumtot = zero
      do i=1,M0
           partsum(i)=zero
           sumtot = sumtot + S(i)**2
           do j=i,M0
             partsum(i)=partsum(i)+S(j)**2
           end do
      end do
c
c
      open (unit=13,file='ruleN.out',status='unknown')
      nretain = 1
      do i=1,min(30,M0)
         eigen = S(i)**2/sumtot
         write (13,434) i,
     $     S(i)**2,eigen,one-real(partsum(i+1)/sumtot),
     $     ruleN3,ruleN1,ruleN2
         if (eigen.gt.ruleN3) nretain = i
      end do
      close (unit=13)
      write (6,*) 'rule N gives ',nretain,
     $       ' eofs to be retained'
      if (ifilter.eq.2) then
         neofs = nretain
      endif
434   format (i4,6f9.4)
c
c
      write (6,*) 'nproxy = ',nproxy
      open (unit=15,file='TRAIN/pca-standard.out',status='unknown')
c
      do j=1,ipc
         write (15,*) j,datave(j),sd(j)
      end do
      do j=1,nproxy
         write (15,*) j,proxave(j),sdprox(j)
      end do
      close (unit=15)
c
      do i=1,ipc
        id1 = i/10
        id2 = i-id1*10
        name2 = 'TRAIN/pc'
     $     //char(48+id1)//char(48+id2)//'.out'
        open (unit=30+i,file=name2,status='unknown')
        do j=1,nyearf
           iyear = itrainmin+j-1
           tpc(j,i)=anew(j,i)*sd(i)+datave(i)
           write (30+i,*) iyear,tpc(j,i)
        end do
        close (unit=30+i)
      end do
c
c     determine the pc-filtered fields
c
      do i=1,nyearf
         do j=1,iabv
            a(i,j)=zero
            do l=1,nfilter
               ll = iwhich(l)
               a(i,j)=a(i,j) + S0(ll)*sigma(j)
     $          *tpc(i,ll)*eof(j,ll)/weight(j)
            end do
         end do
      end do
c
c     write out raw and eof-filtered seasonal/annual temperature
c
      open (unit=31,file='RECONDATA/raw.dat',status='unknown')
      do i=1,iabv
        do j=1,nyearf2
           write (31,*) i,j+itrainmin-1,annual(j,i)
c           if (i.eq.20) write (31,*) j+itrainmin-1,annual(j,i)
        end do
      end do
      close (unit=31)
c
      open (unit=31,file='RECONDATA/filtered.dat',status='unknown')
c
      do i=1,iabv
        do j=1,nyearf
           write (31,*) i,j+itrainmin-1,a(j,i)
c           if (i.eq.20) write (31,*) j+itrainmin-1,a(j,i)
        end do
      end do
      close (unit=31)
c
c     NOW WE TRAIN THE p PROXY RECORDS AGAINST THE FIRST
c     neofs ANNUAL/SEASONAL RESOLUTION INSTRUMENTAL pcs
c
      N = nyears
      ntotal = ipc
      M = ntotal
      neofmax = M
      if (nfilter.gt.neofmax) then
          write (6,*) 'sorry--too many ceofs kept!'
          write (6,*) 'max is: ',neofmax
          stop
      endif
c
      if (N.lt.M) then
         write (6,*) 'SVD is underdetermined!'
         stop
      endif
c
c
c     set specified weights on data 
c
      open (unit=15,file='TRAIN/pca-weights.out',status='unknown')
c
c     downweight proxy weights to adjust for size
c     of proxy vs. instrumental samples
c
c     weights on PCs are proportional to their singular values
c
      ssum = 0.0
      do j=1,ipc
          ssum = ssum+S0(j)
      end do
      ssum = ssum/float(ipc)
c
      do j=1,ipc
         weight0(j)=S0(j)/ssum
         write (15,*) j,weight0(j)
      end do
      do j=1,nproxy
         write (15,*) j,weightprx(j)
      end do
      close (unit=15)
c
      NU = M
      NV = M
      do i=1,M
        partsum(i)=zero
        sum(i)=zero
      end do
      sumtot = 0.d0
      niter = 0
c
c     now we will compute the proxy <-> instr pc transfer matrix
c     one proxy at a time
c
      if (nproxy.lt.nfilter) then
         write (6,*) 'beware: # proxies less than # eofs!!'
      endif
      M2 = nfilter
      NU2 = M2
      NV2 = M2
c
c
      do ii=1,nproxy
c
c     right side --the ``known'' values of proxy ``i'' during training
c     period
c
      do k=1,ntrain
            iyear = itrainmin+k-nlow
            B0(k)=weightprx(ii)*aprox(iyear,ii)
      end do
c
c     matrix of the neofs training period  PCs 
c
      do k=1,ntrain
         do j=1,nfilter
            jj = iwhich(j)
            AA(k,j)=dcmplx(weight0(jj)*anew(k,jj),zero)
         end do
      end do
c
c     now determine the transfer function vector and store
c     in the nproxy x matrix beta
c
      N2 = ntrain
c
      call CSVD(AA,nmax1,mmax,N2,M2,IP2,NU2,NV2,SS,UU,VV)
      do k=1,nfilter
         work0(k)=zero
          do j=1,ntrain
            work0(k)=work0(k)+real(UU(j,k))*B0(j)/SS(k)
          end do
      end do
      do i=1,nfilter
         x0(i)=zero
         do k=1,nfilter
           x0(i) = x0(i) + real(VV(i,k))*work0(k)
         end do
c
         beta(ii,i)=x0(i)
      end do
c
c
      end do
c
c     determine the reconstructed fields by inverting the transfer
c     function
c
      write (6,*) 'reconstructing patterns from ',iearly,' to ',
     $    ilate
c
c
c     loop over desired time reconstruction interval
c
      iymin0 = itrainmin-nlow+1
      iymax0 = ilate-nlow+1
      iymax1 = iyearmax-nlow+1
c
      do it=iearly,ilate
c
c      determine which proxy indicators are available
c
c       if ((mod(it,20).eq.0).or.(it.eq.iearly)) 
c     $      write(6,*) 'finished year ',it,'...'
       iy = it-nlow+1
       mkeep = 0
       idate = iproxmin
       if (iset.eq.1) idate=it
       do i=1,nproxy
         if (istart(i).le.idate) then
            mkeep = mkeep + 1
            ikeep(mkeep)=i
         endif
       end do
       nkeep(iy)=mkeep
c
       if (mkeep.lt.nfilter) then
         write (6,*) 'beware: # proxies less than # eofs!!'
       endif
       N2 = mkeep
       M2 = nfilter
       NU2 = M2
       NV2 = M2
c
       do i=1,mkeep
c
c        B contains the actual weighted values of the proxy indicators
c
         B0(i) = weightprx(ikeep(i))*
     $      aprox(iy,ikeep(i))
c
         do j=1,nfilter
c
c         determine projection of the jth ceof on the ith
c         proxy record
c
          AA(i,j)=dcmplx(beta(ikeep(i),j),zero)
c
         end do
       end do
c
c      now determine the least-squares optimal weights on
c      each of the "neofs" eofs contained in the array "x"
c
       call CSVD(AA,nmax1,mmax,N2,M2,IP2,NU2,NV2,SS,UU,VV)
       do k=1,nfilter
          work0(k)=zero
          do j=1,mkeep
            work0(k)=work0(k)+real(UU(j,k))*B0(j)/SS(k)
          end do
       end do
       do i=1,nfilter
         x0(i)=zero
         do k=1,nfilter
           x0(i) = x0(i) + real(VV(i,k))*work0(k)
         end do
       end do
c
c      record k weights as the k reconstructed pcs (rpcs) for
c      time iy
c
       do k=1,nfilter
          rpc(k)=x0(k)*sd(k)/weight0(k)+datave(k)
       end do
c
c      store the reconstructed cpcs
c
       do k=1,nfilter
          kk = iwhich(k)
          rpckeep(iy,kk)=rpc(k)
       end do
c
c     end time loop
c
      end do
c
c     insure that reconstructed pcs are properly
c     normalized (amplitude determined by fraction
c     of calibration resolved variance of actual pc)
c
      do k=1,nfilter
         kk = iwhich(k)
         pccalib(kk)=zero
         vartrn1 = zero
         vartrn2 = zero
         vartrn3 = zero
         amean1 = zero
         amean2 = zero
         do it=itrainmin,itrainmax
            iy = it-nlow+1
            jj = it-itrainmin+1
            vartrn1 = vartrn1 + tpc(jj,kk)**2
            vartrn2 = vartrn2 + rpckeep(iy,kk)**2
            amean1 = amean1+tpc(jj,kk)
            amean2 = amean2+rpckeep(iy,kk)
         end do
         amean1 = amean1/float(ntrain)
         amean2 = amean2/float(ntrain)
         do it=itrainmin,itrainmax
            iy = it-nlow+1
            jj = it-itrainmin+1
            pccalib(kk)=pccalib(kk)+(tpc(jj,kk)-amean1)*
     $             (rpckeep(iy,kk)-amean2)
         end do
         sd1 = sqrt(vartrn1/float(ntrain))
         sd2 = sqrt(vartrn2/float(ntrain))
         pccalib(kk)=pccalib(kk)/sqrt(vartrn1*vartrn2)
         do it=ireconmin,ireconmax
            iy = it-nlow+1
c            rpckeep(iy,kk)=pccalib(kk)*rpckeep(iy,kk)*sd1/sd2
             rpckeep(iy,kk)=rpckeep(iy,kk)*sd1/sd2
         end do
      end do
c
c     store and print out the reconstructed pcs
c
c
      do i=1,neofs
        ii = iwhich(i)
        id1 = i/10
        id2 = i-id1*10
        name2 = 'RECONDATA/rpc'
     $     //char(48+id1)//char(48+id2)//'.out'
        open (unit=30+i,file=name2,status='unknown')
        close (unit=30+i)
      end do
c
c        
      do i=1,nfilter
        ii = iwhich(i)
        id1 = ii/10
        id2 = ii-id1*10
        name2 = 'RECONDATA/rpc'
     $     //char(48+id1)//char(48+id2)//'.out'
        open (unit=30+ii,file=name2,status='unknown')
      end do
c
      do it=iearly,ilate
       iy = it-nlow+1
       do k=1,nfilter
          kk=iwhich(k)
          write (30+kk,*) it,rpckeep(iy,kk)
       end do
      end do
c
      do i=1,nfilter
        ii = iwhich(i)
        close (unit=30+ii)
      end do
c
c
c     reconstruct calibrated fields
c
c
c     filter out the nfilter highest variance reconstructed
c     pcs if indicated
c
c
      do it=ireconmin,ireconmax
c
      iy = it-nlow+1
c
      do j=1,iabv
         aint(iy,j)=zero
         do l=1,nfilter
           ll = iwhich(l)
           aint(iy,j)=aint(iy,j) + S0(ll)*sigma(j)
     $          *rpckeep(iy,ll)*eof(j,ll)/weight(j) 
         end do
      end do
c
      end do
c
c     write out calibrated seasonal/annual temeprature data
c     construct raw, filtered, and reconstructed nhem, global, 
c     and nino3 averages
c
      call averages(iabv,1,nyearf2,1082,nmax2,annual,
     $    cslat,newlon,newlat,nhemr,globr,ninor,detrr)
      call averages(iabv,1,nyearf,1082,nmax2,a,
     $    cslat,newlon,newlat,nhemt,globt,ninot,detrt)
c
c     calculate fractional variances explained relative to
c     raw seasonal/annual averages
c
      varrawtot = zero
      amsetrntot = zero
      do i=1,iabv   
        varrawgp(i)=zero
        amsetrngp(i)=zero
        do j=1,nyear
         varrawgp(i)=varrawgp(i)+annual(j,i)**2
         amsetrngp(i)=amsetrngp(i)+(annual(j,i)-a(j,i))**2
         varrawtot = varrawtot + varrawgp(i)
         amsetrntot = amsetrntot + amsetrngp(i)
        end do
      end do
c
      varrawglob = zero
      varrawnhem = zero
      varrawnino = zero
      varrawdetr = zero
      amsetrnglob = zero
      amsetrnnhem = zero
      amsetrnnino = zero
      amsetrndetr = zero
      do j=1,nyear
         varrawglob = varrawglob + globr(j)**2
         varrawnhem = varrawnhem + nhemr(j)**2
         varrawnino = varrawnino + ninor(j)**2
         varrawdetr = varrawdetr + detrr(j)**2
         amsetrnglob = amsetrnglob + (globr(j)-globt(j))**2
         amsetrnnhem = amsetrnnhem + (nhemr(j)-nhemt(j))**2
         amsetrnnino = amsetrnnino + (ninor(j)-ninot(j))**2
         amsetrndetr = amsetrndetr + (detrr(j)-detrt(j))**2
      end do
c
      open (unit=9,file='betas-train.out',status='unknown')
      write (9,*) 'globe:        ',one-amsetrnglob/varrawglob
      write (9,*) 'nhem:         ',one-amsetrnnhem/varrawnhem
      write (9,*) '  "   (detr): ',one-amsetrndetr/varrawdetr
      write (9,*) 'nino3:        ',one-amsetrnnino/varrawnino
      write (9,*)
      write (9,*) 'multivariate: ',one-amsetrntot/varrawtot
      write (9,*)
      write (9,*) 'gridpoints:'
      do i=1,iabv
         write (9,*) i,one-amsetrngp(i)/varrawgp(i)
      end do
      close (unit=9)
c
      open (unit=31,file='RECONDATA/nhem-raw.out',status='unknown')
      open (unit=32,file='RECONDATA/glob-raw.out',status='unknown')
      open (unit=33,file='RECONDATA/nino-raw.out',status='unknown')
      open (unit=34,file='RECONDATA/nhem-trn.out',status='unknown')
      open (unit=35,file='RECONDATA/glob-trn.out',status='unknown')
      open (unit=36,file='RECONDATA/nino-trn.out',status='unknown')
      open (unit=37,file='RECONDATA/detr-raw.out',status='unknown')
      open (unit=38,file='RECONDATA/detr-trn.out',status='unknown')
c
      do i=1,nyearf2
         write (31,*) itrainmin+i-1,nhemr(i)
         write (32,*) itrainmin+i-1,globr(i)
         write (37,*) itrainmin+i-1,detrr(i)
      end do
      do i=1,nyearf
         write (34,*) itrainmin+i-1,nhemt(i)
         write (35,*) itrainmin+i-1,globt(i)
         write (33,*) itrainmin+i-1,ninor(i)
         write (36,*) itrainmin+i-1,ninot(i)
         write (38,*) itrainmin+i-1,detrt(i)
      end do
      close (unit=31)
      close (unit=32)
      close (unit=33)
      close (unit=34)
      close (unit=35)
      close (unit=36)
      close (unit=37)
      close (unit=38)
c
c     write out reconstructed seasonal/annual temeprature data
c
      open (unit=31,file='RECONDATA/recon.dat',status='unknown')
c
      do i=1,iabv
        do it=iearly,ilate
           j = it-nlow+1
           write (31,*) i,it,aint(j,i)
c           if (i.eq.20) write (31,*) it,aint(j,i)
        end do
      end do
      close (unit=31)
c
c     reconstruct nhem, global, and nino3 averages
c
      call averages(iabv,iymin0,iymax1,1082,nmax0,aint,
     $    cslat,newlon,newlat,nhemc,globc,ninoc,detrc)
c
c
c     calculate fractional variances explained relative to
c     raw seasonal/annual averages
c
      varrawtot = zero
      amsecaltot = zero
      do i=1,iabv 
        varrawgp(i)=zero
        amsecalgp(i)=zero
        do j=1,nyear
         iy = itrainmin+j-nlow
         varrawgp(i)=varrawgp(i)+annual(j,i)**2
         amsecalgp(i)=amsecalgp(i)+(annual(j,i)-aint(iy,i))**2
         varrawtot = varrawtot + varrawgp(i)
         amsecaltot = amsecaltot + amsecalgp(i)
        end do
      end do
c
c
      do i=1,nfilter
        ii = iwhich(i)
        varrawpc(ii)=zero
        amsecalpc(ii)=zero
        do j=1,nyear
         iy = itrainmin+j-nlow
         varrawpc(ii)=varrawpc(ii)+tpc(j,ii)**2
         amsecalpc(ii)=amsecalpc(ii)+
     $     (tpc(j,ii)-rpckeep(iy,ii))**2
        end do
      end do
c
c
      varrawglob = zero
      varrawnhem = zero
      varrawnino = zero
      varrawdetr = zero
      amsecalglob = zero
      amsecalnhem = zero
      amsecalnino = zero
      amsecaldetr = zero
      do j=1,nyear
         iy = itrainmin+j-nlow
         varrawglob = varrawglob + globr(j)**2
         varrawnhem = varrawnhem + nhemr(j)**2
         varrawnino = varrawnino + ninor(j)**2
         varrawdetr = varrawdetr + detrr(j)**2
         amsecalglob = amsecalglob + (globr(j)-globc(iy))**2
         amsecalnhem = amsecalnhem + (nhemr(j)-nhemc(iy))**2
         amsecalnino = amsecalnino + (ninor(j)-ninoc(iy))**2
         amsecaldetr = amsecaldetr + (detrr(j)-detrc(iy))**2
      end do
c
      open (unit=9,file='betas-recon1.out',status='unknown')
      write (9,*) 'globe:        ',one-amsecalglob/varrawglob
      write (9,*) 'nhem:         ',one-amsecalnhem/varrawnhem
      write (9,*) ' "   (detr):  ',one-amsecaldetr/varrawdetr
      write (9,*) 'nino3:        ',one-amsecalnino/varrawnino
      write (9,*)
      write (9,*) 'multivariate: ',one-amsecaltot/varrawtot
      write (9,*)
      write (9,*) 'pcs:'
      do i=1,nfilter
         ii = iwhich(i)
         write (9,*) ii,one-amsecalpc(ii)/varrawpc(ii)
      end do
      write (9,*)
      write (9,*) 'gridpoints: '
      do i=1,iabv
         write (9,*) i,one-amsecalgp(i)/varrawgp(i)
      end do
      close (unit=9)
c
c     calculate correlations with raw seasonal/annual averages
c
      do i=1,iabv 
        varrawgp(i)=zero
        varcalgp(i)=zero
        corrgp(i)=zero
        amean1 = 0.0
        amean2 = 0.0
        do j=1,nyear
         iy = itrainmin+j-nlow
         amean1 = amean1+annual(j,i)
         amean2 = amean2+aint(iy,i)
        end do  
        amean1 = amean1/float(nyear)
        amean2 = amean2/float(nyear)
        do j=1,nyear
         iy = itrainmin+j-nlow
         varrawgp(i)=varrawgp(i)+(annual(j,i)-amean1)**2
         varcalgp(i)=varcalgp(i)+(aint(iy,i)-amean2)**2
         corrgp(i)=corrgp(i)+(annual(j,i)-amean1)*
     $     (aint(iy,i)-amean2)
        end do
        corrgp(i)=corrgp(i)/sqrt(varcalgp(i)*varrawgp(i))
      end do
c
c
      do i=1,nfilter
        ii = iwhich(i)
        varrawpc(ii)=zero
        varcalpc(ii)=zero
        corrpc(ii)=zero
        amean1 = zero
        amean2 = zero
        do j=1,nyear
           iy = itrainmin+j-nlow
           amean1 = amean1 + tpc(j,ii)
           amean2 = amean2 + rpckeep(iy,ii)
        end do
        do j=1,nyear
         iy = itrainmin+j-nlow
         varrawpc(ii)=varrawpc(ii)+tpc(j,ii)**2
         varcalpc(ii)=varcalpc(ii)+rpckeep(iy,ii)**2
         corrpc(ii)=corrpc(ii)+tpc(j,ii)*rpckeep(iy,ii)
        end do
        corrpc(ii)=corrpc(ii)/sqrt(varcalpc(ii)*varrawpc(ii))
      end do
c
c
      varrawglob = zero
      varcalglob = zero
      varrawnhem = zero
      varcalnhem = zero
      varrawnino = zero
      varcalnino = zero
      varrawdetr = zero
      varcaldetr = zero
      corrglob = zero
      corrnhem = zero
      corrnino = zero
      corrdetr = zero
      do j=1,nyear
         iy = itrainmin+j-nlow
         varrawglob = varrawglob + globr(j)**2
         varcalglob = varcalglob + globc(iy)**2
         varrawnhem = varrawnhem + nhemr(j)**2
         varcalnhem = varcalnhem + nhemc(iy)**2
         varrawnino = varrawnino + ninor(j)**2
         varcalnino = varcalnino + ninoc(iy)**2
         varrawdetr = varrawdetr + detrr(j)**2
         varcaldetr = varcaldetr + detrc(iy)**2
         corrglob = corrglob + globr(j)*globc(iy)
         corrnhem = corrnhem + nhemr(j)*nhemc(iy)
         corrnino = corrnino + ninor(j)*ninoc(iy)
         corrdetr = corrdetr + detrr(j)*detrc(iy)
      end do
      corrglob = corrglob/sqrt(varrawglob*varcalglob)
      corrnhem = corrnhem/sqrt(varrawnhem*varcalnhem)
      corrnino = corrnino/sqrt(varrawnino*varcalnino)
      corrdetr = corrdetr/sqrt(varrawdetr*varcaldetr)
c
c     now calculate correlation with soi over
c     abbreviated period -- read in as "soir(iy)"
c
      nyear2 = iysoi2-iysoi1+1
      asum = zero
      lnino = 0
      open (unit=7,file='winsoi.dat',status='old')
      do i=1,nyear2
         iyear = iysoi1+i-1
         iy = iyear-nlow+1
         read (7,*) idum,soir(iy)
         if ((iyear.ge.1902).and.(iyear.le.1980)) then
           asum = asum+soir(iy)
           lnino = lnino + 1
         endif
      end do
      close (unit=7)
c
      asum = asum/float(lnino)
      do i=1,nyear2
         iyear = iysoi1+i-1
         iy = iyear-nlow+1
         soir(iy) = soir(iy)-asum
      end do
c
c
      varrawnino = zero
      varcalsoi = zero
      corrsoi = zero
      corrsoi = corrsoi/sqrt(varrawnino*varcalsoi)
c
      open (unit=19,file='corrs-cal.out',status='unknown')
      write (19,*) 'win soi/nino3:     ',corrsoi,corrsoi**2
      close (unit=19)
c
      open (unit=9,file='corrs-recon.out',status='unknown')
      write (9,*) 'globe:        ',corrglob,corrglob**2
      write (9,*) 'nhem:         ',corrnhem,corrnhem**2
      write (9,*) ' "   (detr):  ',corrdetr,corrdetr**2
      write (9,*) 'nino3:        ',corrnino,corrnino**2
      write (9,*)
      write (9,*) 'pcs:'
      do i=1,nfilter
         ii = iwhich(i)
         write (9,*) ii,corrpc(ii),corrpc(ii)**2
      end do
      write (9,*)
      write (9,*) 'gridpoints: '
      do i=1,iabv
         write (9,*) i,corrgp(i),corrgp(i)**2
      end do
      close (unit=9)
c
      open (unit=34,file='RECONDATA/nhem-cal.out',status='unknown')
      open (unit=35,file='RECONDATA/glob-cal.out',status='unknown')
      open (unit=36,file='RECONDATA/nino-cal.out',status='unknown')
      open (unit=37,file='RECONDATA/detr-cal.out',status='unknown')
c
      do iyear = ireconmin,ireconmax
         iy = iyear-nlow+1
         write (34,*) iyear,nhemc(iy)
         write (35,*) iyear,globc(iy)
         write (36,*) iyear,ninoc(iy)
         write (37,*) iyear,detrc(iy)
      end do
      close (unit=34)
      close (unit=35)
      close (unit=36)
      close (unit=37)
c
c
c     now compute verification statistics
c
c
      if (iverif.lt.2) then
c
c
c     if indicated, read in 19th century surface temperature data for use in cross-validation
c          (these correspond to the 219 quasi-continuous
c          instrumental gridpoint records available from 1854-1901 as described in Mann et al' 98
c          supplementary information)
c
c
      write (6,*) 'reading in long (back to 1854) gridpoint data...'
c
      open (unit=1,file=
     $ '/p0/tape/big2/DATA/JONESBRIFFA/1854-1993/globe-1854.dat',
     $    status='old')
      open (unit=3,file=
     $ '/p0/tape/big2/DATA/JONESBRIFFA/1854-1993/globe-1854.mask',
     $    status='old')
      imskcount = 0
      do i=1,maxlong
         read (1,*,end=55) idum1,idum2,lon(i),lat(i)
         read (3,*,end=55) imask(i)
         if (imask(i).eq.1) imskcount = imskcount+1
      end do
55    close (unit=1)
      close (unit=3)
      numlong = i-1
c
      do i=7,numlong/10+7
        ix1 = i-6
        id1 = ix1/10
        id2 = ix1-id1*10
        if (iseason.eq.0) then
         name = '/p0/tape/big2/DATA/JONESBRIFFA/1854-1993/glb-long-all'
     $   //char(48+id1)//char(48+id2)//'.int'
        else
         if (iseason.eq.1) then
       name = '/p0/tape/big2/DATA/JONESBRIFFA/1854-1993/glb-long-cold'
     $   //char(48+id1)//char(48+id2)//'.int'
         else
        name = '/p0/tape/big2/DATA/JONESBRIFFA/1854-1993/glb-long-warm'
     $  //char(48+id1)//char(48+id2)//'.int'
         endif
        endif
        open (unit=1,file=name,status='old')
        do iyr=1,nmax0
          iyear = iyr+1854-nlow
          read (1,*,end=34) ayr,(anom0(iyear,k+10*(ix1-1)),k=1,10)
          do k=1,10
             if (anom0(iyear,k+10*(ix1-1)).lt.-15.d0) 
     $          anom0(iyear,k+10*(ix1-1)) = zero
          end do
        end do
34      close (unit=1)
      end do
c
c
      write (6,*) numlong,
     $     ' total gridpoints back to 1854 available'
      write (6,*) imskcount, ' gridpoints in ocean-only mask'
c     check for which gridpoints match a gridpoint in the training
c     1902-19XX set
c
      do i=1,iabv
         index0(i)=0
      end do
      nlong = 0
      do j=1,numlong
         lon1 = lon(j)
         lat1 = lat(j)
c
         do i=1,iabv
c
           lon2 = newlon(i)
           lat2 = newlat(i)
c 
           if (  (abs(lon1-lon2).lt.half).and.
     $            (abs(lat1-lat2).lt.half)  ) then
              nlong = nlong+1
              index0(i)=j
           endif
c
          end do
      end do
c
      write (6,*) 'found ',nlong,' matches...'
      write (6,*) 
c
c     asign verification status to any gridpoint
c     in the matching set--otherwise nullify
c
      open (unit=31,file='RECONDATA/verif.dat',status='unknown')
      open (unit=32,file='RECONDATA/verifmeans.dat',
     $      status='unknown')
c
      do j=1,iabv
         if (index0(j).gt.0) then
            write (32,*) j,mean(j)
         endif
      enddo
      do i=ivermin,iyearmax
         iy = i-nlow+1
         do j = 1,iabv
            if (index0(j).gt.0) then
              aver0(iy,j)=anom0(iy,index0(j))-mean(j)
            else
              aint(iy,j)=outofbounds
              aver0(iy,j)=outofbounds
            endif
         end do
      end do
c
c
      iiii = 0
      do j=1,iabv
         if (index0(j).gt.0) then
            iiii=iiii+1
            do i=ivermin,iyearmax
              iy = i-nlow+1
              write (31,*) j,i,aver0(iy,j)
            end do
         endif
      end do       
      close (unit=31)
      close (unit=32)
c
      endif
c
c     calculate verification series and cross-validation stats
c
      iymin = ivermin-nlow+1
      iymax = ivermax-nlow+1
      iy0 = iyearmax-nlow+1
c
c     resample reconstructed nhem, global averages
c     based on sparse early grid
c
c     do not recalculate nino3 since it doesn't exist
c     in the sparse grid 
c
c     "ninod" is a dummy placeholder
c
      call averages(iabv,iymin,iy0,1082,nmax0,aint,
     $    cslat,newlon,newlat,nhemc,globc,ninod,detrc)
c
c
c     determine verification nhem,global,nino3
c     (except nino3--see above)
c
c     first use all early (land+ocean) data
c
      call averages(iabv,iymin,iy0,1082,nmax0,aver0,
     $    cslat,newlon,newlat,nhemv,globv,ninod,detrv)
c
c
c     calculate fractional variances explained relative to
c     raw seasonal/annual averages
c
      varrawtot = zero
      amsevertot = zero
      do i=1,iabv 
c
c       make sure gridpoint is in verification set
c
        if (index0(i).gt.0) then
c
        varrawgp(i)=zero
        amsevergp(i)=zero
        do iy=iymin,iymax
         j = iy-iymin+1
         series1(j) = aver0(iy,i)
         series2(j) = aint(iy,i)
         varrawgp(i)=varrawgp(i)+series1(j)**2
         amsevergp(i)=amsevergp(i)+(series1(j)-series2(j))**2
         amsevertot = amsevertot + amsevergp(i)
         varrawtot = varrawtot + varrawgp(i)
c
        end do
c
        endif
c
      end do
c
c     read in soi as "ninov"
c
      iysoi2 = ivermax
      asum1 = zero
      asum2 = zero
      nlennino = iysoi2-iysoi1+1
      open (unit=7,file='winsoi.dat',status='old')
      do i=1,nlennino
         iyear = iysoi1+i-1
         iy = iyear-nlow+1
         read (7,*) idum,soi(iy)
         enso1(iy)=soi(iy)
         asum1 = asum1+soi(iy)
         enso2(iy)=ninoc(iy)
         asum2 = asum2+ninoc(iy)
      end do
      asum1 = asum1/float(nlennino)
      asum2 = asum2/float(nlennino)
      do i=1,nlennino
         iyear = iysoi1+i-1
         iy = iyear-nlow+1
         enso1(iy)=enso1(iy)-asum1
         enso2(iy)=enso2(iy)-asum2
      end do
c
      close (unit=7)
c
c
c     calculate correlations with raw seasonal/annual averages
c
      do i=1,iabv 
        varrawgp0(i)=zero
        varcalgp(i)=zero
        corrgp(i)=zero
        amean1 = 0.0
        amean2 = 0.0
        do iy=iymin,iymax
           amean1 = amean1+aver0(iy,i)
           amean2 = amean2+aint(iy,i)
        end do
        amean1 = amean1/float(iymax-iymin+1)
        amean2 = amean2/float(iymax-iymin+1)
        do iy=iymin,iymax
         varrawgp0(i)=varrawgp0(i)+(aver0(iy,i)-amean1)**2
         varcalgp(i)=varcalgp(i)+(aint(iy,i)-amean2)**2
         corrgp(i)=corrgp(i)+(aver0(iy,i)-amean1)
     $         *(aint(iy,i)-amean2)
        end do
        if (varcalgp(i)*varrawgp0(i).gt.zero) then
           corrgp(i)=corrgp(i)/sqrt(varcalgp(i)*varrawgp0(i))
        else
           corrgp(i)=1.0
        endif
      end do
c
c
c
      amean1 = zero
      amean2 = zero
      amean3 = zero
      amean4 = zero
      varverglob = zero
      varvernhem = zero
      varcalglob = zero
      varcalnhem = zero
      corrglob = zero
      corrnhem = zero
      do iy=iymin,iymax
         amean1 = amean1+globv(iy)
         amean2 = amean2+nhemv(iy)
         amean3 = amean3+globc(iy)
         amean4 = amean4+nhemc(iy)
      end do
      amean1 = amean1/float(iymax-iymin+1)
      amean2 = amean2/float(iymax-iymin+1)
      amean3 = amean3/float(iymax-iymin+1)
      amean4 = amean4/float(iymax-iymin+1)
      do iy=iymin,iymax
         varcalglob = varcalglob + (globc(iy)-amean3)**2
         varcalnhem = varcalnhem + (nhemc(iy)-amean4)**2
         varverglob = varverglob + (globv(iy)-amean1)**2
         varvernhem = varvernhem + (nhemv(iy)-amean2)**2
         corrglob = corrglob + (globv(iy)-amean1)
     $     *(globc(iy)-amean3)
         corrnhem = corrnhem + (nhemv(iy)-amean2)
     $     *(nhemc(iy)-amean4)
      end do
      corrglob = corrglob/sqrt(varverglob*varcalglob)
      corrnhem = corrnhem/sqrt(varvernhem*varcalnhem)
c
c
      varenso1 = zero
      varenso2 = zero
      corrvernino = zero
      varrawglob = zero
      varrawnhem = zero
      varrawdetr = zero
      amseverglob = zero
      amsevernhem = zero
      amseverdetr = zero
      do iy=iymin,iymax
         j = iy-iymin+1
         iyear = ivermin+j-1
c
           varrawglob = varrawglob + globv(iy)**2
           varrawnhem = varrawnhem + nhemv(iy)**2
           varrawdetr = varrawdetr + detrv(iy)**2
           amseverglob = amseverglob + (globv(iy)-globc(iy))**2
           amsevernhem = amsevernhem + (nhemv(iy)-nhemc(iy))**2
           if ((iyear.ge.iysoi1).and.(iyear.le.iysoi2)) then
           corrvernino = corrvernino + enso1(iy)*enso2(iy)
           varenso1 = varenso1+enso1(iy)**2
           varenso2 = varenso2+enso2(iy)**2
           endif
           amseverdetr = amseverdetr + (detrv(iy)-detrc(iy))**2
      end do
      corrvernino = corrvernino/sqrt(varenso1*varenso2)
c
      open (unit=9,file='betas-verif1.out',status='unknown')
      write (9,*) 'globe:        ',one-amseverglob/varrawglob
      write (9,*) 'nhem:         ',one-amsevernhem/varrawnhem
      write (9,*)
      write (9,*) 'nino3/soi r   r^2: ',corrvernino,
     $        corrvernino**2
      write (9,*)
      write (9,*) 'multivariate: ',one-amsevertot/varrawtot
      write (9,*)
      write (9,*) 'gridpoints:'
      do i=1,iabv
         write (9,*) i,one-amsevergp(i)/varrawgp(i)
      end do
      close (unit=9)
c
      open (unit=9,file='corrs-verif1.out',status='unknown')
c
      write (9,*) 'globe:        ',corrglob,corrglob**2
      write (9,*) 'nhem:         ',corrnhem,corrnhem**2
      write (9,*)
      write (9,*) 'gridpoints: '
      do i=1,iabv
         write (9,*) i,corrgp(i),corrgp(i)**2
      end do
c
      close (unit=9)
c
c
      open (unit=34,file='RECONDATA/nhem-ver.out',status='unknown')
      open (unit=35,file='RECONDATA/glob-ver.out',status='unknown')
      open (unit=36,file='RECONDATA/nino-ver.out',status='unknown')
      open (unit=37,file='RECONDATA/detr-ver.out',status='unknown')
c
      do iyear = ivermin,iyearmax
         iy = iyear-nlow+1
         j = iyear-ivermin+1
         write (34,*) iyear,nhemv(iy),nhemc(iy)
         write (35,*) iyear,globv(iy),globc(iy)
         write (36,*) iyear,enso1(iy),enso2(iy)
         write (37,*) iyear,detrv(iy),detrc(iy)
      end do
      close (unit=34)
      close (unit=35)
      close (unit=36)
      close (unit=37)
c
c     endif
c
c     now perform verification using long historical
c     records and Quinn index, if indicated
c
c
      if (iverif.ne.1) then
c
c     first, read in actual Quinn chronology.
c
c     scale:
c 
c      1 = M-     4=S
c      2 = M      5=S+
c      3 = M+     6=VS
c
      open (unit=4,file='quinn.dat',status='old')
c
      ii = 0
      do i=1,nmax0
         read (4,*,end=707) iyear,qqq
         iy = iyear-nlow+1
         quinn(iy) = qqq
         if (quinn(iy).gt.zero) then
           ii = ii+1
           quinnsort(ii)=quinn(iy)
         endif
      end do
707   close (unit=4)
      nsort = ii
c
c     now create a *ranked* Quinn INDEX of El Nino events
c     based on deciles of calibration period NINO3
c
c
c     sort the quinn index values in increasing
c     order
c
      do i=1,nsort
         do j=i+1,nsort
            if (quinnsort(j).lt.quinnsort(i)) then
               adum = quinnsort(j)
               quinnsort(j)=quinnsort(i)
               quinnsort(i)=adum
            endif
         end do
      end do
c
c     now determine the 10-iles
c
      do  j=1,10
        jj = j-1
        jindex = 1+jj*nsort/10
        quinnthresh(j)=quinnsort(jindex)
      end do
c
c     now determine the "raw" ranked quinn
c     series from 1902-1980
c 
      do i=1,ntrain
        qrankraw(i)=0.0
        iyear = itrainmin+i-1
        iy = iyear-nlow+1
        do j=1,10
         if (quinn(iy).ge.quinnthresh(j)) 
     $      qrankraw(iy)=float(j)
        end do
      end do
c
c     determine the reconstructed ranked quinn series
c
      do i=ireconmin,ireconmax
         iy = i-nlow+1
         do j=1,10
            if (quinn(iy).ge.quinnthresh(j))
     $        qrankcal(iy) = float(j)
         end do
      end do
c
c     now compute an index of ranked warm events based on
c     the deciles of the positive reconstructed NINO3 events
c
      do i=1,ntrain
         ninosort(i)=ninor(i)
      end do
c
c     first, select out only neutral/warm events in
c     the raw data
c
      ii = 0
      do i=1,ntrain
         if (ninor(i).ge.zero) then
            ii = ii+1
            ninosort(ii)=ninor(i)
         endif
      end do
      nsort = ii
c
c     sort the raw positive NINO3 values in increasing
c     order
c
c
c     sort the raw nino3 index in increasing values
c
      do i=1,nsort
         do j=i+1,nsort
            if (ninosort(j).lt.ninosort(i)) then
               adum = ninosort(j)
               ninosort(j)=ninosort(i)
               ninosort(i)=adum
            endif
         end do
      end do
c
c     now determine the 10-iles
c
      do  j=1,10
        jj = j-1
        jindex = 1+jj*nsort/10
        ninothresh(j)=ninosort(jindex)
      end do
c
c     now determine the associated "raw" pseudo-quinn
c     series from 1902-19xx
c 
      do i=1,ntrain
        pquinnraw(i)=0.0
        do j=1,10
         if (ninor(i).ge.ninothresh(j)) 
     $      pquinnraw(i)=float(j)
        end do
      end do
c
c     determine the reconstructed pseudu-quinn series
c
      do i=ireconmin,ireconmax
         iy = i-nlow+1
         do j=1,10
            if (ninoc(iy).gt.ninothresh(j))
     $        pquinncal(iy) = float(j)
         end do
      end do
c
c     first, compute correlation between ranked Quinn
c     index and raw and reconstructed ranked NINO3 warm
c     event series during 
c     training period
c
      do i=1,ntrain
         iyear = itrainmin+i-1
         iy = iyear-nlow+1
         series1(i) = qrankcal(iy)
         series2(i) = pquinnraw(i)
         series3(i) = pquinncal(iy)
      end do
c
c     calculate calibration congruences
c
      corrraw = 0.0
      corrcal = 0.0
      var1 = 0.0
      var2 = 0.0 
      var3 = 0.0
      do i=1,ntrain
         var1 = var1+series1(i)**2
         var2 = var2+series2(i)**2
         var3 = var3+series3(i)**2
         corrraw = corrraw + series1(i)*series2(i)
         corrcal = corrcal + series1(i)*series3(i)
      end do
      corrraw = corrraw/sqrt(var1*var2)
      corrcal = corrcal/sqrt(var1*var3)
c     
      open (unit=8,file='betas-recon2.out',status='unknown')
      open (unit=9,file='betas-verif2.out',status='unknown')
c
      write (9,*) 'congruences: g  g^2'
      write (8,*) 'congruences: g  g^2'
      write (8,*) '1902-',nhigh,' Quinn/NINO3 raw: ',corrraw,
     $    corrraw**2
      write (8,*) '1902-',nhigh,' Quinn/NINO3 cal: ',corrcal,
     $    corrcal**2
c
c     now calculate verification congruences
c
      corrraw = 0.0
      corrcal = 0.0
      var1 = 0.0
      var2 = 0.0 
      var3 = 0.0
      iquinnmin = 1525
      iyrstart = max(iquinnmin,ireconmin)
c
      do iyear=iyrstart,ivermax
         iy = iyear-nlow+1
         i = iyear-iyrstart+1
         series1(i)=qrankcal(iy)
         series2(i)=pquinncal(iy)
      end do
      ntot = ivermax-iyrstart+1
      do i=1,ntot
         var1 = var1+series1(i)**2
         var2 = var2+series2(i)**2
         corrcal = corrcal + series1(i)*series2(i)
      end do
      corrcal = corrcal/sqrt(var1*var2)
      write (9,*) iyrstart,'-',ivermax,' Quinn/NINO3 cal: ',corrcal,
     $    corrcal**2
c
      open (unit=34,
     $    file='RECONDATA/qrank-raw.dat',status='unknown')
      do i=1,ntrain
         iyear = itrainmin+i-1
         iy = iyear-nlow+1
         write (34,343) iyear,pquinnraw(i),qrankraw(iy)
      end do
      close (unit=34)
c
      open (unit=34,
     $    file='RECONDATA/qrank-cal.dat',status='unknown')
      do i=ireconmin,ireconmax
         iyear = i
         iy = i-nlow+1
         write (34,343) iyear,pquinncal(iy),qrankcal(iy)
      end do
      close (unit=34)
343   format (i6,2f5.1)
c
c     now calculate monte carlo confidences if indicated
c
c     determine the 90,95,99% points of the distribution
c     of g^2 for the quinn index and time-randomized versions
c     of the reconstructed index
c
      if (imonte.eq.1) then
c
c
      niter = itermax
      write (6,*) 'performing ',niter,' bootstrap realizations...'
c
      do iter=1,niter
c
c     produce randomized pquinncal sequence
c
      do i=1,ntot
c
         iyear =  iyrstart+i-1
         iy = iyear-nlow+1
c
8181     iseq(i)=rand(iseed)*ntot+1
         if (iseq(i).gt.ntot) goto 8181
         do j=1,i-1
            if (iseq(j).eq.iseq(i)) goto 8181
         end do
         iy2 = iyrstart+iseq(i)-nlow
         pquinnrand(iy)=pquinncal(iy2)
      end do
c
c     now correlate the randomized pquinncal sequence
c     with the Quinn sequence...
c
c     now calculate verification congruences
c
      corr = 0.0
      var1 = 0.0
      var2 = 0.0 
c
      do iyear=iyrstart,ivermax
         iy = iyear-nlow+1
         i = iyear-iyrstart+1
         series1(i)=qrankcal(iy)
         series2(i)=pquinnrand(iy)
      end do
      ntot = ivermax-iyrstart+1
      do i=1,ntot
         var1 = var1+series1(i)**2
         var2 = var2+series2(i)**2
         corr = corr + series1(i)*series2(i)
      end do
c
c     store the correlation
c
      corrrand(iter) = corr**2/(var1*var2)
c
c
      end do
c
c     now calculate the distribution of the random
c     congruences
c
c     first sort the congruences by value
c
      do i=1,niter
         do j=i+1,niter
            if (corrrand(j).lt.corrrand(i)) then
                adum = corrrand(j)
                corrrand(j)=corrrand(i)
                corrrand(i)=adum
            endif
         end do
      end do
c
c     write out the distribution
c
      write (9,*) 'confidence levels for g^2'
      write (9,*) '99% ',corrrand(99*niter/100)    
      write (9,*) '95% ',corrrand(95*niter/100)    
      write (9,*) '90% ',corrrand(90*niter/100)
      write (9,*) '85% ',corrrand(85*niter/100)
      write (9,*) '80% ',corrrand(80*niter/100)
      write (9,*) '70% ',corrrand(70*niter/100)    
      write (9,*) '50% ',corrrand(50*niter/100)    
c
      endif
c
c
c     now determine MULT verification statistics based
c     on very long historical temperature records
c
c     use 1820-ivermax period during which  10 temperature
c     gridpoint records are available
c
      write (6,*) 'reading in long historical temperature records...'
c
      nhistor = 10
c
443   format (1a24)
      open (unit=1,
     $  file='/p0/tape/big2/DATA/PROXY-ANNUAL/names-longtemp',
     $  status='old')
      open (unit=2,
     $  file='/p0/tape/big2/DATA/PROXY-ANNUAL/temp-1820.loc',
     $  status='old')
      do i=1,nhistor
         read (1,443) nome2(i)
      end do
      close (unit=1)
c
      ihistor = 0
      do j=1,nhistor
         flnm2 = '/p0/tape/big2/DATA/PROXY-ANNUAL/'//nome2(j)
         ihistor = ihistor+1
         read (2,*) idum,alat,alon
         alath(ihistor)=alat
         alonh(ihistor)=alon
         open (unit=29,file=flnm2,status='old')
         iflag = 0
         do i=1,nmax0
            read (29,*,end=99) ayear,adat
            iyr = int(ayear+0.5)
            if (iflag.eq.0) then
                istart(ihistor)=iyr
                iflag = 1
            endif
            if ((iyr.ge.nlow).and.(iyr.le.nhigh)) then
               iyear = iyr-nlow+1
               ahist(iyear,ihistor)=adat
            endif
         end do
99       close (unit=29)
         istart(ihistor)=max(istart(ihistor),ireconmin)
888   end do    
      nhistor = ihistor
c
      write (6,*) 'read in ',nhistor,
     $' long historical temp records'
c
c     check for which long historical "gridpoints" match a gridpoint in 
c     the training/reconstruction set
c
      do i=1,iabv
         imatch(i)=0
      end do
      do j=1,nhistor
         lon1 = alonh(j)
         lat1 = alath(j)
c
         do i=1,iabv
c
           lon2 = newlon(i)
           lat2 = newlat(i)
c 
           if (  (abs(lon1-lon2).lt.half).and.
     $            (abs(lat1-lat2).lt.half)  ) then
              imatch(j)=i
           endif
c
          end do
      end do
c
c
c     now remove 1902-19XX reference mean and calculate
c     verification betas
c
c
      do j=1,nhistor
         amean = zero
         do iyr=itrainmin,itrainmax
           iy =iyr-nlow+1
           amean = amean + ahist(iy,j)
         end do
         amean = amean/float(nraw)
         do iyr=istart(j),irawmax
            iy = iyr-nlow+1
            ahist(iy,j)=ahist(iy,j)-amean
         end do
      end do
c
c
      varrawtot = zero
      varcaltot = zero
      amsevertot = zero
      amsecaltot = zero
      do j=1,nhistor
        varrawgp(j)=zero
        amsevergp(j)=zero
        varcalgp(j)=zero
        amsecalgp(j)=zero
        jj  = imatch(j)
        do iyr=istart(j),ivermax
         iy = iyr-nlow+1
         varrawgp(j)=varrawgp(j)+ahist(iy,j)**2
         amsevergp(j)=amsevergp(j)+(ahist(iy,j)-aint(iy,jj))**2
         if (iyr.ge.iearlyinst) then
          varrawtot = varrawtot + ahist(iy,j)**2
          amsevertot = amsevertot + (ahist(iy,j)-
     $       aint(iy,jj))**2
         endif   
        end do
c
        do iyr=itrainmin,itrainmax
           ii = iyr-itrainmin+1
           iy = iyr-nlow+1 
           varcalgp(j)=varcalgp(j)+ahist(iy,j)**2
           amsecalgp(j)=amsecalgp(j)+(ahist(iy,j)-
     $       annual(ii,jj))**2
           varcaltot = varcaltot + ahist(iy,j)**2
           amsecaltot = amsecaltot + (ahist(iy,j)-
     $       annual(ii,jj))**2
        end do
      end do
c
      write (9,*) 'individual gridpoint betas:'
      write (9,*) 'lat   lon    beta '
      write (8,*) 'individual gridpoint betas:'
      write (8,*) 'lat   lon    beta '
      do j=1,nhistor
         write (9,421) j,istart(j),alath(j),alonh(j),
     $    one-amsevergp(j)/varrawgp(j)
      end do
      do j=1,nhistor
         write (8,422) j,alath(j),alonh(j),
     $    one-amsecalgp(j)/varcalgp(j)
      end do
      write (9,*)
      write (9,*) 'multivariate beta: (start=',iearlyinst,') ',
     $    one-amsevertot/varrawtot
      write (8,*)
      write (8,*) 'multivariate beta: ',
     $    one-amsecaltot/varcaltot
c
      close (unit=8)
      close (unit=9)
c
421   format (2i6,2f8.1,f7.3)
422   format (i6,2f8.1,f7.3)
      endif
c
      stop
      end
c
c
      SUBROUTINE averages(iabv,nyear0,nyear,mmax,nmax,field,
     $    cs,lon,lat,nhem,glob,nino,detr)
c
      parameter (zero=0.d0,outbound=-900.d0)
      integer iabv, nyear,nyear0,nn
      real *8 field(nmax,mmax)
      real *8 cs(mmax),lon(mmax),lat(mmax)
      real *8 nhem(nmax),glob(nmax),nino(nmax),detr(nmax)
      real *8 sum1,sum2,sum5,sum6,slope,ainter
c
c
      do i=1,nyear
         nino(i)=zero
         nhem(i)=zero
         glob(i)=zero
         nin3 = 0
         area = zero
         areanh = zero
c
         do j=1,iabv
           if (field(i,j).gt.outbound) then
c
           if ((lat(j).ge.-5.0).and.(lat(j).le.5.0)
     $     .and.(lon(j).ge.210.0).and.(lon(j).le.270.0) ) then
            nin3=nin3+1
            nino(i) = nino(i) + field(i,j)
           endif
c
           if (lat(j).gt.zero) then
            nhem(i) = nhem(i)+ 
     $          cs(j)*field(i,j)
            areanh = areanh+cs(j)
           endif
c
           glob(i)=glob(i)+cs(j)*field(i,j)
           area = area + cs(j)
c
           endif
         end do
c
         if (nin3.eq.0) then
           nino(i)=-999.0
         else
           nino(i)=nino(i)/float(nin3)
         endif
         nhem(i)=nhem(i)/areanh
         glob(i)=glob(i)/area
      end do
c
c     determine detrended nhem series
c
      sum1 = 0.0
      sum2 = 0.0
      sum5 = 0.0
      sum6 = 0.0
      do i=nyear0,nyear
         sum1 = sum1 + float(i)
         sum2 = sum2 + float(i**2)
         sum5 = sum5 + nhem(i)
         sum6 = sum6 + nhem(i)*float(i)
      end do
      nn = nyear-nyear0+1
c
      slope = (nn*sum6-sum1*sum5)
     $       /(nn*sum2-sum1**2)
      ainter = (sum2*sum5-sum1*sum6)/
     $      (nn*sum2-sum1**2)
c
      do i=1,nyear
         detr(i)=nhem(i)-slope*float(i)
     $      - ainter
      end do
      return
      end
c
c
      SUBROUTINE CSVD (A,MMAX,NMAX,M,N,IP,NU,NV,S,U,V)
C
C   Singular value decomposition of an  M by N  complex matrix  A,
C   where M .GT. N .  The singular values are stored in the vector
C   S. The first NU columns of the M by M unitary matrix U and the
C   first NV columns of the N by N unitary matrix V  that minimize
C   Det(A-USV*) are also computed.
C
C
C         P.A. Businger and G.H. Golub, "Singular Value Decomposition
C         of a Complex Matrix," Communications of the ACM, vol. 12,
C         pp. 564-565, October 1969.
C
C   This algorithm is reprinted by permission, Association for
C   Computing Machinery; copyright 1969.
C
      COMPLEX *16 A(MMAX,NMAX),U(MMAX,MMAX),V(NMAX,NMAX),Q,R
      REAL *8 S(NMAX),B(100000),C(100000),T(100000),
     $   zero,one,ETA,TOL,
     $        EPS
      ETA = 1.2d-7
      TOL = 2.4d-32
      zero = 0.d0
      one = 1.d0
      NP=N+IP
      N1=N+1
C   Householder reduction
      C(1)=zero
      K=1
10    K1=K+1
C   Elimination of A(I,K) , I=K+1,...,M
      Z=zero
      DO 20 I=K,M
20      Z=Z+REAL(A(I,K))**2+dIMAG(A(I,K))**2
      B(K)=zero
      IF (Z .LE. TOL)  GO TO 70
      Z=SQRT(Z)
      B(K)=Z
      W=CdABS(A(K,K))
      Q=dcmplx(one,zero)
      IF (W .NE. zero)  Q=A(K,K)/W
      A(K,K)=Q*(Z+W)
      IF (K .EQ. NP)  GO TO 70
      DO 50 J=K1,NP
        Q=dcmplx(zero,zero)
        DO 30 I=K,M
30        Q=Q+dCONJG(A(I,K))*A(I,J)
        Q=Q/(Z*(Z+W))
        DO 40 I=K,M
40        A(I,J)=A(I,J)-Q*A(I,K)
50      CONTINUE
C   Phase transformation
      Q=-dCONJG(A(K,K))/CdABS(A(K,K))
      DO 60 J=K1,NP
60      A(K,J)=Q*A(K,J)
C   Elimination of A(K,J) , J=K+2,...,N
70    IF (K .EQ. N)  GO TO 140
      Z=zero
      DO 80 J=K1,N
80      Z=Z+REAL(A(K,J))**2+dIMAG(A(K,J))**2
      C(K1)=zero
      IF (Z .LE. TOL)  GO TO 130
      Z=SQRT(Z)
      C(K1)=Z
      W=CdABS(A(K,K1))
      Q=dcmplx(one,zero)
      IF (W .NE. zero)  Q=A(K,K1)/W
      A(K,K1)=Q*(Z+W)
      DO 110 I=K1,M
        Q=dcmplx(zero,zero)
        DO 90 J=K1,N
90        Q=Q+dCONJG(A(K,J))*A(I,J)
        Q=Q/(Z*(Z+W))
        DO 100 J=K1,N
100       A(I,J)=A(I,J)-Q*A(K,J)
110     CONTINUE
C   Phase transformation
      Q=-dCONJG(A(K,K1))/CdABS(A(K,K1))
      DO 120 I=K1,M
120     A(I,K1)=A(I,K1)*Q
130   K=K1
      GO TO 10
C   Tolerance for negligible elements
140   EPS=zero
      DO 150 K=1,N
        S(K)=B(K)
        T(K)=C(K)
150   EPS=DMAX1(EPS,S(K)+T(K))
      EPS=EPS*ETA
C   Initialization of U and V
      IF (NU .EQ. 0)  GO TO 180
      DO 170 J=1,NU
        DO 160 I=1,M
160       U(I,J)=dcmplx(zero,zero)
170     U(J,J)=dcmplx(one,zero)
180   IF (NV .EQ. 0)  GO TO 210
      DO 200 J=1,NV
        DO 190 I=1,N
190       V(I,J)=dcmplx(zero,zero)
200     V(J,J)=dcmplx(one,zero)
C   QR diagonalization
210   DO 380 KK=1,N
        K=N1-KK
C   Test for split
220     DO 230 LL=1,K
          L=K+1-LL
          IF (ABS(T(L)) .LE. EPS)  GO TO 290
          IF (ABS(S(L-1)) .LE. EPS)  GO TO 240
230       CONTINUE
C   Cancellation of B(L)
240     CS=zero
        SN=one
        L1=L-1
        DO 280 I=L,K
          F=SN*T(I)
          T(I)=CS*T(I)
          IF (ABS(F) .LE. EPS)  GO TO 290
          H=S(I)
          W=SQRT(F*F+H*H)
          S(I)=W
          CS=H/W
          SN=-F/W
          IF (NU .EQ. 0)  GO TO 260
          DO 250 J=1,N
            X=REAL(U(J,L1))
            Y=REAL(U(J,I))
            U(J,L1)=dCMPLX(X*CS+Y*SN,0.)
250         U(J,I)=dCMPLX(Y*CS-X*SN,0.)
260       IF (NP .EQ. N)  GO TO 280
          DO 270 J=N1,NP
            Q=A(L1,J)
            R=A(I,J)
            A(L1,J)=Q*CS+R*SN
270         A(I,J)=R*CS-Q*SN
280       CONTINUE
C   Test for convergence
290     W=S(K)
        IF (L .EQ. K)  GO TO 360
C   Origin shift
        X=S(L)
        Y=S(K-1)
        G=T(K-1)
        H=T(K)
        F=((Y-W)*(Y+W)+(G-H)*(G+H))/(2.d0*H*Y)
        G=SQRT(F*F+one)
        IF (F .LT. zero)  G=-G
        F=((X-W)*(X+W)+(Y/(F+G)-H)*H)/X
C   QR step
        CS=one
        SN=one
        L1=L+1
        DO 350 I=L1,K
          G=T(I)
          Y=S(I)
          H=SN*G
          G=CS*G
          W=SQRT(H*H+F*F)
          T(I-1)=W
          CS=F/W
          SN=H/W
          F=X*CS+G*SN
          G=G*CS-X*SN
          H=Y*SN
          Y=Y*CS
          IF (NV .EQ. 0)  GO TO 310
          DO 300 J=1,N
            X=REAL(V(J,I-1))
            W=REAL(V(J,I))
            V(J,I-1)=dCMPLX(X*CS+W*SN,0.)
300         V(J,I)=dCMPLX(W*CS-X*SN,0.)
310       W=SQRT(H*H+F*F)
          S(I-1)=W
          CS=F/W
          SN=H/W
          F=CS*G+SN*Y
          X=CS*Y-SN*G
          IF (NU .EQ. 0)  GO TO 330
          DO 320 J=1,N
            Y=REAL(U(J,I-1))
            W=REAL(U(J,I))
            U(J,I-1)=dCMPLX(Y*CS+W*SN,0.)
320         U(J,I)=dCMPLX(W*CS-Y*SN,0.)
330       IF (N .EQ. NP)  GO TO 350
          DO 340 J=N1,NP
            Q=A(I-1,J)
            R=A(I,J)
            A(I-1,J)=Q*CS+R*SN
340         A(I,J)=R*CS-Q*SN
350       CONTINUE
        T(L)=zero
        T(K)=F
        S(K)=X
        GO TO 220
C   Convergence
360     IF (W .GE. zero)  GO TO 380
        S(K)=-W
        IF (NV .EQ. 0)  GO TO 380
        DO 370 J=1,N
370       V(J,K)=-V(J,K)
380     CONTINUE
C   Sort singular values
      DO 450 K=1,N
        G=-one
        J=K
        DO 390 I=K,N
          IF (S(I) .LE. G)  GO TO 390
          G=S(I)
          J=I
390       CONTINUE
        IF (J .EQ. K)  GO TO 450
        S(J)=S(K)
        S(K)=G
        IF (NV .EQ. 0)  GO TO 410
        DO 400 I=1,N
          Q=V(I,J)
          V(I,J)=V(I,K)
400       V(I,K)=Q
410     IF (NU .EQ. 0)  GO TO 430
        DO 420 I=1,N
          Q=U(I,J)
          U(I,J)=U(I,K)
420       U(I,K)=Q
430     IF (N .EQ. NP)  GO TO 450
        DO 440 I=N1,NP
          Q=A(J,I)
          A(J,I)=A(K,I)
440       A(K,I)=Q
450     CONTINUE
C   Back transformation
      IF (NU .EQ. 0)  GO TO 510
      DO 500 KK=1,N
        K=N1-KK
        IF (B(K) .EQ. zero)  GO TO 500
        Q=-A(K,K)/CdABS(A(K,K))
        DO 460 J=1,NU
460       U(K,J)=Q*U(K,J)
        DO 490 J=1,NU
          Q=dcmplx(zero,zero)
          DO 470 I=K,M
470         Q=Q+dCONJG(A(I,K))*U(I,J)
          Q=Q/(CdABS(A(K,K))*B(K))
          DO 480 I=K,M
480         U(I,J)=U(I,J)-Q*A(I,K)
490       CONTINUE
500     CONTINUE
510   IF (NV .EQ. 0)  GO TO 570
      IF (N .LT. 2)  GO TO 570
      DO 560 KK=2,N
        K=N1-KK
        K1=K+1
        IF (C(K1) .EQ. zero)  GO TO 560
        Q=-dCONJG(A(K,K1))/CdABS(A(K,K1))
        DO 520 J=1,NV
520       V(K1,J)=Q*V(K1,J)
        DO 550 J=1,NV
          Q=dcmplx(zero,zero)
          DO 530 I=K1,N
530         Q=Q+A(K,I)*V(I,J)
          Q=Q/(CdABS(A(K,K1))*C(K1))
          DO 540 I=K1,N
540         V(I,J)=V(I,J)-Q*dCONJG(A(K,I))
550       CONTINUE
560     CONTINUE
570   RETURN
      END