diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90
index 0913c54ebd..a35c2d2a90 100644
--- a/src/core/MOM_PressureForce_FV.F90
+++ b/src/core/MOM_PressureForce_FV.F90
@@ -49,6 +49,7 @@ module MOM_PressureForce_FV
                             !! timing of diagnostic output.
   integer :: MassWghtInterp !< A flag indicating whether and how to use mass weighting in T/S interpolation
   logical :: correction_intxpa !< If true, apply a correction to surface intxpa under ice.
+  logical :: correction_intxpa_5pt ! Use 5 point quadrature to calculate surface intxpa
   logical :: use_inaccurate_pgf_rho_anom !< If true, uses the older and less accurate
                             !! method to calculate density anomalies, as used prior to
                             !! March 2018.
@@ -660,9 +661,14 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
   integer, dimension(2) :: EOSdom_h ! The i-computational domain for the equation of state at tracer points
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
-  integer :: i, j, k
+  integer :: i, j, k, m
+  real :: T5(5), S5(5) ! Temperatures and salinities at five quadrature points [C ~> degC] and [S ~> ppt]
+  real :: p5(5)        ! Pressures at five quadrature points [R L2 T-2 ~> Pa]
+  real :: r5(5)        ! Densities at five quadrature points [R ~> kg m-3]
   real, parameter :: C1_6 = 1.0/6.0    ! A rational constant [nondim]
+  real, parameter :: C1_90 = 1.0/90.0  ! A rational constant [nondim]
   real, parameter :: C1_12 = 1.0/12.0  ! A rational constant [nondim]
+  real :: wt_R       ! A weighting factor [nondim]
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   nkmb=GV%nk_rho_varies
@@ -945,7 +951,37 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
             0.25*((rho_top(i+1,j)+rho_top(i,j))-2.0*rho_ref) * dz_geo_sfc**2) then
           ! The pressure difference is at least half the size of the difference expected by hydrostatic
           ! balance.  This test gets rid of pressure differences that are small, e.g. open ocean.
-          intx_pa_cor(I,j) = C1_12 * (rho_top(i+1,j)-rho_top(i,j)) * dz_geo_sfc
+          if (CS%correction_intxpa_5pt) then
+            !! Use 5 point quadrature to calculate intxpa 
+            T5(1) = T_t(I,j,1) ; T5(5) = T_t(I+1,j,1)
+            S5(1) = S_t(I,j,1) ; S5(5) = S_t(I+1,j,1)
+            ! Pressure input to density EOS should be real pressure not rho_ref, I think
+            p5(1) = pa(I,j,1) - (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref)
+            p5(5) = pa(I+1,j,1) - (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref)
+            do m=2,4
+              wt_R =  0.25*real(m-1)
+              T5(m) = T5(1)+(T5(5)-T5(1))*wt_R !Quadratic: + (T5(5)-T5(1))*B*wt_R*(wt_R-1);
+              S5(m) = S5(1)+(S5(5)-S5(1))*wt_R !+ (S5(5)-S5(1))*B*wt_R*(wt_R-1);
+              p5(m) = p5(1)+(p5(5)-p5(1))*wt_R
+            enddo !m
+            call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)
+            ! add rhoref back in
+            do m=1,5
+              p5(m) = p5(m) + (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref)
+            enddo
+            do m=2,4
+              ! Make pressure curvature a difference from the linear fit of pressure between the two points
+              ! Do this by integrating pressure between each of the 5 points and adding up
+              ! This way integration direction doesn't matter when adding up pressure from previous point
+              p5(m) = p5(m-1) + ((0.25*(p5(5)-p5(1)) + 0.125*(r5(5)+r5(1))*dz_geo_sfc) - &
+                                 0.125*(r5(m)+r5(m-1))*dz_geo_sfc)
+            enddo
+            intx_pa(I,j,1) = C1_90*(7.0*(p5(1)+p5(5)) + 32.0*(p5(2)+p5(4)) + 12.0*p5(3))
+            ! Get correction from difference between this and linear average. This is clunky and repetitive.
+            intx_pa_cor(I,j) = -0.5*(pa(i,j,1) + pa(i+1,j,1)) + intx_pa(I,j,1)
+          else ! Do not use 5-point quadrature.
+            intx_pa_cor(I,j) = C1_12 * (rho_top(i+1,j)-rho_top(i,j)) * dz_geo_sfc
+          endif
         endif
       endif
       intx_pa(I,j,1) = 0.5*(pa(i,j,1) + pa(i+1,j,1)) + intx_pa_cor(I,j)
@@ -960,7 +996,37 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
             0.25*((rho_top(i,j+1)+rho_top(i,j))-2.0*rho_ref) * dz_geo_sfc**2) then
           ! The pressure difference is at least half the size of the difference expected by hydrostatic
           ! balance.  This test gets rid of pressure differences that are small, e.g. open ocean.
-          inty_pa_cor(i,J) = C1_12 * (rho_top(i,j+1)-rho_top(i,j)) * dz_geo_sfc
+           if (CS%correction_intxpa_5pt) then
+            !! Use 5 point quadrature to calculate intxpa
+            T5(1) = T_t(I,j,1) ; T5(5) = T_t(i,j+1,1)
+            S5(1) = S_t(I,j,1) ; S5(5) = S_t(i,j+1,1)
+            ! Pressure input to density EOS should be real pressure not rho_ref, I think
+            p5(1) = pa(i,j,1) - (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref)
+            p5(5) = pa(i,j+1,1) - (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref)
+            do m=2,4
+              wt_R =  0.25*real(m-1)
+              T5(m) = T5(1)+(T5(5)-T5(1))*wt_R !Quadratic: + (T5(5)-T5(1))*B*wt_R*(wt_R-1);
+              S5(m) = S5(1)+(S5(5)-S5(1))*wt_R !+ (S5(5)-S5(1))*B*wt_R*(wt_R-1);
+              p5(m) = p5(1)+(p5(5)-p5(1))*wt_R
+            enddo !m
+            call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)
+            ! add rhoref back in
+            do m=1,5
+              p5(m) = p5(m) + (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref)
+            enddo
+            do m=2,4
+              ! Make pressure curvature a difference from the linear fit of pressure between the two points
+              ! Do this by integrating pressure between each of the 5 points and adding up
+              ! This way integration direction doesn't matter when adding up pressure from previous point
+              p5(m) = p5(m-1) + ((0.25*(p5(5)-p5(1)) + 0.125*(r5(5)+r5(1))*dz_geo_sfc) - &
+                                 0.125*(r5(m)+r5(m-1))*dz_geo_sfc)
+            enddo
+            inty_pa(I,j,1) = C1_90*(7.0*(p5(1)+p5(5)) + 32.0*(p5(2)+p5(4)) + 12.0*p5(3))
+            ! Get correction from difference between this and linear average. This is clunky and repetitive.
+            inty_pa_cor(I,j) = -0.5*(pa(i,j,1) + pa(i,j+1,1)) + inty_pa(I,j,1)
+          else ! Do not use 5-point quadrature.
+            inty_pa_cor(i,J) = C1_12 * (rho_top(i,j+1)-rho_top(i,j)) * dz_geo_sfc
+          endif
         endif
       endif
       inty_pa(i,J,1) = 0.5*(pa(i,j,1) + pa(i,j+1,1)) + inty_pa_cor(i,J)
@@ -1226,6 +1292,9 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
   call get_param(param_file, mdl, "CORRECTION_INTXPA",CS%correction_intxpa, &
                  "If true, use a correction for surface pressure curvature in intx_pa.", &
                  default = .false.)
+  call get_param(param_file, mdl, "CORRECTION_INTXPA_5PT",CS%correction_intxpa_5pt, &
+                 "If true, use 5point quadrature to calculate intxpa. This requires "//&
+                 "CORRECTION_INTXPA = True.",default = .false.)
   call get_param(param_file, mdl, "USE_INACCURATE_PGF_RHO_ANOM", CS%use_inaccurate_pgf_rho_anom, &
                  "If true, use a form of the PGF that uses the reference density "//&
                  "in an inaccurate way. This is not recommended.", default=.false.)