Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add abiotic tracers #441

Merged
merged 64 commits into from
Oct 16, 2023
Merged

Conversation

mnlevy1981
Copy link
Collaborator

Port the abio_dic_dic14_mod.F90 module from POP to MARBL. We should be able to run with the two abiotic tracers without turning on what the base (non-isotopic) MARBL tracers.

Added base_tracers_on and abio_on to the list of MARBL parameters, and then
added ABIO_DIC and ABIO_DI14C to the tracer list if abio_on = .true.

At this point, we initialize the tracers correctly when abio_on = .true. (and
also handle abio_on = .true. and base_tracers_on = .false. correctly), but
there are no base_tracer_on checks in surface_flux_compute() or
interior_tendency_compute().
This supports creating marbl_in with abio tracers when running with 4p2z or
cocco (4p1z)
We'll refer to the abiotic carbon tracer module as abio_dic so abio_on ->
abio_dic_on (and a few other similar renamings). Also changed base_tracers_on
-> base_bio_on; at the same time, I changed references to "ecosys_base" to be
"base_tracers" instead (e.g. tracer indexing type).

Also cleaned up the stand-alone "request *" tests to be consistent in output
when there are none of whatever is being requested... e.g. abio_dic does not
request any forcing fields for the interior tendency computation.
MARBL will abort in either of the following two situations:

1. user tries to run without any tracer packages (i.e. base_bio_on,
   abio_dic_on, and ciso_on all set to .false.)

2. user tries to run with ciso tracers but not base biotic tracers (i.e.
   ciso_on = .true. while base_bio_on = .false.)

I also modified how init defines the ciso tracer indices; even though MARBL
will not run with base_bio_on = .false. and ciso_on = .true., allowing the
tracer_index_type to define the indices for the CISO tracers in that situation
lets the initialization process get to the consistency check and report the
actual error rather than crashing because the DI13C index is -1
marbl_abio_dic_surface_flux_mod contains marbl_abio_dic_surface_flux_compute(),
which currently sets the abio surface fluxes to 0 and then returns (assuming
abio_dic_on = .true.) but will eventually compute the correct surface fluxes.

marbl_abio_dic_interior_tendency_mod contains
marbl_abio_dic_interior_tendency_compute(), which doesn't really do anything
but will return immediately if abio_dic_on = .false.

I also added a couple of return statements so that we can run with abio_dic_on
= .true. and base_bio_on = .false. without running into memory / indexing
issues. This will be useful for testing the next set of commits, when we start
computing surface fluxes and interior tendencies for the abio_dic tracers and
when we start turning off settings / diagnostics when base_bio is turned off.
If lsource_sink is false, we want to return 0s for the tendencies right away;
no need to call abio routines.
@mnlevy1981
Copy link
Collaborator Author

Currently we are using same atmospheric CO2 forcing from GCM for biotic and abiotic tracers; while the majority of runs will likely use the same value for both, we do want to support running both modules with different CO2 forcings (Keith L gave the example of using emissions-driven CO2 for base biotic tracers in a coupled run but something different for abiotic)

First pass at the marbl_abio_dic_surface_flux_compute() so it doesn't just
return 0s. Still need to compare these fluxes to those computed by the POP
module.

Note that I updated the saved state type to add an abiotic surface pH (and
updated the internal metadata to specify the existing saved state variables
were for the biotic tracers).
marbl_abio_dic_interior_tendency_compute() returns tendencies. I also cleaned
up how the CISO module handles C14 decay, because those terms are also needed
in the abio module.
I don't know how the test suite didn't catch this, but the local
bot_flux_to_tend was the wrong type and never actually used.
@mnlevy1981
Copy link
Collaborator Author

If abio_on = .true. and base_biotic_on = .false., then we need autotroph_cnt and zooplankton_cnt set to 0 (may as well set max_grazer_prey_cnt as well, for completeness) so that setup_local_tracers() doesn't try to set tracers that don't exist. This issue is cropping up now because I just moved the call to setup_local_tracers() out of the if (base_biotic_on) block.

Updated some of the scripts to get base biotic diagnostics right when running
in cgs (was creating CaCO3_FLUX_10000m instead of _100m). Also updated the
R14C_atm computation, but there are still differences in the first time step
when comparing MARBL abio to POP; I think something funny is happening with ice
fraction, because STF maps show large differences in the polar regions.
MARBL_generate_diagnostics_file.py needs to pass unit system to the
settings_class and diagnostics_class
Starting with just a handful of the diagnostics provided by POP while I
investigate differences between POP's abio_dic_dic14 module and the new
abio_dic module in MARBL

Also changed the long name of a couple of the existing ECOSYS_ diagnostics to
clarify they are only used by the base biotic module
Had a couple of unused variables, and also wasn't setting the right index for
the ABIO_ALK_SURF diagnostic
With the exception of the derivative terms, all surface flux diagnostics from
the POP abio module are available in MARBL now.
store_diagnostics_ciso_interior() ->
marbl_ciso_diagnostics_interior_tendency_compute()

This change brings the routine in line with other public subroutines (by
beginning with {module}_), and when I copy this function to the abio module it
will also follow convention.
set interior diagnostics to 0 in interior_tendency_compute() so we don't
over-write the ABIO diagnostics; also, pass specific fields for abio interior
tendency diag computation instead of tracer_local and indices.
This was causing an issue in POP when trying to run without MARBL's abio
because the ecosys_diagnostics_operators file still listed ABIO_CO2_SCHMIDT as
a possible MARBL diagnostic even though it wasn't available
if autotroph_cnt = 0, then we don't want to process autotroph_settings(:) when
generating marbl_in
If base_bio_on == .false., autotroph_cnt, zooplankton_cnt, and
max_grazer_prey_cnt should all be 0. Note that I also modified
MARBL_settings_file_class to make sure that base_bio_on being false takes
precendence over value of PFT_defaults when setting autotroph_cnt, etc
If base_bio_on = .false. then PFT_defaults is treated as "None" regardless of
its value (a message is added to the log, similar to the way cocco and 4p2z are
treated as user-specified).

Also realized a better solution to the problem addressed in the previous commit
is to move PFT_defaults to general_parms2 so that its default can depend on
base_bio_on; this means that autotroph_cnt and other _cnt variables can be set
to zero when PFT_defaults = None, which is the new default for base_bio_on =
.false.

Lastly, cleaned up the diagnostics yaml file to note which diagnostics should
only be defined when base_bio_on is true
Moving PFT_defaults from general_parms to general_parms2 affected this script
need to move PFT_defaults to general_parms2
piston_velocity, xkw_ice, schmidt_co2, and pv_co2 can be computed once and
stored in surface_flux_internal rather than being computed in both abio and
base surface_flux_compute() routines.
Added a "defined_if" key to the YAML schema for settings files, and now
MARBL_generate_settings_file.py will only include variables that can be used by
the selected tracer modules (I also re-organized the YAML to put all the
[module]_on variables in a separate section that gets processed before
general_parms).
run_test_suite.sh now calls call_compute_subroutines test with

abio_dic_on = .true.
base_bio_on = .false.

And compares against a baseline generated on cheyenne (with gfortran)
There was a change to the units of Fefree in the last commit, so the netcdf
comparison script could not compare the field to the values in the baseline.
The new baselines for the default test, ciso test, and 4p2z now have proper
units.

Also, the baseline for the abio test only contains the abio fields (the last
commit changed what fields are written when base_biotic_on = .false.) and has
been renamed to explicitly note that it is abio_only rather than both the abio
and base biotic modules.
1. MARBL_settings_file_class.py should assume that the input_file contains
   variable values in the proper units (rather than expecting user to know
   units in YAML file); this is consistent with how the Fortran code reads
   input files
2. MARBL_share.py can't assume particulate_flux_ref_depth will be a float; if
   it is specified in an input file it is stored internally as a str

There were also two issues with running the test suite:
1. results from the netcdf metadata check were not being reported
2. the 4p2z metadata check needs to use the 4p2z settings file to generate the
   internal MARBL_settings object
There is a try / except block to avoid cases when particulate_flux_ref_depth is
not in MARBL_settings.settings_dict (currently just when base_biotic_on is
false) but I had an if statement checking the value of
MARBL_settings.settings_dict['particulate_flux_ref_depth']['attrs']['units']
outside that block.
Need to convert particulate_flux_ref_depth to float regardless of unit system,
otherwise setting it via input_file results in a string value.
Rather than do some kludgy unit clean-up in netcdf_metadata_check.py, the units
specified in the diagnostics YAML file are updated to the correct unit system
in MARBL_diagnostics_class (it's still pretty kludgy).

I also converted print statements to logging statements in
netcdf_metadata_check.py, and used some stack overflow magic to manage to the
log formatting in a new class in MARBL_share.py (it would be nice to propogate
this to other .py scripts in MARBL_tools)
Apparently using f-strings with logging messages is frowned upon.
If num_levels is a dimension of the netcdf variable, that should correspond to
"layer_avg" in the JSON; we don't have a num_iface dimension that would
correspond to "layer_iface", so everything else should have vertical_grid =
"none"

Also tried to clean up some comments.
There's no reason to use a different keyword for dependencies in the settings
file than we do for the tracers in the settings file or the diagnostics file
Inserted missing word in a comment, and removed assumption that any array of
length zero we want to remove from settings file is in the PFT_defaults
category (instead we loop through all categories looking for the variable we
want to remove)
Following code review, a few small changes (mostly updating comments); I did
rename the indices in the saved state type, which required updating a few other
modules as well
Turned some divisions into multiplications, cleaned up local variable names,
and added a new settings variable to control whether we compute the derivative
diagnostics. Also fixed a missing-import bug in MARBL_share.py
Need to turn all tracer modules on to be able to test all possible config flags
Want to use lflux_gas_co2 with both base_biotic and abio_dic tracer modules;
setting default value for lflux_gas_o2 = base_bio_on then cleans up a lot of
the logic in surface_flux_compute(). Also added spaces after several ".not."
statements and took xkw_ice out of the interior share type (it's only needed in
surface_flux_compute)
Remove "in permil" from long name of variables ("permil" is already noted in
the "units" attribute) and dropped "for base biotic tracer fluxes" from the few
variables where I had added it.

Also updated the long name of the abiotic derivative diagnostics to refer to
ABIO_FG_DIC (what the diagnostic is called) instead of STF_ABIO_DIC.
We set surface_fluxes(:,:) and interior_tendencies(:,:) to c0 at the beginning
of their respective compute routines, so the abio and ciso versions of those
routines don't need to set the specific subset to c0 as well.

Note that these subroutines no longer require abio_ind_beg, abio_ind_end, or
the ciso equivalents.
1. renamed lconstructed() -> lallocated()
2. only allocate indices for autotroph / zooplankton diags if base_bio_on is
   true (requires checking base_bio_on before deallocating)
   -- this reverts some of the allocate statements to put them back where they
     are on development (modulo additional whitespace because it is now in a
     new "if (base_bio_on)" block)
There were a few messages being added to the log inside if statements that
should always evaluate to .false. (there will always be at least one tracer,
and currently there is guaranteed to be at least one forcing field requested
for surface_flux_compute() as well)
1. ECOSYS_IFRAC, ECOSYS_XKW, and ECOSYS_ATMPRESS are only available if either
   lflux_gas_co2 or lflux_gas_o2 are true
2. PV_CO2 and SCHMIDT_CO2 are only available if lflux_gas_co2 is true

Adding [1] required adding a mechanism for specifying that only one of a set of
dependencies need to evaluate to true (the default behavior is that all
dependencies must evaluate to true).

I tested this by running with

base_bio_on = .false.
abio_dic_on = .true.
lflux_gas_co2 = .false.

And needed to add an if (lflux_gas_co2) block to
abio_dic_surface_flux_compute() to avoid seg faults.
Rather than computing totalChl_local in setup_local_tracers(), we pass
autotroph_local to compute_PAR and let that subroutine compute the sum of the
Chl tracer values.
The tracer initialization was moved into marbl_init_mod, but the next commit
will refactor it once more to create a marbl_init_tracer_metadata module
This splits marbl_init_mod.F90 up and makes it a little easier to read. I also
renamed all the private functions in both marbl_init_mod.F90 and the new
marbl_init_tracer_metadata_mod.F90 modules so they no longer begin with
"marbl_" -- this keeps with the convention that public subroutines should begin
with marbl_ to avoid namespace conflicts in the GCM.
The abiotic diagnostic long names were not consistent in how they were
capitalized, but now all of them are treated as titles (capitalizing all words
except articles, prepositions, etc. A few base biotic diagnostic long names had
comments suggesting improved names, and there was also a comment about changing
a variable name in the indexing type -- I addressed those as well.

Also changed some logic in how autotroph metadata is set -- rather than making
sure the index of every ciso tracer is positive, I just wrapped all of the
metadata updates in an "if (ciso_on)" block.
@mnlevy1981
Copy link
Collaborator Author

The tests that fail in 5a1f06e pass on my laptop, but fail on cheyenne (with gnu) -- I'll investigate further on that machine, it's not clear what changed in that commit to cause problems.

Not every autotroph will have a Ca13CO3 or Ca14CO3 tracer, even if ciso_on is
true. Kind of shocked these tests passed on my laptop when accessing index 0 of
an array...
We don't have a documentation page about what tracers are available in MARBL,
but we do provide instructions for adding a new tracer and some of the sample
code listed on that page has changed with the addition of the abio tracers (and
subsequent refactoring). Updating the documentation also highlighted a few
issues with comments in the Fortran code, and pointed out an issue with how we
were setting lfull_depth_tavg
I cleaned up the diagnostics calls and also the metadata in
diagnostics_latest.yaml so that we only define the abio_dic surface flux
diagnostics if lflux_gas_co2 = .true. Originally I was just trying to track
down a seg-fault when running with abio_dic_on = .true. and lflux_gas_co2 =
.false. -- it turns out that the xco2 forcing is not requested if lflux_gas_co2
= .false., but the abio module was still trying to add it to the diagnostics.
The fix is to treat abio surface fluxes similar to the base biotic surface
fluxes, and only define some of them if lflux_gas_co2 = .true.
@mnlevy1981 mnlevy1981 marked this pull request as ready for review October 16, 2023 16:03
@mnlevy1981 mnlevy1981 merged commit 6e6b2f7 into marbl-ecosys:development Oct 16, 2023
1 check passed
@mnlevy1981 mnlevy1981 deleted the add_abio branch October 16, 2023 21:59
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

MARBL_settings_class assumes units in input_file match YAML, not unit_system
1 participant