diff --git a/src/core_ocean/Registry.xml b/src/core_ocean/Registry.xml
index b182b7ac8a..ab8ff8be34 100644
--- a/src/core_ocean/Registry.xml
+++ b/src/core_ocean/Registry.xml
@@ -943,6 +943,10 @@
 					description="Dimensionless bottom drag coefficient, $c_{drag}$."
 					possible_values="any positive real, typically 1.0e-3"
 		/>
+		<nml_option name="config_use_implicit_bottom_roughness" type="logical" default_value=".false." units="unitless"
+					description="If true, depends on bottom roughness z0"
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_use_implicit_bottom_drag_variable" type="logical" default_value=".false." units="unitless"
 					description="If true, spatially-variable implicit bottom drag is used on the momentum equation."
 					possible_values=".true. or .false."
@@ -959,6 +963,15 @@
 					description="Dimensionless bottom drag coefficient, $c_{drag}$."
 					possible_values="any positive real, typically 1.0e-3"
 		/>
+		<nml_option name="config_use_topographic_wave_drag" type="logical" default_value=".false." units="unitless"
+					description="If true, topographic wave drag is used on the momentum equation."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_topographic_wave_drag_coeff" type="real" default_value="5.0e-4" units="unitless"
+					description="Dimensionless topographic wave drag coefficient, $c_{topo}$."
+					possible_values="any positive real, typically 5.0e-4"
+		/>
+
 	</nml_record>
 	<nml_record name="wetting_drying" mode="init;forward">
 		<nml_option name="config_use_wetting_drying" type="logical" default_value=".false." units="unitless"
@@ -1283,6 +1296,10 @@
 					description="Disables tendencies on the velocity field from horizontal surface stresses (e.g. wind stress and top drag)."
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_disable_vel_topographic_wave_drag" type="logical" default_value=".false." units="unitless"
+					description="Disables tendencies on the velocity field from topographic wave drag"
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_disable_vel_explicit_bottom_drag" type="logical" default_value=".false." units="unitless"
 					description="Disables tendencies on the velocity field from explicit bottom drag"
 					possible_values=".true. or .false."
@@ -1372,6 +1389,7 @@
 		<package name="tracerBudget" description="This package controls variables that record the tracer budget."/>
 		<package name="gm" description="This package is for GM variables, which are only needed at lower resolutions."/>
 		<package name="tidalPotentialForcingPKG" description="This package controls variables required for tidal potential forcing"/>
+		<package name="topographicWaveDragPKG" description="This package controls variables required for topographic wave drag"/>
 		<package name="gotmPKG" description="This package is for GOTM variables, which are only needed when using GOTM for the vertical turbulence closure."/>
 	</packages>
 
@@ -1958,6 +1976,15 @@
 			<var name="pointCellGlobalID"/>
 		</stream>
 
+		<stream name="topographic_wave_drag"
+				type="input"
+				input_interval="initial_only"
+				filename_template="topographic_wave_drag.nc"
+				runtime_format="single_file"
+				mode="forward;analysis">
+			<var name="topographic_wave_drag"/>
+		</stream>
+
 	</streams>
 
 	<var_struct name="state" time_levs="2">
@@ -2850,6 +2877,10 @@
 		<var name="surfaceFluxAttenuationCoefficientRunoff" type="real" dimensions="nCells Time" units="m"
 			description="The spatially-dependent length scale of exponential decay of river runoff. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
 		/>
+		<var name="topographic_wave_drag" type="real" dimensions="nEdges" units="unitless"
+			description="wave drag coefficient or 1/(rinv) where rinv is the e-folding time used in HyCOM"
+			packages="topographicWaveDragPKG"
+		/>
 		<!-- diagnostic fields for land-ice fluxes -->
 		<var name="landIceFrictionVelocity" type="real" dimensions="nCells Time" units="m s^{-1}"
 			description="The friction velocity $u_*$ under land ice"
diff --git a/src/core_ocean/driver/mpas_ocn_core_interface.F b/src/core_ocean/driver/mpas_ocn_core_interface.F
index ebe58cd70a..cbd8b55c0f 100644
--- a/src/core_ocean/driver/mpas_ocn_core_interface.F
+++ b/src/core_ocean/driver/mpas_ocn_core_interface.F
@@ -126,6 +126,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: tidalForcingActive
       logical, pointer :: tidalPotentialForcingPKGActive
       logical, pointer :: vegetationDragPKGActive
+      logical, pointer :: topographicWaveDragPKGActive
       logical, pointer :: inSituEOSActive
       logical, pointer :: variableShortwaveActive
       logical, pointer :: gmActive
@@ -156,6 +157,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: config_use_frazil_ice_formation
       logical, pointer :: config_use_tidal_forcing
       logical, pointer :: config_use_tidal_potential_forcing
+      logical, pointer :: config_use_topographic_wave_drag
       logical, pointer :: config_use_GM
       logical, pointer :: config_use_Redi
       logical, pointer :: config_use_vegetation_drag
@@ -329,6 +331,15 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
          vegetationDragPKGActive = .true.
       end if
 
+      !
+      ! test for use of topographic wave drag, topographicWaveDragPKG
+      !
+      call mpas_pool_get_package(packagePool, 'topographicWaveDragPKGActive', topographicWaveDragPKGActive)
+      call mpas_pool_get_config(configPool, 'config_use_topographic_wave_drag', config_use_topographic_wave_drag)
+      if (config_use_topographic_wave_drag) then
+         topographicWaveDragPKGActive = .true.
+      end if
+
       !
       ! test for use of gm
       !
diff --git a/src/core_ocean/ocean.cmake b/src/core_ocean/ocean.cmake
index 7a593f55cd..04d7c949fd 100644
--- a/src/core_ocean/ocean.cmake
+++ b/src/core_ocean/ocean.cmake
@@ -53,6 +53,7 @@ list(APPEND RAW_SOURCES
   core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F
   core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F
   core_ocean/shared/mpas_ocn_vel_pressure_grad.F
+  core_ocean/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F
   core_ocean/shared/mpas_ocn_vmix.F
   core_ocean/shared/mpas_ocn_vmix_coefs_redi.F
   core_ocean/shared/mpas_ocn_vmix_cvmix.F
diff --git a/src/core_ocean/shared/Makefile b/src/core_ocean/shared/Makefile
index 52671386b6..112fbfc8a0 100644
--- a/src/core_ocean/shared/Makefile
+++ b/src/core_ocean/shared/Makefile
@@ -67,7 +67,8 @@ OBJS = mpas_ocn_init_routines.o \
 	   mpas_ocn_framework_forcing.o \
 	   mpas_ocn_time_varying_forcing.o \
 	   mpas_ocn_wetting_drying.o \
-	   mpas_ocn_vel_tidal_potential.o
+	   mpas_ocn_vel_tidal_potential.o \
+	   mpas_ocn_vel_forcing_topographic_wave_drag.o
 
 all: $(OBJS)
 
@@ -105,10 +106,12 @@ mpas_ocn_vel_hmix_leith.o: mpas_ocn_constants.o mpas_ocn_config.o
 
 mpas_ocn_vel_hmix_del4.o: mpas_ocn_constants.o mpas_ocn_config.o
 
-mpas_ocn_vel_forcing.o: mpas_ocn_vel_forcing_surface_stress.o mpas_ocn_vel_forcing_explicit_bottom_drag.o mpas_ocn_forcing.o mpas_ocn_constants.o mpas_ocn_config.o
+mpas_ocn_vel_forcing.o: mpas_ocn_vel_forcing_surface_stress.o mpas_ocn_vel_forcing_explicit_bottom_drag.o mpas_ocn_forcing.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_vel_forcing_topographic_wave_drag.o
 
 mpas_ocn_vel_forcing_surface_stress.o: mpas_ocn_forcing.o mpas_ocn_constants.o mpas_ocn_config.o
 
+mpas_ocn_vel_forcing_topographic_wave_drag.o: mpas_ocn_constants.o mpas_ocn_config.o
+
 mpas_ocn_vel_forcing_explicit_bottom_drag.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_forcing.o
 
 mpas_ocn_vel_hadv_coriolis.o: mpas_ocn_constants.o mpas_ocn_config.o
diff --git a/src/core_ocean/shared/mpas_ocn_diagnostics_variables.F b/src/core_ocean/shared/mpas_ocn_diagnostics_variables.F
index 9713067a21..633767a1c2 100644
--- a/src/core_ocean/shared/mpas_ocn_diagnostics_variables.F
+++ b/src/core_ocean/shared/mpas_ocn_diagnostics_variables.F
@@ -80,6 +80,7 @@ module ocn_diagnostics_variables
    real (kind=RKIND), dimension(:), pointer :: landIceFrictionVelocity
    real (kind=RKIND), dimension(:), pointer :: topDrag
    real (kind=RKIND), dimension(:), pointer :: topDragMag
+   real (kind=RKIND), dimension(:), pointer :: topographic_wave_drag
    real (kind=RKIND), dimension(:,:), pointer :: landIceBoundaryLayerTracers
    real (kind=RKIND), dimension(:,:), pointer :: landIceTracerTransferVelocities
 
@@ -282,6 +283,10 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       call mpas_pool_get_array(diagnosticsPool, 'bottomLayerShortwaveTemperatureFlux', bottomLayerShortwaveTemperatureFlux)
       call mpas_pool_get_array(diagnosticsPool, 'indexSurfaceLayerDepth', indexSurfaceLayerDepth)
 
+      if (config_use_topographic_wave_drag) then
+         call mpas_pool_get_array(diagnosticsPool, 'topographic_wave_drag', topographic_wave_drag)
+      end if
+
       !
       ! set pointers for fields related forcing at ocean surface
       !
diff --git a/src/core_ocean/shared/mpas_ocn_vel_forcing.F b/src/core_ocean/shared/mpas_ocn_vel_forcing.F
index ef82187bf6..0a5dc1ffcb 100644
--- a/src/core_ocean/shared/mpas_ocn_vel_forcing.F
+++ b/src/core_ocean/shared/mpas_ocn_vel_forcing.F
@@ -26,6 +26,7 @@ module ocn_vel_forcing
    use ocn_constants
    use ocn_vel_forcing_surface_stress
    use ocn_vel_forcing_explicit_bottom_drag
+   use ocn_vel_forcing_topographic_wave_drag
 
    implicit none
    private
@@ -126,6 +127,10 @@ subroutine ocn_vel_forcing_tend(normalVelocity, sfcFlxAttCoeff, &
                    kineticEnergyCell, layerThicknessEdge, tend, err1)
       err = ior(err, err1)
 
+      call ocn_vel_forcing_topographic_wave_drag_tend(normalVelocity, &
+                   tend, err1)
+
+      err = ior(err, err1)
    !--------------------------------------------------------------------
 
    end subroutine ocn_vel_forcing_tend!}}}
@@ -181,6 +186,14 @@ subroutine ocn_vel_forcing_init(err)!{{{
          err = err2
       endif
 
+      call ocn_vel_forcing_topographic_wave_drag_init(err2)
+      if (err2 /= 0) then
+         call mpas_log_write( &
+            'vel_forcing_init encountered error in topographic wave drag', &
+             MPAS_LOG_ERR)
+         err = err2
+      endif
+
    !--------------------------------------------------------------------
 
    end subroutine ocn_vel_forcing_init!}}}
diff --git a/src/core_ocean/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F b/src/core_ocean/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F
new file mode 100644
index 0000000000..602b51131e
--- /dev/null
+++ b/src/core_ocean/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F
@@ -0,0 +1,192 @@
+! Copyright (c) 2013,  Los Alamos National Security, LLC (LANS)
+! and the University Corporation for Atmospheric Research (UCAR).
+!
+! Unless noted otherwise source code is licensed under the BSD license.
+! Additional copyright and license information can be found in the LICENSE file
+! distributed with this code, or at http://mpas-dev.github.com/license.html
+!
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+!  ocn_vel_forcing_topographic_wave_drag
+!
+!> \brief MPAS ocean topographic wave drag
+!> \author Nairita Pal
+!> \date   Oct 2020
+!> \details
+!>  This module contains the routine for computing
+!>  tendencies from topographic wave drag.
+!
+!-----------------------------------------------------------------------
+
+module ocn_vel_forcing_topographic_wave_drag
+
+   use mpas_derived_types
+   use mpas_pool_routines
+   use mpas_timer
+
+   use ocn_constants
+   use ocn_mesh
+   use ocn_diagnostics_variables
+   use ocn_forcing
+
+   implicit none
+   private
+   save
+
+   !--------------------------------------------------------------------
+   !
+   ! Public parameters
+   !
+   !--------------------------------------------------------------------
+
+   !--------------------------------------------------------------------
+   !
+   ! Public member functions
+   !
+   !--------------------------------------------------------------------
+
+   public :: ocn_vel_forcing_topographic_wave_drag_tend, &
+             ocn_vel_forcing_topographic_wave_drag_init
+
+   !--------------------------------------------------------------------
+   !
+   ! Private module variables
+   !
+   !--------------------------------------------------------------------
+
+   logical :: topographicWaveDragOn
+   real (kind=RKIND) :: topographicWaveDragCoeff
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+!  routine ocn_vel_forcing_topographic_wave_drag_tend
+!
+!> \brief   Computes tendency term from topographic wave drag
+!> \author  Nairita Pal
+!> \date    15 Oct 2020
+!> \details
+!>  This routine computes the topographic wave drag tendency for momentum
+!>  based on current state.
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_vel_forcing_topographic_wave_drag_tend(normalVelocity, & !{{{
+                     tend, err)
+
+      !-----------------------------------------------------------------
+      !
+      ! input variables
+      !
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:,:), intent(in) :: &
+         normalVelocity   !< Input: velocity
+
+      !-----------------------------------------------------------------
+      !
+      ! input/output variables
+      !
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:,:), intent(inout) :: &
+         tend          !< Input/Output: velocity tendency
+
+      !-----------------------------------------------------------------
+      !
+      ! output variables
+      !
+      !-----------------------------------------------------------------
+
+      integer, intent(out) :: err !< Output: error flag
+
+      !-----------------------------------------------------------------
+      !
+      ! local variables
+      !
+      !-----------------------------------------------------------------
+
+      integer :: iEdge, k
+
+      err = 0
+      if ( .not. topographicWaveDragOn ) return
+
+      call mpas_timer_start('vel topographic wave drag')
+
+      !$omp do schedule(runtime) private(k)
+      do iEdge = 1, nEdgesOwned
+        k =  maxLevelEdgeTop(iEdge)
+        ! topographic wave drag term:
+        ! du/dt = ... 1/rinv * u
+        ! rinv is the e-folding time use in HyCOM. Here
+        ! topographic_wave_drag = 1/rinv
+        if (k>0) then
+          tend(k,iEdge) = tend(k,iEdge) - topographicWaveDragCoeff &
+             * topographic_wave_drag(iEdge) * normalVelocity(k,iEdge)
+        endif
+      enddo
+      !$omp end do
+
+      call mpas_timer_stop('vel topographic wave drag')
+
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_vel_forcing_topographic_wave_drag_tend!}}}
+
+!***********************************************************************
+!
+!  routine ocn_vel_forcing_topographic_wave_drag_init
+!
+!> \brief   Initializes ocean topographic wave drag forcing
+!> \author  Nairita Pal
+!> \date    Oct 2020
+!> \details
+!>  This routine initializes quantities related to topographic wave  drag
+!>  in the ocean.
+!
+!-----------------------------------------------------------------------
+
+   subroutine ocn_vel_forcing_topographic_wave_drag_init(err)!{{{
+
+   !--------------------------------------------------------------------
+
+      !-----------------------------------------------------------------
+      !
+      ! call individual init routines for each parameterization
+      !
+      !-----------------------------------------------------------------
+
+      integer, intent(out) :: err !< Output: error flag
+
+      logical, pointer :: config_disable_vel_topographic_wave_drag
+      logical, pointer :: config_use_topographic_wave_drag
+      real (kind=RKIND), pointer :: config_topographic_wave_drag_coeff
+
+      err = 0
+
+      call mpas_pool_get_config(ocnConfigs, 'config_use_topographic_wave_drag', config_use_topographic_wave_drag)
+      call mpas_pool_get_config(ocnConfigs, 'config_topographic_wave_drag_coeff', config_topographic_wave_drag_coeff)
+      call mpas_pool_get_config(ocnConfigs, 'config_disable_vel_topographic_wave_drag', config_disable_vel_topographic_wave_drag)
+
+      topographicWaveDragCoeff = 0.0_RKIND
+
+      if (config_use_topographic_wave_drag) then
+          topographicWaveDragOn = .true.
+          topographicWaveDragCoeff = config_topographic_wave_drag_coeff
+      endif
+
+      if (config_disable_vel_topographic_wave_drag) topographicWaveDragOn = .false.
+
+   !--------------------------------------------------------------------
+
+   end subroutine ocn_vel_forcing_topographic_wave_drag_init!}}}
+
+!***********************************************************************
+
+end module ocn_vel_forcing_topographic_wave_drag
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
diff --git a/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F b/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F
index 57362b26e8..f3d07ebd50 100644
--- a/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F
+++ b/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F
@@ -307,7 +307,9 @@ subroutine ocn_vel_tidal_potential_init(domain,err)!{{{
       type (block_type), pointer :: block_ptr 
       type (mpas_pool_type), pointer :: forcingPool
       type (mpas_pool_type), pointer :: meshPool
+      integer, pointer :: nCells
       integer, pointer :: nTidalPotentialConstituents
+      real (kind=RKIND), dimension(:), pointer :: latCell
 
       ! End preamble
       !-----------------------------------------------------------------
@@ -332,9 +334,17 @@ subroutine ocn_vel_tidal_potential_init(domain,err)!{{{
          !*** eventually replace with allocated module private vars
 
          block_ptr => domain % blocklist
+
          call mpas_pool_get_subpool(block_ptr%structs, 'forcing', &
                                                       forcingPool)
 
+         call mpas_pool_get_subpool(block_ptr%structs, 'mesh', &
+                                                      meshPool)
+
+         call mpas_pool_get_dimension(block_ptr % dimensions, 'nCells', nCells)
+
+         call mpas_pool_get_array(meshPool, 'latCell', latCell)
+
          call mpas_pool_get_array(forcingPool, 'tidalPotentialEta', &
                                                 tidalPotEta)
 
@@ -369,7 +379,7 @@ subroutine ocn_vel_tidal_potential_init(domain,err)!{{{
          call mpas_set_time(refTime, &
                dateTimeString=config_tidal_potential_reference_time)
 
-         do iCell = 1,nCellsAll
+         do iCell = 1,nCells
             tidalPotEta(iCell) = 0.0_RKIND
          end do
 
@@ -428,7 +438,7 @@ subroutine ocn_vel_tidal_potential_init(domain,err)!{{{
                                        tidalConstituentType, &
                                        err)
 
-         do iCell = 1,nCellsAll
+         do iCell = 1,nCells
             lat = latCell(iCell)
             latitudeFunction(iCell,1) = 3.0_RKIND*sin(lat)**2 -1.0_RKIND 
             latitudeFunction(iCell,2) = sin(2.0_RKIND*lat)
diff --git a/src/core_ocean/shared/mpas_ocn_vmix.F b/src/core_ocean/shared/mpas_ocn_vmix.F
index 375d97ed52..01e1de81a8 100644
--- a/src/core_ocean/shared/mpas_ocn_vmix.F
+++ b/src/core_ocean/shared/mpas_ocn_vmix.F
@@ -821,7 +821,16 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forc
          end do
 
          ! average cell-based implicit bottom drag to edges and convert Mannings n to Cd
-         if (config_use_vegetation_drag .AND. config_use_vegetation_manning_equation) then
+         if (config_use_implicit_bottom_roughness) then
+           ! NCOM's formula for bottom drag coefficient
+           ! The original formula is cb(i,j) = max(cbmin, (vonk/log(0.5*depth(i,j)/z0))**2 ))
+           ! where vonk=0.16, z0=0.001, bracketed between [0.0025, 0.1]
+           ! The constant below is 250 = 0.5*0.5/z0 with the extra 0.5 from
+           ! averaging the two bottomDepth cells.
+           implicitCd = max(0.0025_RKIND, &
+                        min(0.1_RKIND, &
+                          0.16_RKIND/log(250.0_RKIND*(bottomDepth(cell1)+bottomDepth(cell2)))**2 ))
+         elseif (config_use_vegetation_drag .AND. config_use_vegetation_manning_equation) then
            implicitCd = gravity*(0.5_RKIND*(vegetationManning(cell1) + vegetationManning(cell2)))**2.0 * &
             (0.5_RKIND * (ssh(cell1) + ssh(cell2) + bottomDepth(cell1) + bottomDepth(cell2)))**(-1.0_RKIND/3.0_RKIND)
          else
@@ -1137,7 +1146,8 @@ subroutine ocn_vmix_implicit(dt, meshPool, statePool, forcingPool, scratchPool,
       if (config_use_implicit_bottom_drag_variable) then
         call ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, dt, kineticEnergyCell, &
           vertViscTopOfEdge, layerThickness, layerThickEdge, normalVelocity, err)
-      else if (config_use_implicit_bottom_drag_variable_mannings) then
+      else if (config_use_implicit_bottom_drag_variable_mannings.or. &
+               config_use_implicit_bottom_roughness) then
         ! update bottomDrag via Cd=g*n^2*h^(-1/3)
         call ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forcingPool, bottomDrag, &
           dt, kineticEnergyCell, &