SET_DTBT using wrong face area

[herrwang0]

Background

Subroutine set_dtbt() is used to estimate barotropic time step size. The subroutine is called in

• barotropic solver initialization
  

  MOM6/src/core/MOM_barotropic.F90

  Line 4907 in 7305528

  call set_dtbt(G, GV, US, CS, gtot_est=gtot_estimate, SSH_add=SSH_extra)

• After a dtbt_reset_period
  

  MOM6/src/core/MOM_dynamics_split_RK2.F90

  Line 646 in 7305528

  if (calc_dtbt) call set_dtbt(G, GV, US, CS%barotropic_CSp, eta, CS%pbce)

In the subroutine, there are currently three options in calculating/estimating face area

1. With optional input argument BT_cont
2. With optional input argument eta, passing to subroutine find_face_areas()
3. With optional input argument add_SSH, passing to subroutine find_face_areas()

And all of them seem to have an issue.

Issues

1. BT_cont

For reasons unknown, BT_cont option is never used in this subroutine.

2. eta

Instead, eta is used in recalculating dtbt. The problem is that eta is almost never used because of this buggy block (NONLINEAR_BT_CONTINUITY is by default False and not used with USE_BT_CONT_TYPE):

MOM6/src/core/MOM_barotropic.F90

Lines 2869 to 2875 in 7305528

	if (use_BT_cont) then
	call BT_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, halo=0)
	elseif (CS%Nonlinear_continuity .and. present(eta)) then
	call find_face_areas(Datu, Datv, G, GV, US, CS, MS, 0, eta=eta)
	else
	call find_face_areas(Datu, Datv, G, GV, US, CS, MS, 0, add_max=add_SSH)
	endif

3. add_SSH

The bug in the if-block above effectively let the model use zero sea surface height (add_SSH is zero by default) to calculate face area, which may underestimate gravity wave speed.

Moreover, without eta, subroutine find_face_areas() will use (Z_ref-bathymetry) to estimate thickness. This does not apply to cases like the Great Lakes, where 1) Within the domain, the reference height is spatially varying 2) portions of the water bodies are above the sea level reference height . (and runtime parameter SSH_EXTRA, which is currently only used to estimate an upper bound of eta during initialization would not help.)

Solutions

1. Use BT_cont when possible
2. Fix the if-block bug related to eta
3. *Make Z_ref a 2D field, rather than a scalar. I put an asterisk because the scalar should be fine for most cases, even with spatially varying reference heights. But I wonder if something weird may come up with OBCs...

The fixes may make the model slower with shorter dtbt.