-
Notifications
You must be signed in to change notification settings - Fork 58
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
Ice shelf melt parameterization fixes #395
Conversation
Codecov Report
@@ Coverage Diff @@
## dev/gfdl #395 +/- ##
============================================
- Coverage 38.08% 38.07% -0.01%
============================================
Files 269 269
Lines 76977 76987 +10
Branches 14207 14209 +2
============================================
Hits 29313 29313
- Misses 42351 42361 +10
Partials 5313 5313
... and 1 file with indirect coverage changes 📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
Hi @claireyung thanks very much for submitting this. I added one very minor comment, but otherwise looks very straightforward. @MJHarrison-GFDL (or maybe @OlgaSergienko) could you review this? |
Both changes, Line 530 and 569, are correct and should be accepted taking into account suggestions by @marshallward and @adcroft |
(1) It turns out there are more places in the ice shelf melt code where terms are out by a factor of the von Karman constant. This occurs every time the parameter ZETA_N or its reciprocal I_ZETA_N is used. It appears that ZETA_N was not written in the code to be the \xi_N (note the different Greek letters) that Holland and Jenkins 1999 use as a stability constant, but rather a scaled version, ZETA_N = \xi_N/k where k is the von Karman constant (VK in MOM6). With \xi_N = 0.052 (HJ99) and k=0.40, this gives ZETA_N=0.13 and resolves a comment on its definition where someone presumably questioned why ZETA_N used to be 1/8 (approximately 0.13) and changed it to 0.052 (see line 265) Just to verify this, this linked document shows all the places where ZETA_N is used and the equivalent equations from HJ99. I have thus reinstated the von Karman constant that I previously suggested removing (though I removed the double division as suggested) and instead redefined the default value of ZETA_N to scale it by the von Karman constant, along with adding an explanation of its meaning. This of course extends the correction of my previous commit to other places where ZETA_N is used. (2) A further issue I found is that in Line 640 there is a line to iterate the boundary salinity that overrides the previous if statement. This means that the false position method in line 635 is never used (or technically, it can be used but is then immediately overwritten). I have committed a change to remove this line. Sorry for the extra additions, but I hope this will improve the accuracy of the code to reflect the original Holland and Jenkins parameterisation. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you very much, @claireyung for this thorough documentation of the derivation of this parameter and the identification of the extra von Karman constant. Moreover, it is clear that missing the update to the buoyancy flux in the iterations is obviously a bug. I am convinced that these changes are correct. Again, thank you very much for this valuable contribution!
However, there is the potential (perhaps certainty) that other people have been doing MOM6 ice-shelf simulations with the old (wrong) value of ZETA_N and the bug you have fixed. We try to avoid making code changes that irreversibly change answers, even if they are incorrect. The way that we deal with this challenge is to set parameters (and revise their default values) via a get_param call, with the parameters stored in a module's control structure (here the ice_shelf_CS).
In this case, there should probably be two new parameters - one to replace the hard-coded parameter ZETA_N, perhaps something like "ICE_SHELF_LINEAR_SHELF_FRAC" to set CS%Zeta_N, which woudl then replace ZETA_N on line 532. A second would be a bug-fix flag to enable the old (wrong) answers to be recovered, perhaps something like "ICE_SHELF_BUOYANCY_FLUX_ITT_BUG" to set CS%buoy_flux_itt_bug, and then line 571 would become if (.not.CS%buoy_flux_itt_bug) wB_flux = wB_flux_new
. About a year after we add this bug-fix flag, we can obsolete it and then remove it from the MOM6 code if there are no objections. I would be happy to work with you on this if there is anything that is unclear about what I am suggesting here.
While we are at it, it would probably make sense to also add runtime parameters to replace VK and RC and the hard-coded 1.0e-4 convergence criterion on about lines 562 and 617 of this same module, but these improvements could sensibly be deferred to a subsequent PR.
Hi @Hallberg-NOAA, thanks for the instructions. I have added parameters for Zeta_N, RC, VK, the convergence criterion and two iteration bug fixes (where the defaults are set to correct the model i.e. change answers). Please let me know if you have any feedback! |
Hi Claire, You are very nearly there, but the automated checking is detecting some minor problems with the dOxygen comments. This is giving us the red "X" on the main list of PRs and on the conversation page. If you click on the "Details" for the failing test, the message you get is:
It should be easy to fix the 3 instances of trailing white space and the two problematic bits of LaTeX syntax. |
…lse position method
…ting false position method bug
11398ce
to
2e4fa1b
Compare
This test is now passing all of our automated code style and self-consistency checks, and it is passing the MOM6-examples pipeline testing at gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/20184. |
Two fixes to the ice shelf melt parameterization in
MOM_ice_shelf.F90
:(1) the removal of an extra von Kármán constant which differs from the Holland and Jenkins (1999) formulation. My working is detailed in this document.
(2) a modification to the it3 loop in
shelf_calc_flux
subroutine, which currently does not iterate becausewB_flux
does not get updated to its new value via the Newton solver method. Instead, the loop either runs 30 times with the same value, or is below the threshold and exits on the first loop. I added a line to update the value ofwB_flux
.