C     PROGRA 'nana.f' to readz500windowanomalies.
c* the rmsd and acc have already been calculated and are read in.
      parameter(im=144,jm=73,m=1,ndx=31,nyr=37)
      dimension mdays(12)
      dimension anom(im,jm),anom1(im,jm),anomp(im,jm)
      dimension jear(nyr*ndx),iday(nyr*ndx),itel(nyr,ndx)
      dimension cosl(jm),rmsd(nyr*ndx,nyr*ndx)
      dimension acc(nyr*ndx,nyr*ndx)
      dimension stra(im,jm,11),strb(im,jm,11)
      dimension ace(11),rmse(11)
      dimension pace(11),prmse(11)
      integer case(300,4),mcase(300,4),iens(11)
      data mdays/31,28,31,30,31,30,31,31,30,31,30,31/
C
      open(81,file='/ptmp/wx51hd/test/z500window.jan1968_2004'
     *,form='unformatted')
c     open(82,file='/ptmp/wx51hd/test/rmsdacc.jan.sh'
      open(82,file='/ptmp/wx51hd/test/rmsdacc.jan.nh'
     *,form='unformatted')
      pi=4.*atan(1.)
      CONV=PI/180.
      DO j=1,jm
         ALAT=90.-FLOAT(j-1)*2.5
         COSL(j)=COS(ALAT*CONV)
      END DO
      ic=51
      lb=45
      le=72
      lb=25
      le=49
      lb=2
      le=29
      fak=1.
      arm=0.
      carm=0.
      iu=81
      nbase=1
      icase=0
      imcase=0
c*
      index=0
      rrmin=1000.
      do n=1,nyr
      do nd=1,ndx
      index=index+1
      jear(index)=n
      iday(index)=nd
      itel(n,nd)=index
      enddo
      enddo
c
      indexs=itel(1,1) 
      indexe=itel(nyr,ndx) 
      read(82)rmsd
      do index=indexs,indexe
      rmsd(index,index)=999.
      acc(index,index)=0.
      do idum=index+1,indexe
      acc(index,idum)=rmsd(idum,index)
      rmsd(idum,index)=rmsd(index,idum)
      acc(idum,index)=acc(index,idum)
      enddo
c     write(6,99)index,idum,rmsd(idum-1,index)
99    format(1h ,2i6,f7.1)
      enddo
c*
      xall=0.
      yall=0.
      ysall=0.
      rmina1=0.
      amaxa1=0.
      amina1=0.
      itot=0
      do 999 index=indexs,indexe,1
      if (iday(index).eq.1) then
      if (index.gt.1)
     *write(6,105)jear(index)+1966,rminav,amaxav*100.,aminav*100.
      write(6,103)
      rminav=0.
      amaxav=0.
      aminav=0.
      endif
      rmin=1000.
      amax=-1000.
      amin=1000.
      do 40 n=1,nyr
      ny=n
      do 30 nd=1,ndx
      if (n.eq.jear(index))goto 30
      itot=itot+1
      idum=itel(n,nd)
      x=rmsd(index,idum)
      xall=xall+x*x
      y=acc(index,idum)
      yall=yall+y
      ysall=ysall+y*y
      if (x.lt.rmin)nmin=ny
      if (x.lt.rmin)ndmin=nd
      if (x.lt.rmin)rmin=x
      if (y.gt.amax)nmax=ny
      if (y.gt.amax)ndmax=nd
      if (y.gt.amax)amax=y
      if (y.lt.amin)namin=ny
      if (y.lt.amin)ndamin=nd
      if (y.lt.amin)amin=y
   30 continue 
c     print *,'  data read in for ',n,ny,m,mdays(m),anom(1,55)
   40 continue
      if (rmin.lt.rrmin)iay=nmin
      if (rmin.lt.rrmin)iby=jear(index)
      if (rmin.lt.rrmin)iad=ndmin
      if (rmin.lt.rrmin)ibd=iday(index)
      if (rmin.lt.rrmin)rrmin=rmin
      if (amin.lt.-0.60)imcase=imcase+1
      if (amax.gt.0.60)icase=icase+1
      if (amax.gt.0.60)then
      case(icase,1)=jear(index)
      case(icase,2)=iday(index)
      case(icase,3)=nmax
      case(icase,4)=ndmax
      write(6,101)jear(index)+1967,iday(index),nmin+1967,ndmin,rmin,
     *nmax+1967,ndmax,amax*100.,namin+1967,ndamin,amin*100
      endif
      if (amin.lt.-0.60)then
      mcase(imcase,1)=jear(index)
      mcase(imcase,2)=iday(index)
      mcase(imcase,3)=namin
      mcase(imcase,4)=ndamin
      if (amax.le.0.60)write(6,101)jear(index)+1967,iday(index),
     *nmin+1967,ndmin,rmin,
     *nmax+1967,ndmax,amax*100.,namin+1967,ndamin,amin*100
      endif
      rminav=rminav+rmin/float(ndx)
      amaxav=amaxav+amax/float(ndx)
      aminav=aminav+amin/float(ndx)
      if (index.eq.indexe)
     *write(6,105)jear(index)+1967,rminav,amaxav*100.,aminav*100.
      if (iday(index).eq.ndx) then
      rmina1=rmina1+rminav/float(nyr)
      amaxa1=amaxa1+amaxav/float(nyr)
      amina1=amina1+aminav/float(nyr)
      endif
999   continue
      write(6,103)
      write(6,105)jear(indexe)+1967,rmina1,amaxa1*100.,amina1*100.
      tot=float(itot)
      xall=sqrt(xall/tot)
      yall=yall/tot
      ysall=sqrt(ysall/tot)
      dof=2.+1./ysall/ysall
      write(6,103)
      write(6,106)icase,imcase,itot,xall,
     *xall/sqrt(2.),yall*10000.,ysall,dof
      write(6,103)
      write(6,100)iby+1967,ibd,iay+1967,iad,rrmin
      write(6,103)
c
c* elaborate on best rms case
c
      indexa=itel(iay,iad)
      indexb=itel(iby,ibd)
      call fill(stra,anom,im,jm,indexa,1)
      call fill(strb,anom,im,jm,indexb,1)
      iby=jear(indexb)
      iay=jear(indexa)
      do idum=1,11
      ibd=iday(indexb-6+idum)
      iad=iday(indexa-6+idum)
      do i=1,im
      do j=1,jm
      anom(i,j)=strb(i,j,idum)
      anom1(i,j)=stra(i,j,idum)
      enddo
      enddo
      x=rmsdif(anom,anom1,cosl,im,jm,fak,lb,le,arm,carm)
      y=ac(anom,anom1,cosl,im,jm,fak,lb,le,sx,sy)
      write(6,100)iby+1967,ibd,iay+1967,iad,x,y*100.,sx,sy
      enddo
      write(6,103)
C
C* reevaluate all analogues
C
      do idum=1,11
      rmse(idum)=0.
      ace(idum)=0.
      prmse(idum)=0.
      pace(idum)=0.
      iens(idum)=0
      enddo
      do icas=1,icase
c     do icas=1,8
      write(6,103)
      iby=case(icas,1)
      ibd=case(icas,2)
      iay=case(icas,3)
      iad=case(icas,4)
      indexa=itel(iay,iad)
      indexb=itel(iby,ibd)
      if (indexa.le.5.or.indexb.le.5) goto 986
      if (indexa.ge.(indexe-5).or.indexb.ge.(indexe-5)) goto 986
      ac0=acc(indexb,indexa)
      ac1=(acc(indexb+1,indexa+1)+acc(indexb-1,indexa-1))/2.
      if (ac1.ge.ac0) goto 986
      call fill(stra,anom,im,jm,indexa,1)
      call fill(strb,anom,im,jm,indexb,1)
      do i=1,im
      do j=1,jm
      anomp(i,j)=stra(i,j,6)
      enddo
      enddo
      iby=jear(indexb)
      iay=jear(indexa)
      do idum=1,11
      ibd=iday(indexb-6+idum)
      ibyy=jear(indexb-6+idum)
      iad=iday(indexa-6+idum)
      iayy=jear(indexa-6+idum)
      if (ibyy.ne.iby.or.iayy.ne.iay) goto 985
      do i=1,im
      do j=1,jm
      anom(i,j)=strb(i,j,idum)
      anom1(i,j)=stra(i,j,idum)
      enddo
      enddo
      x=rmsdif(anom,anom1,cosl,im,jm,fak,lb,le,arm,carm)
      y=ac(anom,anom1,cosl,im,jm,fak,lb,le,sx,sy)
      write(6,100)iby+1967,ibd,iay+1967,iad,x,y*100.,sx,sy,icas,idum
      xp=rmsdif(anom,anomp,cosl,im,jm,fak,lb,le,arm,carm)
      yp=ac(anom,anomp,cosl,im,jm,fak,lb,le,sx,sy)
      ace(idum)=ace(idum)+y
      rmse(idum)=rmse(idum)+x*x
      pace(idum)=pace(idum)+yp
      prmse(idum)=prmse(idum)+xp*xp
      iens(idum)=iens(idum)+1
985   continue
      enddo
986   continue
      enddo
      do idum=1,11
      z=float(iens(idum))
      ace(idum)=ace(idum)/z*100.
      rmse(idum)=sqrt(rmse(idum)/z)
      pace(idum)=pace(idum)/z*100.
      prmse(idum)=sqrt(prmse(idum)/z)
      write(6,108)idum-6,iens(idum),ace(idum),rmse(idum)
     *,pace(idum),prmse(idum)
      enddo
      write(6,103)
C
C* reevaluate all anti-analogues
C
      do idum=1,11
      rmse(idum)=0.
      ace(idum)=0.
      prmse(idum)=0.
      pace(idum)=0.
      iens(idum)=0
      enddo
      do icas=1,imcase
c     do icas=1,8
      write(6,103)
      iby=mcase(icas,1)
      ibd=mcase(icas,2)
      iay=mcase(icas,3)
      iad=mcase(icas,4)
      indexa=itel(iay,iad)
      indexb=itel(iby,ibd)
      if (indexa.le.5.or.indexb.le.5) goto 987
      if (indexa.ge.(indexe-5).or.indexb.ge.(indexe-5)) goto 987
      ac0=-acc(indexb,indexa)
      ac1=-(acc(indexb+1,indexa+1)+acc(indexb-1,indexa-1))/2.
      if (ac1.ge.ac0) goto 987
      call fill(stra,anom,im,jm,indexa,1)
      call fill(strb,anom,im,jm,indexb,-1)
      do i=1,im
      do j=1,jm
      anomp(i,j)=stra(i,j,6)
      enddo
      enddo
      iby=jear(indexb)
      iay=jear(indexa)
      do idum=1,11
      ibd=iday(indexb-6+idum)
      ibyy=jear(indexb-6+idum)
      iad=iday(indexa-6+idum)
      iayy=jear(indexa-6+idum)
      if (ibyy.ne.iby.or.iayy.ne.iay) goto 988
      do i=1,im
      do j=1,jm
      anom(i,j)=strb(i,j,idum)
      anom1(i,j)=stra(i,j,idum)
      enddo
      enddo
      x=rmsdif(anom,anom1,cosl,im,jm,fak,lb,le,arm,carm)
      y=ac(anom,anom1,cosl,im,jm,fak,lb,le,sx,sy)
      write(6,100)iby+1967,ibd,iay+1967,iad,x,y*100.,sx,sy,icas,idum
      xp=rmsdif(anom,anomp,cosl,im,jm,fak,lb,le,arm,carm)
      yp=ac(anom,anomp,cosl,im,jm,fak,lb,le,sx,sy)
      ace(idum)=ace(idum)+y
      rmse(idum)=rmse(idum)+x*x
      pace(idum)=pace(idum)+yp
      prmse(idum)=prmse(idum)+xp*xp
      iens(idum)=iens(idum)+1
988   continue
      enddo
987   continue
      enddo
      do idum=1,11
      z=float(iens(idum))
      ace(idum)=ace(idum)/z*100.
      rmse(idum)=sqrt(rmse(idum)/z)
      pace(idum)=pace(idum)/z*100.
      prmse(idum)=sqrt(prmse(idum)/z)
      write(6,108)idum-6,iens(idum),ace(idum),rmse(idum)
     *,pace(idum),prmse(idum)
      enddo
100   format(1h ,4i5,4f8.1,2i5)
101   format(1h ,4i5,f8.1,2i5,f8.1,2i5,f8.1)
103   format(/)
105   format(1h ,i5,3f8.1,i5)
106   format(1h ,3i8,5f8.1)
108   format(1h ,2i8,4f8.1)
      stop
      end
c*
      function rmsdif(anom,anom1,cosl,im,jm,fak,lb,le,arm,carm)
      dimension anom(im,jm),cosl(jm)
      dimension anom1(im,jm)
      x=0.000001
      y=0.000001
      zz=0.
      do jlat=lb,le
      w=cosl(jlat)
      do i=1,144
      dif=(anom(i,jlat)-anom1(i,jlat))/fak
      difc=anom(i,jlat)/fak
      x=x+(dif*dif)*w
      y=y+(difc*difc)*w
      zz=zz+w
      enddo
      enddo
      z=zz
      if (x.gt.0.) sx=sqrt(x/z)
      if (y.gt.0.) sy=sqrt(y/z)
c     write(6,101)sx,sy,z,lb,le
      rmsdif=sx
      arm=arm+x/z
      carm=carm+y/z
100   format(1h ,2f10.3,i5)
101   format(1h ,3f7.2,3i4)
      return
      end
c*
      subroutine fill(strb,anom,im,jm,indexb,iana)
      dimension strb(im,jm,11),anom(im,jm)
c* iana=1 for regular
c* iana=-1 for antis used as forecast
      fak=float(iana)
      iu=81
      rewind(iu)
      do idum=1,indexb-6
      read(iu)anom 
      enddo
      do idum=1,11
      read(iu)anom
      do i=1,im
      do j=1,jm
      strb(i,j,idum)=anom(i,j)*fak
      enddo 
      enddo 
      enddo 
      return
      end
      function ac(anom,anom1,cosl,im,jm,fak,lb,le,sx,sy)
      dimension anom(im,jm),cosl(jm)
      dimension anom1(im,jm)
      x=0.
      y=0.
      z=0.
      zz=0.
      do jlat=lb,le
      w=cosl(jlat)
      do i=1,144
      ax=anom(i,jlat)/fak
      ay=anom1(i,jlat)/fak
      x=x+(ax*ax)*w
      y=y+(ay*ay)*w
      z=z+(ay*ax)*w
      zz=zz+w
      enddo
      enddo
      if (x.gt.0.) sx=sqrt(x/zz)
      if (y.gt.0.) sy=sqrt(y/zz)
      ac=z/zz/sx/sy
c     write(6,101)sx,sy,ac*100.,lb,le
100   format(1h ,2f10.3,i5)
101   format(1h ,3f7.2,3i4)
      return
      end
