Iterating with NDSolveValue in a Do-loop [duplicate]

Question

I am trying to get a differential solution to converge. There are two separate time periods that I am solving for and the entire solution is periodic. I want to use the output from one as the input for the next and so on and so forth. I have attempted to do this with a Do loop but it does not seem to be working. I don't know if the Do loop is set up to work this way or if I should use another approach. Please let me know if you have any suggestions on how to get this to work.

Here is the code:

<< "NDSolve`FEM`"

ts = .000250; tsl = .000250; (* m *)
ρs = 3980;  ρsl = 958; (* kg/m3 *)
ks = .035;  ksl = .00067; (* kW/m/K *)
cs = .75; csl = 4.22; (* kJ/kg/K *)

ρ = If[0 <= z < ts, ρs, ρsl];
k = If[0 <= z < ts, ks, ksl];
c =  If[0 <= z < ts, cs, csl];
q = 28;

td = 0.060; tg = 0.01;

eqn1 = k*\!\(
\*SubscriptBox[\(∂\), \(z\)]\(T1[z]\)\) + q;
Tbl = 100;
Subscript[Γ1, D] = 
  DirichletCondition[T1[z] == Tbl, z == ts + tsl];

BCr = NDSolveValue[{eqn1 == 0, Subscript[Γ1, D]}, 
   T1, {z, 0, ts + tsl}];

Ti[z_] := \[Piecewise] {
    {BCr[z], 0 <= z <= ts},
    {100, True}
   };

eqn2 = ρ*c*\!\(
\*SubscriptBox[\(∂\), \(t\)]\(T[t, z]\)\) - k*\!\(
\*SubscriptBox[\(∂\), \(z, z\)]\(T[t, z]\)\);
Subscript[Γ, D] = 
  DirichletCondition[T[t, z] == Tbl, z == ts + tsl];
Subscript[Γ, N] = NeumannValue[-q, z == 0];

solnd = NDSolveValue[{eqn2 == Subscript[Γ, N], 
    Subscript[Γ, D], T[0, z] == Ti[z]}, 
   T, {t, 0, td}, {z, 0, ts + tsl}, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "MeshOptions" -> {"MaxCellMeasure" -> (ts + tsl)/1011}}}}];

eqn3 = ρ*c*\!\(
\*SubscriptBox[\(∂\), \(t\)]\(Tg[t, z]\)\) - k*\!\(
\*SubscriptBox[\(∂\), \(z, z\)]\(Tg[t, z]\)\);
TiG[z_?NumericQ] := solnd[td, z];
Tits := solnd[td, ts]
Subscript[Γd, N] = NeumannValue[-q, z == 0];

solng = NDSolveValue[{eqn3 == Subscript[Γd, 
     N], (Tg[t, z] /. z -> ts) == (-300*t) + Tits, 
    Tg[0, z] == TiG[z]}, Tg, {t, 0, tg}, {z, 0, ts}, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "MeshOptions" -> {"MaxCellMeasure" -> (ts)/1011}}}}];

Do[
 Ti[z_?NumericQ] := \[Piecewise] {
    {solng[tg, z], 0 <= z <= ts},
    {100, True}
   };

 solnd = NDSolveValue[{eqn2 == Subscript[Γ, N], 
    Subscript[Γ, D], T[0, z] == Ti[z]}, 
   T, {t, 0, td}, {z, 0, ts + tsl}, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "MeshOptions" -> {"MaxCellMeasure" -> (ts + tsl)/1011}}}}];

 TiG[z_?NumericQ] := solnd[td, z];
 Tits := solnd[td, ts];

 solng = NDSolveValue[{eqn3 == Subscript[Γd, 
     N], (Tg[t, z] /. z -> ts) == (-300*t) + Tits, 
    Tg[0, z] == TiG[z]}, Tg, {t, 0, tg}, {z, 0, ts}, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "MeshOptions" -> {"MaxCellMeasure" -> (ts)/1011}}}}];
 , {i, 10}]

Animate[Plot[solnd[t, z], {z, 0, ts + tsl}, 
  GridLines -> {{ts, ts + .00005}, {0}}, PlotRange -> {99, 111}], {t, 
  0, td}]
Plot[solnd[t, ts], {t, 0, td}, GridLines -> {{0}, {107}}, 
 PlotRange -> Full]

Animate[Plot[solng[t, z], {z, 0, ts}, PlotRange -> {100, 111}], {t, 0,
   tg}]
Plot[solng[t, ts], {t, 0, tg}, PlotRange -> Full]