How to make gaussian pulse move using Boundary Conditions (one spatial dimension)?

Question

I'm trying to solve a system of normalized coupled PDE's to model Raman scattering. Essentially, a gaussian light pulse of 300ps duration enters the raman active medium at ζ = 0 which then drives the molecules and, as a result, gets scattered to a lower frequency (1st Stokes). I'm interested in how the spatial field profile for both the main pulse and stokes pulse evolves with time. I would like to reflect the fact that the pulse is longer than the 8cm crystal and only solve the coupled PDE's within a smaller subdomain. Also, I made some attempt to get the pulse to move through the medium using Dirichlet BC's at ζ = 0 but the frames that are generated show that only the peak is affected and the pulse still remains stationary. I appreciate any help on this as I've been stuck on these issues for awhile now and have looked over Mathematica documentation many times now. Here is the code:

ClearAll
(*Physical quantities and constants*)
eps0 = 8.85*10^-12;
T2 = 25*10^-12;
c = 3*10^8;
λf = 800*10^-9;
λs = 873*10^-9;
gs = 16*10^-11;
ωl = 2*π*c/λf;
ωs = 2*π*c/λs;
tp = 300*10^-12;
pdiam = .003; (*Diameter of pump in [m]]*)
w0 = pdiam/2 ;
n[λ_] := Sqrt[
 2.4069 + .01992/(λ^2 - .03773) - 0.006166*λ^2]
μl = n[.8];
μs = n[.873];
κQ = I*c*μl*μs/(16*π*ωs);
κS = I*ωs/μs;
κL = I*ωl/μl;
Γl = 3.8*10^-4;
Γs = 3.8*10^-4;
vp = 3*10^8/n[.8];
vs = 3*10^8/n[.873];
ZRL = c*tp; 
vpnorm = vp*tp/ZRL;
vsnorm = vs*tp/ZRL;
T = T2/tp;
(*End of Physical Quantities and Constants*)



f[ζ_] := E^(-2 Log[2]*ζ^2) 
pdes = {D[M[τ, ζ], τ] - 
     vsnorm*D[M[τ, ζ], ζ] - 
     vsnorm*(2*ζ/(1 + ζ^2))*M[τ, ζ] + 
     M[τ, ζ]/T == (κQ/T)*1/(1 + ζ^2)*
     EL[τ, ζ]*ES[τ, ζ]\[Conjugate], 
   D[ES[τ, ζ], τ] + Γs*
      ES[τ, ζ] == (1/3)*κS*gs*EL[τ, ζ]^2*
     ZRL*vsnorm*EL[τ, ζ]*
     M[τ, ζ]\[Conjugate], (vpnorm - vsnorm)*
      D[EL[τ, ζ], ζ] + (vpnorm - 
        vsnorm)*(ζ/(1 + ζ^2))*EL[τ, ζ] + 
     D[EL[τ, ζ], τ] + Γl*
      EL[τ, ζ] == (1/3)*κL*gs*EL[τ, ζ]^2*
     ZRL*vpnorm*M[τ, ζ]*ES[τ, ζ]};
initial = {M[0, ζ] == 0, 
   ES[0, ζ] == 0*E^(-2 Log[2]*(ζ/2)^2), 
   EL[0, ζ] == f[ζ] };
bc = DirichletCondition[{M[τ, ζ] == 0, 
    ES[τ, ζ] == 0, 
    EL[τ, ζ] == 
     E^(-2 Log[2]*(ζ - vpnorm*τ)^2)}, ζ == 0 ];
{M, ES, EL} = 
  NDSolveValue[{pdes, initial, bc}, {M, ES, EL}, {τ, 0, 
    2.6*10^-10/tp}, {ζ, 0, .08/ZRL}];
framesTEQ = 
  Table[Plot[{EL[τ, ζ]*EL[τ, ζ]\[Conjugate], 
     ES[τ, ζ]*ES[τ, ζ]\[Conjugate]}, {ζ, 
     0, .08/ZRL}, PlotRange -> {0, 1}], {τ, 
    0, (2.6*10^-10)/tp, (2.6*10^-10)/(20*tp)}];
Manipulate[framesTEQ[[i]], {{i, 3, "time"}, 1, Length[framesTEQ], 2}, 
 SaveDefinitions -> True]