c     Program written 10/11/11 by Vern Kousky
c      modified for PW 9/23/13
c23456789012345678901234567890123456789012345678901234567890123456789012
      parameter (imax=40,jmax=40,maxyr=30)
      dimension x(imax,jmax),y(imax,jmax,365),numd(12)
      dimension yy(imax,jmax,1095),cnt(imax,jmax,365)
      dimension ys365(imax,jmax)
      character path*33
      data numd/31,28,31,30,31,30,31,31,30,31,30,31/
      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//'irwin-daily-clim365-boxes-y1980-2009',
     *access='direct',form='unformatted',recl=imax*jmax*4)
      open (52,file=path//'irwin-daily-clim366-boxes-y1980-2009',
     *access='direct',form='unformatted',recl=imax*jmax*4)
      do 99 i=1,imax
      do 98 j=1,jmax
      do 97 nd=1,365
      y(i,j,nd)=0.
      cnt(i,j,nd)=0.
 97   continue
 98   continue
 99   continue
      lpyr=1980
      nrec=0
      irec=0
      do 999 iyr=1,maxyr
      iyear=1979+iyr
      nd=0
      do 997 imon=1,12
      max=numd(imon)
      if(iyear.eq.lpyr.and.imon.eq.2)max=max+1
      do 998 iday=1,max
      irec=irec+1
      read(11,rec=irec)x
      if(imon.eq.2.and.iday.eq.29)go to 998
      nd=nd+1
      do 1 i=1,imax
      do 2 j=1,jmax
      if(x(i,j).lt.0.)go to 2
      if(x(i,j).gt.310.)go to 2
      cnt(i,j,nd)=cnt(i,j,nd)+1.
      y(i,j,nd)=y(i,j,nd)+x(i,j)
 2    continue
 1    continue
 998  continue
 997  continue
      if(iyear.eq.lpyr)lpyr=lpyr+4
 999  continue
ccc   compute 30-yr climatology for 365-d year
      do 3 i=1,imax
      do 4 j=1,jmax
      do 5 nd=1,365
      y(i,j,nd)=y(i,j,nd)/cnt(i,j,nd)
 5    continue
 4    continue
 3    continue
ccc   create a 3yr file then smooth data with a 31-d running mean filter
      nday=0
      do 6 iyr=1,3
      do 7 nd=1,365
      nday=nday+1
      do 8 i=1,imax
      do 9 j=1,jmax
      yy(i,j,nday)=y(i,j,nd)
 9    continue
 8    continue
 7    continue
 6    continue
ccc   now smooth the data
      nd=0
      ndd=0
      do 12 nday=366,730
      nd=nd+1
      ndd=ndd+1
      min=nday-15
      max=nday+15
      do 10 i=1,imax
      do 11 j=1,jmax
      sum=0.
      do 13 k=min,max
      sum=sum+yy(i,j,k)
 13   continue
      ys365(i,j)=sum/31.
 11   continue
 10   continue 
      write(6,200)nd,yy(20,20,nd),ys365(20,20)
 200  format(i5,2x,3f8.1)
      write(51,rec=nd)ys365
      write(52,rec=ndd)ys365
      if(nd.eq.59)go to 50
      go to 12
 50   ndd=ndd+1
      write(52,rec=ndd)ys365
 12   continue
 990  continue
      stop
      end