We study data-driven prediction of coarse-grained dynamics in multiscale PDE systems. Adopting a closure-free operator-learning viewpoint, we apply a linear coarse-graining map and learn a surrogate evolution operator for the resolved field directly from filtered high-fidelity trajectories. Motivated by the Mori--Zwanzig formalism, we propose a spatiotemporal neural operator mapping a resolved history slab on to a resolved future slab. Spatial mixing uses Fourier convolution, while temporal mixing uses a causal kernel operator with position-attention weights on time lags. This causal temporal operator encodes finite-memory effects in the resolved dynamics while preserving the directionality of the history-to-future map. To improve rollout robustness and suppress nonconservative artifacts, we embed a flux-form inductive bias by parameterizing the windowed update in explicit divergence form. We also provide a data-driven guideline for selecting the memory length via the decorrelation time of a closure-injection diagnostic computed from filtered trajectories. We validate on coarse-grained the viscous Burgers' equation, the Kuramoto--Sivashinsky, and two-dimensional turbulent flows, obtaining stable autoregressive rollouts with improved long-horizon accuracy and statistical fidelity.