c  written by Vern Kousky 8/18/09
c  Program to compute anomaly correlations    
c23456789012345678901234567890123456789012345678901234567890123456789012
      parameter (imax=40,jmax=40,num=840,i1=15,i2=35,j1=12,j2=26)
      parameter (knum=21*15)
      dimension za(imax,jmax),x(knum),y(knum)
      dimension zanum(imax,jmax,num),rcor(num)
      character path*33
      data path/'/cpc/home/wd52vk/samer/satellite/'/
      open (11,file=path//'irwin-SA-daily-anoms-boxes-feb1980-2009'
     *,access='direct',form='unformatted',recl=imax*jmax*4)
      open (51,file='pattern-correlations-irwin-SA-feb1980-2009'
     *,access='direct',form='unformatted',recl=num*4)
      nrec=0
      do 1 n=1,num
      read(11,rec=n)za
      do 2 i=1,imax
      do 3 j=1,jmax
      zanum(i,j,n)=za(i,j)
 3    continue
 2    continue
 1    continue
ccc   compute cosine weighted anomalies
      dlat=1.75
      rad=3.14159/180.
      do 6 n=1,num
      do 4 j=1,jmax
      flat=(j-1)*dlat-54.18
      do 5 i=1,imax
      zanum(i,j,n)=zanum(i,j,n)*cos(flat*rad)
 5    continue
 4    continue
 6    continue
ccc   compute anomaly correlations for selected domain
      do 7 n=1,num
      ii=0
      do 8 j=j1,j2
      do 9 i=i1,i2
      ii=ii+1
      x(ii)=zanum(i,j,n)
 9    continue
 8    continue
      do 10 k=1,num
      jj=0
      do 11 j=j1,j2
      do 12 i=i1,i2
      jj=jj+1
      y(jj)=zanum(i,j,k)
 12   continue
 11   continue
      call cor(x,y,knum,r)
      rcor(k)=r*100.
 10   continue
      write(51,rec=n)rcor
      if(n.gt.140.and.n.le.168)write(6,200)rcor
 200  format(25f5.0)
 7    continue
      stop
      end
ccc
      subroutine cor(x,y,n,r)
ccc   x and y are the input time series
ccc   n is the number of points
ccc   r is the correlation between x and y
      dimension x(n),y(n)
      sum=0.
      tot=0.
      do 1 i=1,n
      sum=sum+x(i)
      tot=tot+y(i)
  1   continue
      avex=sum/float(n)
      avey=tot/float(n)
      sum=0.
      tot=0.
      do 2 i=1,n
      sum=sum+(x(i)-avex)*(x(i)-avex)
      tot=tot+(y(i)-avey)*(y(i)-avey)
  2   continue
      sdvx=sqrt(sum/float(n))
      sdvy=sqrt(tot/float(n))
      lag=0
      nmax=n-lag
      sum=0.
      do 4 i=1,nmax
      sum=sum+(x(i)-avex)*(y(i+lag)-avey)
  4   continue
      if(sdvx.eq.0..or.sdvy.eq.0.) go to 5
      r=sum/(float(nmax)*sdvx*sdvy)
      go to 3
  5   r=0.
  3   continue
      return
      end