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 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