On a Riemannian manifold $(M,g)$ the gradient $\text{grad} f$ of a real-valued function $f$ is the unique vector field satisfying $g(\text{grad} f, X) = \Bbb d f (X)$ for every tangent field $X$. Working in coordinates and letting $X = \partial_j$ we obtain
$$g\left( \sum_i (\text{grad} f)^i \partial_i, \partial_j \right) = (\Bbb d f) (\partial_j)$$
which means that
$$\sum_i (\text{grad} f)^i \ \underbrace{ g( \partial_i, \partial_j) } _{g_{ij}} = \frac {\partial f} {\partial x_j} ,$$
whence, if we multiply the equality by $(g^{ij})$, the matrix inverse to $(g_{ij})$, we get
$$(\text{grad} f)^i = \sum _j g^{ij} \frac {\partial f} {\partial x_j} ,$$
so that finally (and after an interchange of indices, for purely cosmetic reasons)
$$\text{grad} f = \sum _{i,j} g^{ij} \frac {\partial f} {\partial x_i} \partial _j.$$
Concerning the Hessian, remember that the metric $g$ gives birth to the Levi-Civita connection $\nabla$: according to Koszul's formula, $\nabla_X Y$ is the unique tangent field satisfying
$$g(\nabla _X Y, Z) = \frac 1 2 g \Big( X g(Y,Z) + Y g(Z,X) - Z g(X,Y) + g([X,Y],Z) - g([Y,Z],X) - g([Z,X],Y) \Big) .$$
Next, for $T$ a $p$-covariant tensor field, one defines
$$(\nabla T) (X_1, \dots, X_p, X) = (\nabla_X T) (X_1, \dots, X_p) = X \cdot \Big( T(X_1, \dots, X_p) \Big) - \\
T(\nabla_X X_1, X_2, \dots, X_p) - T(X_1, \nabla_X X_2, \dots, X_p) - \dots - T(X_1, X_2, \dots, X_{p-1}, \nabla_X X_p) .$$
This finally allows us to define the Hessian of $f$ as $\nabla^2 f = \nabla (\nabla f) = \nabla \Bbb d f$.
Explicitly,
$$(\nabla^2 f) (Y,X) = (\nabla \Bbb d f) (Y,X) = (\nabla_X \Bbb d f) (Y) = X \cdot (\Bbb d f (Y)) - \Bbb d f (\nabla_X Y) .$$
In coordinates, taking $Y = \partial_i$ and $X = \partial_j$ and remembering that $\nabla_{\partial_i} \partial_j = \sum_k \Gamma_{ij}^k \partial_k$,
$$(\nabla^2 f) (\partial_i, \partial_j) = \partial_j (\partial_i f) - \Bbb d f (\sum_k \Gamma_{ji}^k \partial_k) = \frac {\partial^2 f} {\partial x_j \partial x_i} - \sum_k \Gamma_{ij} ^k \frac {\partial f} {\partial x_k} .$$
Remembering that the Euclidean metric has $\Gamma_{ij}^k = 0$ shows that in the Euclidean case we get what we already know from calculus.
An interesting fact is that the trace of the Hessian (after having raised an index with the help of $(g^{ij})$) is precisely the Beltrami Laplacian $\Delta f$.
The equation $\nabla h = X$ is equivalent to $dh = X^\flat$, where $X^\flat$ is the 1-form dual to $X$ (this is basically the definition of $\nabla$ on a Riemannian manifold). Thus, your question transforms to the following: for which 1-forms $\omega$ there exists a function $h$ with $dh=\omega$?
By the properties of differentiation of forms, $d\omega=d^2h=0$, hence $d\omega=0$ is a necessary condition for the existence of such $h$. For a form $\omega = \sum \omega_i dx_i$ this condition means
$$\frac{\partial\omega_i}{\partial x_j}=\frac{\partial\omega_j}{\partial x_i}$$
So, there are 1-forms with $\omega = dh$ (called exact forms), which we are interested in, and there are forms with $d\omega=0$ (called closed forms). Any exact form is closed, and we are interested in whether the converse is true. This is precisely what the first de Rham cohomology $H^1_{dR}(M)$ measures: it is defined as the space of closed forms modulo the space of exact forms.
For the case of a 2-sphere $M=S^2$, it is known that $H^1_{dR}(S^2)=0$, meaning that closed and exact forms coincide, so for any closed 1-form $\omega$ there is a function $h$ with $dh=\omega$.
For a 2-torus $M=T^2$ we have $H^1_{dR}(T^2)=\mathbb R^2$, meaning that any closed 1-form decomposes as an exact form plus a linear combination of 2 linearly independent non-exact forms. In fact, using a coordinate system $(\theta, \phi)$ with $\theta, \phi \in [0; 2\pi)$ we can construct a basis of the space of non-exact forms as $\{d\theta, d\phi\}$.
Best Answer
Point (1) : Since $g$ is non degenerate, it induces two musical isomorphisms :
In local coordinates $\{x^i\}$ we have $g = g_{ij}\mathrm d x^i \otimes \mathrm d x^j$ and the musicalities are given by :
where $[g^{ij}] = [g_{ij}]^{-1}$. The gradient $\nabla f$ is then : $$ \nabla f = \sharp(\mathrm d f) = \sharp(\partial_i f \mathrm d x^i) = (\partial_i f) g^{ij}\partial_j $$ So here you have the gradient of $f$ in local coordinates : $\nabla f = (\partial_i f) g^{ij}\partial_j$.
Point (2) : This problem is an example of a more generic problem : given $\alpha\in T_p^*M$ and $v\in T_p M$ such that $\alpha(v)>0$, show there exists $g$ on $M$ such that $\sharp(\alpha) = v$. Since a metric defined at a point can be extended to a global metric, all we need to do is show that there exists a scalar product $g_p : T_pM\otimes T_pM\to \mathbb R$ such that $g_p(v,\cdot) = \alpha$. Now, $\alpha(v)>0$ implies $v\ne 0$ and $\alpha \ne 0$. So we can choose coordinates $\{x^i\}$ around $p$ such that $v = \partial_1$. Now, in those coordinates, we are looking for a symmetric positive definite matrix $[g_{ij}]$ such that $$ g_{1j} = \alpha_j $$ Now, $\alpha(v)>0$ implies $\alpha_1>0$. So we are looking for a symmetric positive matrix $[g_{ij}]$ such that $g_{11}>0$ and such $g_{1j} = \alpha_j$. Such a matrix exists.