Use small minimum ustar in mixed layer restratification #251
Conversation
Codecov Report
@@ Coverage Diff @@
## dev/gfdl #251 +/- ##
============================================
- Coverage 37.12% 37.07% -0.05%
============================================
Files 263 263
Lines 73331 72750 -581
Branches 13669 13466 -203
============================================
- Hits 27222 26970 -252
+ Misses 41076 40905 -171
+ Partials 5033 4875 -158
📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more |
|
This PR would address the issue raised in mom-ocean#1168. Some of the discussion there might be relevant for the consideration of this PR. The tiny value of ustar being set in this PR, of order 10^-18 m s-1, is equivalent to a wind stress of order 10^-33 Pa. Along with the fact that this floor on ustar is being used with a max, it seems very unlikely that this would ever do anything except in the very unusual limit where both f and the applied wind forcing (and hence ustar) are 0. In any case, the resultant mixing timescale is likely to be so large that the CFL limits, and not the details of this value, will actually determine the strength of the overturning streamfunction. The only down-side that I see with this suggested correction is that we would have to justify the value of ustar_min. Note that the eddy-growth timescale in the limit of f=0 goes as ~ .135*h_vel/ustar. With a mixed layer depth of 10 m and the proposed tiny ustar gives a timescale of order 10^18 s, which is longer than the age of the universe (~4x10^17 s). In looking at the code, it strikes me that a physically based argument for what arrests the eddies in this unusual non-rotating and unforced limit could be molecular viscosity acting over the vertical scale of the mixed layer. Rather than putting a floor on the turbulent value of ustar, we might consider putting a floor on mom_mixrate that is proportional to pi2 * Kv_molecular / h_vel2. I would suggest coding this max function as something like the following which avoids divisions by potentially zero values: This seems physical enough to me, and it should give us a ceiling on the eddy growth timescale of order (8 / mom_mixrate), or ~ (h_vel**2 / Kv_molecular) with the new limit, or about 10^8 s for a 10 m mixed layer and a molecular viscosity of order 10^-6 m2 s-1. What would you think about this? I would suggest doing this via new runtime parameter (like KV_EDDY_ARREST) to set this Kv_molecular with a default of 0 m2 s-1 to avoid changing answers in any cases that currently work. What would you think of this alternative @angus-g ? |
|
Upon further reflection, this commit seems like it should be safe enough, and it does successfully address a known bug. Moreover, adding this tiny floor on ustar in this calculation should not be in conflict at all with also adding the viscosity based solution to this divide-by-zero problem suggested in the previous comment. Assuming that this passes the pipeline testing (which I am confident that it will), I think that we should merge this PR in. I will also follow up by adding the viscosity-based limit via a new runtime parameter in a separate PR. |
If the forcing ustar is exactly zero, the denominator of the momentum mixing rate term reduces to just the Coriolis parameter. With a grid construction that is symmetric about the equator, this term is also exactly zero, leading to a division by zero. To avoid this, a very small minimum ustar value is introduced, using the Earth's rotation as a velocity timescale, in the same manner as in some of the vertical parameterisations. This should prevent the denominator of the momentum mixing rate from going to zero.
Hallberg-NOAA
force-pushed
the
restrat-neglect
branch
from
6d5c636 to
3515a75
Compare
There was a problem hiding this comment.
This PR looks like a sensible solution to a known problem. Although this could change answers in some cases, in the sense that these cases would avoid divisions by zero, it seems unlikely that there will be adverse changes due to this miniscule floor on ustar for any cases currently in use. This PR has passed TC testing and it has passed pipeline testing at
https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/17600
|
(NOTE: In some circumstances with no background gustiness and wind stresses that go to zero, this commit will change answers!) |
|
There is a rescaling bug in this PR that I just noticed - the expressions for ustar_min have extra factors of US%T_to_s, but CS%omega is already in units of [T-1 ~> s-1], so the effective units of ustar_min are [Z s T-2 ~> m s-1]. Because ustar_min is used as a floor and has tiny (unscaled) values, for some range of values of T_RESCALE_POWER, the answers will be unchanged, but for larger values this bug will change answers. I will put in another pull request shortly that addresses this bug. |
Added the runtime parameters KV_RESTRAT and RESTRAT_USTAR_MIN, to build on the improvements in github.com/NOAA-GFDL/pull/251, and to provide run-time physical parameters to avoid the potential division by zero in the mixed_layer_restrat code noted at github.com/mom-ocean/issues/1168. Once this PR is merged onto the main branch of MOM6, that issue can be closed. By default, these do not change answers in the MOM6-examples test suite, but the default value for RESTRAT_USTAR_MIN was taken from the hard-coded value in PR that PR. The six copies of the eddy growth rate timescale calculations were consolidated into a new internal function, growth_time, with some other related minor refactoring of the code. Also, mixedlayer_restrat_register_restarts now takes a unit_scale_type arguments like many other analogous routines. All answers are bitwise identical, but there are new runtime parameters or comments that lead to changes in the MOM_parameter_doc files. Also clarified in the comments sent to the MOM_parameter_doc files how VISBECK_L_SCALE works as a dimensional scaling factor when it is given a negative value, and rescaled its units when read as though it were always in m.
Added the runtime parameters KV_RESTRAT and RESTRAT_USTAR_MIN, to build on the improvements in github.com/NOAA-GFDL/pull/251, and to provide run-time physical parameters to avoid the potential division by zero in the mixed_layer_restrat code noted at github.com/mom-ocean/issues/1168. Once this PR is merged onto the main branch of MOM6, that issue can be closed. By default, these do not change answers in the MOM6-examples test suite, but the default value for RESTRAT_USTAR_MIN was taken from the hard-coded value in PR that PR. The six copies of the eddy growth rate timescale calculations were consolidated into a new internal function, growth_time, with some other related minor refactoring of the code. Also, mixedlayer_restrat_register_restarts now takes a unit_scale_type arguments like many other analogous routines. All answers are bitwise identical, but there are new runtime parameters or comments that lead to changes in the MOM_parameter_doc files. Also clarified in the comments sent to the MOM_parameter_doc files how VISBECK_L_SCALE works as a dimensional scaling factor when it is given a negative value, and rescaled its units when read as though it were always in m.
Added the runtime parameters KV_RESTRAT and RESTRAT_USTAR_MIN, to build on the improvements in github.com/NOAA-GFDL/pull/251, and to provide run-time physical parameters to avoid the potential division by zero in the mixed_layer_restrat code noted at github.com/mom-ocean/issues/1168. Once this PR is merged onto the main branch of MOM6, that issue can be closed. By default, these do not change answers in the MOM6-examples test suite, but the default value for RESTRAT_USTAR_MIN was taken from the hard-coded value in PR that PR. The six copies of the eddy growth rate timescale calculations were consolidated into a new internal function, growth_time, with some other related minor refactoring of the code. Also, mixedlayer_restrat_register_restarts now takes a unit_scale_type arguments like many other analogous routines. All answers are bitwise identical, but there are new runtime parameters or comments that lead to changes in the MOM_parameter_doc files. Also clarified in the comments sent to the MOM_parameter_doc files how VISBECK_L_SCALE works as a dimensional scaling factor when it is given a negative value, and rescaled its units when read as though it were always in m.
Added the runtime parameters KV_RESTRAT and RESTRAT_USTAR_MIN, to build on the improvements in github.com/NOAA-GFDL/pull/251, and to provide run-time physical parameters to avoid the potential division by zero in the mixed_layer_restrat code noted at github.com/mom-ocean/issues/1168. Once this PR is merged onto the main branch of MOM6, that issue can be closed. By default, these do not change answers in the MOM6-examples test suite, but the default value for RESTRAT_USTAR_MIN was taken from the hard-coded value in PR that PR. The six copies of the eddy growth rate timescale calculations were consolidated into a new internal function, growth_time, with some other related minor refactoring of the code. Also, mixedlayer_restrat_register_restarts now takes a unit_scale_type arguments like many other analogous routines. All answers are bitwise identical, but there are new runtime parameters or comments that lead to changes in the MOM_parameter_doc files. Also clarified in the comments sent to the MOM_parameter_doc files how VISBECK_L_SCALE works as a dimensional scaling factor when it is given a negative value, and rescaled its units when read as though it were always in m.
There are extra US%T_to_s scaling factors in the expressions for ustar_min that were recently introduced with dev/gfdl PR #251; these are duplicative of the scaling factor that is already being applied when the parameter OMEGA is read in. The resulting expressions for ustar_min therefore effectively have units of [Z s T-2 ~> m s-1] when they should have units of [Z T-1 ~> m s-1]. Because ustar_min is a tiny floor on the magnitude of ustar, there is a range of values of T_RESCALE_POWER that will give the same answers as when it is 0, but for large enough values the answers will change, perhaps dramatically. This small commit removes these extra factors. Answers will change for some large values of T_RESCALE_POWER, but they are bitwise identical in the TC testing and in MOM6-examples based regression tests with modest or negative values.
Added the runtime parameters KV_RESTRAT and RESTRAT_USTAR_MIN, to build on the improvements in github.com/NOAA-GFDL/pull/251, and to provide run-time physical parameters to avoid the potential division by zero in the mixed_layer_restrat code noted at github.com/mom-ocean/issues/1168. Once this PR is merged onto the main branch of MOM6, that issue can be closed. By default, these do not change answers in the MOM6-examples test suite, but the default value for RESTRAT_USTAR_MIN was taken from the hard-coded value in PR that PR. The six copies of the eddy growth rate timescale calculations were consolidated into a new internal function, growth_time, with some other related minor refactoring of the code. Also, mixedlayer_restrat_register_restarts now takes a unit_scale_type arguments like many other analogous routines. All answers are bitwise identical, but there are new runtime parameters or comments that lead to changes in the MOM_parameter_doc files. Also clarified in the comments sent to the MOM_parameter_doc files how VISBECK_L_SCALE works as a dimensional scaling factor when it is given a negative value, and rescaled its units when read as though it were always in m.
Added the runtime parameters KV_RESTRAT and RESTRAT_USTAR_MIN, to build on the improvements in github.com//pull/251, and to provide run-time physical parameters to avoid the potential division by zero in the mixed_layer_restrat code noted at github.com/mom-ocean/issues/1168. Once this PR is merged onto the main branch of MOM6, that issue can be closed. By default, these do not change answers in the MOM6-examples test suite, but the default value for RESTRAT_USTAR_MIN was taken from the hard-coded value in PR that PR. The six copies of the eddy growth rate timescale calculations were consolidated into a new internal function, growth_time, with some other related minor refactoring of the code. Also, mixedlayer_restrat_register_restarts now takes a unit_scale_type arguments like many other analogous routines. All answers are bitwise identical, but there are new runtime parameters or comments that lead to changes in the MOM_parameter_doc files. Also clarified in the comments sent to the MOM_parameter_doc files how VISBECK_L_SCALE works as a dimensional scaling factor when it is given a negative value, and rescaled its units when read as though it were always in m.
* Add Leith+E
This commit adds the 2D Leith+E closure, which uses a modified 2D Leith biharmonic viscosity paired with a harmonic backscatter. ('Modified' here is not used in the same sense as 'modified 2D Leith'; it just means that the biharmonic coefficient is modified to account for enstrophy backscatter.) Variables are often named 'leithy' to refer to Leith+E.
The parameterization is controlled by three main entries in user_nl_mom:
1. USE_LEITHY = True
2. LEITH_CK = 1.0
3. LEITH_BI_CONST = 8.0
To use Leith+E you should have LAPLACIAN=True and BIHARMONIC=True. (It doesn't hurt to be explicit and also set LEITH_AH=False, along with any other viscous closures, but this is not required. If USE_LEITHY=True it will not use any of the other schemes. It does use the background value of the biharmonic coefficient as a minimum, but ignores the
background harmonic value.) LEITH_CK is the fraction of energy dissipated by the biharmonic term that gets backscattered by the harmonic term (it's a target; the backscatter rate is not exact.) Recommended values between 0 and 1. LEITH_BI_CONST is Upsilon^6 where Upsilon is the ratio between the grid scale and the dissipation scale for enstrophy.
Values should be greater than or equal to 1; 8 is a good place to start.
The code is sensitive to the background value of Ah; specifically, if Ah is too large, the code is unstable. This is because the backscatter coefficient is proportional to Ah, and if Ah is large then you get large backscatter. If your code is unstable, consider reducing, e.g., `AH_VEL_SCALE`.
* Background Ah
This commit updates the code so that it uses the background Ah as
a minimum. Previously, if `SMAGORINSKY_AH = True`, Leith+E would
use the Smag value of Ah as the minimum, which is incorrect.
* Improve logging
Removed `do_not_log` condition on `USE_LEITHY`
* Fix Leithy Logic
Added one line to fix the fact that the code would only work as
intended if either (i) writing out Ah_h, or (ii) in debug mode.
Also swapped .le. and .lt. for <= and <.
If the forcing ustar is exactly zero, the denominator of the momentum mixing rate term reduces to just the Coriolis parameter. With a grid construction that is symmetric about the equator, this term is also exactly zero, leading to a division by zero.
To avoid this, a very small minimum ustar value is introduced, using the Earth's rotation as a velocity timescale, in the same manner as in some of the vertical parameterisations. This should prevent the denominator of the momentum mixing rate from going to zero. My only question here is whether this is the right velocity scale to be using (and indeed if this is the right approach)?