The following gives a result with a warning:

Clear[x, a, b, d];
data = {{0.1, 1}, {0.5, 0.99245283}, {1, 0.981132075}, {2, 
    0.957629013}, {3, 0.94647978}, {5, 0.892872862}, {8, 
    0.825362586}, {10, 0.794884307}, {15, 0.737820069}, {30, 
    0.688456089}, {50, 0.710852149}, {100, 0.686350949}, {300, 
    0.630188679}, {500, 0.622243133}, {1500, 0.61197978}};
model = Piecewise[{{0, x > a}, {1, 
     x < b}, {(x^(d - 3) - b^(d - 3))/(a^(d - 3) - b^(d - 3)), True}}];
nlm = NonlinearModelFit[data, model, {{a, 1}, {b, 10^6}, {d, 2.8}}, x]

So. The case is that S is the relative saturation, b is the moment in which the saturation is 0 and a is the moment in which the saturation is 1, so that if x(that is the potential in a specific water content)>b saturation is equal to 0, and if the potential in a specific water content (x) is minor than a, so that the saturation is 1.

I think your code has work succesfully, I only changed the way you wrote the conditionals.

I am no expert. It may have something to do with the fact that your model is not differentiable. Sorry I can't help you any further. Here is a code that gives an output, even though with a warning:

data = {{0.1, 1}, {0.5, 0.99245283}, {1, 0.981132075}, {2, 
    0.957629013}, {3, 0.94647978}, {5, 0.892872862}, {8, 
    0.825362586}, {10, 0.794884307}, {15, 0.737820069}, {30, 
    0.688456089}, {50, 0.710852149}, {100, 0.686350949}, {300, 
    0.630188679}, {500, 0.622243133}, {1500, 0.61197978}};
model = Piecewise[{{1, x < a}, {0, 
     x > b}, {(x^(d - 3) - b^(d - 3))/(a^(d - 3) - b^(d - 3)), True}}];
nlm = NonlinearModelFit[data, model, {{a, 1}, {b, 10^6}, {d, 2.8}}, x]

Your SAS code has several syntax errors. "^" should be "**", "b=10^6" should be "b=10E6", and the last "end;" generates an error as it is not needed. Those errors seem a bit more than cut-and-paste issues. Does the code really run in SAS?

Following the good suggestion from Gianluca Gorni you'll find the following result:

Note that the estimate of b is huge and the coefficient -1.02625 x 10^(-11) is essentially zero. That leaves an equation that is linear in the log of x. Given that and that the plot of the fit and the data looks like the following:

this screams for taking the logs of both the dependent and independent variables. From the initial description of the issue of wanting a piecewise fit, one possibility is to have the log of the prediction be zero for $x < a$ and $b + (d-3)\log(x)$ otherwise. Because you probably want the fit to be continuous, you'll want $b + (d-3)\log(x) = 0$ when $x=a$. That means $b=-(d-3)\log(a)$. The associated Mathematica code is

model = Piecewise[{{0, x < a}, {-(d - 3) Log[a] + (d - 3) Log[x], True}}];
logY = Log[data[[All, 2]]];
nlm = NonlinearModelFit[Transpose[{data[[All, 1]], logY}], {model, a > 0.1}, {{a, 0.9}, {d, 2.9}}, x];
nlm["BestFitParameters"]
(* {a -> 0.624483, d -> 2.92726} *)
nlm["BestFit"]
Show[ListLogLogPlot[data, PlotRange -> {{0.1, 1600}, {0.5, 1.1}}, PlotStyle -> {PointSize[0.02], Red}, 
  PlotRangeClipping -> False, Frame -> True],
  LogLogPlot[Exp[nlm[x]], {x, 0.1, 1600}]]

The fit is not great but the residuals are a bit better behaved than not taking logs.