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

Imperfect restarts in Bering domain #1130

Closed
kshedstrom opened this issue Jun 10, 2020 · 15 comments
Closed

Imperfect restarts in Bering domain #1130

kshedstrom opened this issue Jun 10, 2020 · 15 comments

Comments

@kshedstrom
Copy link
Collaborator

By the end of the routine set_BBL_TKE in MOM_set_diffusivity.F90, some values of visc%DKE_BBL don't match between the restart and starting from the beginning. There are mismatches between v2_bbl coming from mismatches in visc%bbl_thick_v, but I haven't figured out where they came from yet. The latter are sometimes only visible in hexadecimal mode in totalview. Surprisingly, bbl_thick_[uv] doesn't report differences in the debug output.

@kshedstrom
Copy link
Collaborator Author

u and v differ coming into step_MOM right from the get go at t=3600.

@adcroft
Copy link
Collaborator

adcroft commented Jun 11, 2020

To clarify, check sums for u,v after restart match but by the end of the first step they don't?

@kshedstrom
Copy link
Collaborator Author

Sorry, u and v do match at time 3600. What I learned is that to do this in debug mode, you have to make your restart file in debug mode too. bbl_thick_v now matches as well, but still not TKE_BBL right at j=4 (southern edge). v2_bbl at j=4 matches for some values of i, not others.

@kshedstrom
Copy link
Collaborator Author

This could be due to that h outside issue we were talking about. What is the best fix for setting h to zero outside on restarts? Or do we fix it by adding OBC code to not use h outside? Both?

@kshedstrom
Copy link
Collaborator Author

Honestly, h outside got set to sensible values in the non-restart case, but was still 1e-10 in the restart case.

@Hallberg-NOAA
Copy link
Collaborator

Hallberg-NOAA commented Jun 12, 2020 via email

@kshedstrom
Copy link
Collaborator Author

Hi Bob and Alistair,

I did have DEFAULT_2018_ANSWERS already set to False.
Fixing the instabilities would be great.

The use of outside h starting at line 1711 of MOM_set_diffusivity.F90 needs fixing. I don't know if that is coming in your pull request or if I should fuss with it tomorrow. I was surprised to find h values of order 10, 20 outside the domain because I thought we had agreed not to set those things away from 1.e-10 or zero or what have you.

@kshedstrom
Copy link
Collaborator Author

PR #1136 pushes worrisome changes down to:

  u-point: mean=  -6.0560599528976952E-02 min=  -1.0000000000000000E+00 max=   1.0000000000000000E+00 u calc_Visbeck_coeffs slope_[xy]
  u-point: c= 112070200 sw= 112099677 se= 111941954 nw= 111986777 ne= 111832860 u calc_Visbeck_coeffs slope_[xy]
  v-point: mean=  -6.9146776804563954E-03 min=  -1.0000000000000000E+00 max=   1.0000000000000000E+00 v calc_Visbeck_coeffs slope_[xy]
  v-point: c= 112102125 sw= 112125928 se= 112059106 nw= 111914938 ne= 111852046 v calc_Visbeck_coeffs slope_[xy]
  u-point: mean=   2.1285611217717305E-05 min=  -1.4786217361362181E-04 max=   6.3638164433869110E-03 u calc_Visbeck_coeffs N2_u, N2_v
  u-point: c=  71412479 u calc_Visbeck_coeffs N2_u, N2_v
  v-point: mean=   2.1196614671105384E-05 min=  -1.4781651039840000E-04 max=   5.8967828118764573E-03 v calc_Visbeck_coeffs N2_u, N2_v
  v-point: c=  71295206 v calc_Visbeck_coeffs N2_u, N2_v
  u-point: mean=   4.7432796016329614E-06 min=   0.0000000000000000E+00 max=   1.9925280330821771E-03 u calc_Visbeck_coeffs SN_[uv]
  u-point: c=   1479060 u calc_Visbeck_coeffs SN_[uv]
  v-point: mean=   4.9393701975289042E-06 min=   0.0000000000000000E+00 max=   1.6428482245135264E-02 v calc_Visbeck_coeffs SN_[uv]

Back to the debugger.

@Hallberg-NOAA
Copy link
Collaborator

Based on this new output, I think that the problem is that calc_isoneutral_slopes is not aware of the faces that have open boundary conditions. For lack anything better to suggest, I think that calc_isoneutral_slopes should return a slope of 0 at OBC velocity points. (For simplicity, I would set the values as is currently done, but then zero them out at the end of the routine.) N2_u and N2_v are also set to be the stratification extrapolated from tracer to velocity points in calc_isoneutral_slopes, so these might need to be dealt with as well, but if the slopes are set to zero at OBC points, it might not matter what the stratification is, since it is the product of the two that are actually being used in the parameterizations.

There may also be a similar issue with calc_slope_functions_using_just_e (in MOM_lateral_mixing_coeffs.F90), and then we probably also want to set CS%SN_[uv] to 0 at OBC points in calc_Visbeck_coeffs.

I suspect that this knot of changes is going to require passing an OBC type into calc_slope_functions, and then on to the other 3 subroutines.

Given that you have the tests that are detecting these issues, @kshedstrom, how would you like to proceed?

@kshedstrom
Copy link
Collaborator Author

I had already started on this and will get back to you.

Speaking of h outside, do you know why T and S are being set outside as well? They too don't match on restart, but I'm not sure they are being used.

@Hallberg-NOAA
Copy link
Collaborator

I believe that the temperatures and salinity values outside of the domain are being used in calculations that are then masked by multiplication by a land-mask field or (in the case of OBCs) should be replaced. This works provided that there are no NaNs outside of the domain, because unfortunately 0*NaN = NaN.

@kshedstrom
Copy link
Collaborator Author

The values of N2_u and N2_v don't match yet. Are they being used? How about the values of slope_x, slope_y one grid outside the domain? Those don't match either. Something leads to differences.

@kshedstrom
Copy link
Collaborator Author

All we need now are halo updates of OBC%tres_x and/or OBC%tres_y. Adding both at once as a vector update causes trouble for dumbbell because it only has one.

@Hallberg-NOAA
Copy link
Collaborator

There are two viable near-term solutions that I could see for the problems you are encountering when doing the halo updates on the open boundary condition tracer reservoirs and similar fields. The first and simplest solution would be to add both of these arrays to do the vector-pair halo update, even if there are only open boundary conditions on X- or Y- faces. This is a little wasteful, but it has the advantage of working even if we move to grids with topologies that mix up the two components, like a cubed-sphere. Also, it can be implemented now in the MOM_open_boundary.F90 without any changes anywhere else.

The second near-term solution would be to extend the code in pass_var to work with variables that are staggered on the cell faces, instead of just the cell centers or corners. The underlying FMS infrastructure supports this, so it should be a minor extension to the MOM6 framework interfaces. We could then do halo updates separately on the two faces if the other does not exist, or together if both are there. I will put this benign extension of the MOM6 framework code on my to-do list, unless someone else steps up and does it first.

A third, longer term option would be to revisit the need to use global arrays to do halo updates on open boundary segment data, but that would be a very involved exercise and would go deeply into the communication code in the FMS infrastructure, and I do not see this happening any time soon.

@kshedstrom
Copy link
Collaborator Author

Closing this now with PR #1159

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

No branches or pull requests

3 participants