anchor.diff:        else
            do s = 1, stations - 1
               if (width(s) .ne. 0.d0) then
  -               forcefactor = rho * vel(s) ** 2.d0 / depth(s) 
  +               forcefactor = rho * vel(s) ** 2 / depth(s) 
        +              * (1.d0 - 2.d0 * depth(s) / width(s))
                  latforce(s) = 0.d0
                  lag(s) = 0.d0
  -               xmagcurvep1 = sign(curvature(s + 1), 1.d0) 
  -               xmagcurve = sign(curvature(s), 1.d0) 
  +               xmagcurvep1 = abs(curvature(s + 1)) 
  +               xmagcurve = abs(curvature(s)) 
                  if ((curvature(s + 1) * curvature(s) .gt. 0.d0 .and.
        +           xmagcurvep1 .gt. xmagcurve .and.
        +           xmagcurvep1 * width(s) .le. 2.d0) .or.
  
positive.code:      +       leftdepth(*), xs(*),
       +       tauwall(*), gaussfactor, xstrt,
       +       xdel, xdepth, gaussian, xdest, tenlambda, xtrmnt
        intrinsic cos, sin, exp
  c
        do s = 1, stnserod
           tauwall(s) = 0.d0
  
negative.code: c    +                    - 0.3606d0)
              rectchan = 0.5d0 * width(s) * depth(s)
              transslope(s) = transfactor * curvature(s) 
              xmagtransslope = abs(transslope(s))
              if (xmagtransslope .le. corner) then
                 Acs(s) = rectchan - width(s) ** 2 * 
       +                  xmagtransslope / 8.d0
                 rightdepth(s) = depth(s) + width(s) * 
       +                         transslope(s) / 2.d0
                 leftdepth(s) = depth(s) - width(s) * transslope(s) / 2.d0
              else
                 Acs(s) = 0.5d0 * depth(s) ** 2 / xmagtransslope
                 if (curvature(s) .lt. 0.d0) then
                    rightdepth(s) = 0.d0
                    leftdepth(s) = 2.d0 * depth(s)
  