anchor.diff:              end if
               vel(s) = flow(s) / depth(s) / width(s)
               shields = radh * slope(s) / 1.65d0 / diam(s) 
  -            if (diam(s) .lt. 0.00015) critshields = 0.080d0	! s*=1.01
  -            if (diam(s) .ge. 0.00015 .and. 
  -     +          diam(s) .lt. 0.00025) critshields = 0.052d0	! s*=2.84
  -            if (diam(s) .ge. 0.00025 .and. 
  -     +          diam(s) .lt. 0.00035) critshields = 0.039d0	! s*=5.22
  -            if (diam(s) .ge. 0.00035 .and. 
  -     +          diam(s) .lt. 0.00045) critshields = 0.035d0	! s*=8.04
  -            if (diam(s) .ge. 0.00045 .and. 
  -     +          diam(s) .lt. 0.00055) critshields = 0.034d0	! s*=11.2
  -            if (diam(s) .ge. 0.00055 .and. 
  -     +          diam(s) .lt. 0.00065) critshields = 0.033d0	! s*=14.8
  -            if (diam(s) .ge. 0.00065 .and. 
  -     +          diam(s) .lt. 0.00075) critshields = 0.034d0	! s*=18.6
  -            if (diam(s) .ge. 0.00075 .and. 
  -     +          diam(s) .lt. 0.00085) critshields = 0.034d0	! s*=22.7
  -            if (diam(s) .ge. 0.00085 .and.
  -     +          diam(s) .lt. 0.00095) critshields = 0.034d0	! s*=27.1
  -            if (diam(s) .ge. 0.00095 .and. 
  -     +          diam(s) .lt. 0.0015) critshields = 0.035d0 	! s*=31.8
  -            if (diam(s) .ge. 0.0015 .and.
  -     +          diam(s) .lt. 0.0025) critshields = 0.045d0 	! s*=89.9
  -            if (diam(s) .ge. 0.0025 .and.
  -     +          diam(s) .lt. 0.0035) critshields = 0.050d0 	! s*=165.
  -            if (diam(s) .ge. 0.0035 .and.
  -     +          diam(s) .lt. 0.0045) critshields = 0.055d0      ! s*=254.
  -            if (diam(s) .ge. 0.0045) critshields = 0.056d0      ! s*>=355
  +            if (diam(s) .lt. 0.00015) then
  +               critshields = 0.080d0 ! s*=1.01
  +            else if (diam(s) .lt. 0.00025) then
  +               critshields = 0.052d0 ! s*=2.84
  +            else if (diam(s) .lt. 0.00035) then
  +               critshields = 0.039d0 ! s*=5.22
  +            else if (diam(s) .lt. 0.00045) then
  +               critshields = 0.035d0 ! s*=8.04
  +            else if (diam(s) .lt. 0.00055) then
  +               critshields = 0.034d0 ! s*=11.2
  +            else if (diam(s) .lt. 0.00065) then
  +               critshields = 0.033d0 ! s*=14.8
  +            else if (diam(s) .lt. 0.00075) then
  +               critshields = 0.034d0 ! s*=18.6
  +            else if (diam(s) .lt. 0.00085) then
  +               critshields = 0.034d0 ! s*=22.7
  +            else if (diam(s) .lt. 0.00095) then
  +               critshields = 0.034d0 ! s*=27.1
  +            else if (diam(s) .lt. 0.0015) then
  +               critshields = 0.035d0 ! s*=31.8
  +            else if (diam(s) .lt. 0.0025) then
  +               critshields = 0.045d0 ! s*=89.9
  +            else if (diam(s) .lt. 0.0035) then
  +               critshields = 0.050d0 ! s*=165.
  +            else if (diam(s) .lt. 0.0045) then
  +               critshields = 0.055d0 ! s*=254.
  +            else
  +               critshields = 0.056d0 ! s*>=355
  +            end if
   c           approximate Engelund diagram for subcritical flow:
  -            if (shields .lt. 1.0 .and. shields .ge. 0.1) then
  +            if (shields .ge. 0.1 .and. shields .lt. 1.0 ) then
                  grainshields = 10.d0 ** (0.74d0 
        +                        * (log10(shields) + 1.03d0) ** 2.d0
        +                        - 1.18d0)
               else if (shields .ge. 1.0 .and. shields .lt. 2.0) then
                  grainshields = 0.4d0 * shields ** 2.d0
  -            else if (shields .ge. 2.0 .or. shields .lt. 0.1) then
  +            else
  +               ! if (shields .ge. 2.0 .or. shields .lt. 0.1) then
                  grainshields = shields
               end if
               corner = 2.d0 * depth(s) / width(s)
  
positive.code:      +                    curvature, dels, Acs, deln,
       +                    rightdepth, leftdepth, depth, vel, 
       +                    delx, dely, transfactor, transslope)
  C
  C  xmagtransslope is the absolute value of the transverse bedslope (transslope)
  C  Acs is cross-sectional area of the half-channel at the inside of the bend
  
negative.code:      +                      lag, latforce, dels, phi, curvature, 
       +                      depth, spreaddelta_x, spreaddelta_y,
       +                      rightdepth, leftdepth, xs, tauwall)
        implicit none
  C
        integer stnserod, stations, s, sp
        real*8 lambda(*), width(*), lag(*), 
       +       latforce(*), dels(*), phi(*), 
       +       curvature(*), depth(*), spreaddelta_x(*), 
       +       spreaddelta_y(*), rightdepth(*),
       +       leftdepth(*), xs(*),
       +       tauwall(*), gaussfactor, xstrt,
       +       xdel, xdepth, gaussian, xdest, tenlambda, xtrmnt
        intrinsic cos, sin, exp
  c
        do s = 1, stnserod
           tauwall(s) = 0.d0
  