NDSolve to solve 2D space 1D time partial integro-differential equations

Question

I have a set of pdes in 2D space and 1D time, one of them is related to the integral in space by Reynolds function. The equations are similar to $$\frac{\partial \phi}{\partial t}=-\frac{\partial \phi_Z}{\partial x}\phi$$ $$ \frac{\partial\phi_Z}{\partial t}=-\int_0^1 \frac{\partial\phi}{\partial x} \frac{\partial \phi}{\partial y} dy $$ $$ Reynolds\equiv -\int_0^1 \frac{\partial\phi}{\partial x} \frac{\partial \phi}{\partial y} dy $$ To calculate the Reynolds stress (averaging in y direction), I used the method from Michael E2's answer on question 175080, which breaks the NDSolve[] apart. That method works well in 1D space, so I modified it like below:

db = 0.01; rho = 0.005; tmax = 0.8;
eqns = {D[ph[t, x, y], t] == -D[PZ[t, x, y], x]*D[ph[t, x, y], y],
         D[PZ[t, x, y], t] == Reynolds[ph[t, x, y], t, x, y]};
incs = {ph[0, x, y] == E^(-((x - 0.5)^2 + (y - 0.5)^2)/0.01), PZ[0, x, y] == 0};
bcs = {ph[t, 0, y] == ph[t, 1, y], ph[t, x, 0] == ph[t, x, 1], 
        PZ[t, 0, y] == PZ[t, 1, y], PZ[t, x, 0] == PZ[t, x, 1]};
sys = {eqns, incs, bcs};
Block[{Reynolds},
Reynolds[pp_, t_ /;t == 0, xn_?NumericQ, yn_?NumericQ] := 0;
  (Uncomment below, it works;)
  (Reynolds[pp_,t_,xn_?NumericQ,yn_?NumericQ]:=0;)
Reynolds[pp_, t_, xA_?MatrixQ, yA_?MatrixQ] := 0;
   (Some integral like below, I use 0 here, not affect the error)
   (ParallelTable[
    Module[{INTfunc = Interpolation[Transpose@{Transpose@{xA, yA}, pp}]}, 
      -NIntegrate[Derivative[1, 0][INTfunc][xx, yi]Derivative[0, 1][INTfunc][xx, yi]
                  , {yi, 0, 1}]], {xx, xA}, {yy, yA}];*)
Clear[foo];
  PrintTemporary@Dynamic@{foo, tmax, Clock[Infinity]};
  (broken down NDSolve call)
InternalInheritedBlock[{MapThread}, {state} = NDSolveProcessEquations[sys, {ph, PZ}, {x, 0, 1}, {y, 0, 1}, {t, 0, tmax}, 
     InterpolationOrder -> All, StepMonitor :> (foo = t), 
     Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> 
                  {"TensorProductGrid", "MinPoints" -> 10}}}];
   Unprotect[MapThread];
   MapThread[f_, data_, 1] /; ! FreeQ[f, Reynolds] := f @@ data;
   Protect[MapThread];
   NDSolveIterate[state, {0, tmax}]; sol = NDSolveProcessSolutions[state]]] // AbsoluteTiming

But the MMA has error like

NDSolve`Iterate::nlnum: The function value {0.,0....0.,<<238>>} is not a list of numbers with dimensions {288} at {t,ph[t,x,y],PZ[t,x,y]} = {8.9003*10^-308,{<<1>>},...}.

I believe the problem comes from the definition of

Reynolds[pp_, t_, xA_?MatrixQ, yA_?MatrixQ] := 0;

This function is not used by the NDSolve`Iterate[] ! Here I want to interpolate pp[xA,yA], so I can do the integral numerically. But seems the argument passed to Reynolds function is Numeric instead of Matrix/Array or Vector. VectorQ is tried:

Reynolds[pp_, t_, xA_?VectorQ, yA_?VectorQ] := 0;

But only NumericQ works.

Reynolds[pp_, t_, xA_?NumericQ, yA_?NumericQ] := 0;

Look like it is not like the 1D space case, where x is a Vector and passed to function. I can't do interpolation with only one numeric value.

Is there any solution for this? For example, force passing Matrix/Array to Reynolds function here? Or other suggestion? Thank you very much!