diff --git a/schemes/beljaars_drag/beljaars_drag.F90 b/schemes/beljaars_drag/beljaars_drag.F90
new file mode 100644
index 00000000..28db5a37
--- /dev/null
+++ b/schemes/beljaars_drag/beljaars_drag.F90
@@ -0,0 +1,105 @@
+! Beljaars Sub-Grid Orographic (SGO) Form Drag Parameterization
+! Returns drag profile and integrated stress associated with subgrid mountains
+! with horizontal length scales nominally below 3 km.
+! Based on Beljaars et al. (2004, QJRMS) https://doi.org/10.1256/qj.03.73.
+!
+! Original author: J. Bacmeister, March 2016, based on TMS
+module beljaars_drag
+  implicit none
+  private
+
+  public :: beljaars_drag_run
+
+contains
+!> \section arg_table_beljaars_drag_run Argument Table
+!! \htmlinclude beljaars_drag_run.html
+  ! Compute Beljaars SGO form drag profile and surface stresses
+  subroutine beljaars_drag_run( &
+    do_beljaars, &
+    ncol, pver, &
+    u, v, delp, zm, sgh30, &
+    gravit, &
+    ! below output:
+    drag, taux, tauy, &
+    errmsg, errflg)
+    use ccpp_kinds, only: kind_phys
+
+    logical,            intent(in)  :: do_beljaars  ! is Beljaars active?
+    integer,            intent(in)  :: ncol
+    integer,            intent(in)  :: pver
+    real(kind_phys),    intent(in)  :: u(:, :)      ! zonal wind [m s-1]
+    real(kind_phys),    intent(in)  :: v(:, :)      ! meridional wind [m s-1]
+    real(kind_phys),    intent(in)  :: delp(:, :)   ! air pressure thickness [Pa]
+    real(kind_phys),    intent(in)  :: zm(:, :)     ! geopotential height wrt surface [m]
+    real(kind_phys),    intent(in)  :: sgh30(:)     ! standard deviation of subgrid orography [m]
+    real(kind_phys),    intent(in)  :: gravit       ! gravitational acceleration [m s-2]
+
+    real(kind_phys),    intent(out) :: drag(:, :)   ! SGO drag profile [s-1]
+    real(kind_phys),    intent(out) :: taux(:)      ! surface zonal wind stress [N m-2]
+    real(kind_phys),    intent(out) :: tauy(:)      ! surface meridional wind stress [N m-2]
+    character(len=*),   intent(out) :: errmsg
+    integer,            intent(out) :: errflg
+
+    ! Local variables
+    integer  :: i, k
+
+    real(kind_phys) :: vmag                         ! velocity magnitude [m s-1]
+
+    real(kind_phys) :: alpha, beta, Cmd, Ccorr, n1, n2, k1, kflt, IH
+    real(kind_phys) :: a1(ncol), a2(ncol)
+
+    errmsg = ''
+    errflg = 0
+
+    if(.not. do_beljaars) then
+      ! if not doing Beljaars, return zero drag and stresses from routine.
+      drag(:,:) = 0._kind_phys
+      taux(:)   = 0._kind_phys
+      tauy(:)   = 0._kind_phys
+      return
+    end if
+
+    alpha =  12.0_kind_phys
+    beta  =  1.0_kind_phys
+    n1    = -1.9_kind_phys
+    n2    = -2.8_kind_phys
+
+    Cmd   = 0.005_kind_phys
+    Ccorr = 0.6_kind_phys * 5.0_kind_phys
+
+    kflt  = 0.00035_kind_phys   ! m-1
+    k1    = 0.003_kind_phys     ! m-1
+    IH    = 0.00102_kind_phys   ! m-1
+
+    a1(1:ncol) = (sgh30(1:ncol) * sgh30(1:ncol)) / (IH * (kflt**n1))
+    a2(1:ncol) = a1(1:ncol) * k1**(n1 - n2)
+
+    do k = 1, pver
+      do i = 1, ncol
+        vmag      = sqrt(u(i, k)**2 + v(i, k)**2)
+        drag(i, k) = -alpha * beta * Cmd * Ccorr * vmag * 2.109_kind_phys * &
+                     exp(-(zm(i, k) / 1500.0_kind_phys) * sqrt(zm(i, k) / 1500.0_kind_phys)) * &
+                     (zm(i, k)**(-1.2_kind_phys)) * a2(i)
+      end do
+    end do
+
+    ! Diagnose effective surface drag in X and Y by integrating in the vertical
+    !
+    ! taux, tauy stresses generated by Beljaars drag here
+    ! uses the pre-vertical diffusion winds (and not the updated winds)
+    ! however, at this point only these winds are available because the diffusion
+    ! solver runs after the orographic drag.
+    ! the actual provisionally updated winds are actually used to recompute taubljx/y
+    ! which are added to the updated residual stress. see the physics scheme
+    ! beljaars_add_updated_residual_stress. (hplin, 5/13/26)
+    taux(:) = 0.0_kind_phys
+    tauy(:) = 0.0_kind_phys
+    do k = 1, pver
+      do i = 1, ncol
+        taux(i) = taux(i) + drag(i, k) * u(i, k) * delp(i, k) / gravit
+        tauy(i) = tauy(i) + drag(i, k) * v(i, k) * delp(i, k) / gravit
+      end do
+    end do
+
+  end subroutine beljaars_drag_run
+end module beljaars_drag
diff --git a/schemes/beljaars_drag/beljaars_drag.meta b/schemes/beljaars_drag/beljaars_drag.meta
new file mode 100644
index 00000000..1ccaad39
--- /dev/null
+++ b/schemes/beljaars_drag/beljaars_drag.meta
@@ -0,0 +1,91 @@
+[ccpp-table-properties]
+  name = beljaars_drag
+  type = scheme
+
+[ccpp-arg-table]
+  name  = beljaars_drag_run
+  type  = scheme
+[ do_beljaars ]
+  standard_name = do_beljaars_drag_in_vertical_diffusion
+  units = flag
+  type = logical
+  dimensions = ()
+  intent = in
+[ ncol ]
+  standard_name = horizontal_loop_extent
+  units = count
+  type = integer
+  dimensions = ()
+  intent = in
+[ pver ]
+  standard_name = vertical_layer_dimension
+  units = count
+  type = integer
+  dimensions = ()
+  intent = in
+[ u ]
+  standard_name = eastward_wind
+  units = m s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ v ]
+  standard_name = northward_wind
+  units = m s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ delp ]
+  standard_name = air_pressure_thickness
+  units = Pa
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ zm ]
+  standard_name = geopotential_height_wrt_surface
+  units = m
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ sgh30 ]
+  standard_name = standard_deviation_of_subgrid_orography_for_turbulent_orographic_form_drag
+  units = m
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = in
+[ gravit ]
+  standard_name = standard_gravitational_acceleration
+  units = m s-2
+  type = real | kind = kind_phys
+  dimensions = ()
+  intent = in
+[ drag ]
+  standard_name = turbulent_orographic_form_drag_coefficent
+  units = s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = out
+[ taux ]
+  standard_name = eastward_beljaars_surface_stress_tbd
+  units = N m-2
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = out
+[ tauy ]
+  standard_name = northward_beljaars_surface_stress_tbd
+  units = N m-2
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = out
+[ errmsg ]
+  standard_name = ccpp_error_message
+  units = none
+  type = character | kind = len=*
+  dimensions = ()
+  intent = out
+[ errflg ]
+  standard_name = ccpp_error_code
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
diff --git a/schemes/beljaars_drag/beljaars_drag_interstitials.F90 b/schemes/beljaars_drag/beljaars_drag_interstitials.F90
new file mode 100644
index 00000000..72e69741
--- /dev/null
+++ b/schemes/beljaars_drag/beljaars_drag_interstitials.F90
@@ -0,0 +1,125 @@
+! Interstitial schemes for Beljaars drag.
+!
+! Put in a separate file than beljaars_drag.F90 because
+! scheme files with multiple schemes cannot have a namelist file attached.
+module beljaars_drag_interstitials
+
+  implicit none
+  private
+
+  ! public CCPP-compliant subroutines
+  public :: beljaars_add_wind_damping_rate_run
+  public :: beljaars_add_updated_residual_stress_run
+contains
+
+  ! Add Beljaars drag to the wind damping rate for vertical diffusion.
+  ! Has to run after vertical_diffusion_wind_damping_rate
+!> \section arg_table_beljaars_add_wind_damping_rate_run Argument Table
+!! \htmlinclude arg_table_beljaars_add_wind_damping_rate_run.html
+  subroutine beljaars_add_wind_damping_rate_run( &
+    ncol, pver, &
+    dragblj, &
+    tau_damp_rate, &
+    errmsg, errflg)
+
+    use ccpp_kinds, only: kind_phys
+
+    ! Input arguments
+    integer,            intent(in)    :: ncol
+    integer,            intent(in)    :: pver
+    real(kind_phys),    intent(in)    :: dragblj(:, :)    ! Drag profile from Beljaars SGO form drag > 0. [s-1]
+
+    ! Input/output arguments
+    real(kind_phys),    intent(inout) :: tau_damp_rate(:,:) ! Rate at which external (surface) stress damps wind speeds [s-1]
+
+    ! Output arguments
+    character(len=*),   intent(out)   :: errmsg  ! error message
+    integer,            intent(out)   :: errflg  ! error flag
+
+    errmsg = ''
+    errflg = 0
+
+    ! Beljaars et al SGO scheme incorporated here. It
+    ! appears as a "3D" tau_damp_rate specification.
+    tau_damp_rate(:ncol,:pver) = tau_damp_rate(:ncol,:pver) + dragblj(:ncol,:pver)
+
+  end subroutine beljaars_add_wind_damping_rate_run
+
+  ! Add Beljaars stress using updated provisional winds to surface stresses used
+  ! for horizontal momentum diffusion - kinetic energy dissipation.
+!> \section arg_table_beljaars_add_updated_residual_stress_run Argument Table
+!! \htmlinclude arg_table_beljaars_add_updated_residual_stress_run.html
+  subroutine beljaars_add_updated_residual_stress_run( &
+    ncol, pver, &
+    gravit, &
+    p, &
+    do_iss, &
+    itaures, &
+    dragblj, &
+    u1, v1, & ! provisionally updated winds
+    ! input/output
+    tauresx, tauresy, &
+    errmsg, errflg)
+
+    use ccpp_kinds, only: kind_phys
+    use coords_1d, only: Coords1D
+
+    ! Input arguments
+    integer,         intent(in)       :: ncol
+    integer,         intent(in)       :: pver
+    real(kind_phys), intent(in)       :: gravit
+    type(coords1d),  intent(in)       :: p                ! Pressure coordinates [Pa]
+    logical,         intent(in)       :: do_iss           ! Flag for implicit surface stress (namelist from vdiff)
+    logical,         intent(in)       :: itaures          ! Flag for updating tauresx tauresy in this subroutine.
+    real(kind_phys), intent(in)       :: dragblj(:, :)    ! Drag profile from Beljaars SGO form drag > 0. [s-1]
+    real(kind_phys), intent(in)       :: u1(:,:)          ! After vertical diffusion u-wind [m s-1]
+    real(kind_phys), intent(in)       :: v1(:,:)          ! After vertical diffusion v-wind [m s-1]
+
+    ! Input/output arguments
+    real(kind_phys), intent(inout)    :: tauresx(:)       ! Partially updated residual surface stress using provisional winds
+    real(kind_phys), intent(inout)    :: tauresy(:)
+    character(len=*),   intent(out)   :: errmsg           ! error message
+    integer,            intent(out)   :: errflg           ! error flag
+
+    integer :: i, k
+    real(kind_phys) :: taubljx(ncol)                      ! recomputed explicit/residual beljaars stress
+    real(kind_phys) :: taubljy(ncol)                      ! recomputed explicit/residual beljaars stress
+
+    errmsg = ''
+    errflg = 0
+
+    if(.not. do_iss) then
+      ! Not added to residual stress if implicit surface stress is off
+      return
+    end if
+
+    if(.not. itaures) then
+      ! Not added to residual stress if residual stress is not requested to be updated.
+      return
+    end if
+
+    do i = 1, ncol
+      ! Add vertically-integrated Beljaars drag to residual stress
+      ! these are calculated using provisionally-updated winds.
+      !
+      ! A previous FIXME (predating CAM6) in beljaars_drag notes:
+      ! "uses 'state' u and v. Should updated u and v's be used?"
+      !
+      ! It may be possible to write back these provisionally updated taubljx/y
+      ! to the model state instead of a local variable (and pbuf in former CAM)
+      ! but it would introduce answer changes, so moving that FIXME to here for now.
+      ! (hplin 5/22/25, 5/13/26)
+      taubljx(i) = 0._kind_phys
+      taubljy(i) = 0._kind_phys
+      do k = 1, pver
+        taubljx(i) = taubljx(i) + (1._kind_phys/gravit)*dragblj(i, k)*u1(i, k)*p%del(i, k)
+        taubljy(i) = taubljy(i) + (1._kind_phys/gravit)*dragblj(i, k)*v1(i, k)*p%del(i, k)
+      end do
+    end do
+
+    tauresx(:ncol) = tauresx(:ncol) + taubljx(:ncol)
+    tauresy(:ncol) = tauresy(:ncol) + taubljy(:ncol)
+
+  end subroutine beljaars_add_updated_residual_stress_run
+
+end module beljaars_drag_interstitials
diff --git a/schemes/beljaars_drag/beljaars_drag_interstitials.meta b/schemes/beljaars_drag/beljaars_drag_interstitials.meta
new file mode 100644
index 00000000..c803c955
--- /dev/null
+++ b/schemes/beljaars_drag/beljaars_drag_interstitials.meta
@@ -0,0 +1,129 @@
+[ccpp-table-properties]
+  name = beljaars_add_wind_damping_rate
+  type = scheme
+
+[ccpp-arg-table]
+  name  = beljaars_add_wind_damping_rate_run
+  type  = scheme
+[ ncol ]
+  standard_name = horizontal_loop_extent
+  units = count
+  type = integer
+  dimensions = ()
+  intent = in
+[ pver ]
+  standard_name = vertical_layer_dimension
+  units = count
+  type = integer
+  dimensions = ()
+  intent = in
+[ dragblj ]
+  standard_name = turbulent_orographic_form_drag_coefficent
+  units = s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ tau_damp_rate ]
+  standard_name = horizontal_wind_damping_rate_due_to_total_surface_drag
+  units = s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = inout
+[ errmsg ]
+  standard_name = ccpp_error_message
+  units = none
+  type = character | kind = len=*
+  dimensions = ()
+  intent = out
+[ errflg ]
+  standard_name = ccpp_error_code
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
+
+[ccpp-table-properties]
+  name = beljaars_add_updated_residual_stress
+  type = scheme
+
+[ccpp-arg-table]
+  name  = beljaars_add_updated_residual_stress_run
+  type  = scheme
+[ ncol ]
+  standard_name = horizontal_loop_extent
+  units = count
+  type = integer
+  dimensions = ()
+  intent = in
+[ pver ]
+  standard_name = vertical_layer_dimension
+  units = count
+  type = integer
+  dimensions = ()
+  intent = in
+[ gravit ]
+  standard_name = standard_gravitational_acceleration
+  units = m s-2
+  type = real | kind = kind_phys
+  dimensions = ()
+  intent = in
+[ p ]
+  standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
+  units = none
+  type = coords1d
+  dimensions = ()
+  intent = in
+[ do_iss ]
+  standard_name = do_implicit_total_surface_stress_in_vertical_diffusion
+  units = flag
+  type = logical
+  dimensions = ()
+  intent = in
+[ itaures ]
+  standard_name = update_residual_stress_at_surface_in_vertical_diffusion
+  units = flag
+  type = logical
+  dimensions = ()
+  intent = in
+[ dragblj ]
+  standard_name = turbulent_orographic_form_drag_coefficent
+  units = s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ u1 ]
+  standard_name = eastward_wind_after_vertical_diffusion
+  units = m s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ v1 ]
+  standard_name = northward_wind_after_vertical_diffusion
+  units = m s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ tauresx ]
+  standard_name = eastward_reserved_stress_at_surface_on_previous_timestep
+  units = N m-2
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = inout
+[ tauresy ]
+  standard_name = northward_reserved_stress_at_surface_on_previous_timestep
+  units = N m-2
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = inout
+[ errmsg ]
+  standard_name = ccpp_error_message
+  units = none
+  type = character | kind = len=*
+  dimensions = ()
+  intent = out
+[ errflg ]
+  standard_name = ccpp_error_code
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
diff --git a/schemes/beljaars_drag/beljaars_drag_namelist.xml b/schemes/beljaars_drag/beljaars_drag_namelist.xml
new file mode 100644
index 00000000..27c5ee33
--- /dev/null
+++ b/schemes/beljaars_drag/beljaars_drag_namelist.xml
@@ -0,0 +1,93 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl"?>
+<entry_id_pg version="2.0">
+
+<!-- Each namelist variable is defined in an <entry> element.  The
+     content of the element is the documentation of how the variable is
+     used.  Other aspects of the variable's definition are expressed as
+     attributes of the <entry> element.  Note that it is an XML requirement
+     that the attribute values are enclosed in quotes.  The attributes are:
+
+     id
+          The variable's name.  *** N.B. *** The name must be lower case.
+          The module convert all namelist variable names to lower case
+          since Fortran is case insensitive.
+
+     type
+          An abbreviation of the fortran declaration for the variable.
+    Valid declarations are:
+
+          char*n
+    integer
+    logical
+    real
+
+    Any of these types may be followed by a comma separated list of
+    integers enclosed in parenthesis to indicate an array.
+
+    The current namelist validation code only distinquishes between
+    string and non-string types.
+
+     input_pathname
+          Only include this attribute to indicate that the variable
+          contains the pathname of an input dataset that resides in the
+          CESM inputdata directory tree.  Note that the variables
+          containing the names of restart files that are used in branch
+          runs don't reside in the inputdata tree and should not be given
+          this attribute.
+
+    The recognized values are "abs" to indicate that an absolute
+          pathname is required, or "rel:var_name" to indicate that the
+          pathname is relative and that the namelist variable "var_name"
+          contains the absolute root directory.
+
+     category
+          A category assigned for organizing the documentation.
+
+     group
+          The namelist group that the variable is declared in.
+
+     valid_values
+          This is an optional attribute that is mainly useful for variables
+          that have only a small number of allowed values. It is a
+          comma-separated list of valid values for the variable.
+
+     desc
+          This is a text description of the variable including its
+          purpose and use.
+
+     values
+          This is a grouping of all the cases where a value can be
+          assigned with the namelist generator has been run.
+
+     value
+          Each <value> tag specifies a case where a value is assigned
+          to this variable when the namelist generator is run. If the
+          tag has no attributes, it is a default value. In general,
+          the namelist generator attempts to find a value with the
+          maximium number of attribute matches (and no non-matches).
+
+     standard_name
+          This is the CCPP Standard Name of the variable
+
+     units
+          This is the CCPP unit specification of the variable (e.g., m s-1).
+                                                                        -->
+
+  <entry id="do_beljaars">
+    <type>logical</type>
+    <category>pbl</category>
+    <group>blj_nl</group>
+    <standard_name>do_beljaars_drag_in_vertical_diffusion</standard_name>
+    <units>flag</units>
+    <desc>
+      Logical switch to turn on the Beljaars SGO form drag scheme.
+    </desc>
+    <values>
+      <value>.false.</value>
+      <value phys_suite="cam6">.true.</value>
+      <value phys_suite="cam7">.true.</value>
+    </values>
+  </entry>
+
+</entry_id_pg>
diff --git a/schemes/sima_diagnostics/beljaars_drag_diagnostics.F90 b/schemes/sima_diagnostics/beljaars_drag_diagnostics.F90
new file mode 100644
index 00000000..a39f21ad
--- /dev/null
+++ b/schemes/sima_diagnostics/beljaars_drag_diagnostics.F90
@@ -0,0 +1,58 @@
+! Diagnostic scheme for Beljaars SGO form drag
+module beljaars_drag_diagnostics
+  implicit none
+  private
+
+  public :: beljaars_drag_diagnostics_init
+  public :: beljaars_drag_diagnostics_run
+
+contains
+
+!> \section arg_table_beljaars_drag_diagnostics_init  Argument Table
+!! \htmlinclude beljaars_drag_diagnostics_init.html
+  subroutine beljaars_drag_diagnostics_init(errmsg, errflg)
+    use cam_history,         only: history_add_field
+    use cam_history_support, only: horiz_only
+
+    character(len=*),   intent(out) :: errmsg
+    integer,            intent(out) :: errflg
+
+    errmsg = ''
+    errflg = 0
+
+    ! 3D field (layer centers)
+    call history_add_field('DRAGBLJ', 'turbulent_orographic_form_drag_coefficent', &
+                           'lev', 'inst', 's-1')
+
+    ! 2D fields (horizontal only)
+    call history_add_field('TAUBLJX', 'eastward_beljaars_surface_stress_tbd', &
+                           horiz_only, 'inst', 'N m-2')
+    call history_add_field('TAUBLJY', 'northward_beljaars_surface_stress_tbd', &
+                           horiz_only, 'inst', 'N m-2')
+
+  end subroutine beljaars_drag_diagnostics_init
+
+!> \section arg_table_beljaars_drag_diagnostics_run  Argument Table
+!! \htmlinclude beljaars_drag_diagnostics_run.html
+  subroutine beljaars_drag_diagnostics_run( &
+    drag, taux, tauy, &
+    errmsg, errflg)
+    use ccpp_kinds,  only: kind_phys
+    use cam_history, only: history_out_field
+
+    real(kind_phys),    intent(in)  :: drag(:, :)   ! SGO drag profile [s-1]
+    real(kind_phys),    intent(in)  :: taux(:)      ! surface zonal wind stress [N m-2]
+    real(kind_phys),    intent(in)  :: tauy(:)      ! surface meridional wind stress [N m-2]
+    character(len=*),   intent(out) :: errmsg
+    integer,            intent(out) :: errflg
+
+    errmsg = ''
+    errflg = 0
+
+    call history_out_field('DRAGBLJ', drag)
+    call history_out_field('TAUBLJX', taux)
+    call history_out_field('TAUBLJY', tauy)
+
+  end subroutine beljaars_drag_diagnostics_run
+
+end module beljaars_drag_diagnostics
diff --git a/schemes/sima_diagnostics/beljaars_drag_diagnostics.meta b/schemes/sima_diagnostics/beljaars_drag_diagnostics.meta
new file mode 100644
index 00000000..d85feab9
--- /dev/null
+++ b/schemes/sima_diagnostics/beljaars_drag_diagnostics.meta
@@ -0,0 +1,53 @@
+[ccpp-table-properties]
+  name = beljaars_drag_diagnostics
+  type = scheme
+
+[ccpp-arg-table]
+  name  = beljaars_drag_diagnostics_init
+  type  = scheme
+[ errmsg ]
+  standard_name = ccpp_error_message
+  units = none
+  type = character | kind = len=*
+  dimensions = ()
+  intent = out
+[ errflg ]
+  standard_name = ccpp_error_code
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
+
+[ccpp-arg-table]
+  name  = beljaars_drag_diagnostics_run
+  type  = scheme
+[ drag ]
+  standard_name = turbulent_orographic_form_drag_coefficent
+  units = s-1
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
+  intent = in
+[ taux ]
+  standard_name = eastward_beljaars_surface_stress_tbd
+  units = N m-2
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = in
+[ tauy ]
+  standard_name = northward_beljaars_surface_stress_tbd
+  units = N m-2
+  type = real | kind = kind_phys
+  dimensions = (horizontal_loop_extent)
+  intent = in
+[ errmsg ]
+  standard_name = ccpp_error_message
+  units = none
+  type = character | kind = len=*
+  dimensions = ()
+  intent = out
+[ errflg ]
+  standard_name = ccpp_error_code
+  units = 1
+  type = integer
+  dimensions = ()
+  intent = out
diff --git a/schemes/vertical_diffusion/diffusion_solver.meta b/schemes/vertical_diffusion/diffusion_solver.meta
index 1063cf92..103167a5 100644
--- a/schemes/vertical_diffusion/diffusion_solver.meta
+++ b/schemes/vertical_diffusion/diffusion_solver.meta
@@ -384,7 +384,7 @@
   dimensions = ()
   intent = in
 [ do_beljaars ]
-  standard_name = do_beljaars_drag_in_vertical_diffusion_tbd
+  standard_name = do_beljaars_drag_in_vertical_diffusion
   units = flag
   type = logical
   dimensions = ()
diff --git a/schemes/vertical_diffusion/diffusion_stubs.F90 b/schemes/vertical_diffusion/diffusion_stubs.F90
index 9f5357d3..cdac0155 100644
--- a/schemes/vertical_diffusion/diffusion_stubs.F90
+++ b/schemes/vertical_diffusion/diffusion_stubs.F90
@@ -9,16 +9,12 @@ module diffusion_stubs
 
   implicit none
   private
-  save
 
   ! CCPP-compliant subroutines
   public :: zero_upper_boundary_condition_init
 
   public :: tms_beljaars_zero_stub_run
   public :: turbulent_mountain_stress_add_drag_coefficient_run
-  public :: beljaars_add_wind_damping_rate_run
-
-  public :: beljaars_add_updated_residual_stress_run
   public :: turbulent_mountain_stress_add_updated_surface_stress_run
 
   public :: vertical_diffusion_not_use_rairv_init
@@ -139,108 +135,6 @@ subroutine turbulent_mountain_stress_add_drag_coefficient_run( &
 
   end subroutine turbulent_mountain_stress_add_drag_coefficient_run
 
-  ! Add Beljaars drag to the wind damping rate for vertical diffusion.
-  ! Has to run after vertical_diffusion_wind_damping_rate
-!> \section arg_table_beljaars_add_wind_damping_rate_run Argument Table
-!! \htmlinclude arg_table_beljaars_add_wind_damping_rate_run.html
-  subroutine beljaars_add_wind_damping_rate_run( &
-    ncol, pver, &
-    dragblj, &
-    tau_damp_rate, &
-    errmsg, errflg)
-
-    ! Input arguments
-    integer,            intent(in)    :: ncol
-    integer,            intent(in)    :: pver
-    real(kind_phys),    intent(in)    :: dragblj(:, :)    ! Drag profile from Beljaars SGO form drag > 0. [s-1]
-
-    ! Input/output arguments
-    real(kind_phys),    intent(inout) :: tau_damp_rate(:,:) ! Rate at which external (surface) stress damps wind speeds [s-1]
-
-    ! Output arguments
-    character(len=512), intent(out)   :: errmsg  ! error message
-    integer,            intent(out)   :: errflg  ! error flag
-
-    errmsg = ''
-    errflg = 0
-
-    ! Beljaars et al SGO scheme incorporated here. It
-    ! appears as a "3D" tau_damp_rate specification.
-    tau_damp_rate(:ncol,:pver) = tau_damp_rate(:ncol,:pver) + dragblj(:ncol,:pver)
-
-  end subroutine beljaars_add_wind_damping_rate_run
-
-  ! Add Beljaars stress using updated provisional winds to surface stresses used
-  ! for horizontal momentum diffusion - kinetic energy dissipation.
-!> \section arg_table_beljaars_add_updated_residual_stress_run Argument Table
-!! \htmlinclude arg_table_beljaars_add_updated_residual_stress_run.html
-  subroutine beljaars_add_updated_residual_stress_run( &
-    ncol, pver, &
-    gravit, &
-    p, &
-    do_iss, &
-    itaures, &
-    dragblj, &
-    u1, v1, & ! provisionally updated winds
-    ! input/output
-    tauresx, tauresy, &
-    errmsg, errflg)
-
-    use coords_1d, only: Coords1D
-
-    ! Input arguments
-    integer,         intent(in)       :: ncol
-    integer,         intent(in)       :: pver
-    real(kind_phys), intent(in)       :: gravit
-    type(coords1d),  intent(in)       :: p                ! Pressure coordinates [Pa]
-    logical,         intent(in)       :: do_iss           ! Flag for implicit surface stress (namelist from vdiff)
-    logical,         intent(in)       :: itaures          ! Flag for updating tauresx tauresy in this subroutine.
-    real(kind_phys), intent(in)       :: dragblj(:, :)    ! Drag profile from Beljaars SGO form drag > 0. [s-1]
-    real(kind_phys), intent(in)       :: u1(:,:)          ! After vertical diffusion u-wind [m s-1]
-    real(kind_phys), intent(in)       :: v1(:,:)          ! After vertical diffusion v-wind [m s-1]
-
-    ! Input/output arguments
-    real(kind_phys), intent(inout)    :: tauresx(:)       ! Partially updated residual surface stress
-    real(kind_phys), intent(inout)    :: tauresy(:)
-    character(len=512), intent(out)   :: errmsg  ! error message
-    integer,            intent(out)   :: errflg  ! error flag
-
-    integer :: i, k
-    real(kind_phys) :: taubljx(ncol)             ! recomputed explicit/residual beljaars stress
-    real(kind_phys) :: taubljy(ncol)             ! recomputed explicit/residual beljaars stress
-
-    errmsg = ''
-    errflg = 0
-
-    if(.not. do_iss) then
-      ! Not added to residual stress if implicit surface stress is off
-      return
-    endif
-
-    if(.not. itaures) then
-      ! Not added to residual stress if residual stress is not requested to be updated.
-      return
-    endif
-
-    do i = 1, ncol
-      ! Add vertically-integrated Beljaars drag to residual stress.
-      ! note: in vertical_diffusion_tend, tautotx has taubljx added to it
-      ! but this is computed separately using UPDATED u, v (and is not written back to pbuf)
-      ! there is a FIXME in beljaars_drag that notes maybe updated u, v could be used there
-      ! but that is too early (before vdiff). Keeping this for now. hplin 5/22/25
-      taubljx(i) = 0._kind_phys
-      taubljy(i) = 0._kind_phys
-      do k = 1, pver
-        taubljx(i) = taubljx(i) + (1._kind_phys/gravit)*dragblj(i, k)*u1(i, k)*p%del(i, k)
-        taubljy(i) = taubljy(i) + (1._kind_phys/gravit)*dragblj(i, k)*v1(i, k)*p%del(i, k)
-      end do
-    enddo
-
-    tauresx(:ncol) = tauresx(:ncol) + taubljx(:ncol)
-    tauresy(:ncol) = tauresy(:ncol) + taubljy(:ncol)
-
-  end subroutine beljaars_add_updated_residual_stress_run
-
   ! Add TMS using updated provisional winds to total surface stress in
   ! horizontal momentum diffusion - kinetic energy dissipation.
 !> \section arg_table_turbulent_mountain_stress_add_updated_surface_stress_run Argument Table
diff --git a/schemes/vertical_diffusion/diffusion_stubs.meta b/schemes/vertical_diffusion/diffusion_stubs.meta
index ff2836e0..83a90f4f 100644
--- a/schemes/vertical_diffusion/diffusion_stubs.meta
+++ b/schemes/vertical_diffusion/diffusion_stubs.meta
@@ -80,7 +80,7 @@
   type = real | kind = kind_phys
   intent = out
 [ do_beljaars ]
-  standard_name = do_beljaars_drag_in_vertical_diffusion_tbd
+  standard_name = do_beljaars_drag_in_vertical_diffusion
   units = flag
   type = logical
   dimensions = ()
@@ -160,136 +160,6 @@
   type = integer
   intent = out
 
-[ccpp-table-properties]
-  name = beljaars_add_wind_damping_rate
-  type = scheme
-
-[ccpp-arg-table]
-  name  = beljaars_add_wind_damping_rate_run
-  type  = scheme
-[ ncol ]
-  standard_name = horizontal_loop_extent
-  units = count
-  type = integer
-  dimensions = ()
-  intent = in
-[ pver ]
-  standard_name = vertical_layer_dimension
-  units = count
-  type = integer
-  dimensions = ()
-  intent = in
-[ dragblj ]
-  standard_name = turbulent_orographic_form_drag_coefficent
-  units = s-1
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
-  intent = in
-[ tau_damp_rate ]
-  standard_name = horizontal_wind_damping_rate_due_to_total_surface_drag
-  units = s-1
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
-  intent = inout
-[ errmsg ]
-  standard_name = ccpp_error_message
-  units = none
-  type = character | kind = len=512
-  dimensions = ()
-  intent = out
-[ errflg ]
-  standard_name = ccpp_error_code
-  units = 1
-  type = integer
-  dimensions = ()
-  intent = out
-
-[ccpp-table-properties]
-  name = beljaars_add_updated_residual_stress
-  type = scheme
-
-[ccpp-arg-table]
-  name  = beljaars_add_updated_residual_stress_run
-  type  = scheme
-[ ncol ]
-  standard_name = horizontal_loop_extent
-  units = count
-  type = integer
-  dimensions = ()
-  intent = in
-[ pver ]
-  standard_name = vertical_layer_dimension
-  units = count
-  type = integer
-  dimensions = ()
-  intent = in
-[ gravit ]
-  standard_name = standard_gravitational_acceleration
-  units = m s-2
-  type = real | kind = kind_phys
-  dimensions = ()
-  intent = in
-[ p ]
-  standard_name = vertical_moist_pressure_coordinates_for_vertical_diffusion
-  units = none
-  type = coords1d
-  dimensions = ()
-  intent = in
-[ do_iss ]
-  standard_name = do_implicit_total_surface_stress_in_vertical_diffusion
-  units = flag
-  type = logical
-  dimensions = ()
-  intent = in
-[ itaures ]
-  standard_name = update_residual_stress_at_surface_in_vertical_diffusion
-  units = flag
-  type = logical
-  dimensions = ()
-  intent = in
-[ dragblj ]
-  standard_name = turbulent_orographic_form_drag_coefficent
-  units = s-1
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
-  intent = in
-[ u1 ]
-  standard_name = eastward_wind_after_vertical_diffusion
-  units = m s-1
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
-  intent = in
-[ v1 ]
-  standard_name = northward_wind_after_vertical_diffusion
-  units = m s-1
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent, vertical_layer_dimension)
-  intent = in
-[ tauresx ]
-  standard_name = eastward_reserved_stress_at_surface_on_previous_timestep
-  units = N m-2
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent)
-  intent = inout
-[ tauresy ]
-  standard_name = northward_reserved_stress_at_surface_on_previous_timestep
-  units = N m-2
-  type = real | kind = kind_phys
-  dimensions = (horizontal_loop_extent)
-  intent = inout
-[ errmsg ]
-  standard_name = ccpp_error_message
-  units = none
-  type = character | kind = len=512
-  dimensions = ()
-  intent = out
-[ errflg ]
-  standard_name = ccpp_error_code
-  units = 1
-  type = integer
-  dimensions = ()
-  intent = out
-
 [ccpp-table-properties]
   name = turbulent_mountain_stress_add_updated_surface_stress
   type = scheme
diff --git a/test/test_suites/suite_beljaars_form_drag.xml b/test/test_suites/suite_beljaars_form_drag.xml
new file mode 100644
index 00000000..aaf73971
--- /dev/null
+++ b/test/test_suites/suite_beljaars_form_drag.xml
@@ -0,0 +1,11 @@
+<?xml version="1.0" encoding="UTF-8"?>
+
+<!-- suite is named beljaars_form_drag
+     as it cannot be the same as the scheme name beljaars_drag -->
+<suite name="beljaars_form_drag" version="1.0">
+  <group name="physics_before_coupler">
+    <!-- Beljaars SGO form drag -->
+    <scheme>beljaars_drag</scheme>
+    <scheme>beljaars_drag_diagnostics</scheme>
+  </group>
+</suite>
diff --git a/to_be_ccppized/error_messages.F90 b/to_be_ccppized/error_messages.F90
index a2a64bca..6bd7d6b8 100644
--- a/to_be_ccppized/error_messages.F90
+++ b/to_be_ccppized/error_messages.F90
@@ -106,7 +106,7 @@ subroutine handle_ncerr( ret, mes, line )
       if ( ret .ne. NF90_NOERR ) then
          if(present(line)) then
             write(iulog,*) mes, line
-         else	
+         else
             write(iulog,*) mes
          end if
          write(iulog,*) nf90_strerror( ret )