Add ice shelf pressure initialisation bug fix (#800)

claireyung · web-flow · commit dd1f3f52f5fe · 2025-01-14T17:06:36.000-05:00
* Add ice shelf pressure initialisation bug fix

This commit adds a new parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX
and associated bug fix. This bug occurs when TRIM_IC_FOR_P_SURF
pressure initialisation is used for ice shelf, whilst in
Boussinesq mode and a nonlinear equation of state. The subroutine
trim_for_ice calls cut_off_column_top. This routine in Boussinesq
mode calls find_depth_of_pressure_in_cell in
MOM_density_integrals.F90. This subsequently calls the function
frac_dp_at_pos which calls the density equation of state based on
given salinity, temperature and depth, which is incorrectly converted
into pressure by z position (which is negative) multiplied by g and
rho0. The bug results in incorrect densities from the equation of state
and therefore an imperfect initialisation of layer thicknesses and sea
surface height due to the displacement of water by the ice.

The bug is fixed by multiplying the pressure by -1.

This commit introduces parameter FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX with
default value False to preserve answers. If true, the ice shelf
initialisation is accurate to a high degree with a nonlinear
equation of state.

Answers should not change, except for the added parameter in
MOM_parameter_doc.

The same functions are called by a unit test in MOM_state_initialization;
in this commit I set the MOM_state_initialization unit test to
use the existing bug (with a false flag).

* Resolve indenting and white space inconsitencies with MOM6 style
for ice shelf pressure initialisation bug fix FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX
diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90
@@ -2001,7 +2001,7 @@ end subroutine diagnose_mass_weight_p
 
 !> Find the depth at which the reconstructed pressure matches P_tgt
 subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_tgt, &
-                       rho_ref, G_e, EOS, US, P_b, z_out, z_tol)
+                       rho_ref, G_e, EOS, US, P_b, z_out, z_tol, frac_dp_bugfix)
   real,                  intent(in)  :: T_t !< Potential temperature at the cell top [C ~> degC]
   real,                  intent(in)  :: T_b !< Potential temperature at the cell bottom [C ~> degC]
   real,                  intent(in)  :: S_t !< Salinity at the cell top [S ~> ppt]
@@ -2020,6 +2020,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t
   real,                  intent(out) :: P_b !< Pressure at the bottom of the cell [R L2 T-2 ~> Pa]
   real,                  intent(out) :: z_out !< Absolute depth at which anomalous pressure = p_tgt [Z ~> m]
   real,                  intent(in)  :: z_tol !< The tolerance in finding z_out [Z ~> m]
+  logical,               intent(in)  :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos
 
   ! Local variables
   real :: dp    ! Pressure thickness of the layer [R L2 T-2 ~> Pa]
@@ -2032,7 +2033,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t
   GxRho = G_e * rho_ref
 
   ! Anomalous pressure difference across whole cell
-  dp = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, 1.0, EOS)
+  dp = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, 1.0, EOS, frac_dp_bugfix)
 
   P_b = P_t + dp ! Anomalous pressure at bottom of cell
 
@@ -2063,7 +2064,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t
      call MOM_error(FATAL, 'find_depth_of_pressure_in_cell completes too many iterations: '//msg)
     endif
     z_out = z_t + ( z_b - z_t ) * F_guess
-    Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS) - ( P_tgt - P_t )
+    Pa = frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, F_guess, EOS, frac_dp_bugfix) - ( P_tgt - P_t )
 
     if (Pa<Pa_left) then
       write(msg,*) Pa_left,Pa,Pa_right,P_t-P_tgt,P_b-P_tgt
@@ -2119,7 +2120,7 @@ end subroutine avg_specific_vol
 
 !> Returns change in anomalous pressure change from top to non-dimensional
 !! position pos between z_t and z_b [R L2 T-2 ~> Pa]
-real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS)
+real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EOS, frac_dp_bugfix)
   real,           intent(in)  :: T_t !< Potential temperature at the cell top [C ~> degC]
   real,           intent(in)  :: T_b !< Potential temperature at the cell bottom [C ~> degC]
   real,           intent(in)  :: S_t !< Salinity at the cell top [S ~> ppt]
@@ -2131,6 +2132,7 @@ real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EO
   real,           intent(in)  :: G_e !< The Earth's gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
   real,           intent(in)  :: pos !< The fractional vertical position, 0 to 1 [nondim]
   type(EOS_type), intent(in)  :: EOS !< Equation of state structure
+  logical,        intent(in)  :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos
 
   ! Local variables
   real, parameter :: C1_90 = 1.0/90.0  ! A rational constant [nondim]
@@ -2150,7 +2152,11 @@ real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EO
     ! Salinity and temperature points are linearly interpolated
     S5(n) = top_weight * S_t + bottom_weight * S_b
     T5(n) = top_weight * T_t + bottom_weight * T_b
-    p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref )
+    if (frac_dp_bugfix) then
+      p5(n) = (-1) * ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref )
+    else
+      p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref )
+    endif !bugfix
   enddo
   call calculate_density(T5, S5, p5, rho5, EOS)
   rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref
diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90
@@ -1205,6 +1205,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read)
                                   ! answers from 2018, while higher values use more robust
                                   ! forms of the same remapping expressions.
   logical :: use_remapping ! If true, remap the initial conditions.
+  logical :: use_frac_dp_bugfix   ! If true, use bugfix. Otherwise, pressure input to EOS is negative.
   type(remapping_CS), pointer :: remap_CS => NULL()
 
   call get_param(PF, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, &
@@ -1227,7 +1228,10 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read)
                  "The tolerance with which to find the depth matching the specified "//&
                  "surface pressure with TRIM_IC_FOR_P_SURF.", &
                  units="m", default=1.0e-5, scale=US%m_to_Z, do_not_log=just_read)
-
+  call get_param(PF, mdl, "FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX", use_frac_dp_bugfix, &
+                 "If true, use bugfix in ice shelf TRIM_IC initialization. "//&
+                 "Otherwise, pressure input to density EOS is negative.", &
+                 default=.false., do_not_log=just_read)
   call get_param(PF, mdl, "TRIMMING_USES_REMAPPING", use_remapping, &
                  'When trimming the column, also remap T and S.', &
                  default=.false., do_not_log=just_read)
@@ -1277,7 +1281,8 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read)
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
     call cut_off_column_top(GV%ke, tv, GV, US, GV%g_Earth, G%bathyT(i,j)+G%Z_ref, min_thickness, &
                tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), &
-               p_surf(i,j), h(i,j,:), remap_CS, z_tol=z_tolerance)
+               p_surf(i,j), h(i,j,:), remap_CS, z_tol=z_tolerance, &
+               frac_dp_bugfix=use_frac_dp_bugfix)
   enddo ; enddo
 
 end subroutine trim_for_ice
@@ -1368,7 +1373,7 @@ end subroutine calc_sfc_displacement
 !> Adjust the layer thicknesses by removing the top of the water column above the
 !! depth where the hydrostatic pressure matches p_surf
 subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, T_t, T_b, &
-                              S, S_t, S_b, p_surf, h, remap_CS, z_tol)
+                              S, S_t, S_b, p_surf, h, remap_CS, z_tol, frac_dp_bugfix)
   integer,               intent(in)    :: nk  !< Number of layers
   type(thermo_var_ptrs), intent(in)    :: tv  !< Thermodynamics structure
   type(verticalGrid_type), intent(in)  :: GV  !< The ocean's vertical grid structure.
@@ -1388,6 +1393,7 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T,
                                                    !! if associated
   real,                  intent(in)    :: z_tol !< The tolerance with which to find the depth
                                                 !! matching the specified pressure [Z ~> m].
+  logical,               intent(in)    :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos
 
   ! Local variables
   real, dimension(nk+1) :: e ! Top and bottom edge positions for reconstructions [Z ~> m]
@@ -1416,7 +1422,8 @@ subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T,
     do k=1,nk
       call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), &
                                           P_t, p_surf, GV%Rho0, G_earth, tv%eqn_of_state, &
-                                          US, P_b, z_out, z_tol=z_tol)
+                                          US, P_b, z_out, z_tol=z_tol, &
+                                          frac_dp_bugfix=frac_dp_bugfix)
       if (z_out>=e(K)) then
         ! Imposed pressure was less that pressure at top of cell
         exit
@@ -3139,7 +3146,8 @@ subroutine MOM_state_init_tests(G, GV, US, tv)
   P_t = 0.
   do k = 1, nk
     call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), P_t, 0.5*P_tot, &
-                                        GV%Rho0, GV%g_Earth, tv%eqn_of_state, US, P_b, z_out, z_tol=z_tol)
+                                        GV%Rho0, GV%g_Earth, tv%eqn_of_state, US, P_b, z_out, z_tol=z_tol, &
+                                        frac_dp_bugfix=.false.)
     write(0,*) k, US%RL2_T2_to_Pa*P_t, US%RL2_T2_to_Pa*P_b, 0.5*US%RL2_T2_to_Pa*P_tot, &
                US%Z_to_m*e(K), US%Z_to_m*e(K+1), US%Z_to_m*z_out
     P_t = P_b
@@ -3158,7 +3166,8 @@ subroutine MOM_state_init_tests(G, GV, US, tv)
   !                             h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff)
   ! endif
   call cut_off_column_top(nk, tv, GV, US, GV%g_Earth, -e(nk+1), GV%Angstrom_H, &
-                          T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS, z_tol=z_tol)
+                          T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS, z_tol=z_tol, &
+                          frac_dp_bugfix=.false.)
   write(0,*) GV%H_to_m*h(:)
   if (associated(remap_CS)) deallocate(remap_CS)
 
