From bee7499cb8b77af9dd01221bf19362f3a3be1e28 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Wed, 1 May 2024 15:14:43 -0600 Subject: [PATCH] Limit KE-conserving correction This commit does two main things. - Limit the magnitude of the multiplicative correction applied to the baroclinic velocity to +25%. This prevents rare occasions where the correction creates very large baroclinic velocities. - Move the diagnostic of KE loss/gain from before the correction to after the correction. Without the limit (above) the correction is exact to machine precision, so there was no point in computing it after the correction. With the new limit it makes sense to compute the diagnostic after the correction. --- src/ALE/MOM_ALE.F90 | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index f2d7bfe02f..6c9b7cca25 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -97,7 +97,7 @@ 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 + 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. @@ -1115,7 +1115,7 @@ 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,dz,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 @@ -1136,13 +1136,6 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u call remapping_core_h(CS%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt, & h_neglect, h_neglect_edge) - 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 - endif - if (variance_option .and. CS%conserve_ke) then ! Conserve ke_u by correcting baroclinic component. ! Assumes total depth doesn't change during remap, and @@ -1161,12 +1154,19 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u 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 = sqrt(ke_c_src / (ke_c_tgt + 1.E-19)) + 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 + 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) @@ -1183,7 +1183,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u ! --- Remap v profiles from the source vertical grid onto the new target grid. - !$OMP parallel do default(shared) private(h1,h2,v_src,dz,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 @@ -1204,13 +1204,6 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u call remapping_core_h(CS%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt, & h_neglect, h_neglect_edge) - 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 * GV%H_to_m * v2h_tot / dt - endif - if (variance_option .and. CS%conserve_ke) then ! Conserve ke_v by correcting baroclinic component. ! Assumes total depth doesn't change during remap, and @@ -1229,12 +1222,19 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u 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 = sqrt(ke_c_src / (ke_c_tgt + 1.E-19)) + 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 * GV%H_to_m * v2h_tot / dt + 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