From d7b110f16ae76308972363dd7403c6d3db4b854a Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Wed, 24 Apr 2024 15:04:06 -0600 Subject: [PATCH 1/4] KE-conserving Remap Correction This commit introduces a method that corrects the remapped velocity so that it conserves KE. The correction is activated by setting `REMAP_VEL_CONSERVE_KE = True` The commit also introduces two new diagnostics: `ale_u2` and `ale_v2` These track the change in depth-integrated KE of the u and v components of velocity before the correction is applied. They can be used even if the remapping correction is not turned on. --- src/ALE/MOM_ALE.F90 | 139 ++++++++++++++++++++++++++++++++++++++++++-- src/core/MOM.F90 | 3 +- 2 files changed, 136 insertions(+), 6 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 77ee1192a2..f2d7bfe02f 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -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 @@ -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 @@ -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.) @@ -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 & + * 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 @@ -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 @@ -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 + !! 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] @@ -1051,6 +1077,15 @@ 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 + ! 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 @@ -1058,6 +1093,16 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u 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 @@ -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,dz,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 @@ -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 (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 + ! 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 = 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%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) @@ -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,dz,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 @@ -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 (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 + ! 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 = 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%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 @@ -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()") diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index b7f8bd3f66..2c3e6b3589 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -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. From bee7499cb8b77af9dd01221bf19362f3a3be1e28 Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Wed, 1 May 2024 15:14:43 -0600 Subject: [PATCH 2/4] 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 From fb2b433213c6cfb6af86478545e217ec23bf0a1d Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Thu, 2 May 2024 09:39:54 -0600 Subject: [PATCH 3/4] Fix dimensional scaling error --- src/ALE/MOM_ALE.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 6c9b7cca25..56dbc95c17 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -1232,7 +1232,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u 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 + dv2h_tot(I,j) = GV%Rho0 * v2h_tot / dt endif if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) then From b182900e7c2422ed527dd5aefaa69f6700a0f91e Mon Sep 17 00:00:00 2001 From: Ian Grooms Date: Mon, 6 May 2024 08:06:04 -0600 Subject: [PATCH 4/4] Correct Units This commit addresses @Hallberg-NOAA's comments on [the PR](https://github.com/NCAR/MOM6/pull/277). Computations of `ale_u2` and `ale_v2` are updated to work correctly in both Boussinesq and non-Boussinesq modes. --- src/ALE/MOM_ALE.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 56dbc95c17..600439d5b2 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -361,15 +361,13 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS) 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 & - * US%s_to_T) + ' this measures the change before the KE-conserving correction is applied.', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2 * US%L_to_Z**2) 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) + ' this measures the change before the KE-conserving correction is applied.', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2 * US%L_to_Z**2) end subroutine ALE_register_diags @@ -1079,12 +1077,13 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u 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 :: ke_c_src, ke_c_tgt ! \int [u_c or v_c]^2 dz on src and tgt grids [H L2 T-2 ~> m3 s-2] real, dimension(SZIB_(G),SZJ_(G)) :: du2h_tot ! The rate of change of vertically integrated ! 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] + real :: I_dt ! 1 / dt [T-1 ~> s-1] logical :: variance_option ! Contains the value of allow_preserve_variance when present, else false logical :: show_call_tree integer :: i, j, k, nz @@ -1096,6 +1095,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u ! Setup related to KE conservation variance_option = .false. if (present(allow_preserve_variance)) variance_option=allow_preserve_variance + if (present(dt)) I_dt = 1.0 / dt if (CS%id_remap_delta_integ_u2>0) du2h_tot(:,:) = 0. if (CS%id_remap_delta_integ_v2>0) dv2h_tot(:,:) = 0. @@ -1164,7 +1164,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u do k=1,nz u2h_tot = u2h_tot + h2(k) * (u_tgt(k)**2) enddo - du2h_tot(I,j) = GV%Rho0 * u2h_tot / dt + du2h_tot(I,j) = GV%H_to_RZ * u2h_tot * I_dt endif if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) & @@ -1232,7 +1232,7 @@ subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u do k=1,nz v2h_tot = v2h_tot + h2(k) * (v_tgt(k)**2) enddo - dv2h_tot(I,j) = GV%Rho0 * v2h_tot / dt + dv2h_tot(I,j) = GV%H_to_RZ * v2h_tot * I_dt endif if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) then