anchor.diff:        if (ksts.eq.1.or.ksts.eq.3) 
        1    open(unit=23,file=stsf,status='unknown')
         if (knor.eq.1) open(unit=24,file=norf,status='unknown')
  -      if (kshr.eq.1) open(unit=25,file=shrf,status='unknown')
  +      if (kshr.ge.1) open(unit=25,file=shrf,status='unknown')
         if (kcou.eq.1) open(unit=26,file=coulf,status='unknown')
   C     Initialize progress indicator
         if (prog.eq.1) then
  
positive.code:               endif
            endif
  C         Resolve stresses onto target fault planes
            if (kshr.eq.1.or.knor.eq.1.or.kcou.eq.1) then
                call coulomb(coul,norml,shear,stress,trgstr,trgdip,trgrak,
       1                     frict)
                if (kshr.eq.1.and.noz.eq.0) then
                    write(25,9992) stlo,stla,stdp*1d-3,shear
                elseif (kshr.eq.1.and.noz.eq.1) then
                    write(25,*) stlo,stla,shear
                endif
                if (knor.eq.1.and.noz.eq.0) then
  
negative.code:      1                          + dcos(str-pi/2.0d0)*dsin(rak)*dcos(dip)
        r(3) =                                         dsin(rak)*dsin(dip)
  C Project shear component on slip vector
        if (kshr.eq.2) then
            shear = dsqrt(s(1)*s(1) + s(2)*s(2) + s(3)*s(3))
        else
            shear = s(1)*r(1) + s(2)*r(2) + s(3)*r(3)
        endif
  C Coulomb stress (recall sign convention: pos = dilation)
        coul = shear + frc*norml
        RETURN
  