Expanding 2D fluxes to 3D

[gmacgilchrist]

At present, 2D surface fluxes are expanded to 3 dimensions so that their divergence can be computed in a manner consistent with 3D fluxes.

The present way of choosing to do this is to determine whether lev_outer is present in the dataset. See here. If it is, the variables are assumed to be 3D and will not be expanded, if it's not, they will be expanded.

However, it is possible for a Dataset to have lev_outer and still the fluxes as 2D - a Dataset only needs one variable to have the depth dimension for it to be present. For example, rsdo is a variable that is commonly posted with the surface fluxes. A "TODO" message related to this is present here.

I need to give some thought to what is a more foolproof way to do this. The code loops through the variable anyway - can we instead just check at each variable whether it is 2 or 3 dimensional already, and expand the 2D ones?

[hdrake]

Another major issue with this is that the existing approach really only works for vertical coordinates for which the zeroth index is at the sea surface:

def expand_surface_to_3d(surfaceflux, z):
    """Expand 2D surface array to 3D array with zeros below surface"""
    return surfaceflux.expand_dims({"lev_outer": z}).where(z == z[0], 0)

Expanding 2D surface fluxes to outcropping rho2 or theta levels, for example, sounds tricky...

[gmacgilchrist]

I'm not completely sure whether this is actually an issue, since this expansion is quasi-redundant - you expand in the vertical dimension only so that you can collapse that dimension again in the next step. So the "realism" of your uppermost level is not important. This is technically correct - the "tendency" due to the surface flux can be thought of as the divergence of that flux across an infinitesimally small surface layer (in Stephen's nomenclature he uses a Dirac delta here). However, I could envision it causing some numerical innacuracy if that layer's thickness was ~zero. I've never been a fan of this step (which I implemented in my original version of these calculations) since it is purely to align the forms WMT calculations from fluxes. I've not yet figured a better or more robust way to do this other than specifying somewhere that this is a surface flux and therefore skipping this step and doing the tendency calculation in a different way. This could be the preferred approach.

[hdrake]

I agree that this is not a fundamental limitation, but I'm pretty sure the .where(z == z[0], 0) part of the current implementation will be wrong for, e.g., rho2 coordinates where the surface fluxes are not in the zeroth index and thus will get masked to zero...

I've not yet figured a better or more robust way to do this other than specifying somewhere that this is a surface flux and therefore skipping this step and doing the tendency calculation in a different way. This could be the preferred approach.

Yes, I agree this is the way forward.

[hdrake]

I need to give some thought to what is a more foolproof way to do this. The code loops through the variable anyway - can we instead just check at each variable whether it is 2 or 3 dimensional already, and expand the 2D ones?

This is what I do in https://github.com/NOAA-GFDL/xwmt/pull/40/commits/da06aea194d162d47538b8304c0c6d732a55a48e, see https://github.com/hdrake/xwmt/blob/da06aea194d162d47538b8304c0c6d732a55a48e/xwmt/wmt.py#L124-L138

Technically this would fail if there was a 2D file that had a single-valued vertical coordinate dimension, but I haven't come across that issue. It would be straightforward to modify the bool in the if-else block anyway.

[hdrake]

At present, 2D surface fluxes are expanded to 3 dimensions so that their divergence can be computed in a manner consistent with 3D fluxes.

Following up on this, I don't think this is actually necessary. The whole reason for expanding the 2D fluxes to 3D is, as far as I can tell, for this step: https://github.com/NOAA-GFDL/xwmt/blob/07b2e6fc95fd5bb715b461656509a8497c21b87a/xwmt/compute.py#L152-L160

However, I am pretty sure Jlam is exactly equal to hlamdot, since $h \dot{\lambda} \equiv \int{z{k}}^{z{k}+h} -\partial{z} \mathcal{F}^{(z)} \text{d}z = -(\mathcal{F{\text{surf}}} - 0) = - \mathcal{F{\text{surf}}} \equiv \mathcal{J}_{\lambda}$ by the divergence theorem.

This means that we should be able to just get rid of this function altogether. However, this still does not address the fact that the surface flux (convergence) needs to be added to the correct vertical coordinate level (trivial for z coordinates, but not obvious otherwise). I think this means we need to find the index of the outcropping layer.

[hdrake]

I think this means we need to find the index of the outcropping layer.

A problem we'll have to deal with at a later date is how to find even find this for a non-monotonic coordinate. I guess just the vertical level whose value is closest to the local value of the 2D surface field (e.g. tos, sos)?

[hdrake]

I think this means we need to find the index of the outcropping layer.

We also would need this for the ocean reference value in the surface mass flux formula $Q(\lambda_{m} - \lambda)$, but my sense is that these terms are almost always negligible anyways.

[hdrake]

PR https://github.com/NOAA-GFDL/xwmt/pull/40 addresses this issue by checking the dimensions of the tendency field: https://github.com/hdrake/xwmt/blob/50e470b42b27a2c5bc6215f7b90c2c750a7bb320/xwmt/wmt.py#L187-L201

For 2D fields, we use the following methods in the WaterMass class to expand the 2D fluxes to 3D (outer levels) and to pick out the outcropping layer: https://github.com/hdrake/xwmt/blob/50e470b42b27a2c5bc6215f7b90c2c750a7bb320/xwmt/wm.py#L154-L232