Generalizing WMT (both 3D and 2D) with restructured classes; no redundant methods/functions

[hdrake]

This PR picks up from https://github.com/NOAA-GFDL/xwmt/pull/35 and aims to simultaneously integrate swmt.py into wmt.py while also addressing various structural limitations discussed in the following issues:

• https://github.com/NOAA-GFDL/xwmt/issues/3: Robust expansion of 2D surface fluxes to 3D within one generalized WaterMassTransformation class method
• https://github.com/NOAA-GFDL/xwmt/issues/6: Uses coordinate information from grid to avoid this ambiguity
• https://github.com/NOAA-GFDL/xwmt/issues/7: separates "thermodynamic" calculations, which are now methods of the parent xwmt.wm.WaterMass class from "wmt" calculations, which are now methods of the child xwmt.wmt.WaterMassTransformations class. Also removes redundant/repeated function/methods definitions, making the package more robust and easier to develop moving forward.
• https://github.com/NOAA-GFDL/xwmt/issues/8: addresses most of these points by adding thermodynamic variables as new variables within the core dataset (WaterMass.ds)
• https://github.com/NOAA-GFDL/xwmt/issues/10: Adds a simple test suite that compares against analytical solutions
• https://github.com/NOAA-GFDL/xwmt/issues/34: integrates surface WMT with fully 3D WMT in consistent way that supports all monotonic generalized vertical coordinates

The commits containing the most important major changes are:

• https://github.com/NOAA-GFDL/xwmt/pull/40/commits/5c6b0806d5007d5f50d06fffc5a0866876684fd5
• https://github.com/NOAA-GFDL/xwmt/pull/40/commits/90ebf498aec39a5b9848ff01637831775c893eba
• https://github.com/NOAA-GFDL/xwmt/pull/40/commits/50e470b42b27a2c5bc6215f7b90c2c750a7bb320
• https://github.com/NOAA-GFDL/xwmt/pull/40/commits/fc602ad163df608928dd7dd531c713f81e6f4805
• https://github.com/NOAA-GFDL/xwmt/pull/40/commits/cc24536e1efd446354aced1ca871bee9f5e30839

Remaining to do list:

• [x] Support calculations with only 2D fields (fluxes and scalars)
• [x] Separate out components of surface salt fluxes (as requested by https://github.com/NOAA-GFDL/xwmt/pull/40#issuecomment-1546430994)
• [x] Support for 3D fluxes on cell interfaces, as in rsdo (3D shortwave penetration)

[hdrake]

@gmacgilchrist, do you have any suggestions for how to replace the xwmt.compute.expand_surface_to_3d function with a more robust approach (that would work for any vertical coordinate)?

[hdrake]

Commit https://github.com/NOAA-GFDL/xwmt/pull/40/commits/32c6b44ff6869db672829c1d49d21284d39dab46 addresses https://github.com/NOAA-GFDL/xwmt/issues/3#issuecomment-1533760980

[hdrake]

@gmacgilchrist and @jetesdal, would be great if you could review when you get a chance. Not urgent, but before our O-div presentation would be great!

[jetesdal]

While testing the new code, I noticed that It is required to name salinity as 'salt':

ValueError: ds must include salinity variable defined by kwarg t_name (default: salt).

However, temperature variable is expected to be thetao. I believe it would be more consistent to use either 'so' and 'thetao' or 'salt' and 'temp' together.

I like the use of the naming convention defined in xwmt/conventions/MOM6.yaml, but I think it would be more appropriate to use 'thetao' and 'so' as variable names. Does that make sense?

[jetesdal]

About the comment in the example notebook (basic_xwmt_tutorial.ipynb) regarding the discrepancy between WMT from 3D boundary forcing and 2D surface fluxes:

There is still a small discrepancy between the 3D surface forcings and the 2D total surface fluxes. I am not totally sure where these errors come from.

@hdrake: I am curious if you have considered shortwave penetration in the calculations. If not, this could potentially be the cause of the mismatch between the 3D and 2D results.

[hdrake]

I like the use of the naming convention defined in xwmt/conventions/MOM6.yaml, but I think it would be more appropriate to use 'thetao' and 'so' as variable names. Does that make sense?

I'm happy with changing the default, but note that the names of the T, S, and thickness variables, which are all that's needed to define a WaterMass, can be changed with keyword arguments:

https://github.com/hdrake/xwmt/blob/eef2c5588113cf06fc152413673ada786b2231f4/xwmt/wm.py#L30-L35

[hdrake]

@hdrake: I am curious if you have considered shortwave penetration in the calculations. If not, this could potentially be the cause of the mismatch between the 3D and 2D results.

That's a great point. I'm pretty sure it is included in the 3D and I guess in the shortwave 2D component, which in the 2D calculation is assumed to (inaccurately) go in at only the surface.

[jetesdal]

I like the use of the naming convention defined in xwmt/conventions/MOM6.yaml, but I think it would be more appropriate to use 'thetao' and 'so' as variable names. Does that make sense?

I'm happy with changing the default, but note that the names of the T, S, and thickness variables, which are all that's needed to define a WaterMass, can be changed with keyword arguments:

https://github.com/hdrake/xwmt/blob/eef2c5588113cf06fc152413673ada786b2231f4/xwmt/wm.py#L30-L35

Ah OK, that's good to know. So t_name,s_name and h_name will overwrite what is listed in budgets_dict? I still think that so for salinity is more the standard in MOM6 (at lease in case of the production runs and CMIP). So this would just be something that will be changed in conventions/MOM6, correct?

[jetesdal]

I am trying to understand why the answers differ in total surface WMT using published version (0.0.1) vs your version (0.03). The differences occur in the lighter densities (sigma2 < 15). This is the same data that is used the example notebook (basic_xwmt_tutorial.ipynb: sim = "wmt_incsurffluxes.natv_rho2_zstr.monthly_daily_hourly.13months") I am confused why this has not been picked up by the test, or are we not running the standard unit test on this Pull Request?

[hdrake]

Ah OK, that's good to know. So t_name,s_name and h_name will overwrite what is listed in budgets_dict?

Yes, actually right now it looks like the "lambda" entry in the budgets_dict is not even read in at any point. The reason I allowed t_name and s_name be provided separately is that I imagine there may be some use cases for the WaterMass class without needed wanted to specify the whole budgets_dict needed for initializing a WaterMassTranformation instance.

I think the more intuitive behavior would instead be that budgets_dict would override the t_name and s_name kwargs, however.

So this would just be something that will be changed in conventions/MOM6, correct?

Following on the above suggestion, we would change it both in the kwargs for WaterMass and in the conventions/MOM6.

[hdrake]

I am trying to understand why the answers differ in total surface WMT using published version (0.0.1) vs your version (0.03). The differences occur in the lighter densities (sigma2 < 15).

Are you comparing zstr for both? The notebook basic_xwmt_tutorial.ipynb was defaulting to rho2 diagnostics to showcase the new, more general, approach.

I am confused why this has not been picked up by the test, or are we not running the standard unit test on this Pull Request?

I was not running the xwmt tests, because I couldn't find how to download the test datasets at the time (I now know they are in https://github.com/NOAA-GFDL/xwmt/blob/main/.github/workflows/ci.yml). I'm not sure why the CI tests are not running on the PR. I am new to CI and testing for python packages (only experience w/ Julia).

[hdrake]

Are you comparing zstr for both? The notebook basic_xwmt_tutorial.ipynb was defaulting to rho2 diagnostics to showcase the new, more general, approach.

@jetesdal, I re-ran the notebook, only changing gridname = 'rho2' --> gridname = 'zstr' at the top, and got something that at least qualitatively looks more much similar to the v0.0.1 results, although it seems there are some small discrepancies still.

[hdrake]

This is the same data that is used the example notebook (basic_xwmt_tutorial.ipynb: sim = "wmt_incsurffluxes.natv_rho2_zstr.monthly_daily_hourly.13months")

Also, you probably caught this, but I just wanted to point out that am leaving out the first month of the dataset with ds.sel(time=slice('1900-02-01 00', '1901-02-01 00')) before my mean('time') computations.

[jetesdal]

Yes, using gridname = 'rho2' instead of gridname = 'zstr' was the main factor for the earlier mismatch. I am just comparing surface WMT, so I am surprised to see that the answers are different betweenrho2 and zstr. In v0.0.1, surface WMT is derived solely from 2d surface fields so it should not matter whether gridname = 'rho2' or gridname = 'zstr', but it does matter for v0.0.3. Is this because you are using 3d thetao/so instead of tos/sos? Is there any control on using the density fields in the dataset (ds['rhopot0']) versus deriving it from T&S? As you said, there are still some minor differences. My suspicion is that this is because of differences in the density calculation (e.g., online versus offline). Could that be the case?

[jetesdal]

I was not running the xwmt tests, because I couldn't find how to download the test datasets at the time (I now know they are in https://github.com/NOAA-GFDL/xwmt/blob/main/.github/workflows/ci.yml). I'm not sure why the CI tests are not running on the PR. I am new to CI and testing for python packages (only experience w/ Julia).

I am not too familiar with how it works. @jkrasting, it appears that the CI tests are not running as intended for the recent Pull Request. I think you set up the testing environment last year, I am hoping you could provide some guidance on how to fix this issue. Do you have any idea why the tests might not be working as expected?

[hdrake]

In v0.0.1, surface WMT is derived solely from 2d surface fields so it should not matter whether gridname = 'rho2' or gridname = 'zstr', but it does matter for v0.0.3. Is this because you are using 3d thetao/so instead of tos/sos?

Yes, I think you're right that this is the problem. I got confused by the way the mapping to 3D worked (see my comments here) and had not realized you use the tos and sos tendcodes to correctly locate the surface fluxes in lambda space. That's certainly more accurate than using the 3D value in the layer that outcrops, which is what I am doing here. I can work on changing that. I might still keep the outcropping version as an option since it would be a useful approximation if tos and sos are missing.

I think there is one minor inconsistency with this approach, which is that you will have zero mass flux, by definition, since you use the same quantities for lm and l` in

return Qm * (lm - l)

[hdrake]

As you said, there are still some minor differences. My suspicion is that this is because of differences in the density calculation (e.g., online versus offline). Could that be the case?

I think probably it's just that even with zstr, the outcropping (top) layer is not going to be the same value as sos and tos. Is that right? I'm not actually sure how the variables are defined. Are 2D surface scalars just the finite-volume average value of the top-most natv layer, or are they somehow interpolated to the surface interface itself?

[hdrake]

In v0.0.1, surface WMT is derived solely from 2d surface fields so it should not matter whether gridname = 'rho2' or gridname = 'zstr', but it does matter for v0.0.3.

@jetesdal, upon further thought I'm a bit confused by this. We know that that the boundary term from the 3D watermass tranformaton calculation in sigma2 space is sensitivity to the diagnostic coordinate. The zstr output may be subject to offline regridding errors that could affect both the 3D and 2D calculations, while the rho2 output has no offline regridding errors by construction. If 2D surface fields are independent of the diagnostic grid (as in v0.0.1), but 3D tendency terms are not, then it is not clear to me how the two can be compared or consistently interpreted within the same framework. This makes me question the validity calculations using zstr output, which I am worried is negatively impacted by offline regridding errors (including for the 2D surface fields).

It actually gives me confidence that the v0.0.3 3D and 2D boundary transformations are fairly self-consistent, regardless of the coordinate choice, despite the fact that the results depend on the chosen vertical coordinate. It seems natural that the rho2 diagnostics would be most reliable for transformations in sigma2/rho2 space.

@jetesdal, did you ever compare estimates of the boundary flux terms from your 2D vs. 3D computations (I guess this would be xwmt v0.0.2)?

---

One caveat here is that, because the rho2 resolution is very low for the lighter water masses shown here, it is possible that zstr is more accurately captured the rho2 structure of the watermass transformations. But I would consider this user error (using discretization bins that do not match the underlying data–see the idealized example notebooks), not a flaw in the approach.

In any case, we should verify that the results for the different diagnostic coordinates converge when we reduce the offline averaging timestep (from monthly to daily to hourly=timestep).

[jkrasting]

@jetesdal - the CI is not running since the target branch of this pull request is dev_wmt3d. The CI is configured to run on pull requests going to main.

[jetesdal]

As you said, there are still some minor differences. My suspicion is that this is because of differences in the density calculation (e.g., online versus offline). Could that be the case?

I think probably it's just that even with zstr, the outcropping (top) layer is not going to be the same value as sos and tos. Is that right? I'm not actually sure how the variables are defined. Are 2D surface scalars just the finite-volume average value of the top-most natv layer, or are they somehow interpolated to the surface interface itself?

I am pretty sure that tos and sos are equal to thetao and so of the topmost z-layer which should be equivalent to the topmost native layer. Maybe @jkrasting or @StephenGriffies can confirm?

[StephenGriffies]

I just discussed this point with @adcroft and @Hallberg-NOAA , and we confirmed the following.

tos and sos are diagnosed from the ocean state passed to the atmosphere. These fields are the values of thetao and so averaged over the top 1m of ocean. This 1m is typically within the top model grid cell (which is 2m of z*), in which case sos and tos are the same as the value of so and thetao of k=1.

Make sense?

[jkrasting]

I'm still confused about this @StephenGriffies. In reading your email, it sounds like tos and thetao.isel(z_l=1) are close, but not identical. Am I interpreting that correctly?

[StephenGriffies]

They are the same so long as the top model layer is 1m or thicker.

Are you seeing a difference in the diagnostics?

[jkrasting]

Wouldn't they be the same if the top layer is exactly 1 m? In the case where the top layer is thicker than 1 m, wouldn't the additional volume of water alter the average for the entire box (k=1)?

Memory is hazy, but I do recall slight differences between the two. I haven't checked recently.

[hdrake]

Memory is hazy, but I do recall slight differences between the two. I haven't checked recently.

I have both variables easily accessible for the Baltic test case and can verify once PPAN is back up; I am pretty sure that setup has dz>1 at the surface.

[StephenGriffies]

@Hallberg-NOAA indicated that the top cell temperature is considered "piecewise constant" so that it is vertically constant throughout the layer.

[StephenGriffies]

Would be good to see @hdrake results to confirm the above.

[jkrasting]

The piecewise constant top layer makes sense now, thanks @StephenGriffies.

[hdrake]

Unless I'm missing something, it doesn't seem like the Sea Surface Temperature (tos) and Potential Temperature (thetao.isel(z_l=0)) have the same values. I tried converting tos to potential temperature and that didn't seem to make a difference.

[StephenGriffies]

I suspect the following explains the difference.

If you are looking at thetao from the 35-level regridded z output, then the top z diagnostic layer is certainly thicker than 1m. So then thetao and tos will certainly differ.

tos captures the upper ocean "SST" to communicate with the atmosphere. MOM6 has decided that SST is to be defined by a layer of thickness Hmix, which OM4.0 has set to 1m. Hmix is important particularly it layers get extremely, so that we do not send the atmosphere a "skin" temperature. In contrast to tos, the thetao output, generally regridded to z layers, is computed as a finite volume average over the z layers, which are generally thicker than Hmix.

[hdrake]

Great point, @StephenGriffies! Sorry for the oversight. The same plot but with natv coordinates confirms they are the same (for this case where dz>1 at the surface)!

[StephenGriffies]

Great! Thanks for confirming. Glad everything is self-consistent.

[jkrasting]

Yes, thanks @hdrake and @StephenGriffies. I believe I was also thinking of the remapped z-level output and not the native layer output. Glad this makes sense all around now.

[jetesdal]

@jetesdal, did you ever compare estimates of the boundary flux terms from your 2D vs. 3D computations (I guess this would be xwmt v0.0.2)?

I put together a comparison between surface WMT estimates using the boundary flux terms and surface WMT from fluxes. See my analysis here which I updated using data sourced directly from the GFDL filesystem. This analysis was done before the xwmt, package and was mainly done to look at the role of shortwave penetration.

When shortwave penetration is accounted for, you can get a pretty close match in WMT between the surface fluxes and the 3D boundary forcing:

Note, that we never implemented WMT calculations from SW penetration (rsdo). Also, note that this utilizes only the 'zstr' output. I always thought that surface fluxes + SW penetration should get you the same results as using the boundary tendency field. As you pointed out, using 'rho2' output might give you slightly different results.

[jetesdal]

@hdrake, upon further review of the output from the example notebook (basic_xwmt_tutorial.ipynb), it seems to me that the current implementation is not entirely accurate, but maybe I am missing something.

In the updated packages you need to isolate the surface fluxes after running wmt.map_transformations:

tmp = G_density[[v for v in G_density.data_vars if "surface_flux" in v]]

It appears that tmp includes only a subset of the water mass transformations related to heat fluxes:

> sorted(tmp.data_vars)
['surface_flux_basal', 'surface_flux_frazil_ice', 'surface_flux_latent', 'surface_flux_longwave', 'surface_flux_mass_transfer', 'surface_flux_sensible', 'surface_flux_shortwave', 'surface_flux_total']

However, the surface_flux_total term is in fact equal to the total surface WMT and matches the sum of the other surface_flux_* terms. I believe that currently all transformations from surface freshwater/salt fluxes (variables 'prlq', 'prsn', 'evs', 'friver', 'ficeberg', 'vprec', 'sfdsi') are being merged in surface_flux_basal.

In v0.0.1, we had a more detailed decomposition of the processes:

Can we still isolate these terms?

[hdrake]

Note, that we never implemented WMT calculations from SW penetration (rsdo). Also, note that this utilizes only the 'zstr' output. I always thought that surface fluxes + SW penetration should get you the same results as using the boundary tendency field. As you pointed out, using 'rho2' output might give you slightly different results.

As of https://github.com/NOAA-GFDL/xwmt/pull/40/commits/bd5b720c4e032571f31a57d505d651e00f91f2d9, the default budget specified by conventions/MOM6.yaml now includes the '3D' penetrative shortwave flux term rsdo instead of the '2D' flux rsntds.

I confirmed that this gets us even closer to the 3D boundary forcing term, although there is still a very small residual error for reasons I do not understand. I also verified that the way I am calculating the divergence of rsdo gives me the same thing as the online divergence diagnosed by rsdoabsorb.

[jetesdal]

I confirmed that this gets us even closer to the 3D boundary forcing term, although there is still a very small residual error for reasons I do not understand. I also verified that the way I am calculating the divergence of rsdo gives me the same thing as the online divergence diagnosed by rsdoabsorb.

I am not sure if I am looking at the right data to reproduce the plot, but I am not seeing rsdo in the Baltic_OM4_025 dataset. Could the mismatch actually be because shortwave penetration is in fact missing in the sum of surface WMT?

[hdrake]

I am not sure if I am looking at the right data to reproduce the plot, but I am not seeing rsdo in the Baltic_OM4_025 dataset.

@jetesdal, are you pointing to /archive/Graeme.Macgilchrist/MOM6-examples/ice_ocean_SIS2/Baltic_OM4_025/wmt_incsurffluxes.natv_rho2_zstr.monthly_daily_hourly.13months? This is a newer run that @gmacgilchrist did that includes rsdo and extends to 13 months to make up for the missing first snapshot (needed for the mass change term in xwmb).

[hdrake]

The most recent commit, https://github.com/NOAA-GFDL/xwmt/pull/40/commits/50e470b42b27a2c5bc6215f7b90c2c750a7bb320, includes some quality-of-life changes enabled by the new dependency package xbudget and a much more comprehensive budgets_dict that details multiple ways of decomposing the tracer budgets.

The most important corollary to this is that I have now sorted out what seemed to be an inconsistency between @jetesdal's 2D surface water mass transformation calculations (which are correct) and the way that I think all of us were interpreting the 3D WMT terms (which was incorrect).

Essentially, the problem was that we were missing a factor of $-\lambda Q_{m}$ on the LHS that comes from the dia-sea surface mass flux boundary condition on the dia-surface velocity. Also, because of the way the 3D boundary_forcing_* terms are diagnosed, they only include $Q{\lambda} + \lambda{m}Q{m}$ and need an extra factor of $-\lambda Q{m}$ to give the full non-advection ocean flux on the RHS.

The budget now closes very accurately, regardless of whether you use the 3D boundary_forcing_* diagnostics or the sum of all the 2D surface fluxes. The 3D diagnostics are still a bit better, and close to machine precision, because they can take advantage of online diagnostic regridding, which the 2D surface fields cannot.

[StephenGriffies]

Glad things are sorted @hdrake . But please clarify something. You say

"The most important corollary to this is that I have now sorted out what seemed to be an inconsistency between @jetesdal's 2D surface water mass transformation calculations (which are correct) and the way that I think all of us were interpreting the 3D WMT terms (which was incorrect)"

Is there any problem with the Groeskamp et al paper, or anything in my book, that needs to be clarified or corrected? I am under the impression that the literature is correct, and that @jetesdal followed the literature...

[hdrake]

No, @StephenGriffies, there is nothing wrong with the Groeskamp et al paper or anything in your book. That said, most of the literature (including those two references) does not discuss the issue of how to use the framework to infer spurious numerical mixing.

The main thing that is confusing is that the three-dimensional (layer-integrated) heat budget diagnostics in MOM6 do not actually line up with the terms that are discussed in the literature, because it lumps all tracer exchange fluxes into one single term, which does not correspond to the non-advective ocean tracer flux term that appears in the literature. This term is just stated, not derived, in Groeskamp; I only could wrap my head around where it comes from--and what that means for computing WMT from the MOM6 diagnostics--after reading several of the chapters in your book.

Another source of confusion is that sometimes calculations are right for the wrong reasons, e.g. because $\Theta{m} - \Theta = 0$ and $S{m} = 0$ by assumption and so omitting those terms would still give quantitatively correct answers but introduce conceptual errors elsewhere in the method.

[StephenGriffies]

Ok, thanks. I will keep these points in mind when revising the book.

[hdrake]

The to-do list for this PR is completed by the most recent commit (https://github.com/NOAA-GFDL/xwmt/pull/40/commits/fc602ad163df608928dd7dd531c713f81e6f4805), which adds support for WMT calculations when only 2D surface information (fluxes and lambda tracer) are available. Important use cases in this limit are presently available CMIP output and space-based remote sensing data products.

Alongside the source code changes is a new example which only uses 2D surface fields and shows the relatively small errors in the surface water mass transformation rates relative to the more accurate 3D calculation (see screenshot below, taken from the final figure in the example notebook).

[hdrake]

@jkrasting, I think this PR is ready to be formally reviewed and hopefully soon merged into main. Could you please assign @gmacgilchrist and @jetesdal (and yourself, if you want) as reviewers?

[jkrasting]

@hdrake - I added @jetesdal and myself as reviewers. I sent @gmacgilchrist an invitation that he must accept before I can add him to the reviewers list.

[Hallberg-NOAA]

Hi Henri,

This figure demonstrates very nice progress with this diagnostic, but the horizontal axis labels on this plot need to be reexamined. They are inconsistent between the two panels, and I suspect that both labels (or their units) are wrong.

• Bob Hallberg

  Dr. Robert Hallberg Oceanographer NOAA Geophysical Fluid Dynamics Laboratory

NOAA GFDL || Phone: (609) 452-6508 Princeton University Forrestal Campus || Cell: (732) 599-0459 201 Forrestal Road || Fax: (609) 987-5063 Princeton, New Jersey 08540-6649 || Email: @.***

On Sun, Aug 13, 2023 at 7:32 PM Henri Drake @.***> wrote:

The to-do list for this PR is completed by the most recent commit (fc602ad https://github.com/NOAA-GFDL/xwmt/commit/fc602ad163df608928dd7dd531c713f81e6f4805), which adds support for WMT calculations when only 2D surface information (fluxes and lambda tracer) are available. Important use cases in this limit are presently available CMIP output and space-based remote sensing data products.

Alongside the source code changes is a new example which only uses 2D surface fields and shows the relatively small errors in the surface water mass transformation rates relative to the more accurate 3D calculation. [image: Screenshot 2023-08-13 at 4 31 26 PM] https://user-images.githubusercontent.com/12971166/260332555-d0d3c32a-c974-41ea-a0e3-f3f5f4a384da.png

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/xwmt/pull/40#issuecomment-1676494805, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABR2BFLT7JVVANLJUWECUCDXVFPWXANCNFSM6AAAAAAXURYFEQ . You are receiving this because you were mentioned.Message ID: @.***>

[hdrake]

Thanks @Hallberg-NOAA, I corrected the figure in https://github.com/NOAA-GFDL/xwmt/pull/40/commits/ba8e41c80429bd2252289e6be12ed1b016599af0 and edited my comment above accordingly.

[hdrake]

A new notebook demonstrates how CMIP6 surface fluxes can be used as a limit case of the updated 3D WMT method. The first figure of the notebook exactly replicates the first figure of @jetesdal's xwmt.swmt tutorial notebook (and apparently with improved performance).

It was also found that xgcm.transform is much slower than xhistogram.xarray.histogram for global binning, so I changed the default back to xhistogram.

[hdrake]

This PR now includes two suites of tests: A) comparing water mass transformations computed from idealized depth-space tendencies in idealized stratification against exact analytical solutions B) replication of hard-coded results for the the full process-based water mass transformations for one day of the Baltic region test experiment

Note that the hard-coded results have changed for two reasons: 1) we now compute water mass transformations instead of water volume transformations (in Boussinesq model, just multiply by rho_ref to convert), 2) the test data file has changed, and now contains daily-averaged output for the date 1990-02-01 whereas it seems the previous test file was for a monthly mean centered on 1903-01-16