Copied dye_example.F90, renamed all public routines. This module does not
actually tie into MARBL at this time.

Adding USE_MARBL_TRACERS = True to override file turns on MARBL tracers. At
this point, we call marbl_instances%init and register all 32 tracers but don't
do anything else (so they are initialized to 0 and there is no source term for
advecting them yet)

CESM will use MARBL via manage_externals, but other systems may need to bring
it in via submodules

I think this is just needed for TravisCI, since it doesn't know to look for
code in pkg/

This looks like it was removed from everywhere else when I merged in the latest
dev/ncar branch

Created new configure_MARBL_tracers() subroutine to be called from
call_tracer_register() [between get_param() calls and register_* calls]

Also uses put_settings() to update the MARBL settings before initialization
(tested by setting ciso_on = .true. via user_nl_marbl)

marbl_instance%StatusLog is written after the call to init and in
marbl_tracers_end (which is now called from tracer_flow_control_end())

Failed a TravisCI test due to missing documentation

Also added a placeholder for parsing the MARBL timing information

Every call to print_marbl_log() is followed by a call to %erase()

Instead of passively advecting 0s, the MARBL_tracers module now correctly
initializes each tracer (but doesn't compute any source-sink terms yet)

Went from 0a806cf to 479f914

call the function from register_MARBL_tracers() rather than
MOM_tracer_flow_control.

MARBL will be brought in to CESM via manage_externals, and we will use
-DUSE_MARBL_TRACERS to build with access to MARBL.

Now that it comes in from manage_externals, we want to ignore it

Can build MOM without pkg/MARBL, but if USE_MARBL_TRACERS is True in the param
file and the build does not include -D_USE_MARBL_TRACERS then the model aborts.

Will allow MARBL diagnostics to be added to history files.

And, from that routine, call marbl_instances%surface_flux_compute().

Note that forcings, surface tracers, and surface fluxes are all zeroed out in
this commit. I'd like to get diagnostics posted in the next commit, and then I
can start updating tracer surface values, saving saved state, and looking into
how to read forcing fields.

Looks like column_physics() is the better place for this call

calls post_data (note that created a temporary data structure to hold both the
diagnostic id and a buffer to fill column-by-column as MARBL runs)

During configuration, set indices for each of the surface flux forcings so that
each forcing can be set to a different value in column_physics(); all are set
to zero except u10_sqr (2.5e5), atmpress (1), and xco2 / xco2_alt_co2 (284.7)

Also, cleaned out remnants of old dye_tracer code

Also updated how a few other forcing fields are passed

MARBL still wants things in cgs, so the conversion factors are private in the
MARBL_tracers module

Also, playing around with some debugging diagnostic output because it appears
that register_restart() isn't actually updating the field values (saved state
is initialized to 0 even when available in saved state)

It needs to be in the tracer_registry stage, not the initialize_tracers stage
so that fields will be updated from the restart file.

computes iron flux from dust and black carbon fluxes, adding some new variables
to parameter file

Also, only call surface_flux_compute() from ocean cells

1. Stale change that I forgot to commit earlier this month about how we apply
scale factor to NDEP forcing

2. Starting to put in calls that will be necessary prior to calling
interior_tendency_compute() [copying saved state, getting forcing fields,
setting up domain, etc]

Still don't actually make the call, but I register interior tendency
diagnostics and hopefully have all the data copies set up properly.

In loop that will eventually call interior_tendency_compute(), update
domain%zw, domain%delta_z, and domain%zt to get depths in m (and then convert
to cm when copying to MARBL structure)

Note that this commit does not support tracer restoring; that is disabled in a
temporary block of code after reading marbl_in.

There's a kludgy work-around where we set KMT to be the bottom-most level that
is more than 1 cm thick (to avoid sediment from accumulating in vanishing
layers) and we still don't read iron sediment flux

Read in (and do unit conversion) for fesedflux and feventflux as part of
initialization. Still need to figure out the vertical remap inside the time
step loop.

Use reintegrate_column() to map fesedflux from WOA z-grid to whatever vertical
levels MOM6 is using at the current time step. This is a very kludgy commit!

1. Loop through all levels (from bottom to top) and move any subsurface
   sediment flux up a level. The way the loop is ordered results in all
   subsurface flux ending up in the bottom-most layer. [On the source grid, so
   vanishing layers are not a concern.]
1. I created a 3D array to store the WOA vertical thicknesses, except I
   modified the thickness of the cell containing the bottom of the column
   (G%bathyT) to only be distance from the cell interface to bathyT, and the
   thickness of all cells below it are 0. This is the kludgy part of the
   commit.
2. apply a dilation factor of sum(thickness)/bathyT to the thicknesses of the
   source grid, so we are effectively always mapping between two columns of the
   same thickness.

1. read_z_edges returns depths that are 0 at surface and positive UP; algorithm
   I was using assumed positive DOWN
2. the logic in determining dz was hard to follow, so I rewrote the if
   statements to make it clearer
3. use dz(:) rather than marbl_instance%domain%delta_z to keep MOM units

This is just for reading nitrogen deposition, and it relies on FMS to handle
time interpolation for climatological forcing.

For now I use the POP naming convention of STF_{tracer} and J_{tracer}

Was trying to convert MOM's variables in units of m to MARBL's desired units of
cm by multiplying by 0.01 instead of 100

Copying the changes to generic_tracers, there are now stubs for the MARBL API
in config_src/external/MARBL so if pkg/MARBL is not available users can still
build MOM. Note that these stubs will trigger cause MOM to abort if they are
run with USE_MARBL_TRACERS=True in the parameters file.

Use a file with _FillValue=-1e34 rather than a file with NaNs

Since we are ignoring vanishingly thin layers at the bottom of the column when
passing data into MARBL, we also need to ignore what MARBL returns from those
levels. Now all tendency and diagnostic values for k>kmt are replaced with the
value from the kmt level. I also replace saved state below kmt with 0s though I
suppose an argument could be made for using the kmt value there as well

Introduced marbl_forcing_type, and fluxes%MARBL_forcings of that type. Also
introduced marbl_forcing_CS to handle all the parameter settings. Both of these
changes are aimed at reducing the footprint of the changes in
config_src/mct_driver needed to run MARBL; that will make it easier to bring
the same changes to config_src/nuopc_driver.

I think the last big change will be to create marbl_ice_ocean_boundary_type and
a function to make it easy to copy data from the new type into
marbl_forcing_type.

All MARBL-related forcing code is ready to be shared with nuopc cap

Passing ice_frac, u10_sqr, and seaice_dust_flux; need another round of updates
to handle the atm_dust_flux (split coarse / fine plus wet / dry) and the black
carbon fluxes (split between hydrophilic and hydrophobic)

Besides cleaning up some formatting, this commit adds {tracer}_zint,
{tracer}_zint_100m, and {tracer}_SURF diagnostics. For MARBL tracers, it also
adds Jint_{tracer} and Jint_100m_{tracer} (where J is the source / sink term
returned by MARBL)

Also renames marbl_forcing_type_main.F90 -> marbl_forcing_type_mod.F90 since I
ended up not needing _main / _aux designations.

Coupler recieves five fluxes used to compute dust_flux and iron_flux for MARBL
(atm fine dust, atm coarse dust, atm black carbon, sea ice dust, and sea ice
black carbon). These five fields are now available in history files.

If present(found), then return found=.false. instead of aborting if the
attribute isn't defined (needed for DEPTH:edges from my initial condition file)

Continued lines that previously extended beyond 132 characters

gfortran imposes line limit of 132 characters, but a script in .testing/ checks
for lines longer than 120...

These diagnostics need to be computed with a length scale in m, since we are
integrating over the column, but I had defined the diagnostic with the H_to_m
conversion. Also, I need to use H_to_m instead of H_to_z (or H_to_RZ)

Introduced Time_riv_flux to force reading the file at the first times step.
Also switched to a newer version of the file, which is already in mmol / m^2 /
s (so no need to convert from nmol / cm^2 / s).

For additive fluxes, e.g. atm_bc_flux = bcphidry + bcphodry, this commit
introduces a temporary array named marbl_work and computes the flux in multiple
steps:

1. atm_bc_flux = 0
2. import bcphidry -> marbl_work
3. atm_bc_flux += marbl_work
4. import bcphodry -> marbl_work
5. atm_bc_flux += marbl_work

as it turns out, state_getimport() is cummulative so this could be accomplished
with

1. atm_bc_flux = 0
2. import bcphidry -> atm_bc_flux
3. import bcphodry -> atm_bc_flux

That's coming in the next commit

Since state_getimport() is cummulative, we call it repeated for the forcing
fields that are the sum of multiple coupler fields

Introduce tracer_forcing_utils_mod to handle common interactions with
time_interp_external (tracking offset if model / data time axes are different,
setting an earliest / latest time to read from the file, etc) and then modify
marbl_forcing_type_mod to use this new code for the river fluxes.

A few lines exceeded the character limit imposed by the "Doxygen and Style"
continuous integration test

currently still using the KMT kludge, so this passes 1/dz in the bottom-most
non-vanishing layer. Next step will be to remove kmt from the driver and use
tracer_vertdiff to compute unit flux.

This commit still uses the KMT kludge, but now computes bot_flux_to_tend(:) by
applying a unit bottom flux to a column of 0s in tracer_vertdiff. This commit
includes a check to ensure that sum(dz(:) * bot_flux_to_tend(:)) == 1, though I
am not sure if we want that in the code base long-term.

This commit also fixes the units of the sfc_flux argument passed to
tracer_vertdiff() -- it should be Rho0 * STF instead of just STF. This requires
an update to MOM_tracer_diabatic (merged to dev/NCAR in 276954f but not yet on
this branch); for testing, I copied the updated version of that file but am not
committing it to make the future merge easier.

Last note: I believe the btm_flux argument to tracer_vertdiff is actually
positive upward, although the comments in MOM_tracer_diabatic.F90 claim it is
negative upward. I'll investigate a little more and open an issue ticket with
GFDL if that turns out to be the case, but in this commit I'm multiplying
bot_flux_to_tend(:) by -m_per_cm instead of just m_per_cm and I suspect it's
due to the sign of btm_flux

Note that this required increasing max_bracket_grow_it (I used 5 instead of 3,
but perhaps 4 would have been sufficient?)

Also moved the conservation check (sum(dz * bot_flux_to_tend) == 1) into MARBL

Functions in marbl_forcing_type_mod.F90 return immediately unless
USE_MARBL_TRACERS=True. I also cleaned up the way that module knows what
directory contains netCDF files to match how it is set up in MARBL_tracers.F90
(and make it easier to switch to DIN_LOC_ROOT when the time comes).

Building without pkg/MARBL was failing because the dummy interface was out of
date

Rather than relying on tracer_vertdiff, compute bot_flux_to_tend such that, for
a parameter BOT_FLUX_MIX_THICKNESS, bot_flux_to_tend = 1/BOT_FLUX_MIX_THICKNESS
for cells entirely contained that close to the ocean floor and then a weighted
value for the cell that is partially within the bottom boundary layer (0
elsewhere).

Relying on tracer_vertdiff led to convergence issues in the CO2 solver which
have not cropped up under this implementation

The previous algorithm had a (bathyT - zw(k)) term in it, and if these values
are close together then we can lose precision in the resulting difference.
However, this term can also be represented as the cummulative sum of dz(:) from
the bottom to layer k and that formulation is much more accurate

This is the first pass at a code cleanup; these indices should be set in
register_MARBL_tracers(), not initialize_MARBL_tracers(). The next commit will
further refactor this code.

There was a bug in the code where CS%tracer_inds would not be set during a
restart. Fixing the bug included a little more clean-up:

1. tracer_inds are set in a new routine rather than inline in
   initialize_MARBL_tracers()
2. This new routine is called from register_MARBL_tracers()
3. It should be possible to turn on BGC tracers in a branch / restart now

If ref_depth is below the bottom of the column, MOM should use _FillValue
rather than whatever value is reported from MARBL