program read_parsivel

c          this program reads parsivel with measured fall velocities

           integer i,j,k,nx,nz,ny
           parameter(nx=940000,nz=32,ny=60)
           parameter(nd=1024,isay=1)

           real vel(nx,nz),dsd(nx,nz)
           real binb(nx,nz,nz),dold(nz),vold(nz)
           real rtot(nx),rrate,rtot2,delt(nz),d_bound(nz)
           real dmm(nz),delta(nz),conave(nz)
           real dtv(ny),tv(ny),velt(nz),time(nx)
           real bina(nx,nz),velo(nz),dvel(nz)
           real dsdcal(nx,nz,nz),dsdc(nx,nz)
           real dsdt(nx,nz)
           real condave(nz),vel1(nz),vel2(nz)
           real vela(nx,nz),velm(nz),binm(nz)

           integer binc(nz,nz),is(nz,nz),iss(nz,nz)
           integer bon(nx,nd),jd(nx),ka(nx)
           integer jday(nx),bin(nx,nd),ma(nx)
           integer kr(nx),mn(nx),isec(nx)
           integer jsec(7),ibin(nx,nz)

           character*40 infile,infile2,out1,out2,out3,out4
           Character*49 out5,out6
           
           open(10,file='fileapu',status='unknown')
           open(12,file='parsivel_diameter.txt',status='unknown')
           open(13,file='parsivel_matrix.txt',status='unknown')


          eps=0.01

          do k=1,nz
          read(12,*)dmm(k),delta(k),velt(k),velo(k)
          enddo

          pi=acos(-1.)

          dcount=0.
          daccept=0.
     
          do kk=1,nz
          read(13,61)(iss(kk,k),k=1,nz)
 61       format(32i2)
          enddo

          print*,'fl'

          do il=1,isay
          read(10,21,err=988,end=988)infile,nsay,out1,out2,out3,
     *    out4,out5,out6
 21       format(a9,1x,i3,1x,a22,1x,a32,1x,a35,1x,a39,1x,a29,1x,a33)

           open(11,file=infile,status='unknown')

           open(20,file=out1,status='unknown')
           open(21,file=out2,status='unknown')
           open(22,file=out3,status='unknown')
           open(23,file=out4,status='unknown')
           open(24,file=out5,status='unknown')
           open(25,file=out6,status='unknown')

          do ii=1,nsay
          read(11,19)infile2
 19       format(a24)
          
          print*,ii,nsay

          open(14,file=infile2,status='unknown')

          k0=1
          kk=1
          do i=1,nx
          read(14,30,err=999,end=999)iyear,mon,iday,kr(k0),mn(k0),
     *    isec(k0),(bin(k0,k),k=1,nd)
c         print*,i
30       format(i4,i2,i2,i2,i2,i2,45x,1024(1x,i3))
          if(mon.eq.1)jday(k0)=iday
          if(mon.eq.2)jday(k0)=iday+31
          if(mon.eq.3)jday(k0)=iday+59
          if(mon.eq.4)jday(k0)=iday+90
          if(mon.eq.5)jday(k0)=iday+120
          if(mon.eq.6)jday(k0)=iday+151
          if(mon.eq.7)jday(k0)=iday+181
          if(mon.eq.8)jday(k0)=iday+212
          if(mon.eq.9)jday(k0)=iday+243
          if(mon.eq.10)jday(k0)=iday+273
          if(mon.eq.11)jday(k0)=iday+304
          if(mon.eq.12)jday(k0)=iday+334

c         jday(k0)=jday(k0)-1

          if(k0.eq.1)jsec(kk)=isec(k0)
          if(k0.ge.2)then
          if(mn(k0).eq.mn(k0-1))then
          kk=kk+1
          jsec(kk)=isec(k0)
          else
          nk=kk
          if(kk.gt.6)nk=6
          kky=2022
          write(20,40)kky,jday(k0-1),kr(k0-1),mn(k0-1),
     *    (jsec(kk),kk=1,nk)
 40       format(4i5,7i5)
          kk=1
          jsec(kk)=isec(k0)
          endif
          endif
          k0=k0+1
999       enddo

          nn=k0-1
          print*,k0

          k1=0
          do i=1,nn
          do k=1,nd
          bon(i,k)=0
          enddo
          if(isec(i).eq.0 .and. isec(i+5).eq.50)then
          k1=k1+1
          jd(k1)=jday(i)
          ka(k1)=kr(i)
          ma(k1)=mn(i)
          do k=1,nd
          bon(k1,k)=bin(i,k)+bin(i+1,k)+bin(i+2,k)+
     *    bin(i+3,k)+bin(i+4,k)+bin(i+5,k)
          enddo
          endif
          enddo


          do i=1,k1
          do k=1,nz
          bina(i,k)=0.
          vela(i,k)=0.
          ibin(i,k)=0
          kb=0
          velave=0.
          do kk=k,nd,nz
          kb=kb+1
          binb(i,kb,k)=0.
          if(iss(kb,k).eq.1)then
          bina(i,k)=bina(i,k)+float(bon(i,kk))
          ibin(i,k)=ibin(i,k)+float(bon(i,kk))
          binb(i,kb,k)=float(bon(i,kk))
         if(bin(i,kk).gt.0)vela(i,k)=vela(i,k)+velo(kb)*float(bon(i,kk))
          endif
          if(iss(kb,k).eq.1)daccept=daccept+float(bon(i,kk))
          dcount=dcount+float(bon(i,kk))
          enddo
          velave=velave/bina(i,k)
          enddo


c         calculation of dsd

          do k=1,nz
          dsdc(i,k)=0.
          dsdt(i,k)=0.
          if(dsd(i,k).ne.0.)then
          dsd(i,k)=10.**(dsd(i,k))
          endif
          enddo

          do k=1,nz
          cs=(180.*(30.-(dmm(k)/2.)))/100.
          do kb=1,nz
          bot=60.*cs*velo(kb)*delta(k)*100. 
          bott=60.*cs*velt(k)*delta(k)*100. 
          dsdc(i,k)=dsdc(i,k)+binb(i,kb,k)*1.e6/bot
          dsdt(i,k)=dsdt(i,k)+binb(i,kb,k)*1.e6/bott
          enddo
          enddo

c         integral rainfall parameters

          rain=0.
          raint=0.
          raind=0.
          top=0.
          wc=0.
          wct=0.
          dbzt=0.
          dbz=0.
          con=0.
          cont=0.
          x4=0.
          x4t=0.
          x3=0.
          x3t=0.
          sig1=0
          sig2=0
          dmax=0.


          do k=1,nz
          cs=(180.*(30.-(dmm(k)/2.)))
          bott=60.*cs*velt(k)
          vol=pi*(dmm(k)**3.)/6.
          top=top+bina(i,k)
          if(bina(i,k).gt.0.)dmax=dmm(k)
          do kb=1,nz
          bot=60.*cs*velo(kb)
          timerain=60.
          con=con+binb(i,kb,k)*1.e6/bot 
          wc=wc+vol*binb(i,kb,k)*1.e3/bot
          dbz=dbz+binb(i,kb,k)*(dmm(k)**6.)*1.e6/bot
          x4=x4+binb(i,kb,k)*(dmm(k)**4.)*1.e6/bot
          x3=x3+binb(i,kb,k)*(dmm(k)**3.)*1.e6/bot
          cont=cont+binb(i,kb,k)*1.e6/bott 
          wct=wct+vol*binb(i,kb,k)*1.e3/bott
          dbzt=dbzt+binb(i,kb,k)*(dmm(k)**6.)*1.e6/bott
          x4t=x4t+binb(i,kb,k)*(dmm(k)**4.)*1.e6/bott
          x3t=x3t+binb(i,kb,k)*(dmm(k)**3.)*1.e6/bott
          raind=raind+binb(i,kb,k)*vol*timerain/cs
          enddo
          enddo
          dbz=10.*alog10(dbz)
          dbzt=10.*alog10(dbzt)
          dmass=x4/x3
          dmasst=x4t/x3t

          do k=1,new
          cs=(180.*(30.-(dmm(k)/2.)))
          vol=pi*(dmm(k)**3.)/6.
          do kb=1,nz
          bot=60.*cs*velo(kb)
          sig=(dmm(k)-dmass)**2.
          sig1=sig1+sig*(dmm(k)**3.)*binb(i,kb,k)*1.e6/bot 
          sig2=sig2+(dmm(k)**3.)*binb(i,kb,k)*1.e6/bot 
          enddo
          enddo
          sig0=(sig1/((dmass**2.)*sig2))**0.5
          sig0=sig0*dmass

          if(raind.ge.0.01 .and. top.ge.10.)then
          write(22,33)iyear,jd(i),ka(i),ma(i),top,con,
     *    wc,raind,dbz,dmass,dmax
          write(23,33)iyear,jd(i),ka(i),ma(i),top,cont,
     *    wct,raind,dbzt,dmasst,dmax
 33       format(4i5,f10.0,f10.3,5f10.3)
          write(24,10)iyear,jd(i),ka(i),ma(i),(dsdc(i,k),k=1,nz)
          write(25,10)iyear,jd(i),ka(i),ma(i),(dsdt(i,k),k=1,nz)
          write(21,11)iyear,jd(i),ka(i),ma(i),(ibin(i,k),k=1,nz)
 11       format(4i5,32i4)
 10       format(4i5,32f10.4)
          endif

          enddo

          enddo

988       enddo

          stop
	  end