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

Jump to bottom

KE-conserving correction to velocity remap #277

Merged
merged 4 commits into from
May 24, 2024
Merged

KE-conserving correction to velocity remap #277

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d7b110f
KE-conserving Remap Correction
iangrooms Apr 24, 2024
bee7499
Limit KE-conserving correction
iangrooms May 1, 2024
fb2b433
Fix dimensional scaling error
iangrooms May 2, 2024
b182900
Correct Units
iangrooms May 6, 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
139 changes: 134 additions & 5 deletions src/ALE/MOM_ALE.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,9 @@ module MOM_ALE
!! values result in the use of more robust and accurate forms of
!! mathematically equivalent expressions.

logical :: conserve_ke !< Apply a correction to the baroclinic velocity after remapping to
!! conserve KE.

logical :: debug !< If true, write verbose checksums for debugging purposes.
logical :: show_call_tree !< For debugging

Expand All @@ -117,6 +120,8 @@ module MOM_ALE
integer :: id_e_preale = -1 !< diagnostic id for interface heights before ALE.
integer :: id_vert_remap_h = -1 !< diagnostic id for layer thicknesses used for remapping
integer :: id_vert_remap_h_tendency = -1 !< diagnostic id for layer thickness tendency due to ALE
integer :: id_remap_delta_integ_u2 = -1 !< Change in depth-integrated rho0*u**2/2
integer :: id_remap_delta_integ_v2 = -1 !< Change in depth-integrated rho0*v**2/2

end type

Expand Down Expand Up @@ -298,6 +303,11 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
if (CS%use_hybgen_unmix) &
call init_hybgen_unmix(CS%hybgen_unmixCS, GV, US, param_file, hybgen_regridCS)

call get_param(param_file, mdl, "REMAP_VEL_CONSERVE_KE", CS%conserve_ke, &
"If true, a correction is applied to the baroclinic component of velocity "//&
"after remapping so that total KE is conserved. KE may not be conserved "//&
"when (CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)", &
default=.false.)
call get_param(param_file, "MOM", "DEBUG", CS%debug, &
"If true, write out verbose debugging data.", &
default=.false., debuggingParam=.true.)
Expand Down Expand Up @@ -341,13 +351,25 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS)

CS%id_dzRegrid = register_diag_field('ocean_model', 'dzRegrid', diag%axesTi, Time, &
'Change in interface height due to ALE regridding', 'm', conversion=GV%H_to_m)
cs%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', diag%axestl, Time, &
CS%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', diag%axestl, Time, &
'layer thicknesses after ALE regridding and remapping', &
thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
cs%id_vert_remap_h_tendency = register_diag_field('ocean_model', &
CS%id_vert_remap_h_tendency = register_diag_field('ocean_model', &
'vert_remap_h_tendency', diag%axestl, Time, &
'Layer thicknesses tendency due to ALE regridding and remapping', &
trim(thickness_units)//" s-1", conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.)
CS%id_remap_delta_integ_u2 = register_diag_field('ocean_model', 'ale_u2', diag%axesCu1, Time, &
'Rate of change in half rho0 times depth integral of squared zonal'//&
' velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//&
' this measures the change before the KE-conserving correction is'//&
' applied.', 'W m-2', conversion=US%R_to_kg_m3*GV%H_to_m*US%L_T_to_m_s**2 &
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the changes suggested below are made, these conversion factors would become conversion=US%R_to_kg_m3*US%Z_to_m*US%L_T_to_m_s**2*US%s_to_T or equivalently conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2. The same comment also applies to the register_diag_field() call for CS%id_remap_delta_integ_v2.

Also, it would be more convenient for verifying the consistency between the conversion and units arguments could be kept entirely on the same line and not split the conversion factor with line continuation. In this case, this could easily be accommodated within the 100 character recommended line length by moving 'applied .` up to the preceding line. Also note that the spaces before and after the
continuation between lines 369 and 370 will lead to a double space that we might not want.

* US%s_to_T)
CS%id_remap_delta_integ_v2 = register_diag_field('ocean_model', 'ale_v2', diag%axesCv1, Time, &
'Rate of change in half rho0 times depth integral of squared meridional'//&
' velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//&
' this measures the change before the KE-conserving correction is'//&
' applied.', 'W m-2', conversion=US%R_to_kg_m3*GV%H_to_m*US%L_T_to_m_s**2 &
* US%s_to_T)

end subroutine ALE_register_diags

Expand Down Expand Up @@ -1020,7 +1042,8 @@ end subroutine ALE_remap_set_h_vel_OBC
!! This routine may be called during initialization of the model at time=0, to
!! remap initial conditions to the model grid. It is also called during a
!! time step to update the state.
subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, debug)
subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, debug, &
dt, allow_preserve_variance)
type(ALE_CS), intent(in) :: CS !< ALE control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
Expand All @@ -1041,6 +1064,9 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
logical, optional, intent(in) :: debug !< If true, show the call tree
real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s]
logical, optional, intent(in) :: allow_preserve_variance !< If true, enables ke-conserving
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this logical parameter necessary?
Looks like CS%conserve_ke is enough to determine if the correction should be applied.

!! correction

! Local variables
real :: h_mask_vel ! A depth below which the thicknesses at a velocity point are masked out [H ~> m or kg m-2]
Expand All @@ -1051,13 +1077,32 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]
real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]
real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
real :: rescale_coef ! Factor that scales the baroclinic velocity to conserve ke [nondim]
real :: u_bt, v_bt ! Depth-averaged velocity components [L T-1 ~> m s-1]
real :: ke_c_src, ke_c_tgt ! \int [u_c or v_c]^2 dz on src and tgt grids [H L T-1 ~> m2 s-1]
real, dimension(SZIB_(G),SZJ_(G)) :: du2h_tot ! The rate of change of vertically integrated
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that du2h_tot and dv2h_tot should have units of [R Z L2 T-3 ~> W m-2], but as currently written they actually have units of R H L2 T-3 ~> W m-2 or kg2 m-5 s-3], as is reflected by the consistent conversion factors in the register_diag_field calls for these diagnostics on lines 365 and 371 above.

This is easily corrected, though, if instead of using GV%Rho0 where du2h_tot and dv2h_tot are calculated at lines 1167 and 1235, we were to use GV%H_to_RZ. In Boussinesq mode, GV%H_to_RZ = GV%Rho0*GV%H_to_m*US%m_to_Z, but in non-Boussinesq mode the thicknesses already include the masses and GV%H_to_RZ is just the (constant) rescaling factor to convert them to units of [kg m-2].

This suggested change has the added benefit of not introducing the Boussinesq reference density (GV%Rho0) into non-Boussinesq calculations. In non-Boussinesq mode, I have been working to avoid the use of GV%H_to_m or GV%H_to_Z whereever possible because they imply division by the arbitrary Boussinesq reference density. In this case the various factors would cancel out mathematically, but multiplying and then diving by 1035 kg m-3 changes answers at roundoff, so they can not be cancelled out after people start using the code without changing answers.

! 0.5 * rho0 * u**2 [R Z L2 T-3 ~> W m-2]
real, dimension(SZI_(G),SZJB_(G)) :: dv2h_tot ! The rate of change of vertically integrated
! 0.5 * rho0 * v**2 [R Z L2 T-3 ~> W m-2]
real :: u2h_tot, v2h_tot ! The vertically integrated u**2 and v**2 [H L2 T-2 ~> m3 s-2 or kg s-2]
logical :: variance_option ! Contains the value of allow_preserve_variance when present, else false
logical :: show_call_tree
integer :: i, j, k, nz

show_call_tree = .false.
if (present(debug)) show_call_tree = debug
if (show_call_tree) call callTree_enter("ALE_remap_velocities()")

! Setup related to KE conservation
variance_option = .false.
if (present(allow_preserve_variance)) variance_option=allow_preserve_variance

if (CS%id_remap_delta_integ_u2>0) du2h_tot(:,:) = 0.
if (CS%id_remap_delta_integ_v2>0) dv2h_tot(:,:) = 0.

if (((CS%id_remap_delta_integ_u2>0) .or. (CS%id_remap_delta_integ_v2>0)) .and. .not.present(dt))&
call MOM_error(FATAL, "ALE KE diagnostics requires passing dt into ALE_remap_velocities")

if (CS%answer_date >= 20190101) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
elseif (GV%Boussinesq) then
Expand All @@ -1070,7 +1115,9 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u

! --- Remap u profiles from the source vertical grid onto the new target grid.

!$OMP parallel do default(shared) private(h1,h2,u_src,h_mask_vel,u_tgt)
!$OMP parallel do default(shared) private(h1,h2,u_src,h_mask_vel,u_tgt, &
!$OMP u_bt,ke_c_src,ke_c_tgt,rescale_coef, &
!$OMP u2h_tot,v2h_tot)
do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (G%mask2dCu(I,j)>0.) then
! Make a 1-d copy of the start and final grids and the source velocity
do k=1,nz
Expand All @@ -1079,9 +1126,47 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
u_src(k) = u(I,j,k)
enddo

if (CS%id_remap_delta_integ_u2>0) then
u2h_tot = 0.
do k=1,nz
u2h_tot = u2h_tot - h1(k) * (u_src(k)**2)
enddo
endif

call remapping_core_h(CS%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt, &
h_neglect, h_neglect_edge)

if (variance_option .and. CS%conserve_ke) then
! Conserve ke_u by correcting baroclinic component.
! Assumes total depth doesn't change during remap, and
! that \int u(z) dz doesn't change during remap.
! First get barotropic component
u_bt = 0.0
do k=1,nz
u_bt = u_bt + h2(k) * u_tgt(k) ! Dimensions [H L T-1]
enddo
u_bt = u_bt / (sum(h2(1:nz)) + h_neglect) ! Dimensions return to [L T-1]
! Next get baroclinic ke = \int (u-u_bt)^2 from source and target
ke_c_src = 0.0
ke_c_tgt = 0.0
do k=1,nz
ke_c_src = ke_c_src + h1(k) * (u_src(k) - u_bt)**2
ke_c_tgt = ke_c_tgt + h2(k) * (u_tgt(k) - u_bt)**2
enddo
! Next rescale baroclinic component on target grid to conserve ke
rescale_coef = min(1.25, sqrt(ke_c_src / (ke_c_tgt + 1.E-19)))
do k=1,nz
u_tgt(k) = u_bt + rescale_coef * (u_tgt(k) - u_bt)
enddo
endif

if (CS%id_remap_delta_integ_u2>0) then
do k=1,nz
u2h_tot = u2h_tot + h2(k) * (u_tgt(k)**2)
enddo
du2h_tot(I,j) = GV%Rho0 * u2h_tot / dt
Copy link

@Hallberg-NOAA Hallberg-NOAA May 5, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please revise this expression from du2h_tot(I,j) = GV%Rho0 * u2h_tot / dt to du2h_tot(I,j) = GV%H_to_RZ * u2h_tot * I_dt with the line if (present(dt)) I_dt = 1.0 / dt added somewhere outside the i- and j- loops for this routine.

Changing out GV%Rho0 for this rescaling factor will place du2h_tot into the documented units of [R Z L2 T-3 ~> W m-2], and it avoids using the Boussinesq reference density in non-Boussinesq cases.

An optimizing compiler will replace the divisions by dt with a multiplication by its reciprocal but without optimizations this won't happen. If we explicitly do this replacement ourselves, we help to ensure that the optimized and non-optimized (debugging) versions of the code are closer to giving the same answers.

endif

if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) &
call mask_near_bottom_vel(u_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz)

Expand All @@ -1091,12 +1176,16 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
enddo !k
endif ; enddo ; enddo

if (CS%id_remap_delta_integ_u2>0) call post_data(CS%id_remap_delta_integ_u2, du2h_tot, CS%diag)

if (show_call_tree) call callTree_waypoint("u remapped (ALE_remap_velocities)")


! --- Remap v profiles from the source vertical grid onto the new target grid.

!$OMP parallel do default(shared) private(h1,h2,v_src,h_mask_vel,v_tgt)
!$OMP parallel do default(shared) private(h1,h2,v_src,h_mask_vel,v_tgt, &
!$OMP v_bt,ke_c_src,ke_c_tgt,rescale_coef, &
!$OMP u2h_tot,v2h_tot)
do J=G%JscB,G%JecB ; do i=G%isc,G%iec ; if (G%mask2dCv(i,J)>0.) then

do k=1,nz
Expand All @@ -1105,9 +1194,47 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
v_src(k) = v(i,J,k)
enddo

if (CS%id_remap_delta_integ_v2>0) then
v2h_tot = 0.
do k=1,nz
v2h_tot = v2h_tot - h1(k) * (v_src(k)**2)
enddo
endif

call remapping_core_h(CS%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt, &
h_neglect, h_neglect_edge)

if (variance_option .and. CS%conserve_ke) then
! Conserve ke_v by correcting baroclinic component.
! Assumes total depth doesn't change during remap, and
! that \int v(z) dz doesn't change during remap.
! First get barotropic component
v_bt = 0.0
do k=1,nz
v_bt = v_bt + h2(k) * v_tgt(k) ! Dimensions [H L T-1]
enddo
v_bt = v_bt / (sum(h2(1:nz)) + h_neglect) ! Dimensions return to [L T-1]
! Next get baroclinic ke = \int (u-u_bt)^2 from source and target
ke_c_src = 0.0
ke_c_tgt = 0.0
do k=1,nz
ke_c_src = ke_c_src + h1(k) * (v_src(k) - v_bt)**2
ke_c_tgt = ke_c_tgt + h2(k) * (v_tgt(k) - v_bt)**2
enddo
! Next rescale baroclinic component on target grid to conserve ke
rescale_coef = min(1.25, sqrt(ke_c_src / (ke_c_tgt + 1.E-19)))
do k=1,nz
v_tgt(k) = v_bt + rescale_coef * (v_tgt(k) - v_bt)
enddo
endif

if (CS%id_remap_delta_integ_v2>0) then
do k=1,nz
v2h_tot = v2h_tot + h2(k) * (v_tgt(k)**2)
enddo
dv2h_tot(I,j) = GV%Rho0 * v2h_tot / dt
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please revise this expression from dv2h_tot(i,J) = GV%Rho0 * v2h_tot / dt to dv2h_tot(i,J) = GV%H_to_RZ * v2h_tot * I_dt, for the reasons discussed above with the corresponding suggested changes to du2h_tot.

endif

if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) then
call mask_near_bottom_vel(v_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz)
endif
Expand All @@ -1118,6 +1245,8 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u
enddo !k
endif ; enddo ; enddo

if (CS%id_remap_delta_integ_v2>0) call post_data(CS%id_remap_delta_integ_v2, dv2h_tot, CS%diag)

if (show_call_tree) call callTree_waypoint("v remapped (ALE_remap_velocities)")
if (show_call_tree) call callTree_leave("ALE_remap_velocities()")

Expand Down
3 changes: 2 additions & 1 deletion src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1641,7 +1641,8 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
endif

! Remap the velocity components.
call ALE_remap_velocities(CS%ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, showCallTree)
call ALE_remap_velocities(CS%ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, showCallTree, &
dtdia, allow_preserve_variance=.true.)

if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.

Expand Down
Loading