Corrections to cell-averaged density computation (mom-ocean#213)

* initial hooks for stochastic EOS modifications

* remove debug statements

* add documentation

* Change ampltiude from 0.39 to sqrt(.39)

* remove global_indexing logic from stoch_eos_init

* switch to using MOM_random and add restart capability

* update random sequence to update each each time-step

* remove tseed0 from MOM_random (leftover from debugging)

* Added necessary submodules and S^2, T^2 diagnostics to MOM_diagnostics

* Added diagnostics for outputting variables related to the stochastic parameterization.

* Diagnostics in MOM_PressureForce_FV updated for stochastic (rather than deterministic) Stanley SGS T variance parameterization.

* Added parentheses for reproducibility.

* Changed diagnostics to account for possible absence of stoch_eos_pattern in MOM_PressureForce_FV,
when deterministic parameterization is on.

* remove mom6_da_hooks and geoKdTree from pkg

* add stochastic compoment to MOM_thickness_diffuse

* fix array size declaration and post_data

* Corrected indexing of loops in MOM_calc_varT

* Changed how parameterization of SGS T variance (deterministic and stochastic) is switched on in PGF and thickness diffusion codes

* Corrected a few typos

* Cleaned up indices, redundant diagnostic, printing

* Fixed diagnostic IDs

* Fixed diagnostics typo

* Corrected indices in calculation of tv%varT

* Minor index fix

* Corrected bug in pressure in Stanley diagnostics

* Fixed whitespace error

* Stoch eos clock (#5)

*Added a clock for the Stanley parameterization

Co-authored-by: jkenigson <jkenigso@gmail.com>

* add halo update to random pattern

* Update MOM_stoch_eos.F90

Fix bug for looping over compute domain (is -> isc etc.)

* Avoid unnessary computations on halo (MOM_stoch_eos) and code clean-up (MOM_thickness_diffuse)

* Removed halo updates before determ param calc

* Update MOM_stoch_eos.F90

Removed unnecessary code

* Bug - indices are transposed

* Changed Stanley stochastic coefficient from exp(X) to exp(aX) (mom-ocean#9)

* Changed Stanley stochastic coefficient from exp(X) to exp(aX)

* Extra spaces removed

* Stoch eos init fix (mom-ocean#10)

* Don't bother calculating tv%varT if stanley_coeff<0

* Missing then added

* Merge Ian Grooms Tvar Discretization (mom-ocean#11)

* Update MOM_stoch_eos.F90

In progress updating stencil for$ | dx \times \nabla T|^2$ calculation

* New discretization of |dx\circ\nablaT|^2

Co-authored-by: Ian Grooms <ian.grooms@colorado.edu>

* Multiplied tvar%SGS by grid cell thickness ratio

* Added limiter for tv%varT

* Stoch eos ncar linear disc (mom-ocean#12)

* Update MOM_stoch_eos.F90

In progress updating stencil for$ | dx \times \nabla T|^2$ calculation

* New discretization of |dx\circ\nablaT|^2

* AR1 timescale land mask

Adds land mask to the computation of the AR1 decorrelation time

* Update dt in call to MOM_stoch_eos_run

The call to `MOM_stoch_eos_run` (which time steps the noise) is from within `step_MOM_dynamics`. `step_MOM_dynamics` advances on time step `dt` (per line 957), but the noise is updated using `dt_thermo`. It seems more appropriate to update the noise using `dt`, since it gets called from within `step_MOM_dynamics`.

* Fixed the units for r_sm_H

* Remove vestigial declarations

The variables `hl`, `Tl`, `mn_T`, `mn_T2`, and `r_sm_H` are no longer used, so I removed their declarations and an OMP private clause

Co-authored-by: Ian Grooms <ian.grooms@colorado.edu>

* Update MOM_thickness_diffuse.F90

Changed index for soft convention

* Update CVMix-src

* Ensure use_varT, etc., initialized

* Don't register stanley diagnostics if scheme is off

* Stanley density second derivs at h pts (mom-ocean#15)

* Change discretization of Stanley correction (drho_dT_dT at h points)

* Limit Stanley noise, shrink limiting value

* Revert t variance discretization

* Reverted variable declarations

* Stanley scheme in mixed_layer_restrat, vert_fill in stoch_eos, code cleanup (mom-ocean#19)

* Test Stanley EOS param in mixed_layer_restrat

* Fix size of TS cov, S var in Stanley calculate_density calls

* Test move stanley scheme initialization

* Added missing openMP directives

* Revert Stanley tvar discretization (mom-ocean#18)

* Perform vertical filling in calculation of T variance

* Variable declaration syntax error, remove scaling from get_param

* Fix call to vert_fill_TS

* Code cleanup, whitespace cleanup

Co-authored-by: Jessica Kenigson <jessicak@cheyenne1.cheyenne.ucar.edu>

* Use Stanley (2020) variance; scheme off at coast

* Comment clean-up

* Remove factor of 0.5 in Tvar

* Don't calculate Stanley diagnostics on halo

* Change start indices in stanley_density_1d

* Stanley param in MOM_isopycnal_slopes (mom-ocean#22)

Stanley param in MOM_isopycnal_slopes and thickness diffuse index fix

* Set eady flag to true if use_stored_slopes is true

* Cleanup, docs, whitespace

* Docs and whitespace

* Docs and whitespace

* Docs and whitespace

* Whitespace cleanup

* Whitespace cleanup

* Clean up whitespace

* Docs cleanup

* use_stanley

* Update MOM_lateral_mixing_coeffs.F90

* Adds link to another TEOS10 module

* Set Stanley off for testing

* Line continuation

Co-authored-by: Phil Pegion <38869668+pjpegion@users.noreply.github.com>
Co-authored-by: Philip Pegion <Philip.Pegion@noaa.gov>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne6.cheyenne.ucar.edu>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne3.cheyenne.ucar.edu>
Co-authored-by: jkenigson <jkenigso@gmail.com>
Co-authored-by: jskenigson <jessica.kenigson@colorado.edu>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne1.cheyenne.ucar.edu>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne5.cheyenne.ucar.edu>
Co-authored-by: Philip Pegion <ppegion@Philips-MacBook-Pro.local>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne4.cheyenne.ucar.edu>

Author: 10 people

.testing/tc2/MOM_input

@@ -297,7 +297,7 @@ BOUND_CORIOLIS = True           !   [Boolean] default = False
                                 ! v-points, and similarly at v-points.  This option would
                                 ! have no effect on the SADOURNY Coriolis scheme if it
                                 ! were possible to use centered difference thickness fluxes.
-PGF_STANLEY_T2_DET_COEFF = 0.5  !   [nondim] default = -1.0
+PGF_STANLEY_T2_DET_COEFF = -1.0 !   [nondim] default = -1.0
                                 ! The coefficient correlating SGS temperature variance with the mean temperature
                                 ! gradient in the deterministic part of the Stanley form of the Brankart
                                 ! correction. Negative values disable the scheme.
@@ -430,7 +430,7 @@ KHTH = 1.0                      !   [m2 s-1] default = 0.0
                                 ! The background horizontal thickness diffusivity.
 KHTH_MAX = 900.0                !   [m2 s-1] default = 0.0
                                 ! The maximum horizontal thickness diffusivity.
-STANLEY_PRM_DET_COEFF = 0.5     !   [nondim] default = -1.0
+STANLEY_PRM_DET_COEFF = -1.0    !   [nondim] default = -1.0
                                 ! The coefficient correlating SGS temperature variance with the mean temperature
                                 ! gradient in the deterministic part of the Stanley parameterization. Negative
                                 ! values disable the scheme.

src/core/MOM.F90

@@ -108,6 +108,7 @@ module MOM
 use MOM_shared_initialization, only : write_ocean_geometry_file
 use MOM_sponge,                only : init_sponge_diags, sponge_CS
 use MOM_state_initialization,  only : MOM_initialize_state
+use MOM_stoch_eos,             only : MOM_stoch_eos_init,MOM_stoch_eos_run,MOM_stoch_eos_CS,mom_calc_varT
 use MOM_sum_output,            only : write_energy, accumulate_net_input
 use MOM_sum_output,            only : MOM_sum_output_init, MOM_sum_output_end
 use MOM_sum_output,            only : sum_output_CS
@@ -242,6 +243,7 @@ module MOM
   logical :: use_ALE_algorithm  !< If true, use the ALE algorithm rather than layered
                     !! isopycnal/stacked shallow water mode. This logical is set by calling the
                     !! function useRegridding() from the MOM_regridding module.
+  type(MOM_stoch_eos_CS) :: stoch_eos_CS !< structure containing random pattern for stoch EOS
   logical :: alternate_first_direction !< If true, alternate whether the x- or y-direction
                     !! updates occur first in directionally split parts of the calculation.
   real    :: first_dir_restart = -1.0 !< A real copy of G%first_direction for use in restart files
@@ -444,6 +446,8 @@ module MOM
 integer :: id_clock_other
 integer :: id_clock_offline_tracer
 integer :: id_clock_unit_tests
+integer :: id_clock_stoch
+integer :: id_clock_varT
 !>@}
 
 contains
@@ -1063,6 +1067,15 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
   showCallTree = callTree_showQuery()
 
   call cpu_clock_begin(id_clock_dynamics)
+  call cpu_clock_begin(id_clock_stoch)
+  if (CS%stoch_eos_CS%use_stoch_eos) call MOM_stoch_eos_run(G,u,v,dt,Time_local,CS%stoch_eos_CS,CS%diag)
+  call cpu_clock_end(id_clock_stoch)
+  call cpu_clock_begin(id_clock_varT)
+  if (CS%stoch_eos_CS%stanley_coeff >= 0.0) then
+    call MOM_calc_varT(G,GV,h,CS%tv,CS%stoch_eos_CS,dt)
+    call pass_var(CS%tv%varT, G%Domain,clock=id_clock_pass,halo=1)
+  endif
+  call cpu_clock_end(id_clock_varT)
 
   if ((CS%t_dyn_rel_adv == 0.0) .and. CS%thickness_diffuse .and. CS%thickness_diffuse_first) then
 
@@ -1229,6 +1242,9 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
   if (IDs%id_u > 0) call post_data(IDs%id_u, u, CS%diag)
   if (IDs%id_v > 0) call post_data(IDs%id_v, v, CS%diag)
   if (IDs%id_h > 0) call post_data(IDs%id_h, h, CS%diag)
+  if (CS%stoch_eos_CS%id_stoch_eos > 0) call post_data(CS%stoch_eos_CS%id_stoch_eos, CS%stoch_eos_CS%pattern, CS%diag)
+  if (CS%stoch_eos_CS%id_stoch_phi > 0) call post_data(CS%stoch_eos_CS%id_stoch_phi, CS%stoch_eos_CS%phi, CS%diag)
+  if (CS%stoch_eos_CS%id_tvar_sgs > 0) call post_data(CS%stoch_eos_CS%id_tvar_sgs, CS%tv%varT, CS%diag)
   call disable_averaging(CS%diag)
   call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
 
@@ -2765,6 +2781,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   call set_visc_init(Time, G, GV, US, param_file, diag, CS%visc, CS%set_visc_CSp, restart_CSp, CS%OBC)
   call thickness_diffuse_init(Time, G, GV, US, param_file, diag, CS%CDp, CS%thickness_diffuse_CSp)
 
+  new_sim = is_new_run(restart_CSp)
+  call MOM_stoch_eos_init(G,Time,param_file,CS%stoch_eos_CS,restart_CSp,diag)
   if (CS%split) then
     allocate(eta(SZI_(G),SZJ_(G)), source=0.0)
     call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, &
@@ -2854,7 +2872,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   endif
 
   ! This subroutine initializes any tracer packages.
-  new_sim = is_new_run(restart_CSp)
   call tracer_flow_control_init(.not.new_sim, Time, G, GV, US, CS%h, param_file, &
              CS%diag, CS%OBC, CS%tracer_flow_CSp, CS%sponge_CSp, &
              CS%ALE_sponge_CSp, CS%tv)
@@ -3065,6 +3082,7 @@ subroutine register_diags(Time, G, GV, US, IDs, diag)
       v_extensive=.true.)
   IDs%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, &
       Time, 'Instantaneous Sea Surface Height', 'm', conversion=US%Z_to_m)
+
 end subroutine register_diags
 
 !> Set up CPU clock IDs for timing various subroutines.
@@ -3097,6 +3115,8 @@ subroutine MOM_timing_init(CS)
   if (CS%offline_tracer_mode) then
     id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=CLOCK_SUBCOMPONENT)
   endif
+  id_clock_stoch = cpu_clock_id('(Stochastic EOS)', grain=CLOCK_MODULE)
+  id_clock_varT = cpu_clock_id('(SGS Temperature Variance)', grain=CLOCK_MODULE)
 
 end subroutine MOM_timing_init
 

src/core/MOM_PressureForce_FV.F90

@@ -58,11 +58,12 @@ module MOM_PressureForce_FV
   integer :: Recon_Scheme   !< Order of the polynomial of the reconstruction of T & S
                             !! for the finite volume pressure gradient calculation.
                             !! By the default (1) is for a piecewise linear method
-  real :: Stanley_T2_det_coeff !< The coefficient correlating SGS temperature variance with
-                            !! the mean temperature gradient in the deterministic part of
-                            !! the Stanley form of the Brankart correction.
+
+  logical :: use_stanley_pgf  !< If true, turn on Stanley parameterization in the PGF
   integer :: id_e_tidal = -1 !< Diagnostic identifier
-  integer :: id_tvar_sgs = -1 !< Diagnostic identifier
+  integer :: id_rho_pgf = -1 !< Diagnostic identifier
+  integer :: id_rho_stanley_pgf = -1 !< Diagnostic identifier
+  integer :: id_p_stanley = -1 !< Diagnostic identifier
   type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< Tides control structure
 end type PressureForce_FV_CS
 
@@ -167,7 +168,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
   if (.not.CS%initialized) call MOM_error(FATAL, &
        "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.")
 
-  if (CS%Stanley_T2_det_coeff>=0.) call MOM_error(FATAL, &
+  if (CS%use_stanley_pgf) call MOM_error(FATAL, &
        "MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet"//&
        "implemented in non-Boussinesq mode.")
 
@@ -473,6 +474,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
     S_t, S_b, T_t, T_b ! Top and bottom edge values for linear reconstructions
                        ! of salinity and temperature within each layer.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
+    rho_pgf, rho_stanley_pgf ! Density [kg m-3] from EOS with and without SGS T variance
+                                        ! in Stanley parameterization.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
+    p_stanley ! Pressure [Pa] estimated with Rho_0
+  real :: rho_stanley_scalar ! Scalar quantity to hold density [kg m-3] in Stanley diagnostics.
+  real :: p_stanley_scalar ! Scalar quantity to hold pressure [Pa] in Stanley diagnostics.
   real :: rho_in_situ(SZI_(G)) ! The in situ density [R ~> kg m-3].
   real :: p_ref(SZI_(G))     !   The pressure used to calculate the coordinate
                              ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).
@@ -487,11 +495,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
   logical :: use_ALE         ! If true, use an ALE pressure reconstruction.
   logical :: use_EOS         ! If true, density is calculated from T & S using an equation of state.
   type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
-  real :: Tl(5)              ! copy and T in local stencil [degC]
-  real :: mn_T               ! mean of T in local stencil [degC]
-  real :: mn_T2              ! mean of T**2 in local stencil [degC2]
-  real :: hl(5)              ! Copy of local stencil of H [H ~> m]
-  real :: r_sm_H             ! Reciprocal of sum of H in local stencil [H-1 ~> m-1]
   real, parameter :: C1_6 = 1.0/6.0
   integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
@@ -511,49 +514,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
   use_ALE = .false.
   if (associated(ALE_CSp)) use_ALE = CS%reconstruct .and. use_EOS
 
-  if (CS%Stanley_T2_det_coeff>=0.) then
-    if (.not. associated(tv%varT)) call safe_alloc_ptr(tv%varT, G%isd, G%ied, G%jsd, G%jed, GV%ke)
-    do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1
-      ! Strictly speaking we should estimate the *horizontal* grid-scale variance
-      ! but neither of the following blocks make a rotation to the horizontal
-      ! and instead work along coordinate.
-
-      ! This block calculates a simple |delta T| along coordinates and does
-      ! not allow vanishing layer thicknesses or layers tracking topography
-      !! SGS variance in i-direction [degC2]
-      !dTdi2 = ( ( G%mask2dCu(I  ,j) * G%IdxCu(I  ,j) * ( tv%T(i+1,j,k) - tv%T(i,j,k) ) &
-      !          + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( tv%T(i,j,k) - tv%T(i-1,j,k) ) &
-      !          ) * G%dxT(i,j) * 0.5 )**2
-      !! SGS variance in j-direction [degC2]
-      !dTdj2 = ( ( G%mask2dCv(i,J  ) * G%IdyCv(i,J  ) * ( tv%T(i,j+1,k) - tv%T(i,j,k) ) &
-      !          + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( tv%T(i,j,k) - tv%T(i,j-1,k) ) &
-      !          ) * G%dyT(i,j) * 0.5 )**2
-      !tv%varT(i,j,k) = CS%Stanley_T2_det_coeff * 0.5 * ( dTdi2 + dTdj2 )
-
-      ! This block does a thickness weighted variance calculation and helps control for
-      ! extreme gradients along layers which are vanished against topography. It is
-      ! still a poor approximation in the interior when coordinates are strongly tilted.
-      hl(1) = h(i,j,k) * G%mask2dT(i,j)
-      hl(2) = h(i-1,j,k) * G%mask2dCu(I-1,j)
-      hl(3) = h(i+1,j,k) * G%mask2dCu(I,j)
-      hl(4) = h(i,j-1,k) * G%mask2dCv(i,J-1)
-      hl(5) = h(i,j+1,k) * G%mask2dCv(i,J)
-      r_sm_H = 1. / ( ( hl(1) + ( ( hl(2) + hl(3) ) + ( hl(4) + hl(5) ) ) ) + GV%H_subroundoff )
-      ! Mean of T
-      Tl(1) = tv%T(i,j,k) ; Tl(2) = tv%T(i-1,j,k) ; Tl(3) = tv%T(i+1,j,k)
-      Tl(4) = tv%T(i,j-1,k) ; Tl(5) = tv%T(i,j+1,k)
-      mn_T = ( hl(1)*Tl(1) + ( ( hl(2)*Tl(2) + hl(3)*Tl(3) ) + ( hl(4)*Tl(4) + hl(5)*Tl(5) ) ) ) * r_sm_H
-      ! Adjust T vectors to have zero mean
-      Tl(:) = Tl(:) - mn_T ; mn_T = 0.
-      ! Variance of T
-      mn_T2 = ( hl(1)*Tl(1)*Tl(1) + ( ( hl(2)*Tl(2)*Tl(2) + hl(3)*Tl(3)*Tl(3) ) &
-                                    + ( hl(4)*Tl(4)*Tl(4) + hl(5)*Tl(5)*Tl(5) ) ) ) * r_sm_H
-      ! Variance should be positive but round-off can violate this. Calculating
-      ! variance directly would fix this but requires more operations.
-      tv%varT(i,j,k) = CS%Stanley_T2_det_coeff * max(0., mn_T2)
-    enddo ; enddo ; enddo
-  endif
-
   h_neglect = GV%H_subroundoff
   dz_neglect = GV%H_subroundoff * GV%H_to_Z
   I_Rho0 = 1.0 / GV%Rho0
@@ -690,7 +650,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
   do k=1,nz
     ! Calculate 4 integrals through the layer that are required in the
     ! subsequent calculation.
-
     if (use_EOS) then
       ! The following routine computes the integrals that are needed to
       ! calculate the pressure gradient force. Linear profiles for T and S are
@@ -701,13 +660,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
         if ( CS%Recon_Scheme == 1 ) then
           call int_density_dz_generic_plm(k, tv,  T_t, T_b, S_t, S_b, e, &
                     rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
-                    G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, &
+                    G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa, intz_dpa, intx_dpa, inty_dpa, &
                     useMassWghtInterp=CS%useMassWghtInterp, &
                     use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=G%Z_ref)
         elseif ( CS%Recon_Scheme == 2 ) then
           call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
                     rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
-                    G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, &
+                    G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa, intz_dpa, intx_dpa, inty_dpa, &
                     useMassWghtInterp=CS%useMassWghtInterp, Z_0p=G%Z_ref)
         endif
       else
@@ -797,8 +756,26 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
     endif
   endif
 
+  if (CS%use_stanley_pgf) then
+    do j=js,je ; do i=is,ie ;
+      p_stanley_scalar=0.0
+      do k=1, nz
+        p_stanley_scalar = p_stanley_scalar + 0.5 * h(i,j,k) * GV%H_to_Pa !Pressure at mid-point of layer
+        call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley_scalar, 0.0, 0.0, 0.0, &
+          rho_stanley_scalar, tv%eqn_of_state)
+        rho_pgf(i,j,k) = rho_stanley_scalar
+        call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley_scalar, tv%varT(i,j,k), 0.0, 0.0, &
+          rho_stanley_scalar, tv%eqn_of_state)
+        rho_stanley_pgf(i,j,k) = rho_stanley_scalar
+        p_stanley(i,j,k) = p_stanley_scalar
+        p_stanley_scalar = p_stanley_scalar + 0.5 * h(i,j,k) * GV%H_to_Pa !Pressure at bottom of layer
+     enddo; enddo; enddo
+   endif
+
   if (CS%id_e_tidal>0) call post_data(CS%id_e_tidal, e_tidal, CS%diag)
-  if (CS%id_tvar_sgs>0) call post_data(CS%id_tvar_sgs, tv%varT, CS%diag)
+  if (CS%id_rho_pgf>0) call post_data(CS%id_rho_pgf, rho_pgf, CS%diag)
+  if (CS%id_rho_stanley_pgf>0) call post_data(CS%id_rho_stanley_pgf, rho_stanley_pgf, CS%diag)
+  if (CS%id_p_stanley>0) call post_data(CS%id_p_stanley, p_stanley, CS%diag)
 
 end subroutine PressureForce_FV_Bouss
 
@@ -859,14 +836,16 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS
                  "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, "PGF_STANLEY_T2_DET_COEFF", CS%Stanley_T2_det_coeff, &
-                 "The coefficient correlating SGS temperature variance with "// &
-                 "the mean temperature gradient in the deterministic part of "// &
-                 "the Stanley form of the Brankart correction. "// &
-                 "Negative values disable the scheme.", units="nondim", default=-1.0)
-  if (CS%Stanley_T2_det_coeff>=0.) then
-    CS%id_tvar_sgs = register_diag_field('ocean_model', 'tvar_sgs_pgf', diag%axesTL, &
-        Time, 'SGS temperature variance used in PGF', 'degC2')
+  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
+    CS%id_rho_pgf = register_diag_field('ocean_model', 'rho_pgf', diag%axesTL, &
+        Time, 'rho in PGF', 'kg m3')
+    CS%id_rho_stanley_pgf = register_diag_field('ocean_model', 'rho_stanley_pgf', diag%axesTL, &
+        Time, 'rho in PGF with Stanley correction', 'kg m3')
+    CS%id_p_stanley = register_diag_field('ocean_model', 'p_stanley', diag%axesTL, &
+        Time, 'p in PGF with Stanley correction', 'Pa')
   endif
   if (CS%tides) then
     CS%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &