[Hallberg-NOAA]

Revised the dimensional rescaling in MARBL_tracer to work in non-Boussinesq mode. This includes using the subroutine thickness_to_dz() to convert layer thicknesses into vertical distances, which in turn requires a new thermo_var_ptrs argument to MARBL_tracers_column_physics(), and dozens of minor changes inside of the MARBL_tracers module.

The dimensions of 3 arrays in the MARBL_tracers_CS were changed from [Z ~> m] to [H ~> m or kg m-2], and another changed in the opposite way to reflect how they are actually being used. The units of two internal variables in MARBL_tracers_column_physics() were similarly changed. As a part of this 2 of related rescaling factors were changed in scale= arguments, and others omitted as they are no longer needed. Relatedly, two conversion factors were added in register_diag_field() calls and two others modified.

In 8 places where G%bathyT was being used in isolation, the code was modified to reflect that this is a depth relative to an arbitrary reference height, but this can be accounted for by including the mean sea-level.

In 4 places calls to the Fortran sum() function were replaced with explicit vertical sums so that the order of arithmetic in these sums is not potentially compiler dependent.

The hard-coded tolerance for the interface heights to match the bathymetry was replaced in 3 places with a variable as a possible prelude to making it a runtime parameter.

A total of 12 GV%H_to_Z or GV%Z_to_H conversion factors and 2 direct references to GV%Rho0 were eliminated. In non-Boussinesq mode, these are problematic because they introduce a dependency on the Boussinesq reference density. In Boussinesq mode, these changes take advantage of the fact that GV%Rho0*GV%H_to_Z = GV%H_to_RZ, while in non-Boussinesq mode GV%H_to_RZ does not depend on GV%Rho0.

None of these changes (apart from possibly the replacement of the sum function for some compilers) will change any answers in Boussinesq mode, but in non-Boussinesq mode there should no longer be any dependency in MARBL_tracer on the value of the Boussinesq reference density, RHO_0. There is a new argument to the publicly visible subroutine MARBL_tracers_column_physics().

[Hallberg-NOAA]

There are several minor issues addressed within this single commit PR along with the main issue of eliminating the dependency on the Boussinesq reference density when the model is run in non-Boussinesq mode. If it would ease the review process on the NCAR side, these minor issues (for example, eliminating the possible complier-dependent answer changes due to the use of the sum() function) could be handled in separate PRs.

[mnlevy1981]

I have merged the most recent dev/ncar into a local copy of this branch (there were small conflicts from #459) and am running the CESM test suite.

@Hallberg-NOAA if it's okay with you, I will push my merge as well as one or two small fixes caught by the test suite to your branch and give you a chance to approve those changes before we merge in

[mnlevy1981]

The k dimension of CS%fesedflux_dz has length CS%fesedflux_nz, not GV%ke; this line was flagged in a CESM test run using DEBUG mode with

Subscript #3 of the array FESEDFLUX_DZ has value 35 which is greater than the upper bound of 34

I think the proper solution is to make h_stretched allocatable. To avoid allocating memory every time we enter MARBL_tracers_column_physics, h_stretched can be added to the control structure for this module.