c     Program written 7/11/12 by Vern Kousky
c     program to combine nstn files for 1948-2009.
c23456789012345678901234567890123456789012345678901234567890123456789012
      parameter (imax=40,jmax=40,ndmax=900,nyrs=30,ndmon=30)
      dimension x(imax,jmax),p(imax,jmax,ndmax),y(ndmax),z(ndmax)
      dimension fdate(imax,jmax,ndmax),kndx(ndmax),xd(imax,jmax)
      character path*33
      data path/'/cpc/home/wd52vk/samer/satellite/'/
      open (11,file=path//'irwin-SA-daily-anoms-boxes-apr1980-2009',
     *access='direct',form='unformatted',recl=imax*jmax*4)
      open (51,file=path//'rank-irwin-SA-daily-anoms-boxes-apr1980-2009'
     *,access='direct',form='unformatted',recl=imax*jmax*4)
      open (52,file=path//'rank-dates-irwin-SA-daily-anoms-boxes-apr1980
     *-2009',access='direct',form='unformatted',recl=imax*jmax*4)
      nrec=0 
      do 1 nd=1,ndmax
      read(11,rec=nd)x
      do 2 i=1,imax
      do 3 j=1,jmax
      p(i,j,nd)=x(i,j)
 3    continue
 2    continue
 1    continue
ccc   Now rank the daily irwin values at each grid point 
      do 4 i=1,imax
      do 5 j=1,jmax
      do 6 nd=1,ndmax
      y(nd)=p(i,j,nd)
 6    continue
      call order(y,ndmax,z,kndx)
      do 7 nd=1,ndmax
      p(i,j,nd)=z(nd)
 7    continue
      do 71 iday=1,ndmax
      nd=0
      do 72 iyr=1,nyrs
      iyear=1979+iyr
      do 73 id=1,ndmon
      nd=nd+1
      if(kndx(iday).ne.nd)go to 73
      fdate(i,j,iday)=iyear*100+id
      go to 71
 73   continue
 72   continue
 71   continue 
 5    continue
 4    continue
      do 8 nd=1,ndmax
      do 9 i=1,imax
      do 10 j=1,jmax
      x(i,j)=p(i,j,nd)
      xd(i,j)=fdate(i,j,nd)
 10   continue
 9    continue
      write(51,rec=nd)x
      write(52,rec=nd)xd
      write(6,200)nd,x(26,21),xd(26,21)
 200  format(i5,2f10.1)
 8    continue
      stop
      end
ccc
      subroutine order(x,kmax,y,kind)
      dimension x(kmax),y(kmax),kind(kmax)
      do 1 k=1,kmax
      y(k)=x(k)
      kind(k)=k
 1    continue
ccc   first find largest value
      maxi=1
 99   continue
      ymax=y(maxi)
      lmax=maxi
      do 2 k=maxi,kmax
      if(y(k).le.ymax)go to 2
      ymax=y(k)
      lmax=k
 2    continue
      xy=y(maxi)
      ky=kind(maxi)
      y(maxi)=ymax
      kind(maxi)=kind(lmax)
      y(lmax)=xy
      kind(lmax)=ky
      if(maxi.eq.kmax)go to 999
      maxi=maxi+1
      go to 99
 999  continue
      return
      end