Retrieving forthcoming step size from NDSolve

Question

I'm solving a system of ODEs with a random component, which changes at each calculation step and depends on the step size. So I'm trying to realize something like this: (1) size of the forthcoming step is retrieved from NDSolve (2) random component is calculated (3) only then the step is evaluated with a corresponding random component. Unfortunately, I don't know how to squeeze between the moment, when NDSolve selects step size and when it evaluates the step. StepMonitor (commented in the code) doesn't allow that, as I understand, since it returns information after the step is evaluated.

ps1 Fixed step method is 5x longer, which means hours of difference, since I need 10^5-10^6 evaluations for a single final plot, so this would be a really undesirable option.

ps2 Using previous step size also doesn't seem good due to abrupt changes at some points

Here is the code:

Remove["Global`*"]

(* PARAMETERS *)
γ = 
 1.76*10.^7*10.^4 (*T^-1s^-1*);    α = 0.06;   θsh = \
0.3;   Ms = 7.*10.^5; (*A/m*)   tFM = 2.*10.^-9; (*m*)
currAMP = 30.*10.^-3 (*A*);   sHM = 
 10.*10.^-6*3.*10.^-9(* m^2, tHM = 3 nm *);
Hamp = 0.03; Hk = 0.1;
a = 0. °; b = 0. °; c = 0. °; Ang = 
 0. °; 
temp = 3; rad = 20.*10.^-9; (*Thermal effect*)
kTL = 2.;
tSt = 0.;  tFin = 1.5*10.0^-8;  tStep = 1.*10.0^-12;

(* CURRENT PULSE *)
tPulse = 5.*10^-9; tRise = 0.5*10^-9;
jcurr[t_] := Piecewise[{{1., -1 <= t <= tPulse}}, 0.];
Cjconst = (3.29375`*^-16) (θsh/Ms/tFM)*currAMP/sHM;
Cj[t_] := Cj[t] = jcurr[t]*Cjconst;

(* FIELD COMPONENTS *)
σf[temp_, rad_] := Sqrt[(
 2.*α*1.38*10.^-23*temp)/(γ*Ms*tFM*3.1415*rad^2*
  tStep)];  σ = σf[temp, rad];  (*stand.dev*)
RandComp[] := 
 RandomVariate[
  NormalDistribution[0., σ]]/Sqrt[3.]; (*random component*)
Happ = ({
    {Cos[c], -Sin[c], 0.},
    {Sin[c], Cos[c], 0.},
    {0., 0., 1.}
   }).({
    {Cos[b], 0., Sin[b]},
    {0., 1., 0.},
    {-Sin[b], 0., Cos[b]}
   }).({
    {1., 0., 0.},
    {0., Cos[a], -Sin[a]},
    {0., Sin[a], Cos[a]}
   }).( {
    {Hamp*Cos[Ang]},
    {Hamp*(-Sin[Ang])},
    {0.}
   }); (*permanent applied component*)
Hrand[t_] := Hrand[t] = {RandComp[], RandComp[], RandComp[]};
H[t_, mz_] := Flatten[Happ + Hrand[t] + {0., -kTL*Cj[t], Hk*mz}]

(* ODE SYSTEM *)
t0 = 0;
NDSolve[( {
   {mX'[t] == γ (-mY[t]*H[t, mZ[t]][[3]] + 
        mZ[t]*H[t, mZ[t]][[
          2]] - α (mY[
             t] (mX[t]*H[t, mZ[t]][[2]] - mY[t]*H[t, mZ[t]][[1]]) - 
           mZ[t] (mZ[t]*H[t, mZ[t]][[1]] - mX[t]*H[t, mZ[t]][[3]])) + 
        Cj[t]*mY[t]*mX[t] + α*Cj[t]*mZ[t]), mX[0] == 0.},
   {mY'[t] == γ (-mZ[t]*H[t, mZ[t]][[1]] + 
        mX[t]*H[t, mZ[t]][[
          3]] - α (mZ[
             t] (mY[t]*H[t, mZ[t]][[3]] - mZ[t]*H[t, mZ[t]][[2]]) - 
           mX[t] (mX[t]*H[t, mZ[t]][[2]] - mY[t]*H[t, mZ[t]][[1]])) - 
        Cj[t] (mX[t]^2 + mZ[t]^2)), mY[0] == 0.},
   {mZ'[t] == γ (-mX[t]*H[t, mZ[t]][[2]] + 
        mY[t]*H[t, mZ[t]][[
          1]] - α (mX[
             t] (mZ[t]*H[t, mZ[t]][[1]] - mX[t]*H[t, mZ[t]][[3]]) - 
           mY[t] (mY[t]*H[t, mZ[t]][[3]] - mZ[t]*H[t, mZ[t]][[2]])) + 
        Cj[t]*mY[t]*mZ[t] - α*Cj[t]*mX[t]), mZ[0] == 1.}
  } ), {mX, mY, mZ}, {t, tSt, tFin}, StartingStepSize -> 10.^-16, 
 Method -> {"Projection", Method -> "Automatic", 
   "Invariants" -> (mX[t]^2 + mY[t]^2 + mZ[t]^2 == 1)}, 
 AccuracyGoal -> 15, PrecisionGoal -> 10000 (*,StepMonitor:>Block[{},
 δt=t-t0;t0=t;
 Print["Current time: ",t];
 Print["Current step: ",δt];
 Print["mX: ",mX[t]];
 Print["mY: ",mY[t]];
 Print["mZ: ",mZ[t],"
 "];
 AbsoluteTiming[Pause[1]]
 ]*)
 ]

Thanks!