Enable stochastic EOS with Bodner

docs/zotero.bib

@@ -1,5 +1,34 @@
 % This file is generated by zotero. Manual edits will be lost!
 
+@article{agarwal2023,
+  title={Impact of Stochastic Ocean Density Corrections on Air-Sea Flux Variability},
+  author={Agarwal, Niraj and Small, R Justin and Bryan, Frank O and Grooms, Ian and Pegion, Philip J},
+  journal={Geophysical Research Letters},
+  volume={50},
+  number={13},
+  pages={e2023GL104248},
+  year={2023},
+}
+
+@article{kenigson2022,
+  title={Parameterizing the impact of unresolved temperature variability on the large-scale density field: 2. Modeling},
+  author={Kenigson, JS and Adcroft, A and Bachman, SD and Castruccio, F and Grooms, I and Pegion, P and Stanley, Z},
+  journal={Journal of Advances in Modeling Earth Systems},
+  volume={14},
+  number={3},
+  pages={e2021MS002844},
+  year={2022},
+}
+
+@article{stanley2020,
+  title={Parameterizing the Impact of Unresolved Temperature Variability on the Large-Scale Density Field: {Part 1.} Theory.},
+  author={Stanley, Z and Grooms, I and Kleiber, W and Bachman, SD and Castruccio, F and Adcroft, A},
+  journal={Journal of Advances in Modeling Earth Systems},
+  volume={12},
+  number={12},
+  pages={e2020MS002185},
+  year={2020},
+}
 
 @article{redi1982,
 	title = {Oceanic {Isopycnal} {Mixing} by {Coordinate} {Rotation}},

src/core/MOM_PressureForce_FV.F90

@@ -2030,8 +2030,6 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, ADp, SAL
   type(tidal_forcing_CS), intent(in), target, optional :: tides_CSp !< Tides control structure
 
   ! Local variables
-  real :: Stanley_coeff    ! Coefficient relating the temperature gradient and sub-gridscale
-                           ! temperature variance [nondim]
   integer :: default_answer_date ! Global answer date
   logical :: use_temperature   ! If true, temperature and salinity are used as state variables.
   logical :: use_EOS           ! If true, density calculated from T & S using an equation of state.
@@ -2043,6 +2041,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, ADp, SAL
 # include "version_variable.h"
   character(len=40)  :: mdl  ! This module's name.
   logical :: use_ALE       ! If true, use the Vertical Lagrangian Remap algorithm
+  logical :: stoch_eos     ! Can't use Stanley param here unless stoch_eos is true
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
 
   isd  = G%isd  ; ied  = G%ied  ; jsd  = G%jsd  ; jed  = G%jed ; nz = GV%ke
@@ -2182,16 +2181,15 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, ADp, SAL
                  "boundary cells is extrapolated, rather than using PCM "//&
                  "in these cells. If true, the same order polynomial is "//&
                  "used as is used for the interior cells.", default=.true.)
+  call get_param(param_file, mdl, "STOCH_EOS", stoch_eos, &
+                 default=.false., do_not_log=.true.)
   call get_param(param_file, mdl, "USE_STANLEY_PGF", CS%use_stanley_pgf, &
                  "If true, turn on Stanley SGS T variance parameterization "// &
                  "in PGF code.", default=.false.)
   if (CS%use_stanley_pgf) then
-    call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, &
-                 "Coefficient correlating the temperature gradient and SGS T variance.", &
-                 units="nondim", default=-1.0, do_not_log=.true.)
-    if (Stanley_coeff < 0.0) call MOM_error(FATAL, &
-                 "STANLEY_COEFF must be set >= 0 if USE_STANLEY_PGF is true.")
-
+    if (.not.stoch_eos) then
+      call MOM_error(FATAL, "PressureForce_FV_init: USE_STANLEY_PGF requires STOCH_EOS")
+    endif
     CS%id_rho_pgf = register_diag_field('ocean_model', 'rho_pgf', diag%axesTL, &
         Time, 'rho in PGF', 'kg m-3', conversion=US%R_to_kg_m3)
     CS%id_rho_stanley_pgf = register_diag_field('ocean_model', 'rho_stanley_pgf', diag%axesTL, &

src/core/MOM_stoch_eos.F90

@@ -4,7 +4,7 @@ module MOM_stoch_eos
 ! This file is part of MOM6. See LICENSE.md for the license.
 use MOM_diag_mediator,    only : register_diag_field, post_data, diag_ctrl
 use MOM_error_handler,    only : MOM_error, FATAL
-use MOM_file_parser,      only : get_param, param_file_type
+use MOM_file_parser,      only : get_param, log_version, param_file_type
 use MOM_grid,             only : ocean_grid_type
 use MOM_hor_index,        only : hor_index_type
 use MOM_isopycnal_slopes, only : vert_fill_TS
@@ -62,27 +62,32 @@ logical function MOM_stoch_eos_init(Time, G, GV, US, param_file, diag, CS, resta
   type(MOM_restart_CS),    pointer       :: restart_CS !< A pointer to the restart control structure.
 
   ! local variables
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
   integer :: i,j
 
   MOM_stoch_eos_init = .false.
 
   CS%seed = 0
 
+  call log_version(param_file, "MOM_stoch_eos", version, "")
   call get_param(param_file, "MOM_stoch_eos", "STOCH_EOS", CS%use_stoch_eos, &
-                 "If true, stochastic perturbations are applied "//&
-                 "to the EOS in the PGF.", default=.false.)
+                 "If true, computes stochastic perturbations that can be applied "//&
+                 "to the EOS in various places.", default=.false.)
   call get_param(param_file, "MOM_stoch_eos", "STANLEY_COEFF", CS%stanley_coeff, &
                  "Coefficient correlating the temperature gradient "//&
                  "and SGS T variance.", units="nondim", default=-1.0)
+  if ((CS%stanley_coeff < 0.0) .and. CS%use_stoch_eos) call MOM_error(FATAL, &
+                 "STANLEY_COEFF must be set >= 0 if STOCH_EOS is true.")
   call get_param(param_file, "MOM_stoch_eos", "STANLEY_A", CS%stanley_a, &
                  "Coefficient a which scales chi in stochastic perturbation of the "//&
                  "SGS T variance.", units="nondim", default=1.0, &
-                 do_not_log=((CS%stanley_coeff<0.0) .or. .not.CS%use_stoch_eos))
+                 do_not_log=.not.CS%use_stoch_eos)
   call get_param(param_file, "MOM_stoch_eos", "KD_SMOOTH", CS%kappa_smooth, &
                  "A diapycnal diffusivity that is used to interpolate "//&
                  "more sensible values of T & S into thin layers.", &
                  units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T, &
-                 do_not_log=(CS%stanley_coeff<0.0))
+                 do_not_log=.not.CS%use_stoch_eos)
 
   ! Don't run anything if STANLEY_COEFF < 0
   if (CS%stanley_coeff >= 0.0) then
@@ -258,4 +263,30 @@ subroutine MOM_calc_varT(G, GV, US, h, tv, CS, dt)
   endif
 end subroutine MOM_calc_varT
 
+!> \namespace mom_stoch_eos
+!!
+!! This module provides the foundation of the Stanley parameterization (\cite stanley2020) for correcting the
+!! computation of density. Density is not a prognostic variable in MOM6; it is computed for various purposes
+!! in various places. The correction to this calculation provided by this module has been implemented
+!! in some places where density is used, but not all.
+!!
+!! To use the correction, first set <code>STOCH_EOS=True</code>. Then, choose the constant c from (25) of
+!! \cite stanley2020. This is controlled using <code>STANLEY_COEFF</code>. Setting a negative value will
+!! result in an error. \cite stanley2020 found a value of 0.2 offline, coarsening from 0.1 to 1 degree
+!! resolution. \cite kenigson2022 proposed a value of 0.5 in a 2/3 degree resolution model.
+!!
+!! Whether the correction is deterministic or stochastic can be controlled using the variable
+!! <code>STANLEY_A</code>. Setting this to 0.0 uses the deterministic version, while a value of 1.0 produces
+!! the stochastic version. Reducing from 1 to 0 smoothly transitions from stochastic to deterministic.
+!!
+!! To turn the correction on in various parts of the code, use
+!! - <code>USE_STANLEY_PGF=True</code> for the pressure gradient force (cf. \cite kenigson2022)
+!! - <code>USE_STANLEY_ISO=True</code> to correct the computation of isopycnal slopes (used in many places)
+!! - <code>USE_STANLEY_GM=True</code> to use the parameterization within GM (cf. \cite agarwal2023)
+!! - <code>USE_STANLEY_ML=True</code> to use the parameterization within the mixed-layer restratification
+!!   parameterization. It applies to both the OM4 and Bodner schemes. (cf. \cite agarwal2023)
+!!
+!! For ensemble simulations, the random number generator seed can be controlled using the parameter
+!! <code>SEED_STOCH_EOS</code>
+
 end module MOM_stoch_eos

src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

@@ -1335,9 +1335,8 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
                                   ! scaled by the resolution function.
   logical :: better_speed_est ! If true, use a more robust estimate of the first
                               ! mode wave speed as the starting point for iterations.
-  real :: Stanley_coeff    ! Coefficient relating the temperature gradient and sub-gridscale
-                           ! temperature variance [nondim]
   logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for calculating the EBT structure
+  logical :: stoch_eos        ! Can't use Stanley param here unless stoch_eos is true
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   character(len=40)  :: mdl = "MOM_lateral_mixing_coeffs" ! This module's name.
@@ -1461,15 +1460,13 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
 
   call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.)
 
+  call get_param(param_file, mdl, "STOCH_EOS", stoch_eos, &
+                 default=.false., do_not_log=.true.)
   call get_param(param_file, mdl, "USE_STANLEY_ISO", CS%use_stanley_iso, &
                  "If true, turn on Stanley SGS T variance parameterization "// &
                  "in isopycnal slope code.", default=.false.)
-  if (CS%use_stanley_iso) then
-    call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, &
-                 "Coefficient correlating the temperature gradient and SGS T variance.", &
-                 units="nondim", default=-1.0, do_not_log=.true.)
-    if (Stanley_coeff < 0.0) call MOM_error(FATAL, &
-                 "STANLEY_COEFF must be set >= 0 if USE_STANLEY_ISO is true.")
+  if (CS%use_Stanley_ISO .and. .not.stoch_eos) then
+    call MOM_error(FATAL, "VarMix_init: USE_STANLEY_ISO requires STOCH_EOS")
   endif
 
   if (CS%Resoln_use_ebt .or. CS%khth_use_ebt_struct .or. CS%kdgl90_use_ebt_struct &

src/parameterizations/lateral/MOM_mixed_layer_restrat.F90

@@ -863,6 +863,9 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
            "To use the Bodner et al., 2023, MLE parameterization, either MLE_USE_PBL_MLD or "// &
            "Bodner_detect_MLD must be True.")
   endif
+  if (CS%use_Stanley_ML .and. .not.GV%Boussinesq) call MOM_error(FATAL, &
+       "MOM_mixedlayer_restrat: The Stanley parameterization is not"//&
+       "available without the Boussinesq approximation.")
 
   if (associated(bflux)) &
     call pass_var(bflux, G%domain, halo=1)
@@ -1644,15 +1647,14 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
   real :: flux_to_kg_per_s ! A unit conversion factor for fluxes. [kg T s-1 H-1 L-2 ~> kg m-3 or 1]
   real :: omega            ! The Earth's rotation rate [T-1 ~> s-1].
   real :: ustar_min_dflt   ! The default value for RESTRAT_USTAR_MIN [Z T-1 ~> m s-1]
-  real :: Stanley_coeff    ! Coefficient relating the temperature gradient and sub-gridscale
-                           ! temperature variance [nondim]
   integer :: default_answer_date  ! The default setting for the various ANSWER_DATE flags
   ! This include declares and sets the variable "version".
   character(len=200) :: inputdir   ! The directory where NetCDF input files
   character(len=240) :: mle_fl_filename ! A file from which chl_a concentrations are to be read.
   character(len=128) :: mle_fl_file ! Data containing MLE front-length scale. Used
                                     ! when reading from file.
   character(len=32)  :: fl_varname ! Name of front-length scale variable in mle_fl_file.
+  logical :: stoch_eos     ! Can't use Stanley param here unless stoch_eos is true
 
 # include "version_variable.h"
   integer :: i, j
@@ -1689,6 +1691,14 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
       "This sets the default value for the various _ANSWER_DATE parameters.", &
       default=99991231, do_not_log=.true.)
   call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
+  call get_param(param_file, mdl, "STOCH_EOS", stoch_eos, &
+                 default=.false., do_not_log=.true.)
+  call get_param(param_file, mdl, "USE_STANLEY_ML", CS%use_Stanley_ML, &
+                 "If true, turn on Stanley SGS T variance parameterization "// &
+                 "in ML restrat code.", default=.false.)
+  if (CS%use_Stanley_ML .and. .not.stoch_eos) then
+    call MOM_error(FATAL, "mixedlayer_restrat_init: USE_STANLEY_ML requires STOCH_EOS")
+  endif
   call openParameterBlock(param_file,'MLE') ! Prepend MLE% to all parameters
   if (GV%nkml==0) then
     call get_param(param_file, mdl, "USE_BODNER23", CS%use_Bodner, &
@@ -1751,9 +1761,6 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
              "Fraction by which to extend the mixed-layer restratification "//&
              "depth used for a smoother stream function at the base of "//&
              "the mixed-layer.", units="nondim", default=0.0)
-    call get_param(param_file, mdl, "USE_STANLEY_TVAR", CS%use_Stanley_ML, &
-             "If true, turn on Stanley SGS T variance parameterization "// &
-             "in ML restrat code.", default=.false.)
     call get_param(param_file, mdl, "USE_CR_GRID", CS%Cr_grid, &
              "If true, read in a spatially varying Cr field.", default=.false.)
     call get_param(param_file, mdl, "USE_MLD_GRID", CS%MLD_grid, &
@@ -1811,17 +1818,6 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
              "geostrophic kinetic energy or 1 plus the square of the "//&
              "grid spacing over the deformation radius, as detailed "//&
              "by Fox-Kemper et al. (2011)", units="nondim", default=0.0)
-    ! These parameters are only used in the OM4-era version of Fox-Kemper
-    call get_param(param_file, mdl, "USE_STANLEY_ML", CS%use_Stanley_ML, &
-                   "If true, turn on Stanley SGS T variance parameterization "// &
-                   "in ML restrat code.", default=.false.)
-    if (CS%use_Stanley_ML) then
-      call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, &
-                   "Coefficient correlating the temperature gradient and SGS T variance.", &
-                   units="nondim", default=-1.0, do_not_log=.true.)
-      if (Stanley_coeff < 0.0) call MOM_error(FATAL, &
-               "STANLEY_COEFF must be set >= 0 if USE_STANLEY_ML is true.")
-    endif
     call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, &
                    'The value the von Karman constant as used for mixed layer viscosity.', &
                    units='nondim', default=0.41)

src/parameterizations/lateral/MOM_thickness_diffuse.F90

@@ -97,9 +97,6 @@ module MOM_thickness_diffuse
                                  !! When this is true, it breaks rotational symmetry.
   logical :: use_GM_work_bug     !< If true, use the incorrect sign for the
                                  !! top-level work tendency on the top layer.
-  real :: Stanley_det_coeff      !< The coefficient correlating SGS temperature variance with the mean
-                                 !! temperature gradient in the deterministic part of the Stanley parameterization.
-                                 !! Negative values disable the scheme. [nondim]
   logical :: read_khth           !< If true, read a file containing the spatially varying horizontal
                                  !! isopycnal height diffusivity
   logical :: use_stanley_gm      !< If true, also use the Stanley parameterization in MOM_thickness_diffuse
@@ -2195,6 +2192,7 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
                                  ! available.
   logical :: use_meke = .false. ! If true, use the MEKE formulation for the thickness diffusivity.
   integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
+  logical :: stoch_eos           ! Can't use Stanley param here unless stoch_eos is true
   integer :: i, j
 
   CS%initialized = .true.
@@ -2322,15 +2320,13 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
                  "streamfunction formulation, expressed as a fraction of planetary "//&
                  "rotation, OMEGA. This should be tiny but non-zero to avoid degeneracy.", &
                  default=1.e-15, units="nondim", scale=US%Z_to_L, do_not_log=.not.CS%use_FGNV_streamfn)
+  call get_param(param_file, mdl, "STOCH_EOS", stoch_eos, &
+                 default=.false., do_not_log=.true.)
   call get_param(param_file, mdl, "USE_STANLEY_GM", CS%use_stanley_gm, &
                  "If true, turn on Stanley SGS T variance parameterization "// &
                  "in GM code.", default=.false.)
-  if (CS%use_stanley_gm) then
-    call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, &
-                 "Coefficient correlating the temperature gradient and SGS T variance.", &
-                 units="nondim", default=-1.0, do_not_log=.true.)
-    if (Stanley_coeff < 0.0) call MOM_error(FATAL, &
-                 "STANLEY_COEFF must be set >= 0 if USE_STANLEY_GM is true.")
+  if (CS%use_Stanley_GM .and. .not.stoch_eos) then
+    call MOM_error(FATAL, "thickness_diffuse_init: USE_STANLEY_GM requires STOCH_EOS")
   endif
   call get_param(param_file, mdl, "OMEGA", omega, &
                  "The rotation rate of the earth.", &