program readdata
       parameter (ntot=71399-1, ny=2021-1978, nday=31+30+31+30+31)
c  nday: March - July
       parameter (ni=81, nj=41)
c ni,nj: number of grid points for US domain

       real td(ni,nj,nday,ny), td7(ni,nj,nday,ny)

       real x(81), y(41)

       character*8 cyr, cmon, cday, cmag, cslon, cslat

       open(11, file='d.year')
       open(12, file='d.month')
       open(13, file='d.day')
       open(14, file='d.mag')

       open(15, file='d.slon')
       open(16, file='d.slat')

cc       open(21,file='td.MAMJJ.1979-2021.txt')

       open(22,file='td.MAMJJ.1979-2021.map.gr',
     1  access='direct',form='unformatted',recl=ni*nj*nday*ny*4)

cc       open(31,file='td.MAMJj.1979-2021.7d-mean.txt')

       open(32,file='td.MAMJj.1979-2021.7d-mean.map.gr',
     1  access='direct',form='unformatted',recl=ni*nj*nday*ny*4)

c 1x1 grid
       do i=1,81
        x(i)=(i-1)*1.0  + 230.0
       end do

       do j=1,41
        y(j)=(j-1)*1.0 + 20.0
       end do

       print *, 'Total Years   ', ny 
       print *, 'Total Days    ', nday 

       read(11,*) cyr
       read(12,*) cmon
       read(13,*) cday
       read(14,*) cmag
       read(15,*) cslon
       read(16,*) cslat

       print *, 'LSR:    ', cyr, '     ',cmon,'    ',cday,'    ',cmag,
     1          '    ',cslon,'    ',cslat

       do k=1,ny
       do n=1,nday
       do i=1,ni
       do j=1,nj
        td(i,j,n,k)=0.0
       end do
       end do
       end do
       end do

       do ii=1,ntot
       read(11,*) iyear 
       read(12,*) imon
       read(13,*) iday
       read(14,*) imag
       read(15,*) slon
       read(16,*) slat

       slon=360.0 + slon

       if(ii.le.10) then
         print *, iyear, imon, iday, imag, slon, slat
       end if

c Year:  1979-2021
c Month: March - July
       if(iyear.ge.1979.and.iyear.le.2021) then
       if(imon.ge.3.and.imon.le.7) then

         iy=iyear - 1978

         if(imon.eq.3) nd=iday
         if(imon.eq.4) nd=iday+31
         if(imon.eq.5) nd=iday+61
         if(imon.eq.6) nd=iday+92
         if(imon.eq.7) nd=iday+122 

c March - July: 31+30+31+30+31 = 153 days

         if(imag.ge.1) then

         x0=slon
         y0=slat

         do i=1,ni
         do j=1,nj
c Td in 5x5 box
          a=x(i)-2.5
          b=x(i)+2.5
          c=y(j)-2.5
          d=y(j)+2.5

          if(x0.ge.a.and.x0.lt.b) then
          if(y0.ge.c.and.y0.lt.d) then
           td(i,j,nd,iy)=td(i,j,nd,iy) + 1.0
          end if
          end if
         end do
         end do

         end if

       end if
       end if

       end do

c 7-day mean: March - July.21

       do iy=1,ny
       do nd=1,nday-10

       do i=1,ni
       do j=1,nj
       td7(i,j,nd,iy)=(td(i,j,nd,iy)+td(i,j,nd+1,iy)+td(i,j,nd+2,iy)
     1                +td(i,j,nd+3,iy)+td(i,j,nd+4,iy)
     &                +td(i,j,nd+5,iy)+td(i,j,nd+6,iy))/7.0

       end do
       end do

       end do

C For July.22-31

       do nd=nday-10+1, nday

       do i=1,ni
       do j=1,nj
       td7(i,j,nd,iy)=-9.99
       end do
       end do

       end do

       end do

c       do iy=1,ny
c       do nd=1,nday
c       write(21,200) iy+1978, nd, td(nd,iy)
c       write(31,200) iy+1978, nd, td7(nd,iy) 
c200      format(I4,4x,I3,4x,F8.2)
c       end do
c       end do

       write(22,rec=1) 
     1  ((((td(i,j,nd,iy),i=1,ni),j=1,nj),nd=1,nday),iy=1,ny)
       write(32,rec=1) 
     1  ((((td7(i,j,nd,iy),i=1,ni),j=1,nj),nd=1,nday),iy=1,ny)

       stop
       end