Add MARBL to MOM6 (#157)

* Update MOM6_DA_hooks submodule

* Create a placeholder MARBL tracer module

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

* Add flag to turn on MARBL tracers

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)

* Add MARBL as submodule

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

* Add doxygen documentation to marbl_instances

* Update submodules to use https not ssh

* Softlink to MARBL source code

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

* Remove diag_to_Z_CS

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

* Move MARBL%init out of register_MARBL_tracers

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

* Read marbl_in from run directory

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

* Add ability to write MARBL log to stdout

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

* Add doxygen documentation for print_marbl_log()

Failed a TravisCI test due to missing documentation

* Call MARBL's shutdown() routine

Also added a placeholder for parsing the MARBL timing information

* Erase MARBL log after printing to stdout

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

* Read in MARBL tracers IC file

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

* Update MARBL to latest development commit

Went from 0a806cf to 479f914

* Change git submodule commit for MOM6_DA_hooks

* Make configure_MARBL_tracers() private

call the function from register_MARBL_tracers() rather than
MOM_tracer_flow_control.

* Removed MARBL as submodule

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

* Add pkg/MARBL to .gitignore

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

* Introduce _USE_MARBL_TRACERS macro

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.

* Register MARBL diagnostics

Will allow MARBL diagnostics to be added to history files.

* call MARBL_tracers_surface_state()

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.

* Move call to surface_flux_compute()

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

* add surface flux diags to history output

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)

* Initialize surface flux forcings better

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)

* Move module memory into CS type

* Provide T & S for surface forcing

Also, cleaned out remnants of old dye_tracer code

* Add saved state for surface fluxes

* Get ice fraction from the coupler

* Get u10_sqr from coupler

Also updated how a few other forcing fields are passed

* Clean up old comment

* Better unit conversion

* Module parameters for unit conversion

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

* Easy clean-up based on feedback from Andrew

* Add missing variable declaration

* Add saved state to restart file

* Add surface flux to tracer_vertdiff call

* If tracers mandatory in restart, so is saved state

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)

* Move call to setup_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.

* Get dust and iron fluxes from coupler

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

* Add option to read NDEP from a file

Also, only call surface_flux_compute() from ocean cells

* Update NDEP scale factor, work to interior_tend

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]

* More support for interior_tendency_compute() call

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

* Update MARBL domain per-column

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)

* Add indices for interior tendency forcings

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

* Copy t & s -> interior_tendency_forcing

* Call marbl_interior_tendency_compute()

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

* Add parameters to read in fe sed / vent fluxes

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.

* Get FESEDFLUX into MARBL

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.

* Had a few issues in last commit

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

* Use time_interp_external() not MOM_read_data()

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

* Add surface_flux and tendency to diagnostic output

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

* Fix bug in units of z-coordinate vars

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

* Remove _USE_MARBL_TRACERS CPP macro

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.

* Update ndep file

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

* Ignore levels below kmt from marbl_instance object

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

* Refactor where MARBL forcings are defined

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.

* Finish refactoring mct driver

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

* Start updating forcings in 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)

* Add several new diagnostics for tracers

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.

* Register and post FLUX_CPL diagnostics

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.

* First pass at adding support for river fluxes

* Fix to read_attribute_str()

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)

* Code clean-up:

Continued lines that previously extended beyond 132 characters

* Update doxygen documentation

* One more round of formatting clean-up

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

* Proper units for _zint and _zint_100m

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)

* More river flux clean-up

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).

* Add more fluxes to nuopc cap

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

* Remove marbl_work from mom_cap_methods

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

* Refactor interacting with time_interp_external

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.

* Fix formatting to make doxygen happy

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

* Pass bot_flux_to_tend to MARBL

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.

* Use bot_flux_to_tend from tracer_vertdiff

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

* Remove KMT kludge

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

* Add bot_flux_to_tend to diagnostic output

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

* Avoid allocating MARBL memory unless needed

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).

* Add bot_flux_to_tend to dummy interface

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

* Split long line that doxygen flagged

* Add BOT_FLUX_MIX_THICKNESS parameter

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

* Reformulate bot_flux_to_tend algorithm

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

* Set default thickness for bot_flux_to_tend to 1m

* Set tracer_inds outside of (.not. restart) block

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.

* Clean up how / where tracer_inds are set

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

* Clean up some comments

* Check ref_depth for 2D diagnostics from MARBL

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

* Initialize negative tracer concentrations to 0

If the model is not initializing from a restart file, then we treat all
negative tracer ICs as 0.

I also swapped a couple of i- and j- loops to run through contiguous memory

* Fix bugs in iron_flux computation

Two major issues:

1. iron_flux was missing a few terms
2. I was keeping the units as kg / m^2 / s when MARBL wanted nmol / cm^2 / s

I also switched the default riv_nut file to one interpolated from JRA

* Need marbl_constants_mod.F90 for non-MARBL builds

I added a reference to molw_Fe in marbl_forcing.F90, so I needed to update the
dummy driver in config_src/externals/MARBL to recognize that use statement

* Update doxygen documentation

Forgot to document molw_Fe in the dummy marbl_constants module

* Clean up how we read FESEDFLUX files

1. Default files now have DEPTH_EDGES variable so dz can be length ke
2. Improve logic for moving flux from below the ocean floor to bottom column
3. use v_extensive=.true. when registering FESEDFLUX diagnostic

Also cleaned up how conversion factors are applied to dust / iron flux

* Don't modify values read from restart!

I had the logic to set negative tracer values to 0 in the wrong place; now it
only applies to tracers read from initial conditions instead of also applying
to tracers read from restart

* Fix bug in accumulating 2D fields at ref_depth

Use post_data(mask) instead of accumulating missing_val - this avoids
inadvertently introducing round-off error when taking the average of an array
full of missing_val (which would then not be missing_val any more)

* Only allocate memory we plan to use

Don't need to allocate memory in the diagnostic type for diagnostics that MOM6
is not including in any history files

* Include (i,j) indices in MARBL errors

When a specific column in MARBL returns an error, MARBL will now print both the
global lat,lon (which it printed previously) and the global (i,j) indices

* NUOPC cap improvements

Fields that were posted from the cap (MARBL forcing fields) were not appearing
correctly in the netCDF output because we needed to enable averaging

* Bugfix in NUOPC cap

I was using (i,j) indices instead of (i-i0,j-j0) when pulling data from
MARBL_IOB

* Fix formatting

A couple of lines failed the line-length check in the CI

* Add KPP Nonlocal Terms to MARBL tracers

This commit mimics the changes made to MOM_CFC_cap.F90 and
pseudo_salt_tracer.F90 to apply KPP Nonlocal terms to that MARBL tracers

* Convert dust_flux to cgs before sending to MARBL

* Move riv_flux to applyTracerBoundaryFluxesInOut

Rather than adding river fluxes to CS%STF, created new array CS%RIV_FLUX that
is used as in_flux_optional argument to applying boundary fluxes

Also added a missing doxygen comment

* Major refactor of MARBL forcing fields

Removed marbl_forcing_type and marbl_ice_ocean_boundary_type, keeping the flat
structure of forcing and ice_ocean_boundary_type, respectively. This greatly
reduces the amount of code in marbl_forcing_type_mod.F90 (which now needs a new
name), and adds a little more overhead to the driver layer but with the benefit
of reusing pointers for things like ice_fraction and u10_sqr which are needed
in the CFC cap as well.

I also started the process of adding multiple ice category support, but there's
more to do for that feature.

* REVERT MCT CAP

From this point forward, MCT will not support the MARBL driver. CESM users must
use NUOPC cap for runs with MARBL

* Rename marbl_forcing_type_mod.F90

There is no longer a forcing type, but this module does handle some forcing
fields for MARBL so marbl_forcing_mod.F90 is a more appropriate name.

* NUOPC cap set to receive all ice cat fields

If CPL_I2O_PER_CAT=TRUE, then the nuopc cap allocates memory in
ice_ocean_boundary_type to store the five fields that POP uses with MCOG:

* sf_afrac
* sf_afracr
* Foxx_swnet_afracr
* Fioi_swpen_ifrac_n
* Si_ifrac_n

Next step will be to work from the forcing_type side, and make sure those
fields get copied into the appropriate arrays if the user requests running with
multiple ice categories.

* More per-category forcing updates

1. ice_ncat is stored as-is, rather than storing ice_ncat+1
2. Memory is now allocated on the forcing type if USE_ICE_CATEGORIES is true
   (Default is false, want default to be true when running with MARBL)

* More ice category cleanup

MARBL_tracers will get ICE_NCAT from parameters file (still need to set default
correctly!) instead of passing it down the calling tree.

* More ice category bug fixes

1. several off-by-one errors due to using ice_ncat+1 instead of ice_ncat
2. needed i0 and j0 when copying from IOB to fluxes
3. NUOPC is case sensitive when getting field, but doesn't abort when case is
wrong

Also added the fields passed to MARBL as MOM6 diagnostics (FRACR_CAT_N and
QSW_CAT_N, for N from 1 to ice_ncat+1) -- should probably switch range so it is
0 to ice_ncat, with category 0 representing open ocean

* don't post diagnostics if not requested

when the _CAT_ diags were not in the diag_table, the run was aborting. Adding a
check to make sure the id > 0 avoids that.

* Remove kludgy threshold

Before tracking down issue with i0 and j0, I had set an artificial threshold on
fracr_cat, where I was treating values below 1e-5 as 0. This commit undoes
that, and only treats negative values as 0 (leaving small positive values
alone)

* Rudimentary MARBL support in solo_driver

All the MARBL-specific forcings the NUOPC cap gets from the mediator are set to
0:

* ice_fraction
* u10_sqr
* dust_flux
* iron_flux

Will need to work out details on how else to populate these fields
(ice_fraction is set to 0 when allocated, so that's probably okay)

* Change dust / iron parameter defaults

Fortran defaults are now what we want for C / G compsets, and MOM_interface
will override them for B compsets

* atm_press = 1 atm when p_surf_full is unavailable

A kludge to set surface pressure to 1 atm when using solo_driver (which
allocates memory for p_surf but not p_surf_full and then does not set p_surf)

* Update IC file and units of RIV_FLUXES

The old MARBL initial condition file had some issues at depth when the MOM
topography was deeper than POP; these cells were all set to 0 because the
restart file we mapped tracers from does not include a land mask. Fix was to
mask out the POP data prior to mapping it, then the lateral fill was applied
correctly.

Additionally, Keith noticed a bug in our call to
applyTracerBoundaryFluxesInOut() -- the in_flux_optional argument should be in
units of conc m, but we were passing CS%RIV_FLUXES with units of conc m/s.
Since the in_flux should be the time-integrated value, we set

CS%RIV_FLUXES = CS%RIV_FLUXES * dt

prior to passing it to applyTracerBoundaryFluxesInOut()

The comments in the declaration of RIV_FLUXES was also updated to account for
the unit change.

* update fesedflux files

Somewhere in transition from testing in cesm2_3_alpha05b to testing in
cesm2_3_beta08, I forgot to commit a change to update the iron sediment flux
forcing files. The previous files were generated in a buggy manner, and these
new files provide better forcing.

* Use data_override for some MARBL forcing

ice fraction, u10_sqr, and the various dust / black carbon fluxes that MOM and
POP receive from the CESM coupler in the NUOPC cap can now be read in from
netCDF file using data_override in solo_driver/

also, added a flag (READ_RIV_FLUXES, default: .true.) to let us turn off
looking for river flux files.

* Support CHL_FROM_FILE=FALSE when using with MARBL

If CHL_FROM_FILE is FALSE, then MARBL_tracers will request total_Chl from MARBL
and MOM_tracer_flow_control::get_chl_from_model() can retrieve it.

* Expand dummy cap for MARBL

The previous commit used more of the MARBL interface, so I expanded the dummy
cap to allow MOM6 to build without the full library (as before, setting
USE_MARBL_TRACERS=TRUE but building the dummy cap will result in a FATAL error)

* Refactor to use MARBL's get_output_for_GCM()

Removed code that relied on interior_tendency_output since that no longer
exists in MARBL

* API for MARBL's get_output_for_GCM() changed

Updated MARBL_tracers_get_output_for_GCM() to account for fact that MARBL
doesn't want tracers passed in as an argument, it just wants to use
self%tracers(). Also did some minor clean-up to
MARBL_tracers_get_output_for_GCM():

1. pass in G and GV so that we can index arrays properly / loop through i,j
   (and skip land cells)
2. the do j= and do i= loops are on the same line
3. abort if the MARBL function returns an error

Lastly, cleaned up the interface to hopefully pass CI again (needed to create
marbl_settings_mod and add dummy get_output_for_GCM() function to interface)

* Use MOM_initialize_tracer_from_Z not tracer_Z_init

I added an optional argument ongrid to MOM_initialize_tracer_from_Z() which
gets passed through to horiz_interp_and_extrap_tracer()

* More updates for reading IC file

Switched CESM default to use IC file written on 1x1 grid (with WOA depths)
instead of using a file on the tx0.66v1 grid. Also changed some of the defaults
set in get_param() calls in MOM_initialize_tracer_from_Z() to match the
defaults elsewhere in the code.

* Update dummy interface for get_output_from_GCM()

In MARBL, this changed from a function to a subroutine but I forgot to make the
corresponding change in the dummy MARBL cap

* Get NDEP from NUOPC instead of reading from file

* Update solo_driver to handle ndep

Last commit broke solo_driver because the API to
convert_marbl_IOB_to_forcings() changed. Also, cleaned up some comments in
MARBL_tracers.F90 and the NUOPC cap.

* Remove NDEP_SCALE_FACTOR from parameters file

I think this was a parameter because I was mimicking shr_stream (which includes
a scale factor as part of the namelist options), but it makes far more sense to
combine it with the existing ndep_conversion variable.

Given the way the parentheses were used, this should be bit-for-bit.

* Add MARBL_TRACERS_INIT_VERTICAL_REMAP_ONLY option

If the initial condition file for the MARBL tracers is already on the MOM grid,
this will skip the horizontal interpolation step. This option is not necessary
for CESM-MOM6 (we want to interpolate ICs from the WOA grid), but is useful for
MOM6-examples, where we are interpolating to a single column grid offline.

* Update MARBL tracer IC file

This one has been updated to set negative values -> 0 and then was run through
the autotroph consistency check (if any of Chl, C, P, Fe, or Si are 0 they
should all be 0)

* Move atm_co2 and atm_alt_co2 to MOM_forcing_type

This commit also moves ATM_CO2_CONST and ATM_ALT_CO2_CONST to marbl_forcing_mod
and uses those values for the new forcing fields as preparation for possibly
getting Sa_co2diag and Sa_co2prog from the mediator.

Unrelated, I changed a handful of instances of "else if" to "elseif" to match
what is done elsewhere in the code (these are all elseif statements I
introduced on this branch earlier in development)

* NUOPC cap can receive CO2 if provided

MOM6 advertises for it, but if it is not available then it gets removed from
importState and memory is deallocated so there is no attempt to get it (i.e.
there is no error condition if the atmosphere does not provide it)

* Support using coupler-provided atm_co2

Added ATM_CO2_OPT and ATM_ALT_CO2_OPT, which default to "const" but also
support "prognostic" or "diagnostic" if the coupler is providing those fields.

Also renamed marbl_forcing_mod -> MARBL_forcing_mod to be consistent with
capitalization in MARBL_tracers.F90

* Replace logical flags with integer

For ATM_CO2_OPT and ATM_ALT_CO2_OPT, we do a string -> integer conversion
instead of trying to track all possible options via logicals. This introduces
some new module-level parameters in MARBL_forcing_mod.F90 (and I noticed a
formatting issue in MARBL_tracers.F90 when checking how other modules handle
parameters)

* NUOPC cap can pass CO2_FLUX to atmosphere

Added Faoo_fco2_ocn to the exportState and modified sfc_state to be able to
pass ocn_co2 to ocean_public, which then gets exported in
mom_cap_methods:mom_export()

In src/tracer/, I needed to set up MARBL's surface_flux_output and use it to
request co2_flux. MOM_tracer_flow_control now calls
MARBL_tracers_surface_state(), which copies the flux from MARBL's control
structure to sfc_state.

I also did a bit of code cleanup, deallocating more arrays in the MARBL control
structure prior to deallocating CS itself.

* Only copy co2 to srf_state if memory was allocated

solo_driver doesn't allocate memory in sfc_state to pass CO2 flux back to the
atmosphere and was seg-faulting without the if (allocated) check

* Code clean-up following review

1. sfc_state%sfc_co2 -> sfc_state%fco2
2. re-order some if states and do loops (want if statements outside do loops as
   much as possible)

* Remove spaces in "end if" and "end do"

The MOM6 style guide explicitly states we should use endif and enddo with no
space

* More code clean-up

1. solo_driver has callTree_leave() call in MARBL forcing override routine
2. allocate() statements use the source= argument as much as possible
3. cleaned up spacing in if statements [a mix of "if(condition) then" and "if
   (condition)then"
4. removed unnecessary use statement in mom_cap_methods, and cleaned up vague
   comment
5. renamed ocean_public%ocn_co2 -> ocean_public%fco2_ocn

* Clean up NUOPC cap

Move more calls into if (cesm_coupled) blocks because they require CESM
forcings from CMEPS

* Updates to use support_mks branch of MARBL

This branch of MARBL requires the unit_system argument in a few places;
sticking with cgs introduces a few round-off level differences due to
restructing some internal MARBL computations.

* Update dummy driver to add unit_system args

Because the API to MARBL changed, I needed to update some code in
config_src/externals/MARBL

* Use MARBL in mks, not cgs

Requires the support_mks branch of MARBL as well as generating marbl_in with
--unit_system mks

* Update default for iron forcing files

MARBL_FESEDFLUX_FILE and MARBL_FEVENTFLUX_FILE were updated to remove the
incorrect 1D horizontal dimensions (and rename the dimensions nx and ny so that
categorize_axes could still figure out how to read the files)

* Updated pkg/MARBL, which changed API

A couple arguments changed in functions called from MARBL_tracers.F90 and I
caught a mistake in the comments of MARBL_forcing_mod.F90.

Also updated the dummy API so the model continues to build when MARBL is not
available.

* Dummy MARBL API needs one more function

The MARBL driver calls marbl_instances%get_conc_flux_units(), which I had
forgotten to add to the dummy cap

* MARBL API changed

replaced surface_flux_output%add_output with add_output_from_GCM(); this commit
uses the new API from MARBL and also updates the dummy API so we can still
build without MARBL

* Don't need to overwrite tracer_restore_vars

A MARBL update (reflected in MOM_interface) puts empty strings in
tracer_restore_vars by default for MOM6, so we don't need to remove the POP
defaults in the Fortran code anymore

* Use time_interp_external for restoring

This is still a work in progress, but I cleaned up how the MARBL tracer
restoring fields are translated into something the MARBL_tracers control
structure can parse and also call init_external_field. One issue I'm having is
that time_interp_external needs a time_type argument, and that's not available
in MARBL_tracers_column_physics. I think the solution will be to move a lot of
this code to MARBL_forcing.F90, but I want to do that in a separate commit.

* First pass at implementing tracer restoring

Uses time_interp_external to temporally interpolate restoring fields; assumes
data is on correct spatial grid, and then does vertical interpolation to go
from data's vertical grid to current MOM6 grid (similar to iron sediment flux).

Still to do: vertical interpolation of restoring time scale

* Clean-up to avoid truncation errors

Was failing some CI tests due to truncation issues

* Switch from interpolate_column to remapping_core_h

interpolate_column() is meant for interpolating from cell interfaces, not cell
centers. MOM_initialize_tracer_from_Z() uses ALE_remap_scalar(), which just
calls remapping_core_h() under the hood, so that's what we want to use for
vertical interpolation of the restoring fields and time scales.

* enable_averaging -> enable_averages

previous merge contained a change to this function name, so my MARBL code
additions were stuck calling the wrong one

* Move river flux code into MARBL_tracers

1. MARBL_forcing_mod should just be for fields that pass through the MOM6 cap
2. River fluxes shouldn't be in the forcing datatype
3. MARBL_tracers_set_forcing() is already set up to handle fields read in with
   time_interp_external

* move post_data calls for river nutrient fluxes

Diags are posted after files are read in set_forcing(). Also, cleaned up where
I apply the dt factor to convert from flux to time-integrated flux because
set_forcing is only called once during the first two time steps so we were
inadvertently applying a factor of dt^2 on the second time step.

This commit is bit-for-bit with 69202f7 in my CESM testing, which wasn't the
case for 5540525

* Code cleanup: doxygen test

Several "Line length exceeded" messages in CI, addressed by adding continuation
lines instead

* More doxygen clean-up

* Missing "<" in one comment

Had a ! instead of !< when documenting a variable in a class

* First pass at adding ABIO

Still need to update d14c forcing, currently hardcoded to use -4 (we want to
read from a file with dimensions time,lat_band and interpolate in time)

* Add MARBL_TRACERS_MAY_REINIT to param file

The default is still false, but in some cases (branching off a run that did not
have MARBL enabled) it would be useful to set as true instead

* Add dummy get_setting() to marbl_interface_class

Now that I call get_setting() to determine which MARBL tracer modules are
enabled, the dummy cap also needs this function

* Skip some processes when not base_bio_on

If base_bio_on is false, we don't want to request any output from the GCM or
deal with the iron flux forcing fields

* Add support for reading d14c forcing from netcdf

Use time_interp_external to read in d14c in three latitude bands; in putting
this together, I also found a bug in tracer_forcing_utils that resulted in
being off by a year when reading constant forcing (river fluxes were
interpolated to Jan 1, 1901, rather than Jan 1, 1900; fixing it also meant
updating the forcing file so there was data to read on Jan 1, 1900, since the
original dataset begins on July 1 of that year).

Also, following the GFDL MOM6 call, I added parentheses around the square term
in "a * b**2" constructs [this was a bit-for-bit change on derecho, but some
machines treat "a * b**2" as "(a*b)*b" instead of "a*(b*b)"]

* Update to support marbl0.46.0

That tag changed how total 3d chlorophyll is passed from MARBL to the GCM

* Update interface to build without MARBL

marbl0.46.0 updated the MARBL interface, so that needs to be reflected in the
config_src/ version

* Shorten line that exceeded max length

doxygen test was failing because I added some whitespace between variable
declaration and inline comment; I broke the comment over two lines to fix

* Fix whitespace

MARBL_tracers.F90 and MARBL_forcing_mod.F90 now comply with whitespace rules
from the MOM6 style guide

* Check abio_dic_on and base_bio_on before posting

There are a few diagnostics that are only defined if base_bio_on=.true. (the
river flux nutrient forcing fields), and one that is only defined if
abio_dic_on=.true. (the d14c forcing); some compilers won't initialize the
diagnostic ids to 0 in the control structure, so we need to either explicitly
initialize all the ids or only call post_data when we know the ids have been
set. This commit does the latter.

* Updates for dimensional scaling test

Currently fails T-scaling test with solo driver, probably fails lots of other
scaling tests as well. This commit

1. Adds debug output to MARBL_tracers.F90
2. Gets dimensions correct in comments of MOM_forcing_type, MARBL_forcing_mod,
   and MARBL_tracers
3. Scales forcings correctly for the MARBL surface_flux_compute() step (at
   least in T); output highlights issues in computing source / sink term from
   interior_tendency_compute()

One of the biggest changes from this commit is the handling of units for the
nitrogen deposition fluxes. It looks like they were coming in as kg/m^2/s,
being converted to mol/L^2/T in fluxes%{nhx_dep,noy_dep}, and then converted to
mmol/m^2/s when copied into MARBL. Now the intermediate stage is mmol/m^3 Z/T;
this is not bit-for-bit with the previous setup because I went from multiplying
by (1000/14) (kg -> mol) and then another 1000 in the third step (mol -> mmol)
to just multiplying by 1e6/14 (kg -> mmol) in the second step.

* More dimensional scaling updates

With solo_driver, the following runs are all bit-for-bit with non-scaled runs:

C_RESCALE_POWER = 10
H_RESCALE_POWER = 10
L_RESCALE_POWER = 10
S_RESCALE_POWER = 10
T_RESCALE_POWER = 10
Z_RESCALE_POWER = 10

* Clean up line-lengths in some comments

Should pass doxygen test again

* pass phys units to convert_marbl_IOB_to_forcings()

The function is meant to help copy fields from the ice_ocean_boundary_type
(which is in physical units in all the caps) to the forcing_type (which wants
scaled units). So the solo_driver should NOT scale the dust, black carbon, or
NDEP inputs from data_override, and instead that scaling should happen in
MARBL_forcing_mod.F90

* scale riv flux

applyTracerBoundaryFluxesInOut expects in_flux_optional in units of conc H, and
we were passing conc m T/s. Since riv_flux_loc is now conc H, I also added a
debug-gated hchksum on it.

* Introduce MARBL_IC_MIN_VAL for testing

The dimensional scaling tests fail if the MARBL tracer concentrations are very
very small (O(1e-300)); this can be avoided by setting the minimum tracer value
to be 1e-100 instead of 0. We don't want to do this for production runs,
though, so the default for this parameter is still 0.

* Fixed a few area correction bugs

Sa_co2prog and Sa_co2diag should not be area corrected (they are states) but
Faoo_fco2_ocn should be (it's a flux)

* No support for global ops yet

When calling marbl_instance%init(), we should tell MARBL that MOM6 doesn't have
the global operators that MARBL expects (global sums / running means) so we get
the appropriate error message when trying to run with ladjust_bury_coeff = True

* Add chksum calls for MARBL forcings

Updated ice_ocn_bnd_type_chksum() in the NUOPC cap, though I don't think this
function is ever called

* MARBL input data is now in INPUTDATA

I had created CESM_INPUTDATA as a parameter to point to my work directory, but
it is no longer necessary because INPUTDATA points to the CESM input data
repository and I've moved necessary files there

* Changes following code review

-- cleaned up a lot of comments and whitespace
-- used source argument in more allocate statements, and deallocated more
   arrays
-- 3D diags now have zl:mean in cell_methods attribute
-- marbl_instances%domain%kmt is set once (during initialization)

* Call MARBL_tracers_stock()

* Only use MARBL for Chl when using base_bio tracers

If MARBL is not configured to provide the base biotic tracers, then it will not
be able to provide chlorophyll. In that case, if CHL_FROM_FILE=False, MOM6
needs to get chlorophyll from the generic tracers.

* tracer_forcing_utils moved into MOM_interpolate

To make these subroutines more accessible, they were moved out of src/tracer/
and made available through MOM_interpolate

* Fix whitespace in comments

* Add some variable descriptions

If variable was described in POP comment, I copied the comment over. Otherwise
I came up with a description on my own.

* Use do loops instead of ':'

time_interp_external() does not update halo regions, so running CESM with
DEBUG=TRUE was triggering some overflows from uninitialized memory. Intead of
copying the entire array, we now loop through (is:ie,js:je) when accessing an
array returned from time_interp_external()

* Add parameter to change restoring time scale name

Most use cases don't include restoring for MARBL tracers, but when that feature
is enabled and the time scale is read from a file the user can specify what
variable to read from the netCDF file (default is I_TAU to match naming
convention in MOM6, but some test cases are based on POP files and will need to
read RTAU)

.gitignore

@@ -4,6 +4,15 @@
 html
 
 
+# Build output
+*.o
+*.mod
+MOM6
+build/
+deps/
+pkg/MARBL
+
+
 # Autoconf output
 aclocal.m4
 autom4te.cache/