program readata
       parameter (ny=64)
       parameter (ix=81, il=41)
c ny: 1955-2018

c MAM: 31 + 30 + 31 = 92 days
       real dt(ix,il,ny,92)
       real dh(ix,il,ny,92)

       real t(ix,il,ny)
       real c(ix,il),sdv(ix,il)

       real s(ix,il)

c 1955-2018: 64 years
c MAM: 92 days

c PPH: Tornado
       open(11,file='td-PPH_1955-2018_MAM.92d.1x1.gr',
     1 form='unformatted',access='direct',recl=ix*il*64*92*4)
c PPH: Hail
       open(12,file='hl-PPH_1955-2018_MAM.92d.1x1.gr',
     1 form='unformatted',access='direct',recl=ix*il*64*92*4)
c Land-sea mask
       open(15,file='us.land.1x1.gr',
     1 form='unformatted',access='direct',recl=ix*il*4)

c Output: 1981-2010 climatology and standrad deviation of event days
c For MAM
       open(21,file='MAM.clm.sdv.1981-2010.eventday.gr',
     1 form='unformatted',access='direct',recl=ix*il*4)


        read(11,rec=1) 
     1  ((((dt(i,j,iy,k),i=1,ix),j=1,il),k=1,92),iy=1,ny)
        read(12,rec=1) 
     1  ((((dh(i,j,iy,k),i=1,ix),j=1,il),k=1,92),iy=1,ny)

        read(15,rec=1) ((s(i,j),i=1,ix),j=1,il)

c PPH*100: % value
        do i=1,ix
        do j=1,il
        do k=1,92
        do iy=1,ny
         dt(i,j,iy,k)=dt(i,j,iy,k)*100.0
         dh(i,j,iy,k)=dh(i,j,iy,k)*100.0
        end do
        end do
        end do
        end do

        do i=1,ix
        do j=1,il
        do iy=1,ny
         t(i,j,iy)=0.0
         do k=1,92
         if(dt(i,j,iy,k).ge.10.or.dh(i,j,iy,k).ge.30) then
         t(i,j,iy)=t(i,j,iy) + 1.0
         end if
         end do
        end do
        end do
        end do

        do i=1,ix
        do j=1,il
        do iy=1,ny
        if(s(i,j).lt.0.3) then
         t(i,j,iy)=-8888.8
        end if
        end do
        end do
        end do

c 1981-2010 climatology
        do i=1,ix
        do j=1,il
         c(i,j)=0.0

c Over land (US)
         if(s(i,j).ge.0.3) then
c 1981-2010
         do iy=27,56
         c(i,j)=c(i,j)+t(i,j,iy)
         end do
          c(i,j)=c(i,j)/30.0
         else
          c(i,j)=-8888.8
         end if
        end do
        end do
          
c sdv
        do i=1,ix
        do j=1,il
         sdv(i,j)=0.0
         if(s(i,j).ge.0.3) then
         do iy=27,56
         sdv(i,j)=sdv(i,j)+(t(i,j,iy)-c(i,j))**2
         end do
          sdv(i,j)=sqrt(sdv(i,j)/30.0)
         else
          sdv(i,j)=-8888.8
         end if
        end do
        end do

        write(21,rec=1) ((c(i,j),i=1,ix),j=1,il)
        write(21,rec=2) ((sdv(i,j),i=1,ix),j=1,il)

c a single number for US
        ct=0.0
        nn=0.0
        a=0.0
        do i=1,ix
        do j=1,il
         if(s(i,j).ge.0.3) then
          ct=ct + c(i,j)
          nn=nn + 1
          if(c(i,j).gt.a) a=c(i,j)
         end if
        end do
        end do

        ct=ct/real(nn)

        print *, ct, nn, a

        
        stop
        end