c sed.f - determines tauc, settling velocity 
c settling velocity based on Dietrich, 1982. Settling Velocity of
c	Natural Particles. Water Resources Research. 18(6):1615-1625.
	real nu, d
	character ans1,ans2,d_type
        character*50 infil
c       open input file
	if (iargc().gt.0) then
	    call getarg (1,infil)
	    lunin=2
	    open (2,file=infil,status='unknown')
	    rewind 2
	end if
	g=980.
c	write (6,"('calculate critical shear stress? (y or n): '$)")
c	read (5,*) ans1
	ans1='y'
c	write (6,"('calculate settling velocity? (y or n): '$)")
c	read (5,*) ans2
	ans2='y'
c 	write (6,"('diameter (cm): '$)")
c 	read (5,*) d
c	write (6,"('sediment density (g/cc): '$)")
c	read (5,*) rhos
	rhos=2.65
c	write (6,"('water density (g/cc): '$)")
c	read (5,*) rho
	rho=1.02
c	if (rho.eq.0.) rho = 1.00
c 	write (6,"('viscosity (cm2/s): '$)")
c 	read (5,*) nu
        nu = 0.013
 	if (ans2.eq.'y') then
c	    write (6,"('corey shape factor:  '$)")
c	    read (5,*) csf
	    csf=.7
c	    write (6,"('Powers roundness scale (0-6; 3.5 is normal;
c     1  	    6.0 is for sphere):  '$)")
c	    read (5,*) p
	    p=3.5
 	end if
	a2=(rhos-rho)/rho
	if (iargc().gt.0) then
c		read in whether phi sizes or cm sizes entered
		read (2,*) d_type
	else
		write (6,*) 'Enter p if using phi sizes, d if using cms'
		read (5,'(a)') d_type
	endif
	if (d_type.eq.'p'.or.d_type.eq.'P') then
		d_type = 'p'
	else 
		d_type = 'd'
	endif
c
c
c       start loop here - ckh
        do 200, k=1,100
c
	if (iargc().gt.0) then
            read (2,*,END=99) d
	    if (d_type.eq.'p') d=phitod(d)
        elseif (d_type.eq.'d') then
       	    write (6,"('diameter (cm), enter 999 to quit: '$)")
	    read (5,*) d
	else
       	    write (6,"('diameter (phi), enter 999 to quit: '$)")
	    read (5,*) d
	    if (d_type.eq.'p') d=phitod(d)
        endif
	if (d.eq.999.) goto 99
        if (d.eq.0.0) goto 99  
c
      	zeta=a2*g*d**3/nu**2
      	zetl=alog10(zeta)
	if (ans1.eq.'y') then
	    if (zetl.le.3.6) tst=.00056*zetl**3+.00115*zetl*zetl
     1   	    -0.0316*zetl+.10455
	    if (zetl.gt.3.6) tst=-.00070*zetl**3+.01204*zetl*zetl
     1   	    -0.0584*zetl+.11918
	    tauc=tst*(rhos-rho)*g*d
	    ustc=sqrt(tauc/rho)
c	    write (6,10) tauc,ustc,zeta,tst
10	    format ('tauc =',f7.2,'dy/cm2   u*c=',f7.2,'cm/s',/,
     1          '    (zeta*=',e12.5'  tau*=',f7.3')')
	end if
	if (ans2.eq.'y') then
c     	    wstl = -.04079*(zetl**2)+.97489*zetl-1.18191
c     	    wset = (10.0**wstl)*(nu/d)
	    a=-3.76715+1.92944*zetl-0.09815*zetl*zetl-0.00575*zetl**3
     1          +0.00056*zetl**4
	    omc=1.-csf
	    b=alog10(1.-omc/0.85)-omc**2.3*tanh(zetl-4.6)+.3*(.5-csf)*
     1          omc*omc*(zetl-4.6)
	    c=(.65-(csf/2.83)*tanh(zetl-4.6))**(1.+(3.5-p)/2.5)
	    wstr=c*10**(a+b)
	    wset=(a2*g*nu*wstr)**(1./3.)
	    write (6,11) wset
11	    format ('settling velocity =',f7.2,'cm/s')
            write (4, *) d, wset
	end if
 200    continue 
 99     stop 
	end

        function phitod (phi)
c
c       returns the diameter size in cm for the input phi size
c
        real phitod, phi, mm
c
        mm = 2.0**(-1.0*phi)
        phitod = mm/10.0
        return
        end