equation solving – Speed up FindRoot

I need to calculate Sv(zm,zs,zh) for a given range of zm and zs with fixed zh. I have written a code where zm is a function of zs given by zmzs(zs, zh, tb, zmmin, zmmax) so that Sv(zm(zs),zs,zh) = Sv(zmzs(zs, zh, tb, zmmin, zmmax),zs,zh) and now I only need to give the range of zs only. To add, zh and tb are given constants, zmmin and zmmax are the min and max of the range in FindRoot in zmzs(zs, zh, tb, zmmin, zmmax).

d = 3;
ag = 10;
pg = 10;
wp = 20;
f(z_, zh_) := 1 - (z/zh)^(d + 1);
tzsint(z_?NumericQ, zs_?NumericQ, zh_?NumericQ) := Module({zr, zsr, zhr}, {zr, zsr, zhr} = Rationalize({z, zs, zh}, 0); zr^d/Sqrt(f(zr, zhr) (zsr^(2 d) - zr^(2 d))))
tzs(zs_?NumericQ, zh_?NumericQ) := Module({zsr, zhr}, {zsr, zhr} = Rationalize({zs, zh}, 0); NIntegrate(tzsint(z, zsr, zhr), {z, 0, zsr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20))
tzmint(z_?NumericQ, zm_?NumericQ, zh_?NumericQ) := Module({zr, zmr, zhr}, {zr, zmr, zhr} = Rationalize({z, zm, zh}, 0); -1/(f(zr, zhr) Sqrt(1 - (zmr^(2 d) f(zr, zhr))/(z^(2 d) f(zmr, zhr)))))
tzm(zm_?NumericQ, zh_?NumericQ) := Module({zmr, zhr}, {zmr, zhr} = Rationalize({zm, zh}, 0); NIntegrate(tzmint(z, zmr, zhr), {z, 0, zhr, zmr}, Method -> PrincipalValue, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20))
Svint1(z_?NumericQ, zm_?NumericQ, zh_?NumericQ) := Module({zr, zmr, zhr}, {zr, zmr, zhr} = Rationalize({z, zm, zh}, 0); zmr^d/(zr^d Sqrt(f(zr, zhr) zmr^(2 d) - f(zmr, zhr) zr^(2 d))))
Svint2(z_?NumericQ, zm_?NumericQ, zs_?NumericQ, zh_?NumericQ) := Module({zr, zmr, zsr, zhr}, {zr, zmr, zsr, zhr} = Rationalize({z, zm, zs, zh}, 0); (((zsr^(2 d) zmr^(2 d) zr)/zhr^(d + 1)) + (zr^d (zsr^(2 d) f(zmr, zhr) - zmr^(2 d))))/((zmr^d Sqrt(zsr^(2 d) - zr^(2 d) ) + zsr^d Sqrt(zmr^(2 d) f(zr, zhr) - zr^(2 d) f(zmr, zhr))) Sqrt((zmr^(2 d) f(zr, zhr) - zr^(2 d) f(zmr, zhr)) (zsr^(2 d) - zr^(2 d)))))
Sv(zm_?NumericQ, zs_?NumericQ, zh_?NumericQ) := Module({zmr, zsr, zhr}, {zmr, zsr, zhr} = Rationalize({zm, zs, zh}, 0); NIntegrate(Svint1(z, zmr, zhr), {z, zsr, zmr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20) + NIntegrate(Svint2(z, zmr, zsr, zhr), {z, 0, zsr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20))
tzmzs(zm_?NumericQ, zs_?NumericQ, zh_?NumericQ, tb_?NumericQ) := Module({zmr, zsr, zhr, tbr}, {zmr, zsr, zhr, tbr} = Rationalize({zm, zs, zh, tb}, 0); NIntegrate(tzmint(z, zmr, zhr), {z, 0, zhr, zmr}, Method -> PrincipalValue, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20) - NIntegrate(tzsint(z, zsr, zhr), {z, 0, zsr}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20) - tbr)
zmzs(zs_?NumericQ, zh_?NumericQ, tb_?NumericQ, zmmin_?NumericQ, zmmax_?NumericQ) := zm /. FindRoot(tzmzs(zm, zs, zh, tb) == 0, {zm, zmmin, zmmax}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20)
zszm(zm_?NumericQ, zh_?NumericQ, tb_?NumericQ, zsmin_?NumericQ, zsmax_?NumericQ) := zs /. FindRoot(tzmzs(zm, zs, zh, tb) == 0, {zs, zsmin, zsmax}, AccuracyGoal -> ag, PrecisionGoal -> pg, WorkingPrecision -> wp, MaxIterations -> 20)

The issue is calculating zmzs(zs, zh, tb, zmmin, zmmax) takes a lot of time, for example I evaluated everything at zs=9.9952388099. zmzs(zs, zh, tb, zmmin, zmmax) takes around 13 secs to calculate, Sv(zmzs(zs, zh, tb, zmmin, zmmax),zs,zh) takes 15 secs to calculate. Using the value of zm=13.127957300691564691 calculated from zmzs(zs, zh, tb, zmmin, zmmax) and plugging it directly to Sv(zm,zs,zh) takes only around 2 secs. This means that Sv itself calculates fairly quick, it is really zmzs(zs, zh, tb, zmmin, zmmax) that takes quite some time.

In the end I want to plot Sv for (zs,1.5,9.9952388097) which will definitely take a lot of time, in fact I tried plotting it and it has already been 4 hrs but so far no result yet!!!

zmzs(9.9952388099, 10, 0.100000100000003174, 10.00099982369714445, 13.1279573)//AbsoluteTiming
{13.1723073, 13.127957300691564691}

Sv(zmzs(9.9952388099, 10, 0.100000100000003174, 10.00099982369714445, 13.1279573), 9.9952388099, 10)//Chop//AbsoluteTiming
{15.0956025, 0.010820232314756221363}

Sv(13.127957300691564691, 9.9952388099, 10)//Chop//AbsoluteTiming
{1.9989049, 0.010820232314756221363}

Is there any way to speed up the code in order to plot out Sv in a shorter amount of time?