Skip to content

Add option to use BT_CONT to calculate dtbt #823

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.

Sign up for GitHub

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

Merged
merged 3 commits into from
Feb 13, 2025
Merged

Add option to use BT_CONT to calculate dtbt #823

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ced594d
Reorder optional input arguments in set_dtbt
herrwang0 Feb 2, 2025
59cd564
Add runtime parameter SET_DTBT_USE_BT_CONT
herrwang0 Jan 29, 2025
bd9e9fa
Minor fixes on setting dtbt in barotropic_init
herrwang0 Feb 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 26 additions & 26 deletions src/core/MOM_barotropic.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2847,26 +2847,26 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,

end subroutine btstep

!> This subroutine automatically determines an optimal value for dtbt based
!! on some state of the ocean.
subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
!> This subroutine automatically determines an optimal value for dtbt based on some state of the ocean. Either pbce or
!! gtot_est is required to calculate gravitational acceleration. Column thickness can be estimated using BT_cont, eta,
!! and SSH_add (default=0), with priority given in that order. The subroutine sets CS%dtbt_max and CS%dtbt.
subroutine set_dtbt(G, GV, US, CS, pbce, gtot_est, BT_cont, eta, SSH_add)
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure
real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta !< The barotropic free surface
!! height anomaly or column mass anomaly [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(in) :: pbce !< The baroclinic pressure
!! anomaly in each layer due to free surface
!! height anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
type(BT_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe
!! the effective open face areas as a
!! function of barotropic flow.
real, optional, intent(in) :: gtot_est !< An estimate of the total gravitational
!! acceleration [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
real, optional, intent(in) :: SSH_add !< An additional contribution to SSH to
!! provide a margin of error when
!! calculating the external wave speed [Z ~> m].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
optional, intent(in) :: pbce !< The baroclinic pressure anomaly in each layer due to free
!! surface height anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
real, optional, intent(in) :: gtot_est !< An estimate of the total gravitational acceleration
!! [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
type(BT_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe the effective open
!! face areas as a function of barotropic flow.
real, dimension(SZI_(G),SZJ_(G)), &
optional, intent(in) :: eta !< The barotropic free surface height anomaly or column mass
!! anomaly [H ~> m or kg m-2].
real, optional, intent(in) :: SSH_add !< An additional contribution to SSH to provide a margin of
!! error when calculating the external wave speed [Z ~> m].

! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: &
Expand Down Expand Up @@ -4424,7 +4424,7 @@ end subroutine bt_mass_source
!> barotropic_init initializes a number of time-invariant fields used in the
!! barotropic calculation and initializes any barotropic fields that have not
!! already been initialized.
subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, &
subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
restart_CS, calc_dtbt, BT_cont, SAL_CSp, HA_CSp)
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
Expand All @@ -4435,9 +4435,6 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G)), &
intent(in) :: eta !< Free surface height or column mass anomaly
!! [Z ~> m] or [H ~> kg m-2].
type(time_type), target, intent(in) :: Time !< The current model time.
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
Expand Down Expand Up @@ -4465,7 +4462,7 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
real :: SSH_extra ! An estimate of how much higher SSH might get, for use
! in calculating the safe external wave speed [Z ~> m].
real :: dtbt_input ! The input value of DTBT, [nondim] if negative or [s] if positive.
real :: dtbt_tmp ! A temporary copy of CS%dtbt read from a restart file [T ~> s]
real :: dtbt_restart ! A temporary copy of CS%dtbt read from a restart file [T ~> s]
real :: wave_drag_scale ! A scaling factor for the barotropic linear wave drag
! piston velocities [nondim].
character(len=200) :: inputdir ! The directory in which to find input files.
Expand Down Expand Up @@ -4978,9 +4975,9 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,

CS%dtbt_fraction = 0.98 ; if (dtbt_input < 0.0) CS%dtbt_fraction = -dtbt_input

dtbt_tmp = -1.0
dtbt_restart = -1.0
if (query_initialized(CS%dtbt, "DTBT", restart_CS)) then
dtbt_tmp = CS%dtbt
dtbt_restart = CS%dtbt
endif

! Estimate the maximum stable barotropic time step.
Expand All @@ -4991,14 +4988,17 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
H_to_Z = GV%H_to_RZ / CS%Rho_BT_lin
do k=1,GV%ke ; gtot_estimate = gtot_estimate + H_to_Z*GV%g_prime(K) ; enddo
endif

! CS%dtbt calculated here by set_dtbt is only used when dtbt is not reset during the run, i.e. DTBT_RESET_PERIOD<0.
call set_dtbt(G, GV, US, CS, gtot_est=gtot_estimate, SSH_add=SSH_extra)

if (dtbt_input > 0.0) then
CS%dtbt = US%s_to_T * dtbt_input
elseif (dtbt_tmp > 0.0) then
CS%dtbt = dtbt_tmp
elseif (dtbt_restart > 0.0) then
CS%dtbt = dtbt_restart
endif
if ((dtbt_tmp > 0.0) .and. (dtbt_input > 0.0)) calc_dtbt = .false.

calc_dtbt = .true. ; if ((dtbt_restart > 0.0) .and. (dtbt_input > 0.0)) calc_dtbt = .false.

call log_param(param_file, mdl, "DTBT as used", CS%dtbt, units="s", unscale=US%T_to_s)
call log_param(param_file, mdl, "estimated maximum DTBT", CS%dtbt_max, units="s", unscale=US%T_to_s)
Expand Down
18 changes: 13 additions & 5 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,8 +154,7 @@ module MOM_dynamics_split_RK2
logical :: split_bottom_stress !< If true, provide the bottom stress
!! calculated by the vertical viscosity to the
!! barotropic solver.
logical :: calc_dtbt !< If true, calculate the barotropic time-step
!! dynamically.
logical :: dtbt_use_bt_cont !< If true, use BT_cont to calculate DTBT.
logical :: store_CAu !< If true, store the Coriolis and advective accelerations at the
!! end of the timestep for use in the next predictor step.
logical :: CAu_pred_stored !< If true, the Coriolis and advective accelerations at the
Expand Down Expand Up @@ -648,7 +647,15 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
endif

call cpu_clock_begin(id_clock_btstep)
if (calc_dtbt) call set_dtbt(G, GV, US, CS%barotropic_CSp, eta, CS%pbce)
if (calc_dtbt) then
if (CS%dtbt_use_bt_cont .and. associated(CS%BT_cont)) then
call set_dtbt(G, GV, US, CS%barotropic_CSp, CS%pbce, BT_cont=CS%BT_cont)
else
! In the following call, eta is only used when NONLINEAR_BT_CONTINUITY is True. Otherwise, dtbt is effectively
! calculated with eta=0. Note that NONLINEAR_BT_CONTINUITY is False if BT_CONT is used, which is the default.
call set_dtbt(G, GV, US, CS%barotropic_CSp, CS%pbce, eta=eta)
endif
endif
if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90")
! This is the predictor step call to btstep.
! The CS%ADp argument here stores the weights for certain integrated diagnostics.
Expand Down Expand Up @@ -1410,7 +1417,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p
"If SPLIT is false and USE_RK2 is true, BEGW can be "//&
"between 0 and 0.5 to damp gravity waves.", &
units="nondim", default=0.0)

call get_param(param_file, mdl, "SET_DTBT_USE_BT_CONT", CS%dtbt_use_bt_cont, &
"If true, use BT_CONT to calculate DTBT if possible.", default=.false.)
call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", CS%split_bottom_stress, &
"If true, provide the bottom stress calculated by the "//&
"vertical viscosity to the barotropic solver.", default=.false.)
Expand Down Expand Up @@ -1545,7 +1553,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p
! Copy eta into an output array.
do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo

call barotropic_init(u, v, h, CS%eta, Time, G, GV, US, param_file, diag, &
call barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, &
CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, CS%SAL_CSp, HA_CSp)

if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. &
Expand Down
18 changes: 13 additions & 5 deletions src/core/MOM_dynamics_split_RK2b.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,8 +157,7 @@ module MOM_dynamics_split_RK2b
logical :: split_bottom_stress !< If true, provide the bottom stress
!! calculated by the vertical viscosity to the
!! barotropic solver.
logical :: calc_dtbt !< If true, calculate the barotropic time-step
!! dynamically.
logical :: dtbt_use_bt_cont !< If true, use BT_cont to calculate DTBT.
logical :: calculate_SAL !< If true, calculate self-attraction and loading.
logical :: use_tides !< If true, tidal forcing is enabled.
logical :: remap_aux !< If true, apply ALE remapping to all of the auxiliary 3-D
Expand Down Expand Up @@ -672,7 +671,15 @@ subroutine step_MOM_dyn_split_RK2b(u_av, v_av, h, tv, visc, Time_local, dt, forc
uh_ptr => uh_in ; vh_ptr => vh_in ; u_ptr => u_inst ; v_ptr => v_inst

call cpu_clock_begin(id_clock_btstep)
if (calc_dtbt) call set_dtbt(G, GV, US, CS%barotropic_CSp, eta, CS%pbce)
if (calc_dtbt) then
if (CS%dtbt_use_bt_cont .and. associated(CS%BT_cont)) then
call set_dtbt(G, GV, US, CS%barotropic_CSp, CS%pbce, BT_cont=CS%BT_cont)
else
! In the following call, eta is only used when NONLINEAR_BT_CONTINUITY is True. Otherwise, dtbt is effectively
! calculated with eta=0. Note that NONLINEAR_BT_CONTINUITY is False if BT_CONT is used, which is the default.
call set_dtbt(G, GV, US, CS%barotropic_CSp, CS%pbce, eta=eta)
endif
endif
if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90")
! This is the predictor step call to btstep.
! The CS%ADp argument here stores the weights for certain integrated diagnostics.
Expand Down Expand Up @@ -1335,7 +1342,8 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US,
"If SPLIT is false and USE_RK2 is true, BEGW can be "//&
"between 0 and 0.5 to damp gravity waves.", &
units="nondim", default=0.0)

call get_param(param_file, mdl, "SET_DTBT_USE_BT_CONT", CS%dtbt_use_bt_cont, &
"If true, use BT_CONT to calculate DTBT if possible.", default=.false.)
call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", CS%split_bottom_stress, &
"If true, provide the bottom stress calculated by the "//&
"vertical viscosity to the barotropic solver.", default=.false.)
Expand Down Expand Up @@ -1461,7 +1469,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US,
! Copy eta into an output array.
do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo

call barotropic_init(u, v, h, CS%eta, Time, G, GV, US, param_file, diag, &
call barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, &
CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, &
CS%SAL_CSp, HA_CSp)

Expand Down