Smoothen the result of FindRoot

I have a set of code for which it involves finding the corresponding c for each a (although I will give a value of a later) and z using the constraint toroot[a,c,z] and then substituting c back into the final expression functionS[a,z]. I also have another function for which there is a change of variable functionSR[l,z] where a->l-0.01.

d = 3;
zh = 1.5;
toroot[a_, c_?NumericQ, z_] := a - NIntegrate[(c z^(d + 1) x^d)/((1 - ((z x)/zh)^(d + 1)) (1 - c^2 (z x)^(2 d)))^(1/2), {x, 0, 1}, MaxRecursion -> 5, PrecisionGoal -> 4, Method -> "LocalAdaptive"]
cz[a_?NumericQ, z_?NumericQ] := c /. FindRoot[toroot[a, c, z], {c, 0.0009, 0.0000001, 10000}, WorkingPrecision -> 5]
intS[a_?NumericQ, z_?NumericQ] := NIntegrate[With[{b = z/zh}, (((-1)/(d - 1)) cz[a, z]^2 z^(2 d)) x^d ((1 - (b x)^(d + 1))/(1 - cz[a, z]^2 (z x)^(2 d)))^(1/2) - ((b^(d + 1) (d + 1))/(2 (d - 1))) x ((1 - cz[a, z]^2 (z x)^(2 d))/(1 - (b x)^(d + 1)))^(1/2) + (b^(d + 1) x)/((1 - (b x)^(d + 1)) (1 - cz[a, z]^2 (z x)^(2 d)))^(1/2)], {x, 0, 1}, MaxRecursion -> 5, PrecisionGoal -> 4, Method -> "LocalAdaptive"]
functionS[a_, z_] = ((-((1 - cz[a, z]^2 z^(2 d)) (1 - (z/zh)^(d + 1)))^(1/2)/(d - 1)) + intS[a, z] + 1)/(z^(d - 1));
functionSR[l_, z_] = Replace[functionS[a, z], a -> (l - 0.01), Infinity];

My problem is when I try to find the minimum of functionS[a,z] and functionSR[l,z] for some a and l, say a=1 and l=1, it gives me an error. I think it is connected to the behavior of c when a=1 or l=1.

In[23]:= FindMinimum[functionS[1, z], {z, 1.2, 1.5}] // 
  Quiet // AbsoluteTiming
FindMinimum[functionSR[1, z], {z, 1.2, 1.5}] // Quiet // AbsoluteTiming

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.697475}. NIntegrate obtained 0.000944548 -0.00149313 I and 0.0006178735732839699` for the integral and error estimates.

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.697475}. NIntegrate obtained 0.000944548 -0.00149313 I and 0.0006178735732839699` for the integral and error estimates.

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.697475}. NIntegrate obtained 0.000949747 -0.00149122 I and 0.000620731102746343` for the integral and error estimates.

During evaluation of In[23]:= General::stop: Further output of NIntegrate::ncvb will be suppressed during this calculation.

During evaluation of In[23]:= FindRoot::reged: The point {1.70561} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= FindRoot::reged: The point {1.70561} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= FindRoot::reged: The point {1.70561} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= General::stop: Further output of FindRoot::reged will be suppressed during this calculation.

During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.436961 -1.38189 I is not a real number at {z} = {1.2}.

During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.436961 -1.38189 I is not a real number at {z} = {1.2}.

Out[23]= {0.760891, FindMinimum[functionS[1, z], {z, 1.2, 1.5}]}

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.699811}. NIntegrate obtained 0.00286247 -0.0000971587 I and 0.0005426332486649041` for the integral and error estimates.

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.699811}. NIntegrate obtained 0.00286247 -0.0000971587 I and 0.0005426332486649041` for the integral and error estimates.

During evaluation of In[23]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 5 recursive bisections in x near {x} = {0.699811}. NIntegrate obtained 0.00286812 -0.0000961916 I and 0.0005442259497809905` for the integral and error estimates.

During evaluation of In[23]:= General::stop: Further output of NIntegrate::ncvb will be suppressed during this calculation.

During evaluation of In[23]:= FindRoot::reged: The point {1.68855} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= FindRoot::reged: The point {1.68855} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= FindRoot::reged: The point {1.68855} is at the edge of the search region {1.0000*10^-7,10000.} in coordinate 1 and the computed search direction points outside the region.

During evaluation of In[23]:= General::stop: Further output of FindRoot::reged will be suppressed during this calculation.

During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.439434 -1.36539 I is not a real number at {z} = {1.2}.

During evaluation of In[23]:= FindMinimum::nrnum: The function value 0.439434 -1.36539 I is not a real number at {z} = {1.2}.

Out[24]= {0.771827, FindMinimum[functionSR[1, z], {z, 1.2, 1.5}]}

For a=0.1, the plot is much more smooth

For a=1, the plot contains more bumps

Is my code badly written to extract c? Are there any changes that can be done? I have read somewhere that Reduce can also be used instead of FindRoot but I am still figuring it out. Also, is using LocalAdaptive as a method for NIntegrate suitable here?

UPDATE:
Please note of the typo, I have corrected it. In the plots before, I wrote c=0.1 and c=1 but it should be a=0.1 and a=1.

The expressions of my problem is given by,

$$a = c z_s^{d+1}int_0^1 dx frac{x^d}{sqrt{(1-(z_s/z_h)^{d+1} x^{d+1})(1-c^2 z_s^{2d} x^{2d})}} tag{1}label{1}$$

begin{align}
S &= frac{1}{4 z_s^{d-1}}Bigg(1 -frac{sqrt{(1-c^2 z_s^{2d})(1-b^{d+1})}}{d-1} – frac{1}{d-1} c^2 z_s^{2d} int^1_0 dx x^d sqrt{frac{(1-(b x)^{d+1})}{(1-c^2(z_s x)^{2d})}}\
& -frac{b^{d+1}(d+1)}{2(d-1)} int^1_0 dx x sqrt{frac{(1-c^2(z_s x)^{2d})}{(1-(b x)^{d+1})}}\
& + b^{d+1}int^1_0 dx frac{x}{sqrt{(1-(b x)^{d+1})(1-c^2(z_s x)^{2d})}}Bigg) tag{2}label{2}
end{align}

where $b=frac{z_s}{z_h}$ and note that $c=c(z_s)$ (c=c[z]) although in the code c=c[a,z], $c$ should only depend on $z_s$ (z) since $a$ will be specified in the end.

Also, maybe there is a better way to design finding $c$. Actually, I can have another constraint where $frac{dS}{dz_s} = 0$ (that is because in the end I need to minimize $S$ with respect to $z_s$) and maybe the derivative of $eqref{1}$ with respect to $z_s$, so that these can be used to find $c$?

toroot[a_, c_?NumericQ, z_] := a - NIntegrate[((1 - xd (z /zh)^(d + 1)) 
   (1 - c^2 xd^(2 d/(d + 1))  z^(2 d)))^(-1/2), {xd, 0, 1}, Method -> {"GlobalAdaptive",
    "SingularityHandler" -> "DoubleExponential"}] (c z^(d + 1))/4

cz[a_?NumericQ, z_?NumericQ] := c /. 
    FindRoot[toroot[a, c, z], {c, .5 z^-3, .6 z^-3/2, 0, z^-3}]

Plot[Check[cz[1, z], Null], {z, 1.42, zh}, AxesLabel -> {z, c}, 
    ImageSize -> Large, LabelStyle -> {15, Bold, Black}]

LogPlot[Check[cz[.1, z], Null], {z, .2, zh}, AxesLabel -> {z, c}, 
    ImageSize -> Large, LabelStyle -> {15, Bold, Black}]