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

Possible zero in the denominator in MOM_mixed_layer_restrat #1168

Closed
gustavo-marques opened this issue Jul 24, 2020 · 5 comments
Closed

Possible zero in the denominator in MOM_mixed_layer_restrat #1168

gustavo-marques opened this issue Jul 24, 2020 · 5 comments
Assignees
Labels

Comments

@gustavo-marques
Copy link
Collaborator

We have a global test case where the corner point is on the Equator. Once we set GUST_CONST=0.0, following the discussion about default parameter changes in one of the dev/calls, this test started failing. The reason being a zero in the denominator of the following calculation:

mom_mixrate = (0.41*9.8696)*u_star**2 / & (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)

See relevant lines of code here

There is no issue when GUST_CONST > 0 because that enforces u_star to be > 0. However, since the corner point is on the Equator absf is problematic and perhaps must be bounded:

absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))

@alperaltuntas has come up with a possible way of avoiding this issue, but it would be good to get @Hallberg-NOAA's option before we send a PR.

@Hallberg-NOAA
Copy link
Collaborator

At first glance this is a good candidate for using Adcroft's rule of reciprocals, namely
if (x==0) then Inv_x = 0.0 else Inv_x = 1./x.
Alternately we could just trap the special case where absf = 0 in the mom_mixrate expression on this line, in which a mathematical cancellation of terms gives
if (absf == 0) mom_mixrate = (0.41*9.8696) * u_star / (4.0*(h_vel+dz_neglect)),
but this would just lead to problems a with an infinite mixing growth timescale a few lines later. The code subsequently applies limits on the restratifying streamfunctions to avoid any CFL violations, so a very large (but not infinite) mixing growth timescale might give plausible answers.

Another approach might be to take the viscous limit on the mixing rate as floor, namely something like
mom_mixrate_min = pi2 * Kv / (h_vel + dz_neglect)**2,
where pi2 ~= 9.8696 and Kv is a (molecular or turbulent?) viscosity. This would probably still give such long timescales that the CFL limits would still be what really determined the streamfunctions in the end, but at least there would be a physical explanation for this timescale.

@gustavo-marques
Copy link
Collaborator Author

@alperaltuntas, can you please post here your initial fix to this issue so we can compare all the alternatives?

@Hallberg-NOAA
Copy link
Collaborator

@gustavo-marques and @alperaltuntas, do you have an update on this issue?

@alperaltuntas
Copy link
Collaborator

@Hallberg-NOAA, we haven't implemented a fix for this yet. IIRC, the solution I temporarily tested was to skip the computations on land cells, although I am not sure which mask2d is an appropriate test here. Again, this issue occurs only when ustar==0.0, so our current temporary fix is to run with small GUST_CONST when we have to( i.e., when any grid cell corners coincide with the Equator). Also, I am happy to test any of the solutions you proposed. if I interpret your first solution proposal correctly, the fix is to (effectively) set mom_mixrate and timescale to zero when ustar==0.0 and absf==0. If so, I can implement and test it.

@sanAkel sanAkel assigned sanAkel and unassigned sanAkel May 24, 2022
Hallberg-NOAA added a commit to Hallberg-NOAA/MOM6 that referenced this issue Dec 4, 2022
  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.
Hallberg-NOAA added a commit to Hallberg-NOAA/MOM6 that referenced this issue Dec 5, 2022
  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.
Hallberg-NOAA added a commit to Hallberg-NOAA/MOM6 that referenced this issue Dec 13, 2022
  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.
Hallberg-NOAA added a commit to Hallberg-NOAA/MOM6 that referenced this issue Dec 17, 2022
  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.
Hallberg-NOAA added a commit to Hallberg-NOAA/MOM6 that referenced this issue Dec 18, 2022
  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.
marshallward pushed a commit to NOAA-GFDL/MOM6 that referenced this issue Dec 19, 2022
  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.
@Hallberg-NOAA
Copy link
Collaborator

This issue has been addressed by NOAA-GFDL#266, which has now been merged onto the main branch of MOM6, and hence this issue can be closed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

4 participants