MOM6: +(*)Dimensional consistency completion

src/core/MOM.F90

@@ -1180,7 +1180,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
   call cpu_clock_begin(id_clock_thermo)
   if (.not.CS%adiabatic) then
     if (CS%debug) then
-      call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2)
+      call uvchksum("Pre-diabatic [uv]", u, v, G%HI, haloshift=2, scale=US%L_T_to_m_s)
       call hchksum(h,"Pre-diabatic h", G%HI, haloshift=1, scale=GV%H_to_m)
       call uvchksum("Pre-diabatic [uv]h", CS%uhtr, CS%vhtr, G%HI, &
                     haloshift=0, scale=GV%H_to_m*US%L_to_m**2)

src/core/MOM_dynamics_split_RK2.F90

@@ -480,7 +480,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
   call disable_averaging(CS%diag)
 
   if (CS%debug) then
-    call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym)
+    call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
   endif
   call vertvisc_coef(up, vp, h, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC)
   call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
@@ -670,7 +670,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
 
   if (CS%debug) then
     call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym)
-    call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym)
+    call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
     call hchksum(h_av, "Predictor avg h", G%HI, haloshift=0, scale=GV%H_to_m)
   ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
     call check_redundant("Predictor up ", up, vp, G)
@@ -860,7 +860,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, &
   if (CS%id_v_BT_accel > 0) call post_data(CS%id_v_BT_accel, CS%v_accel_bt, CS%diag)
 
   if (CS%debug) then
-    call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym, vel_scale=1.0)
+    call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
     call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
     call hchksum(h_av, "Corrector avg h", G%HI, haloshift=1, scale=GV%H_to_m)
  !  call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US)

src/core/MOM_forcing_type.F90

@@ -1284,27 +1284,30 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
   ! surface mass flux maps
 
   handles%id_prcme = register_diag_field('ocean_model', 'PRCmE', diag%axesT1, Time,                  &
-        'Net surface water flux (precip+melt+lrunoff+ice calving-evap)', 'kg m-2 s-1',&
+        'Net surface water flux (precip+melt+lrunoff+ice calving-evap)', 'kg m-2 s-1', &
         standard_name='water_flux_into_sea_water', cmor_field_name='wfo',                            &
         cmor_standard_name='water_flux_into_sea_water',cmor_long_name='Water Flux Into Sea Water')
+        ! This diagnostic is rescaled to MKS units when combined.
 
-  handles%id_evap = register_diag_field('ocean_model', 'evap', diag%axesT1, Time,                         &
-       'Evaporation/condensation at ocean surface (evaporation is negative)', 'kg m-2 s-1',&
-       standard_name='water_evaporation_flux', cmor_field_name='evs',                                     &
-       cmor_standard_name='water_evaporation_flux',                                                       &
+  handles%id_evap = register_diag_field('ocean_model', 'evap', diag%axesT1, Time, &
+       'Evaporation/condensation at ocean surface (evaporation is negative)', &
+       'kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, &
+       standard_name='water_evaporation_flux', cmor_field_name='evs', &
+       cmor_standard_name='water_evaporation_flux', &
        cmor_long_name='Water Evaporation Flux Where Ice Free Ocean over Sea')
 
   ! smg: seaice_melt field requires updates to the sea ice model
   handles%id_seaice_melt = register_diag_field('ocean_model', 'seaice_melt',       &
      diag%axesT1, Time, 'water flux to ocean from snow/sea ice melting(> 0) or formation(< 0)', &
-     'kg m-2 s-1',                                                                 &
+     'kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, &
       standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',     &
       cmor_field_name='fsitherm',                                                  &
       cmor_standard_name='water_flux_into_sea_water_due_to_sea_ice_thermodynamics',&
       cmor_long_name='water flux to ocean from sea ice melt(> 0) or form(< 0)')
 
   handles%id_precip = register_diag_field('ocean_model', 'precip', diag%axesT1, Time, &
         'Liquid + frozen precipitation into ocean', 'kg m-2 s-1')
+        ! This diagnostic is rescaled to MKS units when combined.
 
   handles%id_fprec = register_diag_field('ocean_model', 'fprec', diag%axesT1, Time,     &
         'Frozen precipitation into ocean', &
@@ -1324,32 +1327,39 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles,
         units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T)
 
   handles%id_frunoff = register_diag_field('ocean_model', 'frunoff', diag%axesT1, Time,    &
-        'Frozen runoff (calving) and iceberg melt into ocean', 'kg m-2 s-1',               &
+        'Frozen runoff (calving) and iceberg melt into ocean', &
+        units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, &
         standard_name='water_flux_into_sea_water_from_icebergs',                           &
         cmor_field_name='ficeberg',                                                        &
         cmor_standard_name='water_flux_into_sea_water_from_icebergs',                      &
         cmor_long_name='Water Flux into Seawater from Icebergs')
 
   handles%id_lrunoff = register_diag_field('ocean_model', 'lrunoff', diag%axesT1, Time, &
-        'Liquid runoff (rivers) into ocean', 'kg m-2 s-1',                                    &
+        'Liquid runoff (rivers) into ocean', &
+        units='kg m-2 s-1', conversion=US%R_to_kg_m3*US%Z_to_m*US%s_to_T, &
         standard_name='water_flux_into_sea_water_from_rivers', cmor_field_name='friver',      &
         cmor_standard_name='water_flux_into_sea_water_from_rivers',                           &
         cmor_long_name='Water Flux into Sea Water From Rivers')
 
   handles%id_net_massout = register_diag_field('ocean_model', 'net_massout', diag%axesT1, Time, &
         'Net mass leaving the ocean due to evaporation, seaice formation', 'kg m-2 s-1')
+        ! This diagnostic is rescaled to MKS units when combined.
 
   handles%id_net_massin  = register_diag_field('ocean_model', 'net_massin', diag%axesT1, Time, &
         'Net mass entering ocean due to precip, runoff, ice melt', 'kg m-2 s-1')
+        ! This diagnostic is rescaled to MKS units when combined.
 
   handles%id_massout_flux = register_diag_field('ocean_model', 'massout_flux', diag%axesT1, Time, &
         'Net mass flux of freshwater out of the ocean (used in the boundary flux calculation)', &
          'kg m-2', conversion=diag%GV%H_to_kg_m2)
+        ! This diagnostic is calculated in MKS units.
 
   handles%id_massin_flux  = register_diag_field('ocean_model', 'massin_flux', diag%axesT1, Time, &
         'Net mass flux of freshwater into the ocean (used in boundary flux calculation)', 'kg m-2')
+        ! This diagnostic is calculated in MKS units.
+
   !=========================================================================
-  ! area integrated surface mass transport
+  ! area integrated surface mass transport, all are rescaled to MKS units before area integration.
 
   handles%id_total_prcme = register_scalar_field('ocean_model', 'total_PRCmE', Time, diag,         &
       long_name='Area integrated net surface water flux (precip+melt+liq runoff+ice calving-evap)',&

src/diagnostics/MOM_sum_output.F90

@@ -534,9 +534,10 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
 
       call find_eta(h, tv, G, GV, US, eta)
       do k=1,nz ; do j=js,je ; do i=is,ie
-        tmp1(i,j,k) = (eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j)
+        tmp1(i,j,k) = US%Z_to_m*(eta(i,j,K)-eta(i,j,K+1)) * areaTm(i,j)
       enddo ; enddo ; enddo
-      vol_tot = US%Z_to_m*reproducing_sum(tmp1, sums=vol_lay)
+      vol_tot = reproducing_sum(tmp1, sums=vol_lay)
+      do k=1,nz ; vol_lay(k) = US%m_to_Z * vol_lay(k) ; enddo
     else
       do k=1,nz ; do j=js,je ; do i=is,ie
         tmp1(i,j,k) = H_to_kg_m2 * h(i,j,k) * areaTm(i,j)

src/framework/MOM_spatial_means.F90

@@ -43,7 +43,7 @@ function global_area_mean(var, G, scale)
   do j=js,je ; do i=is,ie
     tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j))
   enddo ; enddo
-  global_area_mean = reproducing_sum(tmpForSumming) * (G%US%m_to_L**2 * G%IareaT_global)
+  global_area_mean = reproducing_sum(tmpForSumming) * G%IareaT_global
 
 end function global_area_mean
 
@@ -182,17 +182,20 @@ end function global_mass_integral
 
 !> Determine the global mean of a field along rows of constant i, returning it
 !! in a 1-d array using the local indexing. This uses reproducing sums.
-subroutine global_i_mean(array, i_mean, G, mask, scale)
+subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale)
   type(ocean_grid_type),            intent(inout) :: G    !< The ocean's grid structure
   real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: array  !< The variable being averaged
   real, dimension(SZJ_(G)),         intent(out)   :: i_mean !< Global mean of array along its i-axis
   real, dimension(SZI_(G),SZJ_(G)), &
-                          optional, intent(in)    :: mask !< An array used for weighting the i-mean
-  real,                   optional, intent(in)    :: scale !< A rescaling factor for the variable
+                          optional, intent(in)    :: mask  !< An array used for weighting the i-mean
+  real,                   optional, intent(in)    :: scale !< A rescaling factor for the output variable
+  real,                   optional, intent(in)    :: tmp_scale !< A rescaling factor for the internal
+                                                           !! calculations that is removed from the output
 
   ! Local variables
   type(EFP_type), allocatable, dimension(:) :: asum, mask_sum
   real :: scalefac  ! A scaling factor for the variable.
+  real :: unscale   ! A factor for undoing any internal rescaling before output.
   real :: mask_sum_r
   integer :: is, ie, js, je, idg_off, jdg_off
   integer :: i, j
@@ -201,6 +204,10 @@ subroutine global_i_mean(array, i_mean, G, mask, scale)
   idg_off = G%idg_offset ; jdg_off = G%jdg_offset
 
   scalefac = 1.0 ; if (present(scale)) scalefac = scale
+  unscale = 1.0
+  if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then
+    scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale
+  endif ; endif
   call reset_EFP_overflow_error()
 
   allocate(asum(G%jsg:G%jeg))
@@ -253,31 +260,40 @@ subroutine global_i_mean(array, i_mean, G, mask, scale)
     enddo
   endif
 
+  if (unscale /= 1.0) then ; do j=js,je ; i_mean(j) = unscale*i_mean(j) ; enddo ; endif
+
   deallocate(asum)
 
 end subroutine global_i_mean
 
 !> Determine the global mean of a field along rows of constant j, returning it
 !! in a 1-d array using the local indexing. This uses reproducing sums.
-subroutine global_j_mean(array, j_mean, G, mask, scale)
+subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale)
   type(ocean_grid_type),            intent(inout) :: G    !< The ocean's grid structure
   real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: array  !< The variable being averaged
   real, dimension(SZI_(G)),         intent(out)   :: j_mean !<  Global mean of array along its j-axis
   real, dimension(SZI_(G),SZJ_(G)), &
-                          optional, intent(in)    :: mask !< An array used for weighting the j-mean
-  real,                   optional, intent(in)    :: scale !< A rescaling factor for the variable
+                          optional, intent(in)    :: mask  !< An array used for weighting the j-mean
+  real,                   optional, intent(in)    :: scale !< A rescaling factor for the output variable
+  real,                   optional, intent(in)    :: tmp_scale !< A rescaling factor for the internal
+                                                           !! calculations that is removed from the output
 
   ! Local variables
   type(EFP_type), allocatable, dimension(:) :: asum, mask_sum
   real :: mask_sum_r
   real :: scalefac  ! A scaling factor for the variable.
+  real :: unscale   ! A factor for undoing any internal rescaling before output.
   integer :: is, ie, js, je, idg_off, jdg_off
   integer :: i, j
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   idg_off = G%idg_offset ; jdg_off = G%jdg_offset
 
   scalefac = 1.0 ; if (present(scale)) scalefac = scale
+  unscale = 1.0
+  if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then
+    scalefac = scalefac * tmp_scale ; unscale = 1.0 / tmp_scale
+  endif ; endif
   call reset_EFP_overflow_error()
 
   allocate(asum(G%isg:G%ieg))
@@ -330,6 +346,8 @@ subroutine global_j_mean(array, j_mean, G, mask, scale)
     enddo
   endif
 
+  if (unscale /= 1.0) then ; do i=is,ie ; j_mean(i) = unscale*j_mean(i) ; enddo ; endif
+
   deallocate(asum)
 
 end subroutine global_j_mean

src/initialization/MOM_fixed_initialization.F90

@@ -159,7 +159,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
   call initialize_grid_rotation_angle(G, PF)
 
 ! Compute global integrals of grid values for later use in scalar diagnostics !
-  call compute_global_grid_integrals(G)
+  call compute_global_grid_integrals(G, US=US)
 
 ! Write out all of the grid data used by this run.
   if (write_geom) call write_ocean_geometry_file(G, PF, output_dir, US=US)

src/initialization/MOM_shared_initialization.F90

@@ -1145,17 +1145,21 @@ end subroutine set_velocity_depth_min
 ! -----------------------------------------------------------------------------
 !> Pre-compute global integrals of grid quantities (like masked ocean area) for
 !! later use in reporting diagnostics
-subroutine compute_global_grid_integrals(G)
-  type(dyn_horgrid_type), intent(inout) :: G  !< The dynamic horizontal grid
+subroutine compute_global_grid_integrals(G, US)
+  type(dyn_horgrid_type),          intent(inout) :: G  !< The dynamic horizontal grid
+  type(unit_scale_type), optional, intent(in)    :: US !< A dimensional unit scaling type
 
   ! Local variables
   real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
+  real :: area_scale  ! A scaling factor for area into MKS units
   integer :: i,j
 
+  area_scale = 1.0 ; if (present(US)) area_scale = US%L_to_m**2
+
   tmpForSumming(:,:) = 0.
   G%areaT_global = 0.0 ; G%IareaT_global = 0.0
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    tmpForSumming(i,j) = G%areaT(i,j) * G%mask2dT(i,j)
+    tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j)
   enddo ; enddo
   G%areaT_global = reproducing_sum(tmpForSumming)
 

src/initialization/MOM_state_initialization.F90

@@ -1889,16 +1889,19 @@ end subroutine set_velocity_depth_max
 
 !> Subroutine to pre-compute global integrals of grid quantities for
 !! later use in reporting diagnostics
-subroutine compute_global_grid_integrals(G)
+subroutine compute_global_grid_integrals(G, US)
   type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
+  type(unit_scale_type), intent(in)    :: US !< A dimensional unit scaling type
   ! Local variables
   real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
+  real :: area_scale
   integer :: i,j
 
+  area_scale = US%L_to_m**2
   tmpForSumming(:,:) = 0.
   G%areaT_global = 0.0 ; G%IareaT_global = 0.0
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    tmpForSumming(i,j) = G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j)
+    tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j)
   enddo ; enddo
   G%areaT_global = reproducing_sum(tmpForSumming)
   G%IareaT_global = 1. / (G%areaT_global)

src/parameterizations/lateral/MOM_MEKE.F90

@@ -188,7 +188,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
         call hchksum(MEKE%GM_src, 'MEKE GM_src', G%HI, scale=US%R_to_kg_m3*US%Z_to_m*US%L_to_m**2*US%s_to_T**3)
       if (associated(MEKE%MEKE)) call hchksum(MEKE%MEKE, 'MEKE MEKE', G%HI, scale=US%L_T_to_m_s**2)
       call uvchksum("MEKE SN_[uv]", SN_u, SN_v, G%HI, scale=US%s_to_T)
-      call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, scale=GV%H_to_m)
+      call uvchksum("MEKE h[uv]", hu, hv, G%HI, haloshift=1, scale=GV%H_to_m*US%L_to_m**2)
     endif
 
     sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -3,6 +3,7 @@ module MOM_set_diffusivity
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use MOM_checksums,           only : hchksum_pair
 use MOM_cpu_clock,           only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
 use MOM_cpu_clock,           only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
 use MOM_diag_mediator,       only : diag_ctrl, time_type
@@ -344,8 +345,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, &
 
   if (CS%useKappaShear) then
     if (CS%debug) then
-      call hchksum(u_h, "before calc_KS u_h",G%HI)
-      call hchksum(v_h, "before calc_KS v_h",G%HI)
+      call hchksum_pair("before calc_KS [uv]_h", u_h, v_h, G%HI, scale=US%L_T_to_m_s)
     endif
     call cpu_clock_begin(id_clock_kappaShear)
     if (CS%Vertex_shear) then

src/parameterizations/vertical/MOM_sponge.F90

@@ -420,7 +420,7 @@ subroutine apply_sponge(h, dt, G, GV, US, ea, eb, CS, Rcv_ml)
         eta_anom(i,j) = e_D(i,j,k) - CS%Ref_eta_im(j,k)
         if (CS%Ref_eta_im(j,K) < -G%bathyT(i,j)) eta_anom(i,j) = 0.0
       enddo ; enddo
-      call global_i_mean(eta_anom(:,:), eta_mean_anom(:,K), G)
+      call global_i_mean(eta_anom(:,:), eta_mean_anom(:,K), G, tmp_scale=US%Z_to_m)
     enddo
 
     if (CS%fldno > 0) allocate(fld_mean_anom(G%isd:G%ied,nz,CS%fldno))

src/tracer/MOM_tracer_advect.F90

@@ -485,8 +485,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
         else
           uhh(I) = uhr(I,j,k)
         endif
-       !ts2(I) = 0.5*(1.0 + uhh(I)/(hprev(i+1,j,k)+h_neglect))
-        CFL(I) = - uhh(I)/(hprev(i+1,j,k)+h_neglect) ! CFL is positive
+       !ts2(I) = 0.5*(1.0 + uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j)))
+        CFL(I) = - uhh(I) / (hprev(i+1,j,k) + h_neglect*G%areaT(i+1,j)) ! CFL is positive
       else
         hup = hprev(i,j,k) - G%areaT(i,j)*min_h
         hlos = MAX(0.0,-uhr(I-1,j,k))
@@ -497,8 +497,8 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
         else
           uhh(I) = uhr(I,j,k)
         endif
-       !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect))
-        CFL(I) = uhh(I)/(hprev(i,j,k)+h_neglect) ! CFL is positive
+       !ts2(I) = 0.5*(1.0 - uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)))
+        CFL(I) = uhh(I) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive
       endif
     enddo
 
@@ -573,7 +573,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
           ! Original implementation of PLM
          !flux_x(I,m) = uhh(I)*(Tr(m)%t(i+1,j,k) - slope_x(i+1,m)*ts2(I))
         endif
-       !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect))
+       !ts2(I) = 0.5*(1.0 - uhh(I)/(hprev(i,j,k)+h_neglect*G%areaT(i,j)))
       enddo ; enddo
     endif ! usePPM
 
@@ -856,8 +856,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
         else
           vhh(i,J) = vhr(i,J,k)
         endif
-       !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k)+h_neglect))
-        CFL(i) = - vhh(i,J) / (hprev(i,j+1,k)+h_neglect) ! CFL is positive
+       !ts2(i) = 0.5*(1.0 + vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1))
+        CFL(i) = - vhh(i,J) / (hprev(i,j+1,k) + h_neglect*G%areaT(i,j+1)) ! CFL is positive
       else
         hup = hprev(i,j,k) - G%areaT(i,j)*min_h
         hlos = MAX(0.0,-vhr(i,J-1,k))
@@ -868,8 +868,8 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
         else
           vhh(i,J) = vhr(i,J,k)
         endif
-       !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k)+h_neglect))
-        CFL(i) = vhh(i,J) / (hprev(i,j,k)+h_neglect) ! CFL is positive
+       !ts2(i) = 0.5*(1.0 - vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)))
+        CFL(i) = vhh(i,J) / (hprev(i,j,k) + h_neglect*G%areaT(i,j)) ! CFL is positive
       endif
     enddo