Merge branch 'dev/ncar' into mnlevy1981/add_MARBL

config_src/external/stochastic_physics/stochastic_physics.F90

@@ -16,7 +16,7 @@ module stochastic_physics
 
 !> Initializes the stochastic physics perturbations.
 subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nx, ny, nz, pert_epbl_in, do_sppt_in, &
-                                       mpiroot, mpicomm, iret)
+                                       do_skeb_in,mpiroot, mpicomm, iret)
   real,    intent(in)    :: delt !< timestep in seconds between calls to run_stochastic_physics_ocn [s]
   integer, intent(in)    :: nx   !< number of gridpoints in the x-direction of the compute grid
   integer, intent(in)    :: ny   !< number of gridpoints in the y-direction of the compute grid
@@ -25,6 +25,7 @@ subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nx, ny, nz, pert_
   real,    intent(in)    :: geoLatT(nx,ny) !< Latitude in degrees
   logical, intent(in)    :: pert_epbl_in !< logical flag, if true generate random pattern for ePBL perturbations
   logical, intent(in)    :: do_sppt_in   !< logical flag, if true generate random pattern for SPPT perturbations
+  logical, intent(in)    :: do_skeb_in   !< logical flag, if true generate random pattern for SKEB perturbations
   integer, intent(in)    :: mpiroot !< root processor
   integer, intent(in)    :: mpicomm !< mpi communicator
   integer, intent(out)   :: iret    !< return code
@@ -38,14 +39,20 @@ subroutine init_stochastic_physics_ocn(delt, geoLonT, geoLatT, nx, ny, nz, pert_
     call MOM_error(WARNING, 'init_stochastic_physics_ocn: do_sppt needs to be false if using the stub')
     iret=-1
   endif
+  if (do_skeb_in) then
+    call MOM_error(WARNING, 'init_stochastic_physics_ocn: do_skeb needs to be false if using the stub')
+    iret=-1
+  endif
 
   ! This stub function does not actually do anything.
   return
 end subroutine init_stochastic_physics_ocn
 
+
 !> Determines the stochastic physics perturbations.
-subroutine run_stochastic_physics_ocn(sppt_wts, t_rp1, t_rp2)
+subroutine run_stochastic_physics_ocn(sppt_wts, skeb_wts, t_rp1, t_rp2)
   real, intent(inout) :: sppt_wts(:,:) !< array containing random weights for SPPT range [0,2]
+  real, intent(inout) :: skeb_wts(:,:) !< array containing random weights for SKEB
   real, intent(inout) :: t_rp1(:,:)    !< array containing random weights for ePBL
                                        !! perturbations (KE generation) range [0,2]
   real, intent(inout) :: t_rp2(:,:)    !< array containing random weights for ePBL

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,23 @@ 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%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%RZ3_T3_to_W_m2 * US%L_to_Z**2)
 
 end subroutine ALE_register_diags
 
@@ -1020,7 +1040,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 +1062,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,13 +1075,34 @@ 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 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
 
   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 (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.
+
+  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,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 (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%H_to_RZ * u2h_tot * I_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)
 
@@ -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
@@ -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%H_to_RZ * v2h_tot * I_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
@@ -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()")
 

src/core/MOM.F90

@@ -64,7 +64,7 @@ module MOM
 use MOM_diabatic_driver,       only : diabatic, diabatic_driver_init, diabatic_CS, extract_diabatic_member
 use MOM_diabatic_driver,       only : adiabatic, adiabatic_driver_init, diabatic_driver_end
 use MOM_diabatic_driver,       only : register_diabatic_restarts
-use MOM_stochastics,           only : stochastics_init, update_stochastics, stochastic_CS
+use MOM_stochastics,           only : stochastics_init, update_stochastics, stochastic_CS, apply_skeb
 use MOM_diagnostics,           only : calculate_diagnostic_fields, MOM_diagnostics_init
 use MOM_diagnostics,           only : register_transport_diags, post_transport_diagnostics
 use MOM_diagnostics,           only : register_surface_diags, write_static_fields
@@ -741,7 +741,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
     endif
   endif
   ! advance the random pattern if stochastic physics is active
-  if (CS%stoch_CS%do_sppt .OR. CS%stoch_CS%pert_epbl) call update_stochastics(CS%stoch_CS)
+  if (CS%stoch_CS%do_sppt .OR. CS%stoch_CS%pert_epbl .OR. CS%stoch_CS%do_skeb) &
+    call update_stochastics(CS%stoch_CS)
 
   if (do_dyn) then
     if (nonblocking_p_surf_update) &
@@ -1151,7 +1152,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
       if (CS%VarMix%use_variable_mixing) &
         call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC)
       call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt_thermo, G, GV, US, &
-                             CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp)
+                             CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, &
+                             CS%stoch_CS)
       call cpu_clock_end(id_clock_thick_diff)
       call pass_var(h, G%Domain, clock=id_clock_pass, halo=max(2,CS%cont_stencil))
       if (showCallTree) call callTree_waypoint("finished thickness_diffuse_first (step_MOM)")
@@ -1223,7 +1225,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
     call step_MOM_dyn_split_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, &
                 p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, &
                 CS%eta_av_bc, G, GV, US, CS%dyn_split_RK2_CSp, calc_dtbt, CS%VarMix, &
-                CS%MEKE, CS%thickness_diffuse_CSp, CS%pbv, waves=waves)
+                CS%MEKE, CS%thickness_diffuse_CSp, CS%pbv, CS%stoch_CS, waves=waves)
     if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_split (step_MOM)")
 
   elseif (CS%do_dynamics) then ! ------------------------------------ not SPLIT
@@ -1237,11 +1239,13 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
     if (CS%use_RK2) then
       call step_MOM_dyn_unsplit_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, &
                p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, &
-               CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, CS%pbv)
+               CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, CS%pbv, &
+               CS%stoch_CS)
     else
       call step_MOM_dyn_unsplit(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, &
                p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, &
-               CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, CS%pbv, Waves=Waves)
+               CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, CS%pbv, &
+               CS%stoch_CS, Waves=Waves)
     endif
     if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)")
 
@@ -1282,7 +1286,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
       if (CS%VarMix%use_variable_mixing) &
         call calc_slope_functions(h, CS%tv, dt, G, GV, US, CS%VarMix, OBC=CS%OBC)
       call thickness_diffuse(h, CS%uhtr, CS%vhtr, CS%tv, dt, G, GV, US, &
-                             CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp)
+                             CS%MEKE, CS%VarMix, CS%CDp, CS%thickness_diffuse_CSp, CS%stoch_CS)
 
       if (CS%debug) call hchksum(h,"Post-thickness_diffuse h", G%HI, haloshift=1, scale=GV%H_to_MKS)
       call cpu_clock_end(id_clock_thick_diff)
@@ -1584,6 +1588,10 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
                   Time_end_thermo, G, GV, US, CS%diabatic_CSp, CS%stoch_CS, CS%OBC, Waves)
     fluxes%fluxes_used = .true.
 
+    if (CS%stoch_CS%do_skeb) then
+       call apply_skeb(CS%G,CS%GV,CS%stoch_CS,CS%u,CS%v,CS%h,CS%tv,dtdia,Time_end_thermo)
+    endif
+
     if (showCallTree) call callTree_waypoint("finished diabatic (step_MOM_thermo)")
 
     ! Regridding/remapping is done here, at end of thermodynamics time step
@@ -1641,7 +1649,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.