Theta functions on an elliptic curve and Serre duality
Let $X$ be an elliptic curve $\mathbb{C}/\Lambda$, where $\Lambda$ is a lattice of real rank $2$ in $\mathbb{C}$. A theta function is a holomorphic section of a line bundle $L$ on $X$ whose transition from $U$ to $U + \ell$ is given by$$f(z + \ell) = e^{a_\ell z + b_\ell} f(z)$$for $f \in \Lambda$, where $a_\ell, b_\ell \in \mathbb{C}$ and $f$ is an entire function on $\mathbb{C}$. This theta function can be regarded as an element of $\text{H}^0(X, L)$ when $\text{H}^0(X, L)$ is the kernel and $\text{H}^1(X, L)$ is the cokernel of the map$$C_{0, 0}^\infty(X, L) \overset{\overline{\partial}}{\to} C_{0, 1}^\infty(X, L),$$where $C_{0, q}^\infty(X, L)$ is the space of all $L$-valued $C^\infty(0, q)$-forms on $X$ for $q = 0, 1$. Thus $\text{H}^1(X, L^{-1})$ is the cokernel of$$C_{0,0}^\infty(X, L^{-1}) \overset{\overline{\partial}}{\to} C_{0, 1}^\infty(X, L^{-1}).$$We introduce a smooth metric $H(z) > 0$ for $L$. The compatibility condition for $H(z)$ is$$H(z + \ell) \left|e^{a_\ell z + b_\ell}\right|^2 = H(z)$$for $\ell \in \Lambda$ so that the pointwise norm square$$H(z + \ell)|f(z + \ell)|^2 = H(z)|f(z)|^2$$is a well-defined function on $X$. One way to explicitly write down such a metric $H(z)$ is to choose $H(z) = e^{-\gamma|z|^2}$ so that the curvature$${{\sqrt{-1}}\over{2\pi}}\partial\overline{\partial}(-\log H(z)) = \gamma {{\sqrt{-1}}\over{2\pi}} \text{d}z \wedge \text{d}\overline{z}$$of the metric $H(z)$ of $L$ is positive when $\gamma$ is a positive constant. The compatibility condition$$H(z + \ell) \left|e^{a_\ell z + b_\ell}\right|^2 = H(z)$$corresponds to the conditions $a\overline{\ell} = a_\ell$ and ${1\over2}|\ell|^2 = b_\ell$. With the use of the metric $H(z)^{-1}$ for $L^{-1}$, the cokernel $\text{H}^1(X, L^{-1})$ of$$C_{0,0}^\infty(X, L^{-1}) \overset{\overline{\partial}}{\to} C_{0,1}^\infty(X, L^{-1})$$can be identified as the kernel of the $\overline{\partial}^*$ of the map $\overline{\partial}$ (with respect to the metric $H(z)^{-1}$ of $L^{-1}$). An element $\text{H}^1(X, L^{-1})$ can now be described by a $C^\infty$ $(0, 1)$-form $g(z)\,\text{d}\overline{z}$ on $\mathbb{C}$ such that:
(i) $g(z + \ell) = e^{-a_\ell z - b_\ell}g(z)$ for $z \in \mathbb{C}$ and $\ell \in L$ and
(ii) $e^{\gamma z \overline{z}} g(z)$ is antiholomorphic, which means that $\overline{\partial}^* g(z) = -H(z)\partial_z(H(z)^{-1}g(z))$ (with $H(z) = e^{-\gamma z\overline{z}}$) is identically zero.
In other words, $g(z)\,\text{d}\overline{z}$ defines a $L^{-1}$-valued $(0, 1)$-form on $X$ so that (i) states the fact that $g(z)\,\text{d}\overline{z}$ is $L^{-1}$-valued and (ii) states the fact that $g(z)\,\overline{z}$ is harmonic with respect to the metric $H(z)^{-1}$ of $L^{-1}$. Note that harmonic means both $\overline{\partial}$-closed and $\overline{\partial}^*$-closed, but $\overline{\partial}$-closedness is automatic for a $(0, 1)$-form on the compact Riemann surface $X$.
The Serre duality comes from the pairing$$\text{H}^0(X, L \otimes K_X) \times \text{H}^1(X, L^{-1}) \ni (f, g\,\text{d}\overline{z}) \mapsto \int_X f(z)\,\text{d}z \wedge g(z)\,\text{d}\overline{z}$$because $K_X$ is trivial so that we can identify $f \in \text{H}^0(X, L)$ with $f\,dz \in \text{H}^0(X, L \otimes K_X)$.
Here is a 'low-brow' approach. One type of the result you are talking about has been written up implicitly in Lang's book Introduction to Arakelov theory. The case for cohomology of the elliptic curves was investigated by Coleman (Chapter 2, section 4). The construction of Green functions using theta functions can be found in the next section. It is not difficult to engineer a "dipole Green function" using the given Green functions.
The exact Green function would depend on your metric; if you use Arakelov metric and identify the elliptic curve as $\mathbb{C}/\langle 1, \tau \rangle$, then the function can be expressed using Dedekind $\eta$ functions:
$$ g(z_1,z_2)=\log ||\theta||(\tau, z_1-z_2+\frac{1+\tau}{2})-\log||\eta||(\tau) $$ where $||\theta||(a+ib, x+iy)=b^{1/4}e^{-\pi y^2/b}|\theta(a+ib, x+iy)|$.
Now using this, by taking gradients of harmonic functions, you can "go reverse" and construct meromorphic forms of desired property. For example $x^2-y^2$ gives rise to the holomorphic 1-form $2zdz$. Classically this can then be used to "prove" one side of the Riemann-Roch inequality. But this may be totally overkill for what you are asking about.