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

Support mks unit system #430

Merged
merged 54 commits into from
Aug 31, 2023
Merged

Conversation

mnlevy1981
Copy link
Collaborator

This PR will handle the cgs -> mks conversion in its entirety, rather than focusing on diagnostics first and parameters later. The approach we have decided on:

  1. introduce optional unit_system argument to marbl_instance%init()

    1. If present, call subroutine in marbl_settings_mod.F90 to update several module variables (unit_system is a string that stores the desired unit system, cm2len converts from cm to proper length unit, and g2mass converts from grams to proper mass unit); defaults are 'cgs', 1._r8, and 1._r8.
    2. When setting default parameter values, cm2len and g2mass will be used to convert from cgs to selected system (e.g. caco3_bury_thres_depth = caco3_bury_thres_depth * cm2len)
  2. Update python script that generates settings file to allow --unit-system argument (so POP will generate marbl_in using cgs while MOM generates marbl_in in mks)

  3. Update diagnostic metadata so units are in proper system

  4. There are several fortran parameters (and some hard-coded "magic numbers" as well) that have units that will change based on the unit system

Fixes #41

Default is cgs, but user can change that via the option unit_system argument to
marbl_instance%init(). This is just a first pass to update the default value of
a few parameters, still have several outstanding tasks:

1. Update python script that generates settings file to allow --unit-system
   argument
2. Update diagnostic metadata so units are in proper system
3. There are several fortran parameters (and some hard-coded "magic numbers" as
   well) that have units that will change based on the unit system
@mnlevy1981 mnlevy1981 linked an issue Jun 7, 2023 that may be closed by this pull request
Rather than have more module variables in marbl_settings_mod, there is a
unit_system_type object on the marbl_interface_class and it is passed to
marbl_settings_mod via init() (to set defaults correctly). I also started
going through surface_flux_compute() to make sure units are consistent (Keith L
and I walked through the first few blocks of code in that routine together)
There was a circular dependency when this data type was in
marbl_interface_private_types.F90
1. CI was failing because of trailing whitespace in marbl_settings_mod.F90
2. Removed duplicate comment from marbl_surface_flux_mod.F90 (changed units
from cm -> L and g -> M, but left old version in as well)
By default, tests are all run with cgs... but updating test.nml makes it easy
to run in mks instead.
Rather than using test.nml, added -u (or --user_system) argument to command
line
Variables that are changed in fortran when running in mks are also changed in
MARBL settings file generated by python script
Use different units depending on whether unit system is cgs or mks
Even in cgs, we want IRON_FLUX diagnostic to have units mmol/m^2/s
Use these definitions for defining diagnostic units. Also use these units for
determining proper units of forcing fields.
Since we are storing conc_units, conc_tend_units, and conc_flux_units, it makes
sense to use them when setting up tracer_metadata
The tracer metadata type contains units, tend_units, and flux_units; no reason
for the driver to construct tend_units and flux_units on its own
We only convert forcing to cgs when unit_system = 'cgs'
NHx_SURFACE_EMIS has proper units whether running in cgs or mks.

Also updated xkw_coeff to have proper units regardless of unit_system, so
piston velocities are computed in L/s rather than cm/s
DOP_remin and POC_REMIN_DIC_zint were still 'mmol/m^3/s' and 'mmol/m^3 cm/s',
respectively... not sure how I missed them in an earlier commit
Instead of just listing tracers, this test now prints tracers and units.

Sample:

-----------------
Requested tracers
-----------------

1. PO4 (units: nmol/cm^3)
2. NO3 (units: nmol/cm^3)
3. SiO3 (units: nmol/cm^3)
4. NH4 (units: nmol/cm^3)
Setting up a dictionary to track unit conversion scale factor as well as old /
new units. Currently done by calling a function which returns the dictionary,
which is a terrible way to do it...
There were several mol / kg -> mmol / m^3 conversions that assumed we had
density in g/cm^3; these have been updated to convert M/L^3 -> kg/m^3 ahead of
the conversion. Lots of small file mods to make sure the unit_system_type can
be passed from the interface to marbl_co2calc_mod.F90
Before comparing a variable in the baseline to the same variable in the new
file, look for the 'units' attribute in both. If the attribute exists in both
files but differ between them, look up the conversion factor in a dictionary.
Also update the baseline for call_compute_subroutines to keep CI passing
interior_tendency_compute() passes unit system to co2calc subroutine; also
found the first few instances of parameters in cgs. More updates to
interior_tendency_mod will be needed.
I think compute_particulate_terms() is now correct, though
compute_dissolved_organic_matter() might still be missing a unit conversion or
two. Cleaning up the particulate terms allowed me to remove cperm from
marbl_constants_mod.F90.

I also changed the behavior of particulate_flux_ref_depth: it is now a real
variable instead of an integer, and it gets converted to model length units
instead of being meters. This let me remove the local variable
particulate_flux_ref_depth_cm from compute_particulate_terms().
The ref_depth attribute of the diagnostic types is now in model units, not
hardcoded to be units, so we don't convert it in marbl_io_mod.F90 before
masking out columns that are shallower than the reference depth.
Typo in comments meant I was missing a unit_scale factor when computing PAR,
which affected lots of diagnostics.
Default shallow depth is given in meters and converted to unit length
unlike other fields, which are mks in forcing files, iron sediment flux is cgs
so it didn't have a conversion factor
Unclear why it's needed here, but it lets us get the scavenging terms right
regardless of unit system
in src/, cleaned up conversion for dust terms (getting dust%rho from nmol/g to
mmol/kg when running in mks takes care of the rest of the terms) and also found
a couple of missing conversion terms for iron. Also corrected unit metadata for
a few parameters, and converted them to cgs / mks as needed.

Those metadata changes are reflected in the yaml files in defaults/, and I also
updated MARBL_settings_file_class.py to account for those differences. Lastly,
I cleaned up the netcdf_comparison.py tool.

The following tests show only a few small differences to track down:

1. Run call_compute_subroutines.py with both cgs and mks, compare the output
2. Use MARBL_generate_settings_file.py to generate settings files in both cgs
and mks, use those with the -s option in [1]
I changed the API for requesting output for the GCM, it now requires a
unit_system argument... but the last commit didn't include the update in
call_compute_subroutines_drv to pass the new argument.
Three new variables:

parm_Fe_scavenge_rate0_per_sec
parm_Lig_scavenge_rate0_per_sec
parm_FeLig_scavenge_rate0_per_sec

This way we are not changing the value of the non-per_sec versions after they
are set.
It's okay for new data to be non-zero when baseline is 0 as long as abs value
of new data is less than atol. I think this is a reasonable change, since we
are effectively treating all values smaller in magnitude than atol as 0.
Added a "TODO", fixed a typo, and removed a "^M"
Denominator is now the column max abs value, rather than pointwise absolute
value... we were getting failures in mks near sign changes in some fields,
where absolute difference was roughly the same throughout the column but values
near zero were three or four orders of magnitude smaller than the rest of the
column.
Run call_compute_subroutines with a single instance using cgs and compare
against the baseline, then use mks for comparing 1 inst, 2 inst, and 5 inst
with each other (also comparing the 1 inst mks to the baseline)
pylint was complaining about too many variables in _variable_check_loose() so
I created _compute_rel_err() instead of creating local variables col_max and
rel_denom in the first function.
Forgot to cd to MARBL_tools and set HIST_ROOT before doing the cgs run vs cgs
baseline comparison
Use the max value over a three-point stencil in depth instead of the maximum
over the entire column. (Some variables may decay exponentially with depth, and
we don't want an overly-large denominator to mask issues in the baseline
comparison)
This may end up getting refactored, but two major updates to marbl_io_mod.F90:

1. Hard-coded fallback method for reading non-isotopic tracers and scaling them
   for CISO
2. Hard-coded constant values for d13c and d14c forcing fields (these should
   probably get moved to the input file, where we also include fields like
   atm_co2 that are constant in POP)

setting d14c in the manner described in (2) above requires reading in latitude
from the input file and that changes the API of a few marbl_io_mod.F90
subroutines called from marbl_call_compute_subroutines_drv.F90.

Lastly, I added a test with CISO to run_test_suite.sh. Once we have completed
the mks conversion, I can add an mks vs cgs check as well.
There are a couple of fluxes computed in ciso_interior_tendency_mod that needed
terms from unit_system instead of assuming nmol/cm^2/s, but most of the heavy
lifting was updating the metadata for the tracers and diagnostics (i.e. they
were computed in mks but the metadata presented them as cgs in the netcdf
output).

Added an mks vs cgs ciso comparison to the test suite.
1. Updated marbl_io_mod to initialize cocco tracers (7% of sp values, except
   coccoCaCO3 = spCaCO3 since sp is no longer calcifier) as well as microzooC
   and mesozooC (both equal to zooC)
2. Added marbl_with_4p2z settings files for both cgs and mks
3. Added a new 4p2z baseline (the cgs run on my laptop)
4. Updated run_test_suite.sh to make sure github tests with 4p2z

Expecting the 4p2z mks test to fail on github, locally netcdf_comparison
returns

Variable: POC_REMIN_DIC
... Max relative error (1.4034448149204683e-11) exceeds 1e-11

Variable: PON_REMIN_NH4
... Max relative error (1.403445286451049e-11) exceeds 1e-11

I suspect the github relative errors will be similar, and I'll need to bump the
tolerance to 2e-11.
Had bad filename when comparing mks 4p2z to baseline (expected failure due to
values exceeding tolerance, got FileNotFound instead)
rtol = 2e-11 instead of 1e-11 lets POC_REMIN_DIC and PON_REMIN_NH4 pass CI
This commit reformatted the file to break offending file over two lines...
@mnlevy1981 mnlevy1981 mentioned this pull request Aug 4, 2023
@mnlevy1981
Copy link
Collaborator Author

@klindsay28 and I have had a couple of conversations about determining appropriate thresholds for netcdf_comparison.py and making sure differences between cgs and mks are really the propogation of round-off level errors due to unit conversion. After merging #436 I compared gnu on derecho with gnu on my Mac (both of these are cgs, my mac has an intel chip):

$ ./netcdf_comparison.py -b ../tests/input_files/baselines/call_compute_subroutines.history.nc -n ../tests/regression_tests/call_compute_subroutines/history_1inst.nc --strict loose
Comparing ../tests/regression_tests/call_compute_subroutines/history_1inst.nc to the baseline ../tests/input_files/baselines/call_compute_subroutines.history.nc

Variable: Fe_scavenge
... Max relative error (3.6994203973617207e-11) exceeds 2e-11

Variable: Fefree
... Max relative error (4.3527441876818616e-11) exceeds 2e-11

Variable: J_Fe
... Max relative error (6.007586938336526e-11) exceeds 2e-11

Variable: Lig_scavenge
... Max relative error (3.699352433955277e-11) exceeds 2e-11

Variable: P_iron_PROD
... Max relative error (3.6994203973617207e-11) exceeds 2e-11

Differences found between files!

So I think the threshold change in 53a5632 is justified (and should be increased further, to 1e-10, so derecho vs baseline passes).

Although it's worth noting:

  1. The threshold was previously raised due to some remin terms failing, and those are not the same variables that are getting flagged in the hardware comparison.
  2. cgs vs mks passes on derecho, even with the old threshold

The important takeaway is that I'm pretty convinced the unit conversion has been done correctly, and it would be nice to modify the test thresholds to reflect that (I want a threshold where derecho vs baseline passes and cgs vs mks passes)

@mnlevy1981 mnlevy1981 marked this pull request as ready for review August 10, 2023 20:44
Keith and I reviewed the fortran changes in src/ and this commit addresses
concerns raised during that review. Most of it is straight forward, I do want
to highlight one change that I think still needs improvement:

Instead of passing unit_system_opt ('cgs' or 'mks') to
marbl_instance%surface_flux_output%add_output(), we call
marbl_instance%get_conc_flux_units and then pass the conc_flux_units directly.
This still feels kludgy, but the way the interface and the
marbl_output_for_GCM_type are defined it was the only option for a quick fix.
I'd like to revisit this interface and come up with a complete redesign of the
way we register fields the GCM would like MARBL to return.
Instead of calling marbl_instances%surface_flux_output%add_output() directly,
the driver should call marbl_instances%add_output_for_GCM() which will then
call add_output() with all the data needed from unit system (currently just
conc_flux_units).
Updates to the tests/driver_src/ directory:

1. rename unit_system -> unit_system_opt (matches what is in libmarbl.a)
2. improve comments throughout the driver
3. fix bug in CISO forcing
4. add true fallback option for reading I/O (if 4p2z fields are in IC file,
   don't scale 3p1z values!)
Removed some old python2-specific code (using str instead of type(u''))
Also cleaned up some comments and the formatting of some new dictionaries.
Mentioned unit_system in a handful of pages; also changed the example for
calling init from marbl_init_drv.F90 to marbl_call_compute_subroutines_drv.F90
and that highlighted a missing status log check... so now call_compute checks
labort_MARBL after calling init()
@mnlevy1981 mnlevy1981 merged commit 951d1cd into marbl-ecosys:development Aug 31, 2023
1 check passed
@mnlevy1981 mnlevy1981 deleted the support_mks branch August 31, 2023 15:13
@mnlevy1981 mnlevy1981 mentioned this pull request Feb 23, 2024
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.

Units for MARBL
1 participant