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.

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

MEKE: No-division implementation of src_btm_drag #9

Merged
merged 2 commits into from
Sep 17, 2024
Merged

MEKE: No-division implementation of src_btm_drag #9

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
18ef041
MEKE: No-division implementation of src_btm_drag
marshallward Sep 6, 2024
174e8a3
MEKE: Remove redundant if-block in Strang split
marshallward Sep 8, 2024
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
95 changes: 56 additions & 39 deletions src/parameterizations/lateral/MOM_MEKE.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
src_GM, & ! The MEKE source/tendency from the thickness mixing (GM) [L2 T-3 ~> W kg-1] (= m2 s-3).
src_mom_lp, & ! The MEKE source/tendency from the Laplacian of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).
src_mom_bh, & ! The MEKE source/tendency from the biharmonic of the resolved flow [L2 T-3 ~> W kg-1] (= m2 s-3).
ldamping_Strang1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1].
damp_rate_s1, & ! The MEKE damping rate computed at the 1st Strang splitting stage [T-1 ~> s-1].
MEKE_current, & ! A copy of MEKE for use in computing the MEKE damping [L2 T-2 ~> m2 s-2].
drag_rate_visc, & ! Near-bottom velocity contribution to bottom drag [H T-1 ~> m s-1 or kg m-2 s-1]
drag_rate, & ! The MEKE spindown timescale due to bottom drag [T-1 ~> s-1].
Expand Down Expand Up @@ -233,11 +233,11 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
real :: cdrag2 ! The square of the drag coefficient times unit conversion factors [H2 L-2 ~> nondim or kg2 m-6]
real :: advFac ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1]
real :: mass_neglect ! A negligible mass [R Z ~> kg m-2].
real :: ldamping ! The MEKE damping rate [T-1 ~> s-1].
real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate)
real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split).
real :: sfac ! A factor needed to compute damping due to Strang splitting [nondim]
real :: Isfac ! Inverse of sfac [nondim]
real :: damp_step ! Size of damping timestep relative to sdt [nondim]
real :: damp_rate ! The MEKE damping rate [T-1 ~> s-1].
real :: damping ! The net damping of a field after sdt_damp [nondim]
logical :: use_drag_rate ! Flag to indicate drag_rate is finite
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
real(kind=real32), dimension(size(MEKE%MEKE),NUM_FEATURES) :: features_array ! The array of features
Expand Down Expand Up @@ -287,7 +287,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h

! With a depth-dependent (and possibly strong) damping, it seems
! advisable to use Strang splitting between the damping and diffusion.
sdt_damp = sdt ; if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.) sdt_damp = 0.5*sdt
damp_step = 1.
if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) damp_step = 0.5
sdt_damp = sdt * damp_step

! Calculate depth integrated mass exchange if doing advection [R Z L2 ~> kg]
if (CS%MEKE_advection_factor>0.) then
Expand Down Expand Up @@ -479,19 +481,29 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
! First stage of Strang splitting
!$OMP parallel do default(shared)
do j=js,je ; do i=is,ie
ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
if (MEKE%MEKE(i,j) < 0.) ldamping = 0.
damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
if (MEKE%MEKE(i,j) < 0.) damp_rate = 0.
! notice that the above line ensures a damping only if MEKE is positive,
! while leaving MEKE unchanged if it is negative
Isfac = 1. / (1. + sdt_damp * ldamping)
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
ldamping_Strang1(i,j) = ldamping
src_GM(i,j) = src_GM(i,j) * Isfac
src_mom_lp(i,j) = src_mom_lp(i,j) * Isfac
src_mom_bh(i,j) = src_mom_bh(i,j) * Isfac
sfac = ( 1.0 + sdt_damp*ldamping )
src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) )

damping = 1. / (1. + sdt_damp * damp_rate)

! NOTE: MEKE%MEKE should use `damping` but we must preserve the existing
! expression for bit reproducibility
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1. + sdt_damp * damp_rate)
MEKE_decay(i,j) = damp_rate * G%mask2dT(i,j)

src_GM(i,j) = src_GM(i,j) * damping
src_mom_lp(i,j) = src_mom_lp(i,j) * damping
src_mom_bh(i,j) = src_mom_bh(i,j) * damping

src_btm_drag(i,j) = - MEKE_current(i,j) * ( &
damp_step * (damp_rate * damping) &
)

! Store the effective damping rate if sdt is split
if (CS%MEKE_KH >= 0. .or. CS%MEKE_K4 >= 0.) &
damp_rate_s1(i,j) = damp_rate * damping
enddo ; enddo

if (CS%kh_flux_enabled .or. CS%MEKE_K4 >= 0.0) then
Expand Down Expand Up @@ -647,33 +659,38 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h

! Second stage of Strang splitting
if (CS%MEKE_KH >= 0.0 .or. CS%MEKE_K4 >= 0.0) then
if (sdt>sdt_damp) then
! Recalculate the drag rate, since MEKE has changed.
if (use_drag_rate) then
!$OMP parallel do default(shared)
do j=js,je ; do i=is,ie
drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
enddo ; enddo
endif
! Recalculate the drag rate, since MEKE has changed.
if (use_drag_rate) then
!$OMP parallel do default(shared)
do j=js,je ; do i=is,ie
ldamping = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
if (MEKE%MEKE(i,j) < 0.) ldamping = 0.
! notice that the above line ensures a damping only if MEKE is positive,
! while leaving MEKE unchanged if it is negative
Isfac = 1. / (1. + sdt_damp*ldamping)
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*ldamping)
MEKE_decay(i,j) = ldamping*G%mask2dT(i,j)
src_GM(i,j) = src_GM(i,j) * Isfac
src_mom_lp(i,j) = src_mom_lp(i,j) * Isfac
src_mom_bh(i,j) = src_mom_bh(i,j) * Isfac
src_adv(i,j) = src_adv(i,j) * Isfac
src_mom_K4(i,j) = src_mom_K4(i,j) * Isfac
sfac = ( 1.0 + sdt_damp*ldamping_Strang1(i,j) ) * ( 1.0 + sdt_damp*ldamping )
src_btm_drag(i,j) = MEKE_current(i,j) * ( (1.0 - sfac) / ( sdt * sfac ) )
drag_rate(i,j) = (GV%H_to_RZ * I_mass(i,j)) * sqrt( drag_rate_visc(i,j)**2 + &
cdrag2 * ( max(0.0, 2.0*bottomFac2(i,j)*MEKE%MEKE(i,j)) + CS%MEKE_Uscale**2 ) )
enddo ; enddo
endif
!$OMP parallel do default(shared)
do j=js,je ; do i=is,ie
damp_rate = CS%MEKE_damping + drag_rate(i,j) * bottomFac2(i,j)
if (MEKE%MEKE(i,j) < 0.) damp_rate = 0.
! notice that the above line ensures a damping only if MEKE is positive,
! while leaving MEKE unchanged if it is negative

damping = 1. / (1. + sdt_damp * damp_rate)

! NOTE: As above, MEKE%MEKE should use `damping` but we must preserve
! the existing expression for bit reproducibility.
MEKE%MEKE(i,j) = MEKE%MEKE(i,j) / (1.0 + sdt_damp*damp_rate)
MEKE_decay(i,j) = damp_rate*G%mask2dT(i,j)

src_GM(i,j) = src_GM(i,j) * damping
src_mom_lp(i,j) = src_mom_lp(i,j) * damping
src_mom_bh(i,j) = src_mom_bh(i,j) * damping
src_adv(i,j) = src_adv(i,j) * damping
src_mom_K4(i,j) = src_mom_K4(i,j) * damping

src_btm_drag(i,j) = -MEKE_current(i,j) * ( &
damp_step * damping * (damp_rate + damp_rate_s1(i,j)) &
)
enddo ; enddo
endif ! MEKE_KH>=0

if (CS%debug) then
Expand Down
Loading