Distribution system operators (DSOs) world-wide foresee a rapid roll-out of distributed energy resources. From the system perspective, their reliable and cost effective integration requires accounting for their physical properties in operating tools used by the DSO. This paper describes a decomposable approach to leverage the dispatch flexibility of thermostatically controlled loads (TCLs) for operating distribution systems with a high penetration level of photovoltaic resources. Each TCL ensemble is modeled using the Markov decision process (MDP). The MDP model is then integrated with a chance constrained optimal power flow that accounts for the uncertainty of photovoltaic resources. Since the integrated optimization model cannot be solved efficiently by existing dynamic programming methods or off-the-shelf solvers, this paper proposes an iterative spatiotemporal dual decomposition algorithm (ST-D2). We demonstrate the merits of the proposed integrated optimization and ST-D2 algorithm on the IEEE 33-bus test system.