+Add RESET_INTXPA_INTEGRAL

  Add RESET_INTXPA_INTEGRAL which resets intxpa and intypa at a trusted cell in
the interior (non-vanished and non-MWIPG-affected), and then integrates both up
and down to update intxpa and intypa for the interfaces above and below.  This
also adds the new runtime parameters RESET_INTXPA_H_NONVANISHED and
RESET_INTXPA_MASS_WEIGHT_CLEARANCE to determine when a cell is trusted.
This option is recommended with MASS_WEIGHT_IN_PRESSURE_GRADIENT_IS
for a quiet zstar ice shelf.  By default, this option is not on and no answers
change, but there are new parameters in some MOM_parameter_doc files.

src/core/MOM_PressureForce_FV.F90

@@ -48,8 +48,14 @@ module MOM_PressureForce_FV
                             !! timing of diagnostic output.
   logical :: useMassWghtInterp !< Use mass weighting in T/S interpolation
   logical :: useMassWghtInterpis !< Use mass weighting in T/S interpolation for top boundary
-  logical :: correction_intxpa ! Use correction to surface intxpa
-  logical :: correction_intxpa_5pt ! Use 5 point quadrature to calculate surface intxpa
+  logical :: correction_intxpa !< Use correction to surface intxpa
+  logical :: correction_intxpa_5pt !< Use 5 point quadrature to calculate surface intxpa
+  logical :: reset_intxpa_integral !< In the interior, reset intxpa at a trusted cell (for ice shelf)
+  real    :: h_nonvanished  !< A minimal layer thickness that indicates that a layer is thick enough
+                            !! to usefully reestimate the pressure integral across the interface
+                            !! below it [H ~> m or kg m-2]
+  real    :: dz_MWIPG       !< A vertical distance between in an interface and the depressed sea surface height
+                            !! in a neighboring cell beyond which the mass-weighting is enabled [Z ~> m]
   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.
@@ -577,14 +583,15 @@ 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, m
+  integer :: i, j, k, m, k2
   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]
+  real :: rho_tr, rho_tl ! Store right and left densities in reset intxpa calculation [R ~> kg m-3]
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   nkmb=GV%nk_rho_varies
@@ -826,10 +833,10 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
           ! 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.
           if (CS%correction_intxpa_5pt) then
-            !! Use 5 point quadrature to calculate intxpa 
+            ! 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
+            ! Pressure input to density EOS is the actual pressure not adjusted for rho_ref.
             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
@@ -869,30 +876,30 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
                                                ((-e(i,j+1,1)+e(i,j,1))*GV%g_Earth)) ) 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.
-           if (CS%correction_intxpa_5pt) then
-            !! Use 5 point quadrature to calculate intxpa
+          if (CS%correction_intxpa_5pt) then
+            ! Use 5 point quadrature to calculate intypa
             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
+            ! Pressure input to density EOS is the actual pressure not adjusted for rho_ref.
             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
+              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))*GV%g_Earth*(e(i,j+1,1)-e(i,j,1)) - &
-                   (0.125*(r5(m)+r5(m-1))*GV%g_Earth*(e(i,j+1,1)-e(I,j,1))))
+              ! 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))*GV%g_Earth*(e(i,j+1,1)-e(i,j,1)) - &
+                      (0.125*(r5(m)+r5(m-1))*GV%g_Earth*(e(i,j+1,1)-e(I,j,1))))
             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.
@@ -932,6 +939,62 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
     enddo ; enddo
   enddo
 
+  ! Having stored the pressure gradient info, we can work out where the first nonvanished layers is
+  ! reset intxpa there, then adjust intxpa above and below using the same increments between interfaces as above.
+  ! note: currently assumes pressure varies quadratically between bottom of first nonvanished, nonMWIPG level
+  ! possibly 5 pt quadrature should be implemented as for the surface
+  if (CS%reset_intxpa_integral) then
+    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      kloop: do k=1,nz-1
+        ! Check if both sides are nonvanished and MWIPG is not activated
+        if ((h(i,j,k)>CS%h_nonvanished) .and. (h(i+1,j,k)>CS%h_nonvanished)) then
+          if (.not. (max(0., e(i+1,j,K+1)-e(i,j,1), e(i,j,K+1)-e(i+1,j,1)) > CS%dz_MWIPG)) then
+            ! Calculate pressure at the bottom of this cell (pa are known)
+            ! then we have a "good estimate" for intxpa (it might have quadratic pressure dependence if sloped)
+            ! now we recalculate intx_pa and PFu at each level working up and then down
+            call calculate_density(T_t(i,j,k+1), S_t(i,j,k+1), pa(i,j,K+1), rho_tl, &
+                                   tv%eqn_of_state, rho_ref=rho_ref)
+            call calculate_density(T_t(i+1,j,k+1), S_t(i+1,j,k+1), pa(i+1,j,K+1), rho_tr, &
+                                   tv%eqn_of_state, rho_ref=rho_ref)
+            correction_x(I,j) = C1_12 * (rho_tr-rho_tl)*GV%g_Earth * (e(i+1,j,K+1)-e(i,j,K+1))
+            intx_pa(i,j,K+1) = 0.5*(pa(i,j,K+1) + pa(i+1,j,K+1)) + correction_x(I,j)
+            do k2=1,k
+              intx_pa(I,j,K-K2+1) = intx_pa(I,j,(K-K2+2)) - intx_dpa(i,j,k-k2+1)
+            enddo
+            do k2=k+2,nz
+              intx_pa(I,j,K2) = intx_pa(I,j,K2-1) + intx_dpa(i,j,k2-1)
+            enddo
+            exit kloop
+        endif ; endif
+      enddo kloop
+    enddo ; enddo
+
+    do J=Jsq,Jeq+1 ; do i=is,ie+1
+      kloop2: do k=1,nz-1
+        ! Check if both sides are nonvanished and MWIPG is not activated
+        if ((h(i,j,k)>CS%h_nonvanished) .and. (h(i,j+1,k)>CS%h_nonvanished)) then
+          if (.not. (max(0., e(i,j+1,K+1)-e(i,j,1), e(i,j,K+1)-e(i,j+1,1)) > CS%dz_MWIPG)) then
+            ! calculate pressure at the bottom of this cell (pa are known)
+            ! then we have a "good estimate" for intxpa (it might have quadratic pressure dependence if sloped)
+            ! now we recalculate intx_pa and PFu at each level working up and then down
+            call calculate_density(T_t(i,j,k+1), S_t(i,j,k+1), pa(i,j,K+1), rho_tl, &
+                                   tv%eqn_of_state, rho_ref=rho_ref)
+            call calculate_density(T_t(i,j+1,k+1), S_t(i,j+1,k+1), pa(i,j+1,K+1), rho_tr, &
+                                   tv%eqn_of_state, rho_ref=rho_ref)
+            correction_y(i,J) = C1_12 * (rho_tr-rho_tl) * GV%g_Earth * (e(i,j+1,K+1)-e(i,j,K+1))
+            inty_pa(i,J,K+1) = 0.5*(pa(i,j,K+1) + pa(i,j+1,K+1)) + correction_y(i,J)
+            do k2=1,k
+              inty_pa(i,J,K-K2+1) = inty_pa(i,J,(K-K2+2)) - inty_dpa(i,J,k-k2+1)
+            enddo
+            do k2=k+2,nz
+              inty_pa(i,J,K2) = inty_pa(i,J,K2-1) + inty_dpa(i,J,k2-1)
+            enddo
+            exit kloop2
+        endif ; endif
+      enddo kloop2
+    enddo ; enddo
+  endif ! intx_pa and inty_pa are now reset and should be correct
+
   ! Compute pressure gradient in x direction
   !$OMP parallel do default(shared)
   do k=1,nz ; do j=js,je ; do I=Isq,Ieq
@@ -1140,9 +1203,20 @@ 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, &
+  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, "RESET_INTXPA_INTEGRAL", CS%reset_intxpa_integral, &
+                 "If true, reset INTXPA to match pressures at first nonvanished cell. "//&
+                 "Includes pressure correction. ", default = .false.)
+  call get_param(param_file, mdl, "RESET_INTXPA_H_NONVANISHED", CS%h_nonvanished, &
+                 "A minimal layer thickness that indicates that a layer is thick enough to usefully "//&
+                 "reestimate the pressure integral across the interface below.", &
+                 default=1.0e-6, units="m", scale=GV%m_to_H, do_not_log=.not.CS%reset_intxpa_integral)
+  call get_param(param_file, mdl, "RESET_INTXPA_MASS_WEIGHT_CLEARANCE", CS%dz_MWIPG, &
+                 "A vertical distance between in an interface and the depressed sea surface height "//&
+                 "in a neighboring cell beyond which the mass-weighting is enabled.", &
+                 default=1.0e-8, units="m", scale=US%m_to_Z, do_not_log=.not.CS%reset_intxpa_integral)
   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.)