This commit adds the option to apply a linear decay in the
neutral diffusion fluxes within a transition zone defined
by the boundary layer depths of adjacent columns. This option
is controlled by a new parameter NDIFF_TAPERING, which is only
available when NDIFF_INTERIOR_ONLY=True. By default
NDIFF_TAPERING=False and answers are bitwise identical.

Simplifies and reduces the code by adding hbd to the neutral
diffusion contril structure. This avoid the need to "extract"
hbl multiple times. Answers are bitwise indenticals.

This commit adds a vertical dimension to the tracer diffusivities
(Kh_u and Kh_v) and associated coefficiets (coef_x and coef_y).
The following diagnostics were changed from 2D (lat/lon) to 3D
(lat/lon/depth): KhTr_u, KhTr_v, and KhTr_h.

To preserve old answers, the values of all modified arrays are
depth independent by default. The option to apply the equivalent
barotropic structure as the vertical structure of the tracer
diffusivity is also introduced and this can be controlled
via a new parameter: KHTR_USE_EBT_STRUCT (default is false).

Following up on the previous commit, where a vertical dimension
was added to the tracer diffusivities, this commit modifies the
HBD module to work with this change. To do so, parameter khtr_u
(diffusivity times the time step) is calculated at cell centers
and then remapped onto the HBD vertical grid. All unit tests
in this module were updated to conform with this change.

This commit also makes the default value of HBD_DEBUG equal to
the value set for DEBUG.

This commit modifies the neutral diffusion module to work with
3D diffusivities. When the diffusivities are depth dependent
(KHTR_USE_EBT_STRUCT=True), a new array (Coef_h, with values at
tracer points and at vertical interfaces) with a four-point average
between Coef_x and Coef_y is introduced. This array is then used to
calculate zonal and meridional neutral fluxes via optional arguments
and using an existing four-point average (vertical interfaces of two
tracer cells) inside subroutine neutral_surface_flux. The same approach
is already used when tapering the neutral diffusive fluxes. In this case,
however, the unit of the output from neutral_surface_flux (Flx) is
modified because the flux of the tracer between pairs of neutral layers is
multiplied by the average of Coef_h. To avoid double counting Coef_h,
the code block for updating the tracer concentration from divergence
of neutral diffusive flux components also had to be modified
for when KHTR_USE_EBT_STRUCT=True. Similar for
diagnostics trans_x_2d and trans_y_2d.

This commit also makes the default value of NDIFF_DEBUG equal to
the value set for DEBUG.

Follow Marshall Ward suggestion and rearrange the code
to be closer to what the compilers will do (or we hope
they would do).
This commit aims to potentially enhance performance.
Answers are bit-wise identical.