
         Program particle

         Parameter(nx=4150000,ny=100,nz=5,nn=8)
         Parameter(nu=1,nm=8,nc=90000)

         Integer ksave(nx),msave(nx)
         Integer ip05(ny,nm)
         Integer jwx(nx),ivis,kwx,mwx
         Real am1(nm),bm1(nm),av(nm),bv(nm),ac(nm),bc(nm)
         Real ar(nm),br(nm),ax1(nm),bx1(nm)
         Real am2(nm),bm2(nm),ax2(nm),bx2(nm)
         Real am3(nm),bm3(nm),ax3(nm),bx3(nm)
         Real am4(nm),bm4(nm),ax4(nm),bx4(nm)
         Real dmax(nx,nm),fac(nm),arat(nm),dmx(ny,nm)

         Real bscat(ny),bscatbm05(ny,nm)
         Real bscathw05(ny,nm)
         Real bscatmh05(ny,nm)
         Real bscatkc05(ny,nm)

         Real dmm(nx),velave(ny),cor(ny),f(nz),fth(nz)
         Real vsave(nx),vel05(ny,nm)
         Real hsave(nx),hvel05(ny,nm)

         Real vmed05(ny,nm),vsv05(ny,nm,nc),vsort05(nc)

         Real dln(ny),x1ln(ny),ddln(ny),vln(ny)
         Real drln05(ny),rln05(ny),db1ln(ny),b1ln(ny)
         Real x2ln(ny),b2ln(ny)
         Real xmh(nn),xkc(nn)
         Real x1ln05(ny),x2ln05(ny),xln(ny)
         Real x3ln05(ny),x4ln05(ny)
         Real dbln05(ny),b1ln05(ny),b2ln05(ny)
         Real b3ln05(ny),b4ln05(ny)
         Real dln05(ny),xln05(ny)
         Real ddln05(ny)
         Real vln05(ny)
         Real cln05(ny),dcln05(ny)
         Real twx(nx),vis(nx),airden(nx)
         Real rey05(ny,nm),xbm05(ny,nm),xmbm05(ny,nm)
         Real bkc05(ny,nm),bmh05(ny,nm)
         Real xmmh05(ny,nm)
         Real xmkc05(ny,nm)
         Real xmhw05(ny,nm)
         Real xhw05(ny,nm)

         character*60 infile,infile2,infile3
         Character*60 outfile,outfile2,outfile3
         Character*60 outf01,outf02,outf03,outf04
         Character*60 outf11,outf12,outf13,outf14


         data fth/0.41,0.18,0.1,0.21,0.05/

         data fac/1.0,1.1,1.15,1.2,1.25,1.3,1.35,1.4/

         data xmh/3.8233,-1.5211,0.30065,-0.06104,
     *    0.13074,-0.073429,0.016006,-0.0012483/

         data xkc/3.8816,-1.4579,0.27749,-0.41521,0.57683,
     *   -0.29220,0.06467,-0.0053405/ 

         Open(10,file='filevelmedsnow0',status='unknown')

         open(14,file='AIO_weather.impacts2023',status='unknown')

         open(15,file='filemedmass',status='unknown')

         open(16,file='pip_backscat_ice_no_vol_rgans.txt',
     *   status='unknown')

         open(17,file='filemedbscat',status='unknown')

         open(24,file='PIP003_medmass_012_plot.impact2022',
     *   status='unknown')

         read(14,*)
         k0=1
         do i=1,nx
         read(14,*,err=988,end=988)iy,jwx(k0),kwx,mwx,wd,wsp,tc,wb,
     *   rh,p1,airden(k0),vis(k0)
         twx(k0)=float(kwx)+float(mwx)/60.
         k0=k0+1
 988     enddo
         k0=k0-1

         print*,k0

         do k=1,ny
         read(16,*)dmid60,bscat(k)
         do kk=1,nm
         dmx(k,kk)=dmid60*fac(kk)
         enddo
         enddo
         
         do k=1,nm
         arat(k)=(1./fac(k))**2.
         enddo

         delta0=5.83
         c0=0.6
         delta00=8.0
         c00=0.35
         gg=9.81
         pi=acos(-1.)
         cons00=(delta0**2.)*(c0**0.5)/4.
         cons01=(delta00**2.)*(c00**0.5)/4.

         nw=8

         do ll=1,nu

         read(10,30)infile,outfile,outfile2,outfile3
 30      format(a44,1x,a48,1x,a49,1x,a44)

         read(15,37)outf01,outf02,outf03,outf04
 37      format(a52,1x,a52,1x,a52,1x,a52)

         read(17,36)outf11,outf12,outf13,outf14
 36      format(a53,1x,a53,1x,a53,1x,a53)

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

         open(20,file=outfile,status='unknown')
         open(21,file=outfile2,status='unknown')
         open(22,file=outfile3,status='unknown')

         open(40,file=outf01,status='unknown')
         open(41,file=outf02,status='unknown')
         open(42,file=outf03,status='unknown')
         open(43,file=outf04,status='unknown')

         open(50,file=outf11,status='unknown')
         open(51,file=outf12,status='unknown')
         open(52,file=outf13,status='unknown')
         open(53,file=outf14,status='unknown')

        
         k1=0
         k2=0
         do k=1,nx
         read(13,*,err=999,end=999)ipa,iyear,jd,ka,ma,
     *    dpip,vvel,hvel,ksay,kksay
c        if(mon.eq.10)jday=iday+273
c        if(mon.eq.11)jday=iday+304
c        if(mon.eq.12)jday=iday+334
c        if(mon.eq.1)jday=iday
c        if(mon.eq.2)jday=iday+31
c        if(mon.eq.3)jday=iday+59
c        if(mon.eq.4)jday=iday+90
         k2=k2+1
         if(vvel.ge.0.0)then
         k1=k1+1
         ksave(k1)=ka
         msave(k1)=ma
         dmm(k1)=dpip
         do kl=1,nm
         dmax(k1,kl)=fac(kl)*dpip
         enddo
         vsave(k1)=vvel
         hsave(k1)=abs(hvel)
         endif
 999     enddo

         print*,k1,k2

         if(k1.eq.0)go to 333

c        5-min averaging


         kh=0
         do kk=1,24
         mn=0
         mn1=mn+5
         do il=1,12
         isay=0
         do km=1,ny
         do kl=1,nm
         xbm05(km,kl)=0.
         xhw05(km,kl)=0.
         bmh05(km,kl)=0.
         bkc05(km,kl)=0.
         rey05(km,kl)=0.
         vel05(km,kl)=0.
         vmed05(km,kl)=0.
         do kc=1,nc
         vsv05(km,kl,kc)=0.
         enddo
         hvel05(km,kl)=0.
         ip05(km,kl)=0
         xmbm05(km,kl)=0.
         xmmh05(km,kl)=0.
         xmkc05(km,kl)=0.
         xmhw05(km,kl)=0.
         bscatbm05(km,kl)=0.
         bscatmh05(km,kl)=0.
         bscatkc05(km,kl)=0.
         bscathw05(km,kl)=0.
         enddo
         enddo
         fmin=500.
         do ml=1,nz
         f(ml)=0.
         enddo
         ihabit=0
         do j=1,k1
         if(ksave(j).eq.kh .and. 
     *   msave(j).ge.mn .and. msave(j).lt.mn1)then
         isay=isay+1
         f(1)=f(1)+abs(9.65-10.3*exp(-0.6*dmm(j))-vsave(j))
         f(2)=f(2)+abs(1.3*(dmm(j)**0.94)-vsave(j))
         f(3)=f(3)+abs(0.82*(dmm(j)**0.24)-vsave(j))
         f(4)=f(4)+abs(1.03*(dmm(j)**0.71)-vsave(j))
         f(5)=f(5)+abs(0.74*(dmm(j)**0.35)-vsave(j))
         do kl=1,nm
         a1=0.0
         a2=0.2
         do km=1,ny
         if(dmax(j,kl).ge.a1 .and.dmax(j,kl).lt.a2)then
         ip05(km,kl)=ip05(km,kl)+1
         vel05(km,kl)=vel05(km,kl)+vsave(j)
         vsv05(km,kl,ip05(km,kl))=vsave(j)
         hvel05(km,kl)=hvel05(km,kl)+hsave(j)
         endif
         a1=a1+0.2
         a2=a2+0.2
         enddo
         enddo
         endif
         enddo

         if(isay.gt.0)then
         tsd=float(kh)+float(mn)/60.
         tsd1=float(kh)+float(mn1)/60.
         do mt=1,k0
         if(jd.eq.jwx(mt) .and. twx(mt).ge.tsd
     *   .and. twx(mt).lt.tsd1 .and. airden(mt).gt.0.)then
         densv=airden(mt)
         vissv=vis(mt)
         endif
         enddo

         do ml=1,nz
         f(ml)=f(ml)/float(isay)
         if(f(ml).lt.fmin)then
         khabit=ml
         if(f(ml).gt.0.6)khabit=6
         fmin=f(ml)
         endif
         enddo


         do km=1,ny
         do kl=1,nm
         if(ip05(km,kl).gt.0)then
         if(ip05(km,kl).eq.1)vmed05(km,kl)=vel05(km,kl)
         if(ip05(km,kl).eq.2)vmed05(km,kl)=vel05(km,kl)/2.
         if(ip05(km,kl).gt.3)then
         lsay=ip05(km,kl)
         do kc=1,lsay
         vsort05(kc)=vsv05(km,kl,kc)
         enddo
         call sort(vsort05,nc,lsay)
         md=float(lsay)/2.
         vmed05(km,kl)=vsort05(md)
         endif
         vel05(km,kl)=vel05(km,kl)/float(ip05(km,kl))
         hvel05(km,kl)=hvel05(km,kl)/float(ip05(km,kl))
         rey05(km,kl)=vmed05(km,kl)*dmx(km,kl)*1.e-3*densv/vissv

         xlrey=0.
         if(rey05(km,kl).gt.0)then
         xlrey=alog10(rey05(km,kl))
         endif
         xlmh=0.
         xlkc=0.
         do iy=1,nn
         tyy=float(iy)
         xlmh=xlmh+xmh(iy)*(xlrey**tyy)
         xlkc=xlkc+xkc(iy)*(xlrey**tyy)
         enddo
         xlmh=10.**(xlmh)
         xlkc=10.**(xlkc)
         bmh05(km,kl)=xlmh
         bkc05(km,kl)=xlkc

         xbm05(km,kl)=(cons00*((((4*rey05(km,kl))**0.5)
     *   /delta0)+1)**2 - 1.)**2.

         xhw05(km,kl)=(cons01*((((4*rey05(km,kl))**0.5)
     *   /delta00)+1)**2 - 1.)**2.

c        mass calculation

         xmbm05(km,kl)=1.e3*pi*(vissv**2.)*xbm05(km,kl)*
     *    (arat(kl)**0.25)/(8.*gg*densv)
         xmhw05(km,kl)=1.e3*pi*(vissv**2.)*xhw05(km,kl)*
     *    (arat(kl)**0.5)/(8.*gg*densv)
         xmmh05(km,kl)=1.e3*pi*(vissv**2.)*xlmh*
     *    (arat(kl)**0.25)/(8.*gg*densv)
         xmkc05(km,kl)=1.e3*pi*(vissv**2.)*xlkc*
     *    (arat(kl)**0.25)/(8.*gg*densv)


         volbm05=xmbm05(km,kl)*1.e-3/1.e3
         volmh05=xmmh05(km,kl)*1.e-3/1.e3
         volkc05=xmkc05(km,kl)*1.e-3/1.e3
         volhw05=xmhw05(km,kl)*1.e-3/1.e3

         bscatbm05(km,kl)=(volbm05**2.)*bscat(km)
         bscatmh05(km,kl)=(volmh05**2.)*bscat(km)
         bscatkc05(km,kl)=(volkc05**2.)*bscat(km)
         bscathw05(km,kl)=(volhw05**2.)*bscat(km)
         endif
         enddo
         enddo

c        velocity - diameter (mean)

         do kl=1,nm
         jsay=0
         av(kl)=0.
         bv(kl)=0.
         do km=1,ny
         if(ip05(km,kl).ge.10)then
         jsay=jsay+1
         vln05(jsay)=alog(vel05(km,kl))
         ddln05(jsay)=alog(dmx(km,kl))
         endif
         enddo
         if(jsay.ge.nw)then
         call oreg(ddln05,vln05,ny,jsay,a,b,r,sd)
         av(kl)=exp(a)
         bv(kl)=b
         endif
         enddo

c        velocity - diameter (median)

         do kl=1,nm
         jsay=0
         ac(kl)=0.
         bc(kl)=0.
         do km=1,ny
         if(ip05(km,kl).ge.10)then
         jsay=jsay+1
         cln05(jsay)=alog(vmed05(km,kl))
         dcln05(jsay)=alog(dmx(km,kl))
         endif
         enddo
         if(jsay.ge.nw)then
         call oreg(dcln05,cln05,ny,jsay,a,b,r,sd)
         ac(kl)=exp(a)
         bc(kl)=b
         endif
         enddo

c        reynolds - diameter

         do kl=1,nm
         jsay=0
         ar(kl)=0.
         br(kl)=0.
         do km=1,ny
         if(ip05(km,kl).ge.10)then
         jsay=jsay+1
         rln05(jsay)=alog(rey05(km,kl))
         drln05(jsay)=alog(dmx(km,kl))
         endif
         enddo
         if(jsay.ge.nw)then
         call oreg(drln05,rln05,ny,jsay,a,b,r,sd)
         ar(kl)=exp(a)
         br(kl)=b
         endif
         enddo

c        best - diameter

         do kl=1,nm
         jsay=0
         ax1(kl)=0.
         bx1(kl)=0.
         ax2(kl)=0.
         bx2(kl)=0.
         ax3(kl)=0.
         bx3(kl)=0.
         ax4(kl)=0.
         bx4(kl)=0.
         do km=1,ny
         if(ip05(km,kl).ge.10)then
         jsay=jsay+1
         b1ln05(jsay)=alog(xbm05(km,kl))
         b2ln05(jsay)=alog(xhw05(km,kl))
         b3ln05(jsay)=alog(bkc05(km,kl))
         b4ln05(jsay)=alog(bmh05(km,kl))
         dbln05(jsay)=alog(dmx(km,kl))
         endif
         enddo
         if(jsay.ge.nw)then
         call oreg(dbln05,b1ln05,ny,jsay,a,b,r,sd)
         ax1(kl)=exp(a)
         bx1(kl)=b
         call oreg(dbln05,b2ln05,ny,jsay,a,b,r,sd)
         ax2(kl)=exp(a)
         bx2(kl)=b
         call oreg(dbln05,b3ln05,ny,jsay,a,b,r,sd)
         ax3(kl)=exp(a)
         bx3(kl)=b
         call oreg(dbln05,b4ln05,ny,jsay,a,b,r,sd)
         ax4(kl)=exp(a)
         bx4(kl)=b
         endif
         enddo

c        mass - diameter

         do kl=1,nm
         jsay=0
         am1(kl)=0.
         bm1(kl)=0.
         am2(kl)=0.
         bm2(kl)=0.
         am3(kl)=0.
         bm3(kl)=0.
         am4(kl)=0.
         bm4(kl)=0.
         do km=1,ny
         if(ip05(km,kl).ge.10)then
         jsay=jsay+1
         x1ln05(jsay)=alog(xmbm05(km,kl))
         x2ln05(jsay)=alog(xmhw05(km,kl))
         x3ln05(jsay)=alog(xmkc05(km,kl))
         x4ln05(jsay)=alog(xmmh05(km,kl))
         dln05(jsay)=alog(dmx(km,kl))
         endif
         enddo
         if(jsay.ge.nw)then
         call oreg(dln05,x1ln05,ny,jsay,a,b,r,sd)
         am1(kl)=exp(a)
         bm1(kl)=b
         call oreg(dln05,x2ln05,ny,jsay,a,b,r,sd)
         am2(kl)=exp(a)
         bm2(kl)=b
         call oreg(dln05,x3ln05,ny,jsay,a,b,r,sd)
         am3(kl)=exp(a)
         bm3(kl)=b
         call oreg(dln05,x4ln05,ny,jsay,a,b,r,sd)
         am4(kl)=exp(a)
         bm4(kl)=b
         endif
         enddo

         do kl=1,nm

         write(20,38)iyear,jd,kh,mn,(vmed05(km,kl),km=1,ny)
         write(21,38)iyear,jd,kh,mn,(hvel05(km,kl),km=1,ny)
         write(22,19)iyear,jd,kh,mn,isay,(ip05(km,kl),km=1,ny)

         write(40,35)iyear,jd,kh,mn,(xmbm05(km,kl),km=1,ny)
         write(41,35)iyear,jd,kh,mn,(xmhw05(km,kl),km=1,ny)
         write(42,35)iyear,jd,kh,mn,(xmkc05(km,kl),km=1,ny) 
         write(43,35)iyear,jd,kh,mn,(xmmh05(km,kl),km=1,ny) 

         write(50,35)iyear,jd,kh,mn,(bscatbm05(km,kl),km=1,ny)
         write(51,35)iyear,jd,kh,mn,(bscathw05(km,kl),km=1,ny)
         write(52,35)iyear,jd,kh,mn,(bscatkc05(km,kl),km=1,ny)
         write(53,35)iyear,jd,kh,mn,(bscatmh05(km,kl),km=1,ny)

         write(24,43)iyear,jd,kh,mn,jsay,khabit,fmin,ac(kl),bc(kl),
     *   ar(kl),br(kl),am1(kl),bm1(kl),am2(kl),bm2(kl),
     *   am3(kl),bm3(kl),am4(kl),bm4(kl),ax1(kl),
     *   bx1(kl),ax2(kl),bx2(kl),ax3(kl),bx3(kl),ax4(kl),bx4(kl)
         enddo

 35      format(4i6,120e14.6)
 38      format(4i6,120f15.3)
 19      format(5i6,120i5)
 43      format(6i5,f10.3,11(e14.6,f10.2,e14.6,f10.2))
         endif
         mn=mn+5
         mn1=mn1+5
         enddo
         kh=kh+1
         enddo


         print*,'05'

 333     continue

         enddo

         stop
         end
