Add generic tracer budget diagnostics

[andrew-c-ross]

Work in progress: adds code for budget diagnostics developed by Enhui Lao, Fan Yang, and Mathieu Poupon. This code has been used in the past and confirmed to close, but I need to re-check this re-merged version. This also needs to be paired with a change to MOM6 (NOAA-CEFI-Regional-Ocean-Modeling/MOM6#12) to get all of the diagnostics (I believe one won't break the other, however).

Some of the tracer-specific budget diagnostics snuck in back when we split COBALT off into the CEFI version (for example, here) so I only needed to add these changes that apply to all generic tracers.

[yichengt900]

@andrew-c-ross,
Thank you! The new logic looks excellent. I have a quick question: I noticed that the calculation of g_tracer%vdiffuseh_impl seems identical to that of g_tracer%vdiffuse_impl:

ocean_BGC/generic_tracers/generic_tracer_utils.F90

Lines 3509 to 3510 in 3c965b7

	g_tracer%vdiffuse_impl(i,j,k) = h_old(i,j,k)*g_tracer_com%grid_tmask(i,j,k) * &
	(g_tracer%field(i,j,k,tau) - g_tracer%vdiffuse_impl(i,j,k))/ dt

I also confirmed this in my 1-D case. Am I perhaps overlooking something?

[andrew-c-ross]

Yeah the computation looks the same. The only difference I see is the units in the metadata (mol/kg m s-1 for vdiffuseh and mole/m^2/s for vdiffuse). I think mol/kg m s-1 should be right in Boussinesq mode.

[yichengt900]

Thanks, @andrew-c-ross. I also noticed the unit difference in the metadata, and it seems like vdiffuseh is the one consistent with the computation under Boussinesq mode. My question is, should we keep just one of them for now, or should we keep both?

[andrew-c-ross]

We've been using vdiffuseh for the budget, but it would be easy to switch. Rather than adding a new diagnostic vdiffuseh how about just correcting the units for vdiffuse?

[yichengt900]

Sounds good to me! It's more straightforward and requires fewer code changes.

[yichengt900]

Hi @andrew-c-ross,

Thanks for opening this PR. Overall, it looks good to me, and I tested the new diagnostics using a 1D case, confirming that it at least closes the O2 budget (see figure below). I have just a few minor comments on the code itself.

BATS 1D O2 budget

lhs = ['o2h_tendency', 'o2_advection_xy', 'o2h_tendency_vert_remap']
rhs = ['o2_vdiffuseh_impl', 'jo2']

[andrew-c-ross]

@yichengt900 is the g_tracer_diag routine even necessary? It is called by generic_tracer_diag, but I can't find anything that ever calls generic_tracer_diag. This code is duplicated from g_tracer_send_diag, which does get called by MOM.

[yichengt900]

@andrew-c-ross, good catch. I checked the master branch, and it seems that generic_tracer_diag isn't being used by any other programs. I think it's safe to revert the changes in g_tracer_diag and eventually remove both subroutines in a separate PR.

[andrew-c-ross]

I need to look more closely at the units before this gets merged. There are a few others that are off (some of which existed before this PR). jo2 and jo2c have the same units metadata, for example.

[andrew-c-ross]

Ok, I think this should be good, but would appreciate a second look.

• Concentration terms have units mol kg-1 s-1
• Content term (viffuse) has units mol kg-1 m s-1
• COBALT tendencies (jo2, etc) have units mol m-2 s-1 because they are multiplied by rho_dzt when writing the diagnostic. They need to be divided by rho when summing up the content budget terms afterwards.

I tried to standardize things a little around using "mol" instead of "mole" and using exponents consistently instead of a mix of exponents and division signs. On the MOM6 side, the metadata for the terms are still "mol/kg m s-1", which makes me twitch.

[yichengt900]

@andrew-c-ross, thank you for performing the additional checks on the units in the metadata to ensure greater consistency. I've reviewed both the codes and the diagnostic outputs (e.g. the O2 budget-related diagnostics) and confirmed that the units in the attributes are now consistent.

One minor point to note: some variable descriptions might be a bit confusing for users who don’t check the source code. For example, jo2c and jo2(multiplied by rho_dzt) have the same description in their long names:

	float jo2c(time, zl, yh, xh) ;
		jo2c:_FillValue = 1.e+20f ;
		jo2c:missing_value = 1.e+20f ;
		jo2c:units = "mol kg-1 s-1" ;
		jo2c:long_name = "O2 source concentration" ;
		jo2c:cell_methods = "area:mean zl:mean yh:mean xh:mean time: mean" ;
		jo2c:cell_measures = "volume: volcello area: areacello" ;
		jo2c:time_avg_info = "average_T1,average_T2,average_DT" ;
	float jo2(time, zl, yh, xh) ;
		jo2:_FillValue = 1.e+20f ;
		jo2:missing_value = 1.e+20f ;
		jo2:units = "mol m-2 s-1" ;
		jo2:long_name = "O2 source concentration" ;
		jo2:cell_methods = "area:mean zl:mean yh:mean xh:mean time: mean" ;
		jo2:cell_measures = "volume: volcello area: areacello" ;
		jo2:time_avg_info = "average_T1,average_T2,average_DT" ;

This is a minor issue, and I think it can be addressed in a future diagnostic naming cleanup PR. I'm approving this PR for now.

[charliestock]

Hello Andrew, YC - sorry to be running late on this. Could we spend some time on it tomorrow in the doc&dev and I'll have a thorough review by the end of the week?

[yichengt900]

@andrew-c-ross , Thanks for the additional commits to fix naming issues and add attribution. I think this PR is ready to go and I will merge it soon.