c CORREL
c Program statistically correlates the values of one grid with 
c another by sweeping across both grids with a window of user-determined
c size.  Output is a grid of correlation coefficients. Input window is 
c given as number of grid points.  A window size of 5 or 7 to start is 
c recommended.
c Coded by Tom Hildenbrand.  Modified by Tien Grauch
      dimension f1(1100,21),f2(1100,21),x(1100),af1(21,21),af2(21,21)
      dimension sumg(21),sumt(21)
      character ifil1*50,ifil2*50,ofil*50,id*56,pgm*8
      data dval/0.1701412e+39/
      write(*,'(a,$)')' Enter 1st grid to be correlated(ncol<1101): '
      read(5,10)ifil1
10    format(a50)
	write(*,'(a,$)')' Enter 2nd grid to be correlated (ncol<1101): '
	read(5,10)ifil2
	open(10,file=ifil1,status='old',form='unformatted',
     1 share='denywr')
	open(11,file=ifil2,status='old',form='unformatted',
     1 share='denywr')
      read(10)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      read(11)id,pgm,ncol1,nrow1,nz,xo1,dx1,yo1,dy1
      if(abs(xo-xo1).gt.0.001.or.abs(dx-dx1).gt.0.001) go to 999
      if(abs(yo-yo1).gt.0.001.or.abs(dy-dy1).gt.0.001) go to 999
      if(ncol.ne.ncol1.or.nrow.ne.nrow1) go to 999
      write(6,819)
819   format(' Enter no. of grid pts for window length (must be odd;',
     1  'max=21): ',$)
      read(5,*)nwind
      write(*,'(a,$)')' Enter output correlated file: '
      read(5,10)ofil
	open(12,file=ofil,status='unknown',form='unformatted')
      write(*,'(a,$)')' Enter title: '
      read(5,10)id
      pgm='CORREL**'
      write(12)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      nbeg=(nwind-0.999999)/2.+1.
      nrend=nrow-nbeg+1
      nbegm=nbeg-1
      nwindm=nwind-1
      nwinds=nwind*nwind
      do 20 i=1,nwindm
      read(10)dum,(f1(j,i),j=1,ncol)
20    read(11)dum,(f2(j,i),j=1,ncol)
      do 21 j=1,ncol
21    x(j)=dval
      do 22 i=1,nbegm
22    write(12)dum,(x(j),j=1,ncol)
      irb=nwindm
      do 210 n=nbeg,nrend
      jcr=nbeg
      irb=irb+1
      if(irb.gt.nwind)irb=1
      read(10)dum,(f1(j,irb),j=1,ncol)
      read(11)dum,(f2(j,irb),j=1,ncol)
      do 25 j=1,ncol
25    x(j)=dval
      jr=1
      j1=1
      j2=nwind
      jcb=0
30    do 40 i=1,nwind
      do 40 j=j1,j2
      if(f1(j,i).lt.1.e+30.and.f2(j,i).lt.1.e+30) go to 40
      if(j.gt.jcb)jcb=j
      if(j.eq.j2) go to 50
40    continue
      if(jcb.eq.0) go to 60
50    jr=jcb+1
      j1=j2+1
      j2=jr+nwindm
      if(j2.gt.ncol) go to 200
      jcr=jr+nbegm
      jcb=0
      go to 30
60    j1=jr
      jsa=0
      tsumg=0.
      tsumt=0.
      do 75 j=j1,j2
      jsa=jsa+1
      sumg(jsa)=0.
      sumt(jsa)=0.
      do 70 i=1,nwind
      sumg(jsa)=sumg(jsa)+f1(j,i)
70    sumt(jsa)=sumt(jsa)+f2(j,i)
      tsumg=tsumg+sumg(jsa)
75    tsumt=tsumt+sumt(jsa)
      xtsumg=tsumg/nwinds
      xtsumt=tsumt/nwinds
      do 80 i=1,nwind
      jsa=0
      do 80 j=j1,j2
      jsa=jsa+1
      af1(jsa,i)=f1(j,i)-xtsumg
80    af2(jsa,i)=f2(j,i)-xtsumt
      tsumg2=0.
      tsumt2=0.
      tsumgt=0.
      do 90 i=1,nwind
      do 90 j=1,nwind
      tsumg2=tsumg2+af1(j,i)*af1(j,i)
      tsumt2=tsumt2+af2(j,i)*af2(j,i)
90    tsumgt=tsumgt+af1(j,i)*af2(j,i)
      dc=sqrt((tsumg2)*(tsumt2))
      if(dc.eq.0.0) go to 888
      x(jcr)=tsumgt/dc
      go to 890
888   x(jcr)=0.
890   js=0
100   j2=j2+1
      if(j2.gt.ncol) go to 200
      j1=j1+1
      jcr=jcr+1
      js=js+1
      if(js.gt.nwind)js=1
      tsumg=tsumg-sumg(js)
      tsumt=tsumt-sumt(js)
      sumg(js)=0.
      sumt(js)=0.
      do 110 i=1,nwind
      if(f1(j2,i).ge.1.e+30.or.f2(j2,i).ge.1.e+30) go to 140
      sumg(js)=sumg(js)+f1(j2,i)
110   sumt(js)=sumt(js)+f2(j2,i)
      tsumg=tsumg+sumg(js)
      tsumt=tsumt+sumt(js)
      xtsumg=tsumg/nwinds
      xtsumt=tsumt/nwinds
      do 120 i=1,nwind
      jsa=0
      do 120 j=j1,j2
      jsa=jsa+1
      af1(jsa,i)=f1(j,i)-xtsumg
120   af2(jsa,i)=f2(j,i)-xtsumt
      tsumg2=0.
      tsumt2=0.
      tsumgt=0.
      do 130 i=1,nwind
      do 130 j=1,nwind
      tsumg2=tsumg2+af1(j,i)*af1(j,i)
      tsumt2=tsumt2+af2(j,i)*af2(j,i)
130   tsumgt=tsumgt+af1(j,i)*af2(j,i)
      dc=sqrt((tsumg2)*(tsumt2))
      if(dc.eq.0.0) go to 891
      x(jcr)=tsumgt/dc
      go to 100
891   x(jcr)=0.0
      go to 100
140   j1=j2+1
      j2=nwindm+j1
      if(j2.gt.ncol) go to 200
      jcb=0
      jr=j1
      jcr=jr+nbegm
      go to 30
200   write(12)dum,(x(jw),jw=1,ncol)
210   continue
      do 218 j=1,ncol
218   x(j)=dval
      do 220 i=nrend+1,nrow
220   write(12)dum,(x(j),j=1,ncol)
	close(12)
999   close(10)
      close(11)
      stop
      end


