c    67  1         2         3         4         5         6         7 2
c* Calculation of 2.5-d covar EOFs of monthly base data anomalies (384
c* months, 1979-2010).  The grid is centered on 1.25E-358.75E by
c* 88.75S-88.75N.  EOFs are computed globally.
c
      parameter(im=144,jm=72,ntm=384,nmo=30,nsm=10368)
      real a2(im,jm),da(nsm,ntm),eof(nsm,nmo),eva(nmo)
      real pva(nmo),ets(ntm,nmo),cj(jm)
      integer j1(3),j2(3)
      character path*24,fm*11,st*3
      data path/'/user2/tsmith/stemp/sst/'/
      data fm/'unformatted'/, st/'old'/
      data j1/3,13,23/, j2/14,24,34/
      do j=1,jm
       rlat=0.017453293*(2.5*real(j)-91.25)
       cj(j)=sqrt(cos(rlat))
      end do
c
      open(51,form=fm,file='data/gpcp.mon.eof2.5.ev')
      open(52,form=fm,file='data/gpcp.mon.eof2.5.ts')
c
c* Read in the data and compute annual averages.
      open(21,form=fm,status=st,file='data/gpcp2.5.anom.1979.2010.dat')
      do nt=1,ntm
c       call motan(21,a2)
       read(21) a2
c       call smobmp(a2)
       kn=0
       do 2 j=1,jm
       do 2 i=1,im
        kn=kn+1
        da(kn,nt)=cj(j)*a2(i,j)
        if(a2(i,j).lt.-990.)then
         print*, 'spval: ',nt
         stop 'spval'
        endif
 2     continue
      end do
      close(21)
c
c* Remove the mean.
      rntm=real(ntm)
      do ns=1,nsm
       avg=0.0
       do nt=1,ntm
        avg=avg+da(ns,nt)/(rntm)
       end do
       do nt=1,ntm
        da(ns,nt)=da(ns,nt)-avg
       end do
      end do
      print*, 'Got the anomalies'
c* Get the EOFs
      call ceof(da,nsm,ntm,nmo,eof,eva,pva,ets)
c* Write out the results.
c      write(6,91) (nint(pva(nm)),nm=1,nmo)
      write(6,93) (pva(nm),nm=1,nmo)
      do nm=2,nmo
       pva(nm)=pva(nm-1)+pva(nm)
      end do
      write(6,93) (pva(nm),nm=1,nmo)
c
c      sum=0.0
c      do nm=1,nmo
c       sum=sum+eva(nm)
c      end do
c      scale=100.0/sum
c      write(6,93) (scale*eva(nm),nm=1,nmo)
c
c* Compute the scaled time series.
      do nt=1,ntm
       do nm=1,nmo
        su1=0.0
        su2=0.0
        do ns=1,nsm
         su1=su1+eof(ns,nm)*da(ns,nt)
         su2=su2+eof(ns,nm)*eof(ns,nm)
        end do
        ets(nt,nm)=su1/su2
       end do
      end do
c
      do nm=1,nmo
       do 4 j=1,jm
       do 4 i=1,im
        a2(i,j)=-999.9
 4     continue
       kn=0
       do 6 j=1,jm
       do 6 i=1,im
        kn=kn+1
        a2(i,j)=eof(kn,nm)/cj(j)
 6     continue
       write(51) a2
      end do
c
c* Write out the scaled time series.
      do nt=1,ntm
       write(52) (ets(nt,nm),nm=1,nmo)
      end do
      stop
 91   format(10(10i4,/))
 93   format(10(10f7.2,/))
      end
c
c* To compute annual averages from monthly 5-deg data.  Read data
c* from unit iu.
      subroutine motan(iu,aa)
      parameter(im=72,jm=36)
      real cg(im,jm),ag(im,jm),aa(im,jm)
      do 2 j=1,jm
      do 2 i=1,im
       cg(i,j)=0.0
       ag(i,j)=0.0
 2    continue
      do mon=1,12
       read(iu) aa
       do 4 j=1,jm
       do 4 i=1,im
        if(aa(i,j).gt.-990.)then
         cg(i,j)=cg(i,j)+1.0
         ag(i,j)=ag(i,j)+aa(i,j)
        endif
 4     continue
      end do
      do 6 j=1,jm
      do 6 i=1,im
       aa(i,j)=-999.9
       if(cg(i,j).ge.7.0) aa(i,j)=ag(i,j)/cg(i,j)
 6    continue
      return
      end
c
c* Spatial binomial smoothing.
      subroutine smobmp(a2)
      parameter(im=72,jm=36)
      real a2(im,jm),wk(im,jm),wgt(3,3),cj(jm)
      data cj/0.04,0.13,0.22,0.30,0.38,0.46,0.54,0.61,0.68,0.74,
     2        0.79,0.84,0.89,0.92,0.95,0.98,0.99,1.00,1.00,0.99,
     3        0.98,0.95,0.92,0.89,0.84,0.79,0.74,0.68,0.61,0.54,
     4        0.46,0.38,0.30,0.22,0.13,0.04/
      wgt(1,1)=0.0625
      wgt(1,3)=0.0625
      wgt(3,1)=0.0625
      wgt(3,3)=0.0625
      wgt(2,1)=0.125
      wgt(1,2)=0.125
      wgt(3,2)=0.125
      wgt(2,3)=0.125
      wgt(2,2)=0.25
      do 1 j=1,jm
      do 1 i=1,im
       wk(i,j)=a2(i,j)
 1    continue
      do 3 j0=2,jm-1
      do 3 i0=1,im
       if(wk(i0,j0).gt.-990.)then
        s1=0.0
        s2=0.0
        kj=0
        do 2 jj=j0-1,j0+1
        kj=kj+1
        ki=0
        do 2 ii=i0-1,i0+1
         ki=ki+1
         id=ii
         if(id.lt.1) id=ii+im
         if(id.gt.im)id=ii-im
         if(wk(id,jj).gt.-990.)then
          ww=wgt(ki,kj)*cj(jj)
          s1=s1+wk(id,jj)*ww
          s2=s2+ww
         endif
 2      continue
        a2(i0,j0)=s1/s2
       endif
 3    continue
      return
      end
c
c* EOF calculation.  Note: ntmx is the maximum time length set.
c  Inputs   da(nsm,ntm)...data array;
c           nsm, ntm......space and time dimensions, resp;
c           nmo...........number of modes to keep
c Outputs   ef1(nsm,nmo)..loading patterns (eigenvectors)
c           eva(nmo)......eigenvalues (ordered largest-smallest)
c           pva(nmo)......percent variance that each mode accounts for
c
      subroutine ceof(da,nsm,ntm,nmo,ef1,eva,pva,ets)
      parameter(ntmx=400,nrr=ntmx*(ntmx+1)/2)
      real da(nsm,*),ef1(nsm,nmo),eva(nmo),pva(nmo)
      real ets(ntm,nmo),vt(ntmx),aa(ntmx,ntmx),rr(nrr)
c* Calculate the matrix to pass into AHOUSE.
      do i=1,nrr
       rr(i)=0.0
      end do
      rns=real(nsm)
      l=0
      do 1 jc=1,ntm
      do 1 ic=1,jc
       l=l+1
       do i=1,nsm
        d1=da(i,jc)
        d2=da(i,ic)
        rr(l)=rr(l)+d1*d2
       end do
       rr(l) = rr(l)/rns
 1    continue
c
      trace=0.0
      l=0
      do 2 i=1,ntm
      do 2 j=1,i 
       l=l+1
       if(i.eq.j) trace=trace+rr(l)
 2    continue
      nv=10
      eps=1.E-06
      tol=1.E-30
c* REWT writes the matrix in the form that AHOUSE will know.
      call rewt(rr,ntm)
      call ahouse(eps,tol,ntm,ntmx,rr,vt,nev,aa)
      sum=0.0
c* Compute spatial eigenvalue: lamda=e*ijm/nmax,
c* calculate the % of varaince explained by each modes.
      cnv = real(nsm)/float(ntm)
      do l=1,nmo
       pva(l)=100.0*vt(l)/trace
       vt(l) = cnv*vt(l)
       eva(l)=vt(l)
      end do
      do n = 1,ntm
       sum = 0.0
       do l = 1,ntm
        sum = sum + aa(n,l)*aa(n,l)
       end do
       cvt = sqrt(sum)
       do l = 1,ntm
        aa(n,l) = aa(n,l)/cvt
       end do
      end do
c* Project back to get the eof in space.
      do 24 m=1,nmo
       do i=1,nsm
        ef1(i,m)=0.0
       end do
       do ind=1,ntm
        do i=1,nsm 
         ef1(i,m)=ef1(i,m)+aa(m,ind)*da(i,ind)
        end do
       end do
c  Normalize spatial loadings with amplitudes and s.d..
       sum = 0.0
       do i = 1,nsm
        sum = sum + ef1(i,m)*ef1(i,m)
       end do
       cnorm = sqrt(sum/eva(m))
       do i = 1,nsm
         ef1(i,m) = ef1(i,m)/cnorm
       end do
c  Compute time coefficients.
       do ind = 1,ntm
        ets(ind,m) = 0.0
        do i = 1,nsm
         ets(ind,m) = ets(ind,m) + ef1(i,m)*da(i,ind)
        end do
       end do
 24   continue
      return
      end
c
c* This routine is from Kingtse Mo (1-16-94), who got it from U.C.L.A.
C  EGNV IS EIGENFUNCTION, NEV IS OUTPUT
C   EPS AND TOL SHOULD BE CHOSEN FOR THE MACHINE BEING USED.
C   IF THE NUMBER OF SIGNIFIGANT DECIMAL DIGITS REPRESENTABLE IN THE
C   MACHINE IS SD, THEN EPS=1.E-SD. IF THE SMALLEST NUMBER REPRESENTABL
C   IN THE MACHINE IS TINY, THEN TOL=TINY/EPS.
C   THE SYMMETRIC MATRIX OF ORDER N THAT IS TO BE DECOMPOSED IS
C   STORED IN ARRAY C IN UPPER TRIANGULAR FORM.
C   EIGENVALUES ARE FOUND IN ORDER OF DECREASING SIZE. EV(1) IS THE
C   LARGEST, ETC. ALL EIGENVALUES GE EVMIN*EV(1) ARE FOUND. IT IS
C   FOOLISH TO SET EVMIN GE 1.
C   THE EIGENVECTORS ARE WRITTEN ON LOGICAL UNIT IU.
C   THE EIGENVALUES ARE STORED IN ARRAY EV.
C   THE NUMBER OF EIGENVALUES AND EIGENVECTORS FOUND IS NEV.
c* The dimension nmax is the maximum length.
c
      subroutine ahouse(eps,tol,n,num,c,ev,nev,egnv)
      parameter(nmax=400)
      dimension a(nmax),b(nmax),ev(*),p(nmax),ta(nmax),tb(nmax),
     1          v(nmax),w(nmax),y(nmax),x(nmax),egnv(num,num),c(*)
      nck=n
      umeps=1.-eps
      jskip=0
      kskip=0
      nm1=n-1
      i=1
      idm1=0
      p(1)=0.
      v(1)=0.
      w(1)=0.
    4 ip1=i+1
      nmi=n-i
      kj=idm1
      j=i
    5 jp1=j+1
      vj=v(j)
      k=j
      lj=n-j+1
      jd=kj+1
      if(i.eq.1)go to 6
      pj=p(j)
      wj=w(j)
    6 kj=kj+1
      ckj=c(kj)
      if(i.eq.1)go to 7
      if(kskip.eq.1)go to 7
      dc=-(pj*w(k)+wj*p(k))
      ckj=dc+ckj
      c(kj)=ckj
    7 if(j.gt.i)go to 14
      if(k.gt.j)go to 8
      a(i)=ckj
      k=k+1
      go to 6
    8 y(k)=0.
      v(k)=ckj
      k=k+1
      if(k.le.n)go to 6
      jskip=0
      sum=dot(v(jp1),v(jp1),lj-1)
      if(sum.le.tol)go to 10
      s=sqrt(sum)
      csd=v(jp1)
      if(csd.lt.0.)s=-s
      v(jp1)=csd+s
      c(jd+1)=v(jp1)
      h=sum+csd*s
      b(i)=-s
      go to 12
   10 b(i)=0.
      jskip=1
   12 idm1=kj
      j=jp1
      go to 5
   14 if(jskip.eq.0)go to 15
      k=k+1
      if(k.le.n)go to 6
      j=jp1
      if(j.le.n)go to 5
      go to 215
   15 y(k)=y(k)+ckj*vj
      k=k+1
      if(k.le.n)go to 6
      if(j.eq.n)go to 17
      lj1=lj-1
      jant=jd
      do 1500 iant=1,lj1
      jant=jant+1
      x(iant)=c(jant)
 1500 continue
      y(j)=y(j)+dot(x,v(jp1),lj-1)
      j=jp1
      go to 5
   17 sp=dot(v(ip1),y(ip1),nmi)/(h+h)
      do 21 j=ip1,n
      w(j)=v(j)
   21 p(j)=(y(j)-sp*v(j))/h
  215 kskip=jskip
      i=ip1
      if(i.le.nm1)go to 4
      a(n)=ckj
      b(nm1)=-b(nm1)
      b(n)=0.
      u=abs(a(1))+abs(b(1))
      do 22 i=2,n
   22 u=amax1(u,abs(a(i))+abs(b(i))+abs(b(i-1)))
      bd=u
      rbd=1./u
      do 23 i=1,n
      w(i)=b(i)
      b(i)=(b(i)/u)**2
      a(i)=a(i)/u
      v(i)=0.
   23 ev(i)=-1.
      u=1.
      ik=1
      ndim=kj
 1000 k=ik
      el=ev(k)
   24 elam=.5*(u+el)
      du=(4.*abs(elam)+rbd)*eps
      if(abs(u-el).le.du)go to 42
      iag=0
      q=a(1)-elam
      if(q.ge.0.)iag=iag+1
      do 38 i=2,n
      if(q.eq.0.)z=abs(w(i-1)/bd)/eps
      if(q.ne.0.)z=b(i-1)/q
      q=a(i)-elam-z
      if(q.ge.0.)iag=iag+1
   38 continue
      if(iag.ge.k)go to 39
      u=elam
      go to 24
   39 if(iag.eq.k)go to 41
      m=k+1
      do 40 mm=m,iag
   40 ev(mm)=elam
   41 el=elam
      go to 24
   42 continue
      elam=bd*elam
      ev(k)=elam
      if(ik.eq.1)go to 44
      if(elam.ge.ev(ik-1))ev(ik)=umeps*ev(ik-1)
   44 i=ik
      ii=1
      do 49 j=1,n
   49 y(j)=1.
   50 do 51 k=1,n
      p(k)=0.0
      tb(k)=w(k)
   51 ta(k)=bd*a(k)-ev(i)
      l=n-1
      do 57 j=1,l
      if(abs(ta(j)).lt.abs(w(j)))go to 53
      if(ta(j).eq.0.)ta(j)=eps
      f=w(j)/ta(j)
      go to 55
   53 f=ta(j)/w(j)
      ta(j)=w(j)
      t=ta(j+1)
      ta(j+1)=tb(j)
      tb(j)=t
      p(j)=tb(j+1)
      tb(j+1)=0.0
      t=y(j)
      y(j)=y(j+1)
      y(j+1)=t
   55 tb(j+1)=tb(j+1)-f*p(j)
      ta(j+1)=ta(j+1)-f*tb(j)
      y(j+1)=y(j+1)-f*y(j)
   57 continue
      if(ta(n).eq.0.)ta(n)=eps
      if(ta(n-1).eq.0.)ta(n-1)=eps
      y(n)=y(n)/ta(n)
      y(n-1)=(y(n-1)-y(n)*tb(n-1))/ta(n-1)
      l=n-2
      do 62 j=1,l
      k=n-j-1
      if(ta(k).eq.0.)ta(k)=eps
   62 y(k)=(y(k)-y(k+1)*tb(k)-y(k+2)*p(k))/ta(k)
      ay=abs(y(1))
      do 63 j=2,n
   63 ay=amax1(ay,abs(y(j)))
      do 64 j=1,n
   64 y(j)=y(j)/ay
      ii=ii+1
      if(ii.le.2)go to 50
      id=ndim-2
      do 68 j=1,l
      id=id-j-2
      m=n-j
      h=w(m-1)
      if(h.eq.0.)go to 68
      jp1=j+1
      jant=id
      do 6400 iant=1,jp1
      jant=jant+1
      x(iant)=c(jant)
 6400 continue
      t=dot(x,y(m),jp1)/(h*x(1))
      kj=id
      do 67 k=m,n
      kj=kj+1
   67 y(k)=y(k)+t*c(kj)
   68 continue
      t=dot(y,y,n)
      xnorm=sqrt(t)
      do 70 j=1,n
   70 y(j)=y(j)/xnorm
      do 213 j=1,n
  213 egnv(ik,j)=y(j)
      ik=ik+1
      if(ik.le.nck)go to 1000
      nev=i
      return
      end
c
      function dot(z,zz,n)
      dimension z(*),zz(*)
      s=0.
      do 1 i=1,n
    1  s=s+z(i)*zz(i)
      dot=s
      return
      end
c
c* A(*) has dimension nmax*(nmax+1)/2, others dimension nmax
      subroutine rewt(a,n)
      parameter(nmax=400)
      dimension egnv(nmax,nmax),a(*),irowen(nmax),irowst(nmax)
      data irowst(1),irowen(1)/1,1/
      ncoef=n*(n+1)/2
      do 30 i=2,n
       irowst(i)=irowen(i-1)+1
       irowen(i)=irowst(i)+(i-1)
 30   continue
      do 31 j=1,ncoef
       k=1
  210  if(j.gt.irowen(k)) go to 200
       ivrb=j-irowst(k)+1
       egnv(k,ivrb)=a(j)
       go to 31
  200  k=k+1
       go to 210
 31   continue
c* convert from 2-d mode to upper right triang mode;
c* load back to r
      l=1
      do 32 j=1,n
      do 33 k=j,n
       a(l)=egnv(k,j)
   33  l=l+1
   32 continue
      return
      end