program vaverage
       real d(20000,18),f(5000,18)
       integer dmax,nums,numso,i,j,k,ivm,vint(20),idrx,icat(8),
     &   hal(50)
       real r34,r50,r64,vavx,vmax,r1,r2,dr,dr1,rr,vm,rm,v0,r0,dm,
     &     rm1,alpha,alpha0,ravx1,pi,area,filter,ike,ike2,ikear,fcor
       irn=760013
       ip=0
       itc0=0
       itc1=0
       scal=0. !nonzero for random perturbations (suggest scal=1)
 8     i=1
       k=1
       numso=0
       vmax=0.
       dmax=1
       pi=4.*atan(1.)
       dx=60.
c       filter=4.  !(filter scale in grid points)
       filter=4.  !(filter scale in grid points)
       idrx=101
       dinc=1./8.
c       area=filter*filter*dx*dx/16.
       area=filter*filter*dx*dx/4.  !area is homogeneous (filter/2)**2
       alpha=0.5
       do ivm=1,8
        icat(ivm)=0
       enddo
       do i=1,50
        hal(i)=0
       enddo
 12    read(9,*,end=20)(d(i,j),j=1,18)
       nums=int(d(i,1)/100.)
       dm=0.
       do m=5,8
        dm=max(d(i,m),0.)
       enddo
c       print*,dm,nums,numso
       if (dm.le.0.) goto 12
       if ((nums.ne.numso).and.(numso.gt.0)) then !new storm
        k=k+1  !storm number 
c        print*,k-1,nums,vmax,dmax,d(dmax,2)
        do j=1,18
         f(k-1,j)=d(dmax,j)
c         print*,dmax,d(dmax,j)
         d(1,j)=d(i,j)
        enddo
        write(11,*)float(i),vmax
        vmax=0.
        dmax=1
        i=1
       endif
       if (d(i,2).gt.vmax) then  !max intensity so far
        vmax=d(i,2)
        dmax=i
       endif
c       print*,k,nums,vmax,dmax,d(i,2)
       numso=nums
       i=i+1
       goto 12
c
 20    do j=1,18
        f(k,j)=d(dmax,j)
       enddo 
       kmax=k
c       print*,kmax
       do ivm=1,15
        vint(ivm)=0
       enddo
       do i=1,kmax
        r34=0.
        r50=0.
        r64=0.
        im34=0
        im50=0
        im64=0
        do m=5,8
         if (f(i,m).gt.r34) then
          r34=f(i,m)
          r50=f(i,m+4)
          r64=f(i,m+8)
         endif
        enddo
        write(26,93)i,f(i,2),f(i,4),r64,r50,r34
 93     format(i5,5f7.0)
c fill in RMW
        factor=(1 + scal*(ran(irn)-0.5))
        write(60,*)ip,factor
        if (r64.gt.0.) then
         v0=64.
         r0=r64*factor
        else if (r50.gt.0.) then
         v0=50.
         r0=r50*factor
        else if (r34.gt.0.) then
         v0=34.
         r0=r34*factor
        endif
C Compute rm
        vm=f(i,2)
        fcor=1.458e-4*sin(3.14159*f(i,17)/180.)
        vma=vm*0.5144
        r0a=r0*1852.
        v0a=v0*0.5144
        rm1=f(i,4)
        f(i,4)=( (vma - (vma**2 - (v0a**2 - 0.25*(r0a*fcor)**2))**0.5)/
     &     (v0a/r0a - 0.5*fcor) )/1852.
        rm=f(i,4)
        write(25,*)f(i,1),vm,rm,rm1,r0,v0
        write(6,*)vm,rm,rm1
        vavx=0.
        avvav=0.
        r1=0.
        ike=0.
        ike2=0.
        ikear=0.
        do m=1,idrx*int(1./dinc)
         r2=sqrt(r1**2 + 4.*area/pi)
         dr=r2-r1
         if (m.eq.1) write(52,*)dr,r1,r2
         r00=0.5*(r1+r2)
         vav=0.
         rav=0.
C Compute velocity profiles and filter it
         do n=1,idrx
          dr1=(r2-r1)/(float(idrx-1)) 
          rr=r1 + (n-0.5)*dr1
          if (rr.le.r0) then
           vma=vm*0.5144
           rma=rm*1852.
           ra=rr*1852.
           v=(2.*(rma*vma+0.5*fcor*rma*rma)/(ra + rma**2/ra)
     &        - 0.5*fcor*ra)/0.5144
          else if ((rr.lt.r50).and.(r64.gt.0.)) then
           alpha=alog(50./64.)/alog(r64/r50)
           v=64.*((r64/rr)**(alpha))
          else if ((rr.lt.r34).and.(r50.gt.0.)) then
           alpha=alog(34./50.)/alog(r50/r34)
           v=50.*((r50/rr)**(alpha))
          else if (rr.gt.r34) then
           alpha=0.5
           v=34.*((r34/rr)**(alpha))
          endif
          if ((m.eq.1).and.(vm.ge.140.)) 
     &     write(46,*)f(i,1),v,rr
          vav=vav + v*rr*dr1
          rav=rav + rr*dr1
          if ((r00.lt.300.).and.(n.eq.(idrx-1)/2)) then
           v0=v
          endif
         enddo
         if (vav/rav.gt.vavx) then
          vavx=vav/rav
          ravx1=r00
          drx1=dr
         endif
         avvav=vavx  !avvav + vav/rav/float(idrx)
         r1=r1 + dr*dinc
        enddo
c        if (vavx.ge.34) itc1=itc1+1
c        if (vavx.ge.26) itc0=itc0+1
        write(27,91)i,f(i,2),f(i,3),vavx,rm,ravx1,drx1
 91     format(i5,8f8.1)
        ivm=int(avvav/10.) 
        if (ivm.gt.15) ivm=15
        vint(ivm)=vint(ivm) + 1
c Increment counts in different SS-scale bins
        if (avvav.lt.34) then
         icat(1)=icat(1) + 1
        else if (avvav.lt.50) then
         icat(2)=icat(2) + 1
        else if (avvav.lt.64) then
         icat(3)=icat(3) + 1
        else if (avvav.lt.83) then
         icat(4)=icat(4) + 1
        else if (avvav.lt.96) then
         icat(5)=icat(5) + 1
        else if (avvav.lt.113) then
         icat(6)=icat(6) + 1
        else if (avvav.lt.137) then
         icat(7)=icat(7) + 1
        else 
         icat(8)=icat(8) + 1
        endif
       enddo
       write(58,94)(vint(ivm),ivm=2,15)
 94    format(14i4)
       write(57,99)(icat(ivm),ivm=1,8)
 99    format(8i4)
       ip=ip+1
       rewind(9)
       if (ip.lt.1) goto 8
c       if (ip.lt.100) goto 8  !(set to number of random samples)
       print*,itc0,itc1
       stop
       end