c  written by Vern Kousky 11/14/12
c  Program to compute correlation between irwin (sebr) with irwin SA  
c23456789012345678901234567890123456789012345678901234567890123456789012
      parameter (imax=40,jmax=40,ndays=31,nyrs=30)
      dimension z(imax,jmax),rcor(imax,jmax),x(ndays),y(ndays)
      dimension slp(imax,jmax,ndays)
      character path*33
      data path/'/cpc/home/wd52vk/samer/satellite/'/
      open (11,file=path//'irwin-SA-daily-ave-1.75x1.75deg-boxes-y1980-2
     *009',access='direct',form='unformatted',recl=imax*jmax*4)
      open (51,file=path//'cor-x27y19-irwin-boxes-jan1980-2009',
     *access='direct',form='unformatted',recl=imax*jmax*4)
      iref=27
      jref=19
      mind=1
      maxd=31
      nrec=0
      lpyr=1980
      do 1 iyr=1,nyrs
      iyear=1979+iyr
      nday=0
      do 2 iday=mind,maxd
      nday=nday+1
      read(11,rec=iday)z
      do 3 i=1,imax
      do 4 j=1,jmax
      slp(i,j,nday)=z(i,j)
 4    continue
 3    continue
 2    continue
ccc   compute correlations
ccc   prepare series at reference point
      do 5 iday=1,nday
      x(iday)=slp(iref,jref,iday)
 5    continue
      do 6 i=1,imax
      do 7 j=1,jmax
      do 8 iday=1,nday
      y(iday)=slp(i,j,iday)
 8    continue
      call cor(x,y,nday,r)
      rcor(i,j)=r*100.
 7    continue
 6    continue
      write(6,200)iyear,lpyr,mind,maxd
 200  format(4i10)
      nrec=nrec+1
      write(51,rec=nrec)rcor
      mind=mind+365
      if(iyear.eq.lpyr)mind=mind+1
      maxd=mind+30
      if(iyear.eq.lpyr)lpyr=lpyr+4
 1    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
ccc   avey is the mean value of y
      dimension x(n),y(n)
      sum=0.
      tot=0.
      cnt=0.
      do 1 i=1,n
      if(x(i).lt.0..or.x(i).gt.330.)go to 1
      if(y(i).lt.0..or.y(i).gt.330.)go to 1
      cnt=cnt+1.
      sum=sum+x(i)
      tot=tot+y(i)
  1   continue
      avex=sum/cnt
      avey=tot/cnt
      sum=0.
      tot=0.
      cnt=0.
      do 2 i=1,n
      if(x(i).lt.0..or.x(i).gt.330.)go to 2
      if(y(i).lt.0..or.y(i).gt.330.)go to 2
      cnt=cnt+1.
      sum=sum+(x(i)-avex)*(x(i)-avex)
      tot=tot+(y(i)-avey)*(y(i)-avey)
  2   continue
      sdvx=sqrt(sum/cnt)
      sdvy=sqrt(tot/cnt)
      lag=0
      nmax=n-lag
      sum=0.
      cnt=0.
      do 4 i=1,nmax
      if(x(i).lt.0..or.x(i).gt.330.)go to 4
      if(y(i).lt.0..or.y(i).gt.330.)go to 4
      cnt=cnt+1.
      sum=sum+(x(i)-avex)*(y(i+lag)-avey)
  4   continue
      if(sdvx.eq.0..or.sdvy.eq.0.) go to 5
      r=sum/(cnt*sdvx*sdvy)
      go to 3
  5   r=0.
  3   continue
      return
      end