c* written originally for 102 cds in US at 64 time levels.; but now for a grid
c* has a switch for space time reversal, check isw switch
c* another switch (iter=0, or iter=20) gives either eot, or eof
      parameter(im=144,jm=73,scale=100.,isw=0)  !isw=1/0 normal/reverse
      parameter(nss=540,ntt=54,modemax=4,iter=20)     !choices
      parameter(ns=isw*nss+(1-isw)*ntt,nt=isw*ntt+(1-isw)*nss)
c* ntt is the number of time levels (years)
c* nss is the number of locations in space 
c* modemax<=ns
      dimension help(im,jm),cosl(jm),helpr(36,15)
      dimension anom(ns,nt),anomr(ns,nt),anomp(ns,nt),t4(ns)
      dimension cor(ns,ns),e1(ns),unev(ns),t3(nss),var0(ns),dev(ns)
      dimension alph(modemax,nt),e(modemax,ns)
      dimension ipoint(ns),ipointforce(ns)
      dimension ecorversion(modemax,ns),t2(nss),t1(nss),t12(nss)
      dimension isi(nss),isj(nss)! full grid indexing
      dimension isir(nss),isjr(nss)! reduced grid indexing
c
      open(12,file='hgt3.2',form='unformatted')! input 144X73 grids
c
      pi=4.*atan(1.)
      CONV=PI/180.
      DO j=1,73
         ALAT=90.-FLOAT(j-1)*2.5
         cosl(j)=COS(ALAT*CONV)! for weighting on lat/lon grid
      END DO                                           
c* 44 and 45 are for output; 45 for grads
c     call assign("assign -N ieee -F f77 u:45")
      open(44,file='motime.out',form='unformatted')
      open(45,file='momap.out',form='unformatted')
      open(46,file='sdev.out',form='unformatted')
c
c set up a NH grid >= 20 N, every 5/10 degrees lat/lon   
      is=0
      do i=1,im,4
      do j=29,1,-2
      is=is+1
      isi(is)=i
      isj(is)=j
      isir(is)=1+(i-1)/4
      isjr(is)=1+(j-1)/2
      enddo
      enddo
c
      do idum=1,6
      read(12)i1,i2
      enddo
      do ntime=1,ntt
      read(12)help
      write(6,267)ntime,(help(i,17),i=1,41,4)
267   format(1h ,i5,20f7.2)
       is=0
       do i=1,im,4
       do j=29,1,-2
       is=is+1
       t3(is)=help(i,j)
       enddo
       enddo
      do is=1,nss
c* ts-switch-normal:next line
      if (isw.eq.1)anom(is,ntime)=t3(is)/scale! smaller numbers are more accurate on ibm
c* in case of time space reversal:
c* ts-switch: next line
      if (isw.eq.0)anom(ntime,is)=t3(is)/scale
c* ts-switch-end
        enddo
      enddo
        anomp=anom
c
c* set unexplained variance, initial
      call sdcalc(anomp,unev,cosl,isj,nss,jm,ns,nt,isw)
      var0=unev
c
      call corprep(anomp,cor,cosl,isj,nss,jm,ns,nt,isw)
c
      do 1000 nmode=1,modemax
c     do is=1,10
c     write(6,101)is,(cor(is,isv),isv=1,10)
c     enddo
c
c     write(6,111)(unev(is),is=1,ns)
      call maxcor(cor,unev,var0,dev,cosl,isj,nss,jm,ns,isb,isw)! normal switch
c     isb=ipointforce(nmode)
c     if (nmode.eq.1)isb=30
      if (nmode.le.1)then   !grads unev and dev out!
          t4=sqrt(unev)*scale
          do is=1,ns
          helpr(isir(is),isjr(is))=t4(is)/100.
          enddo
          if (isw.eq.1)write(46)helpr
          write(6,111)(t4(is),is=1,ns)
          write(6,103)
          do is=1,ns
          helpr(isir(is),isjr(is))=dev(is)*4.
          enddo
          if (isw.eq.1)write(46)helpr
          write(6,111)(dev(is)*100.,is=1,ns)
          write(6,103)
          endif
      ipoint(nmode)=isb
      if (isw.eq.1) then
      i=isi(isb)
      alon=float(i-1)*2.5
      if (alon.gt.180.)alon=alon-360.
      j=isj(isb)
      alat=90.-float(j-1)*2.5
      else
      alat=1948+(isb-1)
      alon=1949+(isb-1)
      endif
      write(6,100)nmode,isb,alat,alon
      do is=1,ns
      e1(is)=cor(isb,is)
      enddo
c
      do ntime=1,nt
      alph(nmode,ntime)=anomp(isb,ntime)
      do is=1,ns
      fak=sqrt(unev(is))/sqrt(unev(isb))
      e(nmode,is)=e1(is)*fak
      ecorversion(nmode,is)=e1(is)
c*    anomp(is,ntime)=anomp(is,ntime)-alph(nmode,ntime)*e(nmode,is)
      enddo
      enddo
C*    do is=1,ns
C*    unev(is)=unev(is)*(1.-cor(isb,is)*cor(isb,is))
C*    enddo
      write(6,111)(e(nmode,is),is=1,ns)
      write(6,103)
      write(6,111)(alph(nmode,is),is=1,nt)
c**************** extra for iteration ********************
C* check on iter
c            e(nmode,:)=anomp(:,17+nmode)
             do idum=1,iter
         y=0.
         do is=1,ns
         if (isw.eq.1)cosw=cosl(isj(is))
         if (isw.eq.0)cosw=1.
         y=y+e(nmode,is)*e(nmode,is)*cosw
         enddo
        do ntime=1,nt
        x=0.
        do is=1,ns
        if (isw.eq.1)cosw=cosl(isj(is))
        if (isw.eq.0)cosw=1.
        x=x+e(nmode,is)*anomp(is,ntime)*cosw
        enddo
        alph(nmode,ntime)=x/y
        enddo
          y=0.
          do ntime=1,nt
          if (isw.eq.0)cosw=cosl(isj(ntime))
          if (isw.eq.1)cosw=1.
          y=y+alph(nmode,ntime)*alph(nmode,ntime)*cosw
          enddo
          do is=1,ns
          x=0.
          do ntime=1,nt
          if (isw.eq.0)cosw=cosl(isj(ntime))
          if (isw.eq.1)cosw=1.
          x=x+anomp(is,ntime)*alph(nmode,ntime)*cosw
          enddo
          e(nmode,is)=x/y
          enddo
         write(6,103)
         write(6,111)(e(nmode,is),is=1,ns)
         write(6,103)
         write(6,111)(alph(nmode,is),is=1,nt)
             enddo
      do is=1,ns
      fak=sqrt(unev(is))/sqrt(unev(isb))
      ecorversion(nmode,is)=e(nmode,is)/fak
      enddo
        do is=1,ns
        anomp(is,:)=anomp(is,:)-alph(nmode,:)*e(nmode,is)
        enddo
        if (iter.eq.0)anomp(isb,:)=anom(isb,:)*0.0001
        call sdcalc(anomp,unev,cosl,isj,nss,jm,ns,nt,isw)
        call corprep(anomp,cor,cosl,isj,nss,jm,ns,nt,isw)
c*******end extra for iteration ********************       
c*      call partial(cor,e1,ns,isb)
1000  continue
C* we now have the e's and the alphas
C**********************************************
C* what follows is checking and reconstructing. Nothing essential.
c*************************************************
c* in ts-switch mode write out alphas (in that case really the ortho e's)
c     do nmode=1,modemax
c     do is=1,nss
c     t3(is)=alph(nmode,is)
c     enddo
c     write(44)t3
c     enddo
c* for a grads plot of the modes
      do nmode=1,4
      do is=1,nss
c     helpr(isir(is),isjr(is))=e(nmode,is)
      if (isw.eq.1)helpr(isir(is),isjr(is))=ecorversion(nmode,is)
      if (isw.eq.0)helpr(isir(is),isjr(is))=alph(nmode,is)*scale!reverse switch
      enddo
      write(45)helpr
      enddo
c* end 5 lines ts-switch 
C**********************************************
C* make sure that inner products of alpha and original data yields the e's
c*************************************************
      do nmode=1,modemax
      do is=1,ns
      x=0.
      y=0.
      do ntime=1,nt
      x=x+anom(is,ntime)*alph(nmode,ntime)
      y=y+alph(nmode,ntime)*alph(nmode,ntime)
      enddo
      e1(is)=x/y
      enddo
c     write(6,111)(e1(is),is=1,ns)
      enddo
C**********************************************
C* Here we build up from scratch (anomr) and reduce from the input data, one eof at a time
c*************************************************
c
      evold=0.
      do max=1,modemax
        do ntime=1,nt
        do is=1,ns
        anomr(is,ntime)=0.
        anomp(is,ntime)=anom(is,ntime)
        enddo
        enddo
      do ntime=1,nt
      do nmode=1,max
      do is=1,ns
      anomr(is,ntime)=anomr(is,ntime)+alph(nmode,ntime)*e(nmode,is)
      anomp(is,ntime)=anomp(is,ntime)-alph(nmode,ntime)*e(nmode,is)
      enddo
      enddo
      enddo
      z=0.
      x=0.
      y=0.
      u=0.
      do ntime=1,nt
      do is=1,ns
      cosw=cosl(isj(is))!normal switch lat lon grid
      if (isw.eq.0)cosw=cosl(isj(ntime))! reverse switch
      z=z+cosw
      a=anom(is,ntime)
      x=x+a*a*cosw
      y=y+anomr(is,ntime)*anomr(is,ntime)*cosw
      u=u+anomp(is,ntime)*anomp(is,ntime)*cosw
      enddo
      enddo
c     z=float(nt*ns)
      vx=x/z
      vy=y/z
      vu=u/z
      ev=(vy)/vx*100.
      ev1=(vx-vu)/vx*100.
      isb=ipoint(max)
      if (isw.eq.1) then
      i=isi(isb)
      alon=float(i-1)*2.5
      if (alon.gt.180.)alon=alon-360.
      j=isj(isb)
      alat=90.-float(j-1)*2.5
      else
      alat=1949+(isb-1)
      alon=1949+(isb-1)
      endif
      evmode=ev-evold
      evold=ev
      write(6,112)sqrt(vx)*scale,sqrt(vy)*scale,sqrt(vu)*scale,ev,ev1,
     *evmode,alat,alon,month,max
c     write(6,111)vx,vy,vu,vy+vu,vx-vy-vu
      enddo
c
C*********************************************
C* Here we check orthogonality of the e and the alph
c*************************************************
      do m1=1,5
      do m2=1,5
      if (isw.eq.1)
     *call inprod1(e,ns,cosl,isj,jm,nss,modemax,m1,m2,anorm)
      if (isw.eq.0)call inprod3(e,modemax,m1,m2,anorm,ns)
      t3(m2)=anorm
      enddo
      write(6,101)m1,(t3(m2),m2=1,5)
      enddo
c
      do m1=1,4
      do m2=1,4
      if (isw.eq.0)
     *call inprod1(alph,nt,cosl,isj,jm,nss,modemax,m1,m2,anorm)
      if (isw.eq.1)call inprod3(alph,modemax,m1,m2,anorm,nt)
      t3(m2)=anorm
      enddo
      write(6,101)m1,(t3(m2),m2=1,4)
      enddo
      do m1=1,4
      m2=m1
      if (isw.eq.0)
     *call inprod1(alph,nt,cosl,isj,jm,nss,modemax,m1,m2,anorm)
      if (isw.eq.1)call inprod3(alph,modemax,m1,m2,anorm,nt)
      t3(m2)=anorm
      enddo
c     bring variance into the time series
c       do nmode=1,modemax
c       fak=sqrt(t3(nmode))
c       alph(nmode,:)=alph(nmode,:)/fak
c       do jear=1,ntt
c       e(nmode,jear)=e(nmode,jear)*fak
c       enddo
c       enddo
c     bring variance into the time series
      do ntime=1,ntt
      if (isw.eq.1)then
      x1=alph(1,ntime)*scale
      x2=alph(2,ntime)*scale
      x3=alph(3,ntime)*scale
      x4=alph(4,ntime)*scale
      endif
      if (isw.eq.0)then
      x1=e(1,ntime)
      x2=e(2,ntime)
      x3=e(3,ntime)
      x4=e(4,ntime)
      endif
      write(44)x1
      write(44)x2
      write(44)x3
      write(44)x4
      write(6,101)ntt,x1,x2,x3,x4
      enddo
c
      do m1=1,4
      do m2=1,4
      if (isw.eq.0)
     *call inprod1(alph,nt,cosl,isj,jm,nss,modemax,m1,m2,anorm)
      if (isw.eq.1)call inprod3(alph,modemax,m1,m2,anorm,nt)
      t3(m2)=anorm
      enddo
      write(6,101)m1,(t3(m2),m2=1,4)
      enddo
100   format(1h ,2i6,3f7.1)
101   format(1h ,i6,10f7.2)
111   format(1h ,15f5.1)
112   format(1h ,8f7.2,2i5)
102   format(1h ,i6,10f7.2)
103   format(/)
120   format(2i8)
      stop
      end
c
      subroutine partial(cor,e1,n,isb)
      dimension cor(n,n),e1(n)
      do is=1,n
      do isv=is+1,n
      if (is.eq.isb.or.isv.eq.isb)then
      cor(is,isv)=0.
      goto 999
      endif
      cor(is,isv)=(cor(is,isv)-e1(is)*e1(isv))
     */sqrt(1.-e1(is)*e1(is))
     */sqrt(1.-e1(isv)*e1(isv))
999   cor(isv,is)=cor(is,isv)
      enddo
      enddo
      return
      end
c
      subroutine inprod1(e,n,cosl,isj,jm,nss,modemax,m1,m2,anorm)
c* inner product in space among eofs of length n
      dimension e(modemax,n)
      dimension cosl(jm),isj(nss) 
      x=0.
      do is=1,n
      cosw=cosl(isj(is))
      x=x+e(m1,is)*e(m2,is)*cosw
      enddo
c     write(6,100)x
      anorm=x
100   format(1h ,'ip1= ',f10.6)
      return
      end
c
      subroutine inprod2(e,anom,alph,n,ntime,m1,modemax,nt,anorm)
c* inner product in space among one eof and obs anomaly in jear ntime, both of length n
      dimension e(modemax,n),anom(n,nt)
      dimension alph(modemax,nt)
      x=0.
      y=0.
      do is=1,n
      x=x+e(m1,is)*anom(is,ntime)
      y=y+e(m1,is)*e(m1,is)
      enddo
      write(6,100)x/y,alph(m1,ntime),x/y/alph(m1,ntime)
      anorm=x/y
100   format(1h ,'ip2= ',3f10.6)
      return
      end
c
      subroutine inprod3(alph,modemax,m1,m2,anorm,nt)
c* inner product in space among timeseries of length nt
      dimension alph(modemax,nt)
      x=0.
      do ntime=1,nt
      x=x+alph(m1,ntime)*alph(m2,ntime)
      enddo
c     write(6,100)x
      anorm=x
100   format(1h ,'ip3= ',f10.6)
      return
      end
c
      subroutine maxcor(cor,unev,var0,dev,cosl,isj,nss,jm,n,isb,isw)
      dimension cor(n,n),unev(n),var0(n),dev(n)
      dimension cosl(jm),isj(nss) 
      amax=-1000.
      exvar=0.
      ww=0.
      do is=1,n
      cosw=cosl(isj(is))
      if (isw.eq.0)cosw=1.
      exvar=exvar+(var0(is)-unev(is))*cosw
      ww=ww+var0(is)*cosw
        ev=0.
        w=0.
        do isv=1,n
        cosw=cosl(isj(isv))
        if (isw.eq.0)cosw=1.
        ev=ev+cor(is,isv)*cor(is,isv)*unev(isv)*cosw
        w=w+unev(isv)*cosw
        enddo
       ev=ev/w
       dev(is)=ev
       if (ev.gt.amax)isb=is
       if (ev.gt.amax)wkeep=w
       if (ev.gt.amax)amax=ev
      enddo
      exvar=exvar/ww*100.
      amax=amax*100.
      write(6,102)isb,amax,amax*wkeep/ww,exvar
102   format(1h ,1i6,10f7.2)
      return
      end
c
      subroutine corprep(anomp,cor,cosl,isj,nss,jm,ns,nt,isw)
      dimension anomp(ns,nt),cor(ns,ns)
      dimension cosl(jm),isj(nss) 
      cor=0.0000000000001
      do is=1,ns
      do isv=1,ns
      w=0.
      do ntime=1,nt
      a=anomp(is,ntime)
      b=anomp(isv,ntime)
      cosw=1.
      if (isw.eq.0)cosw=cosl(isj(ntime))  ! reversed, ntime is space
      w=w+cosw
      cor(is,isv)=cor(is,isv)+a*b*cosw
      enddo
      cor(is,isv)=cor(is,isv)/w
      enddo
      enddo
      do is=1,ns
      do isv=is+1,ns
      cor(is,isv)=cor(is,isv)/sqrt(cor(is,is)*cor(isv,isv))
      cor(isv,is)=cor(is,isv)
      enddo
      cor(is,is)=1.
      enddo
      return
      end
c*
      subroutine sdcalc(anomp,sd,cosl,isj,nss,jm,ns,nt,isw)
      dimension anomp(ns,nt),sd(ns)
      dimension cosl(jm),isj(nss) 
      sd=0.
      do is=1,ns
      w=0.
      do ntime=1,nt
      cosw=1.
      if (isw.eq.0)cosw=cosl(isj(ntime))  ! reversed, ntime is space
      a=anomp(is,ntime)
      sd(is)=sd(is)+a*a*cosw
      w=w+cosw
      enddo
      sd(is)=sd(is)/w
      enddo
      return
      end
