I think the problem you’re having is that you’re solving for the optimal $p$ as a standard multivariable optimisation problem. There is a noticeable difference between functional optimisation compared with multivariate optimisation, but this difference is often ignored in economics in favour of heuristic methods.
In the problem you stated, we need to find the optimal $\textit{function}$ $p:(0,\infty)\to \mathbb{R}$ along with an optimal value of a $\textit{variable}$ $K$. Your first question likely stems from this difference.
The Lagrangian is $$\mathcal{L}=\mathbb{E}[S(p(\theta),\theta)-cD(p(\theta),\theta)+\lambda(\theta)(K-D(p(\theta),\theta))]-rK$$ The reason we have $\lambda(\theta)$ inside the expectation is that we have the constraint $D(p(\theta),\theta)\leq K$ holding for all $\theta$. If $\theta$ were drawn from some finite set $\Theta$, we would have the Lagrangian as $$\mathcal{L}=\sum_{\theta'\in \Theta}[S(p(\theta'),\theta')-cD(p(\theta'),\theta')+\lambda(\theta')(K-D(p(\theta'),\theta'))]\mathbb{P}(\theta=\theta')-rK$$ which may be more familiar. The integral simply generalises this to the case when $\Theta$ is a continuum. It is also intuitive that the constraint should be within the expectation operator as it is a stochastic constraint. We know that for any $\theta$ realisation the constraint will bind, but since we don't know what value $\theta$ will take we consider an expectation over the possible realisations of the constraint.
Now we can move on to the optimisation problem.
Consider the variable optimisation problem, i.e., the choice of $K$. Taking the FOC of $\mathcal{L}$ w.r.t. $K$, we find (assuming regularity conditions on $S$ and $D$ allowing us to differentiate under the integral) that $$\begin{align*} \frac{\partial \mathcal{L}}{\partial K} &=\frac{\partial }{\partial K}[\mathbb{E}[S(p(\theta),\theta)-cD(p(\theta),\theta)+\lambda(\theta)(K-D(p(\theta),\theta))]-rK] \\
&= \mathbb{E}\left[\frac{\partial }{\partial K}(S(p(\theta),\theta)-cD(p(\theta),\theta)+\lambda(\theta)(K-D(p(\theta),\theta)))\right]-r \\
&= \mathbb{E}[\lambda(\theta)]-r \end{align*}$$
Second, for the functional optimisation problem, we should really take a variational approach. This answer by tdm gives a good overview of the methodology: FOCs for the Dixit-Stiglitz aggregator: derivative of an integral w.r.t. integrand at one point . Suppose we are given an optimal function choice $p^*$. Consider perturbing this optimal choice by any other $C^1$ function $g$ to $p^*+\varepsilon g$. Define $$\mathcal{L}(f)=\mathbb{E}[S(f(\theta),\theta)-cD(f(\theta),\theta)+\lambda(\theta)(K-D(f(\theta),\theta))]-rK$$ Since $p^*$ is optimal, we must have that $$\frac{\partial}{\partial \varepsilon} \mathcal{L}(p^*+\varepsilon g)\Bigg|_{\varepsilon=0}=0$$ Evaluation of this gives $$\begin{align*} 0 &= \frac{\partial}{\partial \varepsilon} \mathcal{L}(p^*+\varepsilon g)\Bigg|_{\varepsilon=0} \\
&=\frac{\partial}{\partial \varepsilon} [\mathbb{E}[S(p^*(\theta)+\varepsilon g(\theta),\theta)-cD(p^*(\theta)+\varepsilon g(\theta),\theta)+\lambda(\theta)(K-D(p^*(\theta)+\varepsilon g(\theta),\theta))]-rK]\Bigg|_{\varepsilon=0} \\
&=\frac{\partial}{\partial \varepsilon} [\mathbb{E}[S(p^*(\theta)+\varepsilon g(\theta),\theta)-cD(p^*(\theta)+\varepsilon g(\theta),\theta)+\lambda(\theta)(K-D(p^*(\theta)+\varepsilon g(\theta),\theta))]]\Bigg|_{\varepsilon=0}\\
&=\mathbb{E}\left[\frac{\partial}{\partial \varepsilon}(S(p^*(\theta)+\varepsilon g(\theta),\theta)-cD(p^*(\theta)+\varepsilon g(\theta),\theta)+\lambda(\theta)(K-D(p^*(\theta)+\varepsilon g(\theta),\theta)))\right]\Bigg|_{\varepsilon=0} \\
&=\mathbb{E}\left[\frac{\partial S}{\partial p}(p^*(\theta)+\varepsilon g(\theta),\theta)g(\theta)-(c+\lambda(\theta))\frac{\partial D}{\partial p}(p^*(\theta)+\varepsilon g(\theta),\theta)g(\theta) \right]\Bigg|_{\varepsilon=0} \\
&= \mathbb{E}\left[\left(\frac{\partial S}{\partial p}(p^*(\theta)+\varepsilon g(\theta),\theta)-(c+\lambda(\theta))\frac{\partial D}{\partial p}(p^*(\theta)+\varepsilon g(\theta),\theta)\right)g(\theta) \right]\Bigg|_{\varepsilon=0} \\
&= \mathbb{E}\left[\left(\frac{\partial S}{\partial p}(p^*(\theta),\theta)-(c+\lambda(\theta))\frac{\partial D}{\partial p}(p^*(\theta),\theta)\right)g(\theta) \right] \end{align*}$$
The above is equal to zero for any $g\in C^1$. Hence, assuming $S$ and $D$ are continuous (and that the optimal $p^*$ is continuous), by the Fundamental Lemma of the Calculus of Variations, we must have $$\frac{\partial S}{\partial p}(p^*(\theta),\theta)-(c+\lambda(\theta))\frac{\partial D}{\partial p}(p^*(\theta),\theta)=0$$
What I meant in the comments about this being pointwise is that we can get the exact same result if we ignore the the expectation operator in $\mathcal{L}$ and take the derivative w.r.t. $p(\theta)$ for a fixed $\theta$. Generally this is acceptable if derivatives of $p(\theta)$ are not included in $\mathcal{L}$. If the first derivative is included, we can use the Euler-Lagrange equation to get a first-order condition.
