
            Program zswehist

            Parameter(nx=30000,nm=8,nn=4,na=4)
            Parameter(ns=60,nk=300,nz=100,nd=6)

            Integer imatrix(ns,ns)
            Integer jmatrix(ns,ns)
            Integer lmatrix(ns,ns)
            Integer itots(nx),is0(na)
            Integer msay(24),nsay(22)

            Real con(nx),sdepth(nx),iwc(nx,nm,nn)
            Real time(nx),sdb(na,nx)
            Real zdb(na,nx),swe(na,nx),zmm(na,nx)
            Real zdbave(ns,ns,nn,nn),pzdb(ns,ns)
            Real pmatrix(ns,ns),plmatrix(ns,ns)
            Real pjmatrix(ns,ns)
            Real zdbsv(nx),sdbsv(nx)
            Real swenw(nx),swenwv(nx)
            Real swenwsort(nx)
            Real psrat(ns,ns),sratave(ns,ns)
            Real*8 dm1,dm2,xnw1,xnw2


            Character*80 infile1,infile2,infile3,infile4
            Character*80 infile5,infile0
            Character*80 infile01,infile02,infile03

            Character*80 outfile,outfile2,outfile3
            Character*80,outfile4,outfile5,outfile6
            Character*80 outfile7,outfile8,outfile9
            Character*80 outfile0,outevent,outfile1

            open(10,file='filezswehist_lowtrack.impacts2022',
     *      status='unknown')

            read(10,22)infile0
 22         format(a38)
            read(10,30)outfile0
            read(10,30)outfile1
            read(10,30)outfile2
            read(10,30)outfile3
 30         format(a32)

            
            open(11,file=infile0,status='unknown')
            open(20,file=outfile0,status='unknown')
            open(21,file=outfile1,status='unknown')
            open(22,file=outfile2,status='unknown')
            open(23,file=outfile3,status='unknown')

            
            k1=0
            k2=0
            k3=0
            k4=0
            do k=1,nx
            read(11,*,err=999,end=999)iyr,jd,ka,ma,
     *      ifinal,swe0,zmm0,zdb0,srat0
            if(ifinal.eq.1)then
            k1=k1+1
            is0(1)=is0(1)+1
            zdb(1,k1)=zdb0
            zmm(1,k1)=zmm0
            swe(1,k1)=swe0
            sdb(1,k1)=10.*alog10(swe0)
            endif
            if(ifinal.eq.2)then
            k2=k2+1
            is0(2)=is0(2)+1
            zdb(2,k2)=zdb0
            zmm(2,k2)=zmm0
            swe(2,k2)=swe0
            sdb(2,k2)=10.*alog10(swe0)
            endif
            if(ifinal.eq.3)then
            k3=k3+1
            is0(3)=is0(3)+1
            zdb(3,k3)=zdb0
            zmm(3,k3)=zmm0
            swe(3,k3)=swe0
            sdb(3,k3)=10.*alog10(swe0)
            endif
            if(ifinal.eq.4)then
            k4=k4+1
            is0(4)=is0(4)+1
            zdb(4,k4)=zdb0
            zmm(4,k4)=zmm0
            swe(4,k4)=swe0
            sdb(4,k4)=10.*alog10(swe0)
            endif
 999        enddo

            print*,k1,k2,k3,k4
            

            do ii=1,na

            is=is0(ii)

            print*,is

            zmin=40.
            zmax=-20.

            do ll=1,is
            zdbsv(ll)=zdb(ii,ll)
            enddo
            is02=float(is)*0.02
            is98=float(is)*0.98
            call sort(zdbsv,nx,is)
            if(zdbsv(is02).ne.0. .and. zdbsv(is02).lt.zmin)
     *      zmin=zdbsv(is02)
            if(zdbsv(is98).gt.zmax)zmax=zdbsv(is98)

            deltaz=1.
            deltas=0.5

            z1=-20.
            z2=z1+deltaz
            do kw=1,ns
            s1=-20.
            s2=s1+deltas
            do kq=1,ns
            imatrix(kw,kq)=0
            sratave(kw,kq)=0.
            do i=1,is
            if(zdb(ii,i).ge.z1 .and. zdb(ii,i).lt.z2)then
            if(sdb(ii,i).ge.s1 .and. sdb(ii,i).lt.s2)then
            imatrix(kw,kq)=imatrix(kw,kq)+1
c           sratave(kw,kq)=sratave(kw,kq)+srat(ii,i)
            endif
            endif
            enddo

            s1=s1+deltas
            s2=s2+deltas
            enddo
            z1=z1+deltaz
            z2=z2+deltaz
            enddo


            isum=0
            do kw=1,ns
            do kq=1,ns
            pmatrix(kw,kq)=0.
            isum=isum+imatrix(kw,kq)
            psrat(kw,kq)=0.
            enddo
            enddo
            do kw=1,ns
            do kq=1,ns
            pmatrix(kw,kq)=float(imatrix(kw,kq))/float(isum)
            if(imatrix(kw,kq).gt.0.)then
c           psrat(kw,kq)=sratave(kw,kq)/float(imatrix(kw,kq))
            endif
            enddo
            write(ii+19,46)(pmatrix(kw,kq),kq=1,ns)
c           write(21,46)(psrat(kw,kq),kq=1,ns)
 46         format(60e15.8)
            enddo
 

            enddo

            stop
            end
            
    

            
