+*Fix 3-equation ice-ocean flux iteration

  Fix the 3-equation iteration for the buoyancy flux between the ocean and an
overlying ice-shelf when ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX is true and
SHELF_3EQ_GAMMA it false.  This code now uses proper bounding of the
self-consistent solution, avoiding further amplifying the fluxes in the cases
when the differences between the diffusivities of heat and salt to make the
buoyancy flux destabilizing for finite turbulent mixing.  Both the
false-position iterations and the (appropriately chosen) Newton's method
iterations have been extensively examined and determined to be working correctly
via print statements that have subsequently been removed for efficiency.

  Previously, the code to determine the 3-equation solution for the buoyancy
flux between the ocean and an ice shelf had been skipping iteration altogether
or doing un-bounded Newton's method iterations with a sign error in part of the
derivative, including taking the square root of negative numbers, leading to the
issue described at #945.  That issue has
now been corrected and can be closed once this commit has been merged into
the dev/gfdl branch of MOM6.

  This commit also changes the names of the runtime parameters to correct the
ice shelf flux iteration bugs from ICE_SHELF_BUOYANCY_FLUX_ITT_BUG and
ICE_SHELF_SALT_FLUX_ITT_BUG to ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX and
ICE_SHELF_SALT_FLUX_ITT_BUGFIX to avoid confusion with other ..._BUG parameters
where `true` is to turn the bugs on, whereas here `true` fixes them.  The old
names are retained via `old_name` arguments to the `get_param()` calls, so no
existing configurations will be disrupted by these changes.

  Additionally, an expression to determine a scaling factor to limit ice-shelf
bottom slopes in `calc_shelf_driving_stress()` was refactored to avoid the
possibility of division by zero.

  This commit will change (and correct) answers for cases with
ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX set to true, but as these would often fail
with a NaN from taking the square root of a negative value, it is very unlikely
that any such configurations are actively being used, and there seems little
point in retaining the previous answers.  No answers are changed in cases that
do not use an active ice shelf.

Author: Hallberg-NOAA

src/ice_shelf/MOM_ice_shelf.F90

@@ -197,8 +197,8 @@ module MOM_ice_shelf
                                          !! divided by the von Karman constant VK. Was 1/8.
   real :: Vk                             !< Von Karman's constant - dimensionless
   real :: Rc                             !< critical flux Richardson number.
-  logical :: buoy_flux_itt_bug           !< If true, fixes buoyancy iteration bug
-  logical :: salt_flux_itt_bug           !< If true, fixes salt iteration bug
+  logical :: buoy_flux_itt_bugfix        !< If true, fixes buoyancy iteration bug
+  logical :: salt_flux_itt_bugfix        !< If true, fixes salt iteration bug
   real :: buoy_flux_itt_threshold        !< Buoyancy iteration threshold for convergence
 
   !>@{ Diagnostic handles
@@ -298,7 +298,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
   real :: ZETA_N   !< This is the stability constant xi_N = 0.052 from Holland & Jenkins '99
                    !! divided by the von Karman constant VK. Was 1/8. [nondim]
   real :: RC       !< critical flux Richardson number.
-  real :: I_ZETA_N !< The inverse of ZETA_N [nondim].
+  real :: I_2Zeta_N !< Half the inverse of Zeta_N [nondim].
   real :: I_LF     !< The inverse of the latent heat of fusion [Q-1 ~> kg J-1].
   real :: I_VK     !< The inverse of the Von Karman constant [nondim].
   real :: PR, SC   !< The Prandtl number and Schmidt number [nondim].
@@ -327,30 +327,38 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
   real :: dS_ustar ! The difference between the salinity at the ice-ocean interface and the ocean
                    ! boundary layer salinity times the friction velocity [S Z T-1 ~> ppt m s-1]
   real :: ustar_h  ! The friction velocity in the water below the ice shelf [Z T-1 ~> m s-1]
-  real :: Gam_turb ! [nondim]
+  real :: Gam_turb ! A relative turbluent diffusivity [nondim]
   real :: Gam_mol_t, Gam_mol_s ! Relative coefficients of molecular diffusivities [nondim]
   real :: RhoCp     ! A typical ocean density times the heat capacity of water [Q R C-1 ~> J m-3 degC-1]
-  real :: ln_neut
+  real :: ln_neut   ! The log of the ratio of the neutral boundary layer thickness to the molecular
+                    ! boundary layer thickness if it is greater than 1 or 0 otherwise [nondim]
   real :: mass_exch ! A mass exchange rate [R Z T-1 ~> kg m-2 s-1]
   real :: Sb_min, Sb_max ! Minimum and maximum boundary salinities [S ~> ppt]
   real :: dS_min, dS_max ! Minimum and maximum salinity changes [S ~> ppt]
   ! Variables used in iterating for wB_flux.
-  real :: wB_flux_new, dDwB_dwB_in
-  real :: I_Gam_T, I_Gam_S
+  real :: wB_flux_next ! The next interation's guess for wB_flux [Z2 T-3 ~> m2 s-2]
+  real :: wB_flux_new  ! An updated value of wB_flux when Gam_turb is based on wB_flux [Z2 T-3 ~> m2 s-2]
+  real :: wB_flux_max  ! The upper bound on wB_flux [Z2 T-3 ~> m2 s-2]
+  real :: wB_flux_min  ! The lower bound on wB_flux [Z2 T-3 ~> m2 s-2]
+  real :: dDwB_dwB_in  ! The slope of the change in wB_flux between iterations with wB_flux [nondim]
+  real :: DwB_max      ! The change in wB_flux when it is wB_flux_max [Z2 T-3 ~> m2 s-2]
+  real :: DwB_min      ! The change in wB_flux when it is wB_flux_min [Z2 T-3 ~> m2 s-2]
+  real :: I_Gam_T, I_Gam_S  ! Terms that vary inversely with Gam_mol_T or Gam_mol_S and Gam_turb [nondim]
   real :: dG_dwB   ! The derivative of Gam_turb with wB [T3 Z-2 ~> s3 m-2]
   real :: taux2, tauy2 ! The squared surface stresses [R2 L2 Z2 T-4 ~> Pa2].
   real :: u2_av, v2_av ! The ice-area weighted average squared ocean velocities [L2 T-2 ~> m2 s-2]
-  real :: asu1, asu2   ! Ocean areas covered by ice shelves at neighboring u-
-  real :: asv1, asv2   ! and v-points [L2 ~> m2].
+  real :: asu1, asu2   ! Ocean areas covered by ice shelves at neighboring u-points [L2 ~> m2]
+  real :: asv1, asv2   ! Ocean areas covered by ice shelves at neighboring v-points [L2 ~> m2]
   real :: I_au, I_av   ! The Adcroft reciprocals of the ice shelf areas at adjacent points [L-2 ~> m-2]
   real :: Irho0        ! The inverse of the mean density times a unit conversion factor [R-1 L Z-1 ~> m3 kg-1]
   logical :: Sb_min_set, Sb_max_set
+  logical :: root_found
   logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
   logical :: coupled_GL     ! If true, the grounding line position is determined based on
                             ! coupled ice-ocean dynamics.
 
   real, parameter :: c2_3 = 2.0/3.0
-  character(len=160) :: mesg  ! The text of an error message
+  character(len=320) :: mesg  ! The text of an error message
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
   integer :: i, j, is, ie, js, je, ied, jed, it1, it3
   real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation [Z L2 ~> m3]
@@ -395,7 +403,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
   ZETA_N = CS%Zeta_N
   VK = CS%Vk
   RC = CS%Rc
-  I_ZETA_N = 1.0 / ZETA_N
+  I_2Zeta_N = 0.5 / CS%Zeta_N
   I_LF = 1.0 / CS%Lat_fusion
   SC = CS%kv_molec/CS%kd_molec_salt
   PR = CS%kv_molec/CS%kd_molec_temp
@@ -556,68 +564,150 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
             dT_ustar = (ISS%tfreeze(i,j) - sfc_state%sst(i,j)) * ustar_h
             dS_ustar = (Sbdry(i,j) - sfc_state%sss(i,j)) * ustar_h
 
-            ! First, determine the buoyancy flux assuming no effects of stability
-            ! on the turbulence.  Following H & J '99, this limit also applies
-            ! when the buoyancy flux is destabilizing.
-
-            if (CS%const_gamma) then ! if using a constant gamma_T
-              ! note the different form, here I_Gam_T is NOT 1/Gam_T!
+            if (CS%const_gamma) then
+              ! If using a constant gamma_T, there are no effects of the buoyancy flux on the turbulence.
               I_Gam_T = CS%Gamma_T_3EQ
               I_Gam_S = CS%Gamma_S_3EQ
-            else
-              Gam_turb = I_VK * (ln_neut + (0.5 * I_ZETA_N - 1.0))
+              wT_flux = dT_ustar * CS%Gamma_T_3EQ
+              wB_flux = dB_dS * (dS_ustar * CS%Gamma_S_3EQ) + dB_dT * wT_flux
+            else  ! gamma_T and gamma_S are a function of the buoyancy flux.
+              ! Find the root where wB_flux is consistent with the values of gamma with that flux.
+
+              ! First, determine the buoyancy flux assuming no effects of stability
+              ! on the turbulence.  Following H & J '99, this limit also applies
+              ! when the buoyancy flux is destabilizing.
+              Gam_turb = I_VK * (ln_neut + (I_2Zeta_N - 1.0))
               I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb)
               I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb)
-            endif
+              wB_flux = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * (dT_ustar * I_Gam_T)
 
-            wT_flux = dT_ustar * I_Gam_T
-            wB_flux = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * wT_flux
+              if (wB_flux < 0.0) then
+                ! The buoyancy flux is stabilizing and will reduce the turbulent
+                ! fluxes, and iteration is required.
 
-            if (wB_flux < 0.0) then
-              ! The buoyancy flux is stabilizing and will reduce the turbulent
-              ! fluxes, and iteration is required.
-              n_star_term = (ZETA_N * hBL_neut * VK) / (RC * ustar_h**3)
-              do it3 = 1,30
-               ! n_star <= 1.0 is the ratio of working boundary layer thickness
-               ! to the neutral thickness.
-               ! hBL = n_star*hBL_neut ; hSub = 1/8*n_star*hBL
+                ! n_star <= 1.0 is the ratio of working boundary layer thickness
+                ! to the neutral thickness.
+                ! hBL = n_star*hBL_neut ; hSub = 1/8*n_star*hBL
+                n_star_term = (ZETA_N * hBL_neut * VK) / (RC * ustar_h**3)
 
                 I_n_star = sqrt(1.0 - n_star_term * wB_flux)
-                dIns_dwB = 0.5 * n_star_term / I_n_star
                 if (hBL_neut_h_molec > I_n_star**2) then
-                  Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + &
-                                    (0.5*I_ZETA_N*I_n_star - 1.0))
-                  dG_dwB =  I_VK * ( -2.0 / I_n_star + (0.5 * I_ZETA_N)) * dIns_dwB
-                else
-                  !   The layer dominated by molecular viscosity is smaller than
-                  ! the assumed boundary layer.  This should be rare!
-                  Gam_turb = I_VK * (0.5 * I_ZETA_N*I_n_star - 1.0)
-                  dG_dwB = I_VK * (0.5 * I_ZETA_N) * dIns_dwB
+                  Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0))
+                else !   The layer dominated by molecular viscosity is smaller than the boundary layer.
+                  Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0)
                 endif
-
-                if (CS%const_gamma) then ! if using a constant gamma_T
-                  ! note the different form, here I_Gam_T is NOT 1/Gam_T!
-                  I_Gam_T = CS%Gamma_T_3EQ
-                  I_Gam_S = CS%Gamma_S_3EQ
-                else
-                  I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb)
-                  I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb)
+                I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb)
+                I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb)
+
+                if (CS%buoy_flux_itt_bugfix) then
+                  wB_flux_new = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * (dT_ustar * I_Gam_T)
+                  root_found = (abs(wB_flux_new - wB_flux) < &
+                                CS%buoy_flux_itt_threshold*(abs(wB_flux_new) + abs(wB_flux)))
+                  ! Do not update the flux if its maagnitude would be increased by the otherwise
+                  ! stabilizing buoyancy fluxes.  This can happen when the buoyancy flux
+                  ! is stabilizing when one of the heat or salt fluxes are destabilizing due
+                  ! to their different molecular properties.
+                  if (wB_flux_new <= wB_flux) root_found = .true.
+                  ! With the line above, the root is always bracketed.
+                else  ! The bug is never to update the buoyancy flux.
+                  root_found = .true.
                 endif
 
-                wT_flux = dT_ustar * I_Gam_T
-                wB_flux_new = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * wT_flux
+                if (.not.root_found) then
+                  wB_flux_max = 0.0 ; DwB_max = wB_flux
+                  wB_flux_min = wB_flux ; DwB_min = wB_flux_new - wB_flux
+                  ! The test about 10 lines above ensures that that (DwB_max <= 0) and (DwB_min >= 0),
+
+                  if (wB_flux_min*n_star_term <= (1.0 - hBL_neut_h_molec)) then
+                    ! The solution is in the limit where abs(wB_flux) is small and
+                    ! Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0))
+                  elseif (wB_flux_max*n_star_term >= (1.0 - hBL_neut_h_molec)) then
+                    ! The solution is in the limit where abs(wB_flux) is large and
+                    ! Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0)
+                  else
+                    ! The derivative of Gam_turb with wB_flux has a discontinuous change within the
+                    ! bracketed range of values.  For a first guess, take this discontinous slope
+                    ! value, as Newton's method and the false position method do not converge well
+                    ! when this discontinuity is between a guess and the solution.
+
+                    ! Set wB_flux so that I_n_star**2 = hBL_neut_h_molec
+                    ! I_n_star = sqrt(1.0 - n_star_term * wB_flux)
+                    wB_flux = (1.0 - hBL_neut_h_molec) /  n_star_term
+                    I_n_star = sqrt(hBL_neut_h_molec)
+                    Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0)
+                    I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb)
+                    I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb)
+
+                    wB_flux_new = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * (dT_ustar * I_Gam_T)
+                    if (abs(wB_flux_new - wB_flux) <= CS%buoy_flux_itt_threshold*(abs(wB_flux_new) + abs(wB_flux))) then
+                      ! The root has been found to within the tolerance at the kink.  This should be very rare.
+                      root_found = .true.
+                    elseif (wB_flux_new > wB_flux) then
+                      ! The solution is in the limit where abs(wB_flux) is small and
+                      ! Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0))
+                      wB_flux_min = wB_flux ; DwB_min = wB_flux_new - wB_flux
+                    else
+                      ! The solution is in the limt where abs(wB_flux) is large and
+                      ! Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0)
+                      wB_flux_max = wB_flux ; DwB_max = wB_flux_new - wB_flux
+                    endif
+                  endif
+                endif
 
-                ! Find the root where wB_flux_new = wB_flux.
-                if (abs(wB_flux_new - wB_flux) < CS%buoy_flux_itt_threshold*(abs(wB_flux_new) + abs(wB_flux))) exit
+                if (.not.root_found) then
+                  ! Use the false position for the next guess.
+                  wB_flux = wB_flux_min + (wB_flux_max-wB_flux_min) * (DwB_min / (DwB_min - DwB_max))
+
+                  do it3 = 1,30
+                  ! Iterate using Newton's method with bounds or the false position method to find the root.
+
+                    I_n_star = sqrt(1.0 - n_star_term * wB_flux)
+                    dIns_dwB = -0.5 * n_star_term / I_n_star
+                    if (hBL_neut_h_molec > I_n_star**2) then
+                      Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0))
+                      dG_dwB =  I_VK * (( -2.0 / I_n_star + I_2Zeta_N) * dIns_dwB)
+                    else
+                      !   The layer dominated by molecular viscosity is smaller than the boundary layer.
+                      Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0)
+                      dG_dwB = I_VK * (I_2Zeta_N * dIns_dwB)
+                    endif
+                    I_Gam_T = 1.0 / (Gam_mol_t + Gam_turb)
+                    I_Gam_S = 1.0 / (Gam_mol_s + Gam_turb)
+                    wB_flux_new = dB_dS * (dS_ustar * I_Gam_S) + dB_dT * (dT_ustar * I_Gam_T)
+
+                    ! Test for convergence to within tolerance at the point where wB_flux_new = wB_flux.
+                    if (abs(wB_flux_new - wB_flux) <= CS%buoy_flux_itt_threshold*(abs(wB_flux_new) + abs(wB_flux))) &
+                      root_found = .true.
+                    if (root_found) exit
+
+                    dDwB_dwB_in = -dG_dwB * (dB_dS * (dS_ustar * I_Gam_S**2) + &
+                                             dB_dT * (dT_ustar * I_Gam_T**2)) - 1.0
+                    if ((dDwB_dwB_in >= 0.0) .or. &
+                        ( wB_flux - wB_flux_new >= abs(dDwB_dwB_in)*(wB_flux_max - wB_flux)) .or. &
+                        ( wB_flux - wB_flux_new <= abs(dDwB_dwB_in)*(wB_flux_min - wB_flux)) ) then
+                      ! Use the False position method to determine the guess for the next iteration when
+                      ! Newton's method would go out of bounds
+                      wB_flux_next = wB_flux_min + (wB_flux_max-wB_flux_min) * (DwB_min / (DwB_min - DwB_max))
+                    else
+                      ! Use Newton's method for the next guess.
+                      wB_flux_next = wB_flux - (wB_flux_new - wB_flux) / dDwB_dwB_in
+                    endif
+
+                    ! Reset one of the bounds inward.
+                    if (wB_flux_new - wB_flux > 0) then
+                      wB_flux_min = wB_flux ; DwB_min = wB_flux_new - wB_flux
+                    else
+                      wB_flux_max = wB_flux ; DwB_max = wB_flux_new - wB_flux
+                    endif
+
+                    ! Update wB_flux
+                    wB_flux = wB_flux_next
+                  enddo ! it3
+                endif
 
-                dDwB_dwB_in = dG_dwB * (dB_dS * (dS_ustar * I_Gam_S**2) + &
-                                        dB_dT * (dT_ustar * I_Gam_T**2)) - 1.0
-                ! This is Newton's method without any bounds.  Should bounds be needed?
-                wB_flux_new = wB_flux - (wB_flux_new - wB_flux) / dDwB_dwB_in
-                ! Update wB_flux
-                if (CS%buoy_flux_itt_bug) wB_flux = wB_flux_new
-              enddo !it3
-            endif
+              endif  ! End of test for first guess of wB_flux < 0.
+              wT_flux = dT_ustar * I_Gam_T
+            endif  ! End of test for CS%const_gamma
 
             ISS%tflux_ocn(i,j)  = RhoCp * wT_flux
             exch_vel_t(i,j) = ustar_h * I_Gam_T
@@ -688,7 +778,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
                 Sbdry(i,j) = Sbdry_it
               endif ! Sb_min_set
 
-              if (.not.CS%salt_flux_itt_bug) Sbdry(i,j) = Sbdry_it
+              if (.not.CS%salt_flux_itt_bugfix) Sbdry(i,j) = Sbdry_it
 
             endif ! CS%find_salt_root
 
@@ -1686,10 +1776,10 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init,
   call get_param(param_file, mdl, "ICE_SHELF_RC", CS%Rc, &
                  "Critical flux Richardson number for ice melt ", &
                  units="nondim", default=0.20)
-  call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_BUG", CS%buoy_flux_itt_bug, &
-                 "Bug fix of buoyancy iteration", default=.true.)
-  call get_param(param_file, mdl, "ICE_SHELF_SALT_FLUX_ITT_BUG", CS%salt_flux_itt_bug, &
-                 "Bug fix of salt iteration", default=.true.)
+  call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX", CS%buoy_flux_itt_bugfix, &
+                 "Bug fix of buoyancy iteration", default=.true., old_name="ICE_SHELF_BUOYANCY_FLUX_ITT_BUG")
+  call get_param(param_file, mdl, "ICE_SHELF_SALT_FLUX_ITT_BUGFIX", CS%salt_flux_itt_bugfix, &
+                 "Bug fix of salt iteration", default=.true., old_name="ICE_SHELF_SALT_FLUX_ITT_BUG")
   call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_THRESHOLD", CS%buoy_flux_itt_threshold, &
                  "Convergence criterion of Newton's method for ice shelf "//&
                  "buoyancy iteration.", units="nondim", default=1.0e-4)

src/ice_shelf/MOM_ice_shelf_dynamics.F90

@@ -2540,7 +2540,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD)
         endif
 
         if (CS%max_surface_slope>0) then
-          scale = min(CS%max_surface_slope/sqrt((sx**2)+(sy**2)),1.0)
+          scale = CS%max_surface_slope / max( sqrt((sx**2) + (sy**2)), CS%max_surface_slope )
           sx = scale*sx; sy = scale*sy
         endif