      parameter(im=720,jm=360,scale=100.,isw=0)  !isw=1/0 normal/reverse
c* increase ntt by 1 when January data is in: ntt=59 for 1948-2006
      parameter(nss=27951,ntt=59,modemax=30)!choices ntt
      parameter(ns=isw*nss+(1-isw)*ntt,nt=isw*ntt+(1-isw)*nss)
      parameter(iter=100,mfinal=10)! most recent month of data; jan=1 choice iter=0(eot)/20(eof)
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),zmask(im,jm)
      dimension cosl(jm),helpr(231,121),amp(modemax,12,ntt)
      dimension amphelp(modemax,12),vec(12)
      dimension anom(ns,nt),anomr(ns,nt),anomp(ns,nt),t4(ns),t5(ns)
      dimension cor(ns,ns),e1(ns),unev(ns),t3(nss),var0(ns),dev(ns)
      dimension alph(modemax,nt),e(modemax,ns),ipoint(ns),t6(nss)
      dimension ecorversion(modemax,ns),t2(nss),t1(nss),t3mask(nss)
      dimension isi(nss),isj(nss)! full grid indexing
      dimension isir(nss),isjr(nss)! reduced grid indexing
c
      open(11,file='globalw_march_71.00',form='unformatted')! input choice
      open(12,file='global_land_sea.05x05',form='unformatted',
     &     access='direct',recl=720*360*4)! landsea mask
      open(15,file='soilpc.out_mar',form='unformatted')! choice
      open(50,file='soileof.out_mar',form='unformatted')
      open(60,file='alpha_soilNA.out_mar',form='unformatted')
c
      undef=-9.999E+20
c
c* grid is 0.25E,-.25W, at 0.5 degree spacings in E-W (as is input precip)
c     longitude = 0.25+(i-1)*0.5 or i = (lon + 0.25)/0.5, where lon in in degree E
c* grid is -89.75, 89.75, 0.5 degree increment in S-N
c     latitude = -89.75 +(j-1)*0.5 or j = ( lat + 90.25 )/ 0.5
      pi=4.*atan(1.)
      CONV=PI/180.
      DO j=1,jm
         ALAT=-89.75+FLOAT(j-1)*0.5
         cosl(j)=COS(ALAT*CONV)
      END DO                                           
c* 44 and 45 are for output; 45 for grads
      open(45,file='pra24.out',form='unformatted')
      open(46,file='sdev.out',form='unformatted')
c
      amp=-9999999.
      read(12,rec=1)zmask
c set up a grid from 12.25 to 72.25N every 0.5 degrees, and 167.25W to 52.25W.
c set up a grid from j = 205 to 325 every 0.5 degrees, and i = 386 to i = 616.
      iw=386
      ie=616
      jn=325
      js=205
      is=0
      do i=iw,ie
      do j=js,jn
      is=is+1
      isi(is)=i
      isj(is)=j
      isir(is)=1+(i-iw)
      isjr(is)=1+(j-js)
      t3mask(is)=zmask(i,j)
      enddo
      enddo
c
      do month=mfinal,mfinal!month of choice
      mf1=month
      anom=0.
      njearfinal=ntt
      do ntime=1,ntt
      write(6,500)month,ntime+1947
500   format(1h ,2i5)
       read(11)help
       t3=-9999999999.
       is=0
       do i=iw,ie
       do j=js,jn
       a=help(i,j)
       is=is+1
       t3(is)=a
       enddo
       enddo
        do is=1,nss
      if (isw.eq.1)anom(is,ntime)=t3(is)/scale! 
      if (isw.eq.0)anom(ntime,is)=t3(is)/scale
        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)
c     write(6,1111)(isi(is),is=1,nss)
c         write(6,103)
c     write(6,111)(cosl(isj(is)),is=1,nss)
c         write(6,103)
c     write(6,1111)(isj(is),is=1,nss)
c         write(6,103)
      call maxcor(cor,unev,var0,dev,cosl,isj,nss,jm,ns,isb,isw)! normale switch
      if (nmode.le.2)then   !grads unev and dev out!
          t4=sqrt(unev)*scale
c         write(6,111)(t4(is)*10.,is=1,ns)
c         write(6,103)
c         write(6,111)(dev(is)*100.,is=1,ns)
c         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=1948+(isb-1)
      endif
c     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)
citer anomp(is,ntime)=anomp(is,ntime)-alph(nmode,ntime)*e(nmode,is)
      enddo
      enddo
citer do is=1,ns
citer unev(is)=unev(is)*(1.-cor(isb,is)*cor(isb,is))
citer enddo
c     write(6,111)(e(nmode,is),is=1,ns)
c     write(6,103)
c     write(6,111)(alph(nmode,is),is=1,nt)
c**************** extra for iteration ********************
             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
c        write(6,103)
c        write(6,111)(e(nmode,is),is=1,ns)
c        write(6,103)
c        write(6,111)(alph(nmode,is),is=1,nt)
             enddo
        do is=1,ns
        anomp(is,:)=anomp(is,:)-alph(nmode,:)*e(nmode,is)
        enddo
        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
      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
      do is=1,nss
      cosw=cosl(isj(is))
      z=z+cosw
      enddo
      wkeep=z
c
      x=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
      enddo
      enddo
      akeep=x
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)
      akeep=x
      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=1948+(isb-1)
      alon=1948+(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**redo for just the last year to check for poor projection problems
      evold=0.
c     do max=1,modemax
      do max=4,4
      do ispecial=1,njearfinal! over time
        do ntime=1,nt
        do is=ispecial,ispecial
        anomr(is,ntime)=0.
        anomp(is,ntime)=anom(is,ntime)
        enddo
        enddo
      do ntime=1,nt
      do nmode=1,max
      do is=ispecial,ispecial
      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=ispecial,ispecial
      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=1948+(isb-1)
      alon=1948+(isb-1)
      endif
      evmode=ev-evold
      evmode=-99999999.
      evold=ev
      if (max.eq.1.or.max.eq.5.or.max.eq.10.or.max.eq.20)
     *write(6,112)sqrt(vx)*scale,sqrt(vy)*scale,sqrt(vu)*scale,ev,ev1,
     *evmode,alat,alon,month,max,ispecial
c     write(6,111)vx,vy,vu,vy+vu,vx-vy-vu
      enddo! mode loop
      enddo! itime loop
C*********************************************
c**redo for just the last 2 year to check for poor projection problems
      evold=0.
      do max=1,modemax
      do ispecial=njearfinal-1,njearfinal! over time
        do ntime=1,nt
        do is=ispecial,ispecial
        anomr(is,ntime)=0.
        anomp(is,ntime)=anom(is,ntime)
        enddo
        enddo
      do ntime=1,nt
      do nmode=1,max
      do is=ispecial,ispecial
      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
c     if (month.eq.mf1)then
c     if (isw.eq.0)t3(:)=anomr(ispecial,:)
c     if (isw.eq.1)t3(:)=anomr(:,ispecial)
c     call regridr(t3,helpr,isir,isjr,nss)
c     if (isw.eq.0)t3(:)=anomp(ispecial,:)
c     if (isw.eq.1)t3(:)=anomp(:,ispecial)
c     call regridr(t3,helpr,isir,isjr,nss)
c     endif
      z=0.
      x=0.
      y=0.
      u=0.
      do ntime=1,nt
      do is=ispecial,ispecial
      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=1948+(isb-1)
      alon=1948+(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,ispecial
c     write(6,111)vx,vy,vu,vy+vu,vx-vy-vu
      enddo! mode loop
      enddo! itime loop
C*********************************************
C* Here we check orthogonality of the e and the alph
c*************************************************
      do m1=1,modemax
      do m2=1,modemax
      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,modemax)
      enddo
c
      do m1=1,modemax
      do m2=1,modemax
      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,modemax)
      enddo
c
      do m1=1,modemax
      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(m1)=anorm
      enddo
c
c     write(14)alph
            if (month.eq.mf1)then
c* for a grads plot of the 1st 4 eof (raw)
      do idum=1,4
      do is=1,nss
      helpr(isir(is),isjr(is))=alph(idum,is)*scale
      if (t3mask(is).lt.1.)helpr(isir(is),isjr(is))=undef
      enddo
      write(50)helpr
      enddo
             endif
            if (month.eq.mf1)then
      x=0.
      y=0.
      do j=1,njearfinal
      x1=e(1,j)
      x2=e(2,j)
      x=x1*x2
      y=y+x
      write(6,700)j,x,y,x1,x2
700   format(1h ,'orthocheck:',i5,4f7.2)
      x3=e(3,j)
      x4=e(4,j)
      write(15)x1
      write(15)x2
      write(15)x3
      write(15)x4
      enddo
             endif
c
c     bring variance into the time series
      do nmode=1,modemax
      simvar=0.
      fak=sqrt(t3(nmode))
      alph(nmode,:)=alph(nmode,:)/fak
      do jear=1,njearfinal
      e(nmode,jear)=e(nmode,jear)*fak
      a=e(nmode,jear)
      simvar=simvar+a*a
      enddo
      t4(nmode)=simvar
      enddo
      t5=t4
      do nmode=2,modemax
      t5(nmode)=t5(nmode)+t5(nmode-1)
      enddo
      akeep=t5(modemax)
      do nmode=1,modemax
      t3(nmode)=t4(nmode)*100./akeep
      enddo
      t6=t3
      do nmode=2,modemax
      t6(nmode)=t6(nmode)+t6(nmode-1)
      enddo
C* write out eof amplitudes
      do jear=1,njearfinal
c     write(12,358)jear,month,(e(nmode,jear),nmode=1,modemax)
      write(6,357)jear+47,month,(e(nmode,jear),nmode=1,modemax)
      amp(:,month,jear)=e(:,jear)
      enddo
      t4=t4/(wkeep*float(njearfinal))
      t5=t5/(wkeep*float(njearfinal))
      write(6,357)jear,month,(t3(nmode),nmode=1,modemax)
      write(6,357)jear,month,(t6(nmode),nmode=1,modemax)
      t4=sqrt(t4)*100.
      write(6,357)jear,month,(t4(nmode),nmode=1,modemax)
      t5=sqrt(t5)*100.
      write(6,357)jear,month,(t5(nmode),nmode=1,modemax)
c     write(13)alph
      do nmode=1,5
      write(6,359)nmode,(e(nmode,jear),jear=1,njearfinal)
      enddo
      do m1=1,modemax
      do m2=1,modemax
      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,modemax)
      enddo
c
      do m1=1,modemax
      do m2=1,modemax
      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,modemax)
      enddo
c
      do idum=1,modemax
      do is=1,nss
      helpr(isir(is),isjr(is))=alph(idum,is)*scale
      enddo
      write(60)helpr! write out for cas2
      enddo
c*
      enddo! month loop
      do nmode=1,2
      do jear=2,ntt
      do m=1,12
      j=jear
      if (m.ge.mfinal)j=jear-1
      vec(m)=amp(nmode,m,j)
c* element 1-(mfinal-1) present year; mfinal-12 last year
      enddo
      write(6,360)nmode,jear,vec,mfinal
      enddo
      enddo
      do jear=2,ntt
      do m=1,12
      j=jear
      if (m.ge.mfinal)j=jear-1
      amphelp(:,m)=amp(:,m,j)
      enddo
c     write(17)amphelp! for 30 modes, the last 12 seasons
c* the 12 elements are laid out like:
C* current year m=1(djf) thru m=mfinal-1
c* last year m=mfinal thru m=12(ndj)
      enddo
100   format(1h ,2i6,3f7.1)
101   format(1h ,i6,10f7.2)
111   format(1h ,19f4.1)
1111  format(1h ,19i4)
112   format(1h ,6f7.2,2f7.1,4i5)
102   format(1h ,i6,10f7.2)
103   format(/)
120   format(2i8)
357   format(1h ,2i3,1x,30f3.0)
358   format(2i3,30f5.1)
359   format(1h ,'time series of mode:',i5,2(/,25f5.1))
360   format(1h ,2i5,12f5.1,'final month:',i5)
      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.
c     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
c*
      subroutine rms(grid00,grid24,cosl,error,nrec)
      dimension grid00(144,73),grid24(144,73),cosl(73)
      x=0.
      y=0.
      do j=2,72
      w=sqrt(cosl(j))
      do i=1,144
      a=(grid00(i,j)-grid24(i,j))*w
      x=x+a*a
      y=y+w*w
c     if (nrec.eq.1)write(6,101)i,j,grid00(i,j),grid24(i,j),a,x,y,w,nrec
      enddo
      enddo
      error=x/y
      write(6,100)cosl(37),grid00(5,39),grid24(5,39),x,y,error
     *,sqrt(error)
100   format(1h ,7f10.2)
101   format(1h ,2i5,6f10.2,i5)
      return
      end
c*
      subroutine regridr(t3,helpr,isir,isjr,nss)
      dimension helpr(231,121),t3(nss)
      dimension isir(nss),isjr(nss)! reduced grid indexing
      do is=1,nss
      helpr(isir(is),isjr(is))=t3(is)
      enddo
      return
      end
