Make neutral diffusion work with 3D diffusivities

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.