c  written by Vern Kousky 8/18/09
c  Program to compute anomaly correlations    
c23456789012345678901234567890123456789012345678901234567890123456789012
      parameter (imax=40,jmax=40,num=930,i1=15,i2=35,j1=12,j2=26)
      parameter (knum=21*15)
      dimension za(imax,jmax),x(knum),y(imax,jmax),yref(knum)
      dimension zanum(imax,jmax,num),idate(num),rcor(num)
      character path*33
      data path/'/cpc/home/wd52vk/samer/satellite/'/
      open (11,file=path//'irwin-SA-daily-anoms-boxes-jan1980-2009'
     *,access='direct',form='unformatted',recl=imax*jmax*4)
      open (12,file=path//'irwin-comp-23cases-sebr-wet'
     *,access='direct',form='unformatted',recl=imax*jmax*4)
      open (51,file='pattern-correlations-comp-sebr-wet-jan1980-2009'
     *,access='direct',form='unformatted',recl=num*4)
      n=0
      do 99 iyr=1,30
      iyear=1979+iyr
      imon=1
      do 98 iday=1,31
      n=n+1
      idate(n)=iyear*10000+imon*100+iday
 98   continue
 99   continue
      read(12,rec=1)y
      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 20 j=1,jmax
      flat=(j-1)*dlat-54.18
      do 21 i=1,imax
      y(i,j)=y(i,j)*cos(flat*rad)
 21   continue
 20   continue
      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
      jj=0
      do 11 j=j1,j2
      do 12 i=i1,i2
      jj=jj+1
      yref(jj)=y(i,j)
 12   continue
 11   continue
      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
      call cor(x,yref,knum,r)
      rcor(n)=r*100.
 10   continue
      write(6,200)idate(n),rcor(n)
 200  format(i10,2x,f5.0)
 7    continue
      write(51,rec=1)rcor
      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