+(*)Fix numerous issues with MOM_stoch_eos

Corrected several bugs that would prevent cases with STOCH_EOS = True or that set a non-negative value of STANLEY_COEF from running at all. There was also extensive refactoring of the MOM_stoch_eos.F90 code. Specifically this commit makes the following changes: - The new routine stoch_eos_register_restarts was added to register the restart field associated with the MOM_stoch_EOS module before the restarts are read and before restart_registry_lock is called, which was causing any test case with STOCH_EOS=True or a positive value of STANLEY_COEF to have a fatal error. - Added a missing dimensional rescaling factor in the get_param call for KD_SMOOTH in MOM_stoch_eos_init, which would have caused cases exercising the MOM_stoch_eos options to fail the dimensional consistency tests. - MOM_stoch_eos_init was changed from a subroutine into a function that returns a logical value indicating whether this routine is used further. The order of the arguments was modified to match the order used in every other init call. - The new routine post_stoch_eos_diags was added to write diagnostics associated with the stoch_eos module. - The MOM_stoch_eos_CS type was made opaque. - The unused diag argument was removed from MOM_stoch_eos_run. - Unit arguments were added to the get_param calls for STANLEY_COEFF and STANLEY_A. - Four arrays in the MOM_stoch_eos_CS type that had been declared to optionally use static memory allocation were modified to be simple allocatables, as they are not used in the vast majority of MOM6 cases, and there is no reason to always assign memory to them. - The register_restart_field call for "stoch_eos_pattern" was revised to use the newer, more direct form rather than working via a vardesc type. - Return statements were added to several of the MOM_stoch_eos routines in the cases where they are not supposed to do anything. - The comments describing several real variables in the MOM_stoch_eos module were added or modified to describe their units using the standard format. - The module use statements at the start of the MOM_stoch_eos module were updated to reflect these changes. There were also parallel changes in MOM.F90: - Added use_stochastic_EOS element to the MOM_control_struct to indicate whether the stoch_eos calls are to be used. - MOM_stoch_eos_init is only called if temperature and salinities are state variables, as it makes no sense to call it otherwise. - Calls to stoch_eos routines were updates to reflect the new interfaces. - The contents of the MOM_stoch_eos_CS type are no longer used in MOM.F90. All answers are bitwise identical in cases that do not use dimensional rescaling, but answers will change (be corrected) in some cases that do use dimensional consistency tests. Several public interfaces to MOM_stoch_eos routines were altered, and there are changes to multiple MOM_parameter_doc files due to the new units, and due to the fact that stoch_EOS parameters are no longer being logged in cases where they are meaningless.
NOAA-GFDL · Dec 8, 2022 · 47cc149 · 47cc149
1 parent d01c42a
commit 47cc149
Show file tree

Hide file tree

Showing 2 changed files with 132 additions and 85 deletions.
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
@@ -114,7 +114,8 @@ 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_stoch_eos,             only : MOM_stoch_eos_init, MOM_stoch_eos_run, MOM_stoch_eos_CS
+use MOM_stoch_eos,             only : stoch_EOS_register_restarts, post_stoch_EOS_diags, 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
@@ -288,6 +289,7 @@ module MOM
   logical :: thickness_diffuse_first !< If true, diffuse thickness before dynamics.
   logical :: mixedlayer_restrat      !< If true, use submesoscale mixed layer restratifying scheme.
   logical :: useMEKE                 !< If true, call the MEKE parameterization.
+  logical :: use_stochastic_EOS      !< If true, use the stochastic EOS parameterizations.
   logical :: useWaves                !< If true, update Stokes drift
   logical :: use_p_surf_in_EOS       !< If true, always include the surface pressure contributions
                                      !! in equation of state calculations.
@@ -1079,12 +1081,12 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
 
   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)
+  if (CS%use_stochastic_EOS) call MOM_stoch_eos_run(G, u, v, dt, Time_local, CS%stoch_eos_CS)
   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)
+  if (CS%use_stochastic_EOS) then
+    call MOM_calc_varT(G, GV, h, CS%tv, CS%stoch_eos_CS, dt)
+    if (associated(CS%tv%varT)) call pass_var(CS%tv%varT, G%Domain, clock=id_clock_pass, halo=1)
   endif
   call cpu_clock_end(id_clock_varT)
 
@@ -1297,9 +1299,7 @@ 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)
+  if (CS%use_stochastic_EOS) call post_stoch_EOS_diags(CS%stoch_eos_CS, CS%tv, CS%diag)
   call disable_averaging(CS%diag)
   call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
 
@@ -2681,6 +2681,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     call waves_register_restarts(waves_CSp, HI, GV, param_file, restart_CSp)
   endif
 
+  if (use_temperature) then
+    call stoch_EOS_register_restarts(HI, param_file, CS%stoch_eos_CS, restart_CSp)
+  endif
+
   call callTree_waypoint("restart registration complete (initialize_MOM)")
   call restart_registry_lock(restart_CSp)
 
@@ -2966,7 +2970,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     call interface_filter_init(Time, G, GV, US, param_file, diag, CS%CDp, CS%interface_filter_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 (use_temperature) then
+    CS%use_stochastic_EOS = MOM_stoch_eos_init(Time, G, US, param_file, diag, CS%stoch_eos_CS, restart_CSp)
+  else
+    CS%use_stochastic_EOS = .false.
+  endif
 
   if (CS%use_porbar) &
     call porous_barriers_init(Time, US, param_file, diag, CS%por_bar_CS)

diff --git a/src/core/MOM_stoch_eos.F90 b/src/core/MOM_stoch_eos.F90
@@ -2,46 +2,44 @@
 module MOM_stoch_eos
 
 ! This file is part of MOM6. See LICENSE.md for the license.
-use MOM_grid,            only : ocean_grid_type
-use MOM_hor_index,       only : hor_index_type
-use MOM_file_parser,     only : get_param, param_file_type
-use MOM_random,          only : PRNG,random_2d_constructor,random_2d_norm
-use MOM_time_manager,    only : time_type
-use MOM_io,              only : vardesc, var_desc
-use MOM_restart,         only : MOM_restart_CS,is_new_run
-use MOM_diag_mediator,   only : register_diag_field,post_data,diag_ctrl,safe_alloc_ptr
-use MOM_variables,       only : thermo_var_ptrs
-use MOM_verticalGrid,    only : verticalGrid_type
-use MOM_restart,         only : register_restart_field
-use MOM_isopycnal_slopes,only : vert_fill_TS
-!use random_numbers_mod, only : getRandomNumbers,initializeRandomNumberStream,randomNumberStream
+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_grid,             only : ocean_grid_type
+use MOM_hor_index,        only : hor_index_type
+use MOM_isopycnal_slopes, only : vert_fill_TS
+use MOM_random,           only : PRNG, random_2d_constructor, random_2d_norm
+use MOM_restart,          only : MOM_restart_CS, register_restart_field, is_new_run, query_initialized
+use MOM_time_manager,     only : time_type
+use MOM_unit_scaling,     only : unit_scale_type
+use MOM_variables,        only : thermo_var_ptrs
+use MOM_verticalGrid,     only : verticalGrid_type
+!use random_numbers_mod,  only : getRandomNumbers, initializeRandomNumberStream, randomNumberStream
 
 implicit none; private
 #include <MOM_memory.h>
 
 public MOM_stoch_eos_init
 public MOM_stoch_eos_run
+public stoch_EOS_register_restarts
+public post_stoch_EOS_diags
 public MOM_calc_varT
 
 !> Describes parameters of the stochastic component of the EOS
 !! correction, described in Stanley et al. JAMES 2020.
-type, public :: MOM_stoch_eos_CS
-  real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: l2_inv
-                                    !< One over sum of the T cell side side lengths squared
-  real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: rgauss
-                                    !< nondimensional random Gaussian
-  real        :: tfac=0.27          !< Nondimensional decorrelation time factor, ~1/3.7
-  real        :: amplitude=0.624499 !< Nondimensional std dev of Gaussian
+type, public :: MOM_stoch_eos_CS ; private
+  real, allocatable :: l2_inv(:,:)  !< One over sum of the T cell side side lengths squared [L-2 ~> m-2]
+  real, allocatable :: rgauss(:,:)  !< nondimensional random Gaussian [nondim]
+  real        :: tfac=0.27          !< Nondimensional decorrelation time factor, ~1/3.7 [nondim]
+  real        :: amplitude=0.624499 !< Nondimensional standard deviation of Gaussian [nondim]
   integer     :: seed               !< PRNG seed
   type(PRNG)  ::  rn_CS             !< PRNG control structure
-  real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: pattern
-                          !< Random pattern for stochastic EOS [nondim]
-  real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: phi
-                          !< temporal correlation stochastic EOS [nondim]
+  real, allocatable :: pattern(:,:) !< Random pattern for stochastic EOS [nondim]
+  real, allocatable :: phi(:,:)     !< temporal correlation stochastic EOS [nondim]
   logical :: use_stoch_eos!< If true, use the stochastic equation of state (Stanley et al. 2020)
   real :: stanley_coeff   !< Coefficient correlating the temperature gradient
-                          !! and SGS T variance; if <0, turn off scheme in all codes
-  real :: stanley_a       !< a in exp(aX) in stochastic coefficient
+                          !! and SGS T variance [nondim]; if <0, turn off scheme in all codes
+  real :: stanley_a       !< a in exp(aX) in stochastic coefficient [nondim]
   real :: kappa_smooth    !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1]
 
   !>@{ Diagnostic IDs
@@ -52,61 +50,64 @@ module MOM_stoch_eos
 
 contains
 
-!> Initializes MOM_stoch_eos module.
-subroutine MOM_stoch_eos_init(G, Time, param_file, CS, restart_CS, diag)
-  type(param_file_type),  intent(in)    :: param_file   !< structure indicating parameter file to parse
-  type(ocean_grid_type),  intent(in)    :: G            !< The ocean's grid structure.
-  type(time_type),        intent(in)    :: Time         !< Time for stochastic process
-  type(MOM_stoch_eos_CS), intent(inout) :: CS           !< Stochastic control structure
-  type(MOM_restart_CS),   pointer       :: restart_CS   !< A pointer to the restart control structure.
-  type(diag_ctrl),        target, intent(inout) :: diag !< to control diagnostics
+!> Initializes MOM_stoch_eos module, returning a logical indicating whether this module will be used.
+logical function MOM_stoch_eos_init(Time, G, US, param_file, diag, CS, restart_CS)
+  type(time_type),         intent(in)    :: Time       !< Time for stochastic process
+  type(ocean_grid_type),   intent(in)    :: G          !< The ocean's grid structure.
+  type(unit_scale_type),   intent(in)    :: US         !< A dimensional unit scaling type
+  type(param_file_type),   intent(in)    :: param_file !< structure indicating parameter file to parse
+  type(diag_ctrl), target, intent(inout) :: diag       !< Structure used to control diagnostics
+  type(MOM_stoch_eos_CS),  intent(inout) :: CS         !< Stochastic control structure
+  type(MOM_restart_CS),    pointer       :: restart_CS !< A pointer to the restart control structure.
 
   ! local variables
   integer :: i,j
-  type(vardesc) :: vd
-  CS%seed=0
-  ! contants
-  !pi=2*acos(0.0)
+
+  MOM_stoch_eos_init = .false.
+
+  CS%seed = 0
+
   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.)
   call get_param(param_file, "MOM_stoch_eos", "STANLEY_COEFF", CS%stanley_coeff, &
                  "Coefficient correlating the temperature gradient "//&
-                 "and SGS T variance.", default=-1.0)
+                 "and SGS T variance.", units="nondim", default=-1.0)
   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.", default=1.0)
+                 "SGS T variance.", units="nondim", default=1.0, &
+                 do_not_log=((CS%stanley_coeff<0.0) .or. .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)
+                 units="m2 s-1", default=1.0e-6, scale=US%m_to_Z**2*US%T_to_s, &
+                 do_not_log=(CS%stanley_coeff<0.0))
 
-  !don't run anything if STANLEY_COEFF < 0
+  ! Don't run anything if STANLEY_COEFF < 0
   if (CS%stanley_coeff >= 0.0) then
+    if (.not.allocated(CS%pattern)) call MOM_error(FATAL, &
+        "MOM_stoch_eos_CS%pattern is not allocated when it should be, suggesting that "//&
+        "stoch_EOS_register_restarts() has not been called before MOM_stoch_eos_init().")
 
-    ALLOC_(CS%pattern(G%isd:G%ied,G%jsd:G%jed)) ; CS%pattern(:,:) = 0.0
-    vd = var_desc("stoch_eos_pattern","nondim","Random pattern for stoch EOS",'h','1')
-    call register_restart_field(CS%pattern, vd, .false., restart_CS)
-    ALLOC_(CS%phi(G%isd:G%ied,G%jsd:G%jed)) ; CS%phi(:,:) = 0.0
-    ALLOC_(CS%l2_inv(G%isd:G%ied,G%jsd:G%jed))
-    ALLOC_(CS%rgauss(G%isd:G%ied,G%jsd:G%jed))
+    allocate(CS%phi(G%isd:G%ied,G%jsd:G%jed), source=0.0)
+    allocate(CS%l2_inv(G%isd:G%ied,G%jsd:G%jed), source=0.0)
+    allocate(CS%rgauss(G%isd:G%ied,G%jsd:G%jed), source=0.0)
     call get_param(param_file, "MOM_stoch_eos", "SEED_STOCH_EOS", CS%seed, &
                  "Specfied seed for random number sequence ", default=0)
     call random_2d_constructor(CS%rn_CS, G%HI, Time, CS%seed)
     call random_2d_norm(CS%rn_CS, G%HI, CS%rgauss)
-    ! fill array with approximation of grid area needed for decorrelation
-    ! time-scale calculation
+    ! fill array with approximation of grid area needed for decorrelation time-scale calculation
     do j=G%jsc,G%jec
       do i=G%isc,G%iec
-        CS%l2_inv(i,j)=1.0/(G%dxT(i,j)**2+G%dyT(i,j)**2)
+        CS%l2_inv(i,j) = 1.0/(G%dxT(i,j)**2+G%dyT(i,j)**2)
       enddo
     enddo
-    if (is_new_run(restart_CS)) then
-      do j=G%jsc,G%jec
-        do i=G%isc,G%iec
-          CS%pattern(i,j)=CS%amplitude*CS%rgauss(i,j)
-        enddo
-      enddo
+
+    if (.not.query_initialized(CS%pattern, "stoch_eos_pattern", restart_CS) .or. &
+        is_new_run(restart_CS)) then
+      do j=G%jsc,G%jec ; do i=G%isc,G%iec
+        CS%pattern(i,j) = CS%amplitude*CS%rgauss(i,j)
+      enddo ; enddo
     endif
 
     !register diagnostics
@@ -120,10 +121,32 @@ subroutine MOM_stoch_eos_init(G, Time, param_file, CS, restart_CS, diag)
     endif
   endif
 
-end subroutine MOM_stoch_eos_init
+  ! This module is only used if explicitly enabled or a positive correlation coefficient is set.
+  MOM_stoch_eos_init = CS%use_stoch_eos .or. (CS%stanley_coeff >= 0.0)
+
+end function MOM_stoch_eos_init
+
+!> Register fields related to the stoch_EOS module for resarts
+subroutine stoch_EOS_register_restarts(HI, param_file, CS, restart_CS)
+  type(hor_index_type),    intent(in)    :: HI         !< Horizontal index structure
+  type(param_file_type),   intent(in)    :: param_file !< structure indicating parameter file to parse
+  type(MOM_stoch_eos_CS),  intent(inout) :: CS         !< Stochastic control structure
+  type(MOM_restart_CS),    pointer       :: restart_CS !< A pointer to the restart control structure.
+
+  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, do_not_log=.true.)
+
+  if (CS%stanley_coeff >= 0.0) then
+    allocate(CS%pattern(HI%isd:HI%ied,HI%jsd:HI%jed), source=0.0)
+    call register_restart_field(CS%pattern, "stoch_eos_pattern", .false., restart_CS, &
+                                "Random pattern for stoch EOS", "nondim")
+  endif
+
+end subroutine stoch_EOS_register_restarts
 
 !> Generates a pattern in space and time for the ocean stochastic equation of state
-subroutine MOM_stoch_eos_run(G, u, v, delt, Time, CS, diag)
+subroutine MOM_stoch_eos_run(G, u, v, delt, Time, CS)
   type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                            intent(in)    :: u    !< The zonal velocity [L T-1 ~> m s-1].
@@ -132,29 +155,43 @@ subroutine MOM_stoch_eos_run(G, u, v, delt, Time, CS, diag)
   real,                    intent(in)    :: delt !< Time step size for AR1 process [T ~> s].
   type(time_type),         intent(in)    :: Time !< Time for stochastic process
   type(MOM_stoch_eos_CS),  intent(inout) :: CS   !< Stochastic control structure
-  type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics
 
   ! local variables
-  integer ::  i,j
-  integer :: yr,mo,dy,hr,mn,sc
-  real    :: phi,ubar,vbar
+  real    :: ubar, vbar ! Averaged velocities [L T-1 ~> m s-1]
+  real    :: phi        ! A temporal correlation factor [nondim]
+  integer :: i, j
+
+  ! Return without doing anything if this capability is not enabled.
+  if (.not.CS%use_stoch_eos) return
 
   call random_2d_constructor(CS%rn_CS, G%HI, Time, CS%seed)
   call random_2d_norm(CS%rn_CS, G%HI, CS%rgauss)
 
   ! advance AR(1)
   do j=G%jsc,G%jec
     do i=G%isc,G%iec
-      ubar=0.5*(u(I,j,1)*G%mask2dCu(I,j)+u(I-1,j,1)*G%mask2dCu(I-1,j))
-      vbar=0.5*(v(i,J,1)*G%mask2dCv(i,J)+v(i,J-1,1)*G%mask2dCv(i,J-1))
-      phi=exp(-delt*CS%tfac*sqrt((ubar**2+vbar**2)*CS%l2_inv(i,j)))
-      CS%pattern(i,j)=phi*CS%pattern(i,j) + CS%amplitude*sqrt(1-phi**2)*CS%rgauss(i,j)
-      CS%phi(i,j)=phi
+      ubar = 0.5*(u(I,j,1)*G%mask2dCu(I,j)+u(I-1,j,1)*G%mask2dCu(I-1,j))
+      vbar = 0.5*(v(i,J,1)*G%mask2dCv(i,J)+v(i,J-1,1)*G%mask2dCv(i,J-1))
+      phi = exp(-delt*CS%tfac*sqrt((ubar**2+vbar**2)*CS%l2_inv(i,j)))
+      CS%pattern(i,j) = phi*CS%pattern(i,j) + CS%amplitude*sqrt(1-phi**2)*CS%rgauss(i,j)
+      CS%phi(i,j) = phi
     enddo
   enddo
 
 end subroutine MOM_stoch_eos_run
 
+!> Write out any diagnostics related to this module.
+subroutine post_stoch_EOS_diags(CS, tv, diag)
+  type(MOM_stoch_eos_CS), intent(in) :: CS  !< Stochastic control structure
+  type(thermo_var_ptrs),  intent(in) :: tv  !< Thermodynamics structure
+  type(diag_ctrl),        intent(inout) :: diag !< Structure to control diagnostics
+
+  if (CS%id_stoch_eos > 0) call post_data(CS%id_stoch_eos, CS%pattern, diag)
+  if (CS%id_stoch_phi > 0) call post_data(CS%id_stoch_phi, CS%phi, diag)
+  if (CS%id_tvar_sgs > 0) call post_data(CS%id_tvar_sgs, tv%varT, diag)
+
+end subroutine post_stoch_EOS_diags
+
 !> Computes a parameterization of the SGS temperature variance
 subroutine MOM_calc_varT(G, GV, h, tv, CS, dt)
   type(ocean_grid_type),   intent(in)   :: G   !< The ocean's grid structure.
@@ -171,15 +208,17 @@ subroutine MOM_calc_varT(G, GV, h, tv, CS, dt)
                   !! in massless layers filled vertically by diffusion.
     S             !> The filled salinity [S ~> ppt], with the values in
                   !! in massless layers filled vertically by diffusion.
-  integer :: i, j, k
   real :: hl(5)              !> Copy of local stencil of H [H ~> m]
   real :: dTdi2, dTdj2       !> Differences in T variance [C2 ~> degC2]
+  integer :: i, j, k
+
+  ! Nothing happens if a negative correlation coefficient is set.
+  if (CS%stanley_coeff < 0.0) return
 
   ! 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.
-  if (.not. associated(tv%varT)) call safe_alloc_ptr(tv%varT, G%isd, G%ied, G%jsd, G%jed, GV%ke)
-
+  if (.not. associated(tv%varT)) allocate(tv%varT(G%isd:G%ied, G%jsd:G%jed, GV%ke), source=0.0)
   call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, halo_here=1, larger_h_denom=.true.)
 
   do k=1,G%ke
@@ -193,12 +232,12 @@ subroutine MOM_calc_varT(G, GV, h, tv, CS, dt)
 
         ! SGS variance in i-direction [C2 ~> degC2]
         dTdi2 = ( ( G%mask2dCu(I  ,j) * G%IdxCu(I  ,j) * ( T(i+1,j,k) - T(i,j,k) ) &
-              + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( T(i,j,k) - T(i-1,j,k) ) &
-              ) * G%dxT(i,j) * 0.5 )**2
+                  + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( T(i,j,k) - T(i-1,j,k) ) &
+                ) * G%dxT(i,j) * 0.5 )**2
         ! SGS variance in j-direction [C2 ~> degC2]
         dTdj2 = ( ( G%mask2dCv(i,J  ) * G%IdyCv(i,J  ) * ( T(i,j+1,k) - T(i,j,k) ) &
-              + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( T(i,j,k) - T(i,j-1,k) ) &
-              ) * G%dyT(i,j) * 0.5 )**2
+                  + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( T(i,j,k) - T(i,j-1,k) ) &
+                ) * G%dyT(i,j) * 0.5 )**2
         tv%varT(i,j,k) = CS%stanley_coeff * ( dTdi2 + dTdj2 )
         ! Turn off scheme near land
         tv%varT(i,j,k) = tv%varT(i,j,k) * (minval(hl) / (maxval(hl) + GV%H_subroundoff))
@@ -210,7 +249,7 @@ subroutine MOM_calc_varT(G, GV, h, tv, CS, dt)
     do k=1,G%ke
       do j=G%jsc,G%jec
         do i=G%isc,G%iec
-          tv%varT(i,j,k) = exp (CS%stanley_a * CS%pattern(i,j)) * tv%varT(i,j,k)
+          tv%varT(i,j,k) = exp(CS%stanley_a * CS%pattern(i,j)) * tv%varT(i,j,k)
         enddo
       enddo
     enddo