Stanley scheme in mixed_layer_restrat, vert_fill in stoch_eos, code c…

…leanup (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>

src/core/MOM.F90

@@ -1017,9 +1017,9 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
   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)
+    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	
+  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

src/core/MOM_PressureForce_FV.F90

@@ -490,11 +490,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 [degC]
-  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

src/core/MOM_stoch_eos.F90

@@ -12,6 +12,7 @@ module MOM_stoch_eos
 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
 
 implicit none
@@ -40,6 +41,7 @@ module MOM_stoch_eos
   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 
+  real :: kappa_smooth    !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1]
   integer :: id_stoch_eos  = -1, id_stoch_phi  = -1, id_tvar_sgs = -1
 end type MOM_stoch_eos_CS
 
@@ -58,16 +60,20 @@ subroutine MOM_stoch_eos_init(G,Time,param_file,stoch_eos_CS,restart_CS,diag)
   seed=0
   ! contants
   !pi=2*acos(0.0)
-  call get_param(param_file, "MOM", "STOCH_EOS", stoch_eos_CS%use_stoch_eos, &
+  call get_param(param_file, "MOM_stoch_eos", "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", "STANLEY_COEFF", stoch_eos_CS%stanley_coeff, &
+  call get_param(param_file, "MOM_stoch_eos", "STANLEY_COEFF", stoch_eos_CS%stanley_coeff, &
                  "Coefficient correlating the temperature gradient "//&
                  "and SGS T variance.", default=-1.0)
-  call get_param(param_file, "MOM", "STANLEY_A", stoch_eos_CS%stanley_a, &
+  call get_param(param_file, "MOM_stoch_eos", "STANLEY_A", stoch_eos_CS%stanley_a, &
                  "Coefficient a which scales chi in stochastic perturbation of the "//&
                  "SGS T variance.", default=1.0)
-
+  call get_param(param_file, "MOM_stoch_eos", "KD_SMOOTH", stoch_eos_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)
+
   !don't run anything if STANLEY_COEFF < 0
   if (stoch_eos_CS%stanley_coeff >= 0.0) then
 
@@ -77,7 +83,7 @@ subroutine MOM_stoch_eos_init(G,Time,param_file,stoch_eos_CS,restart_CS,diag)
     ALLOC_(stoch_eos_CS%phi(G%isd:G%ied,G%jsd:G%jed)) ; stoch_eos_CS%phi(:,:) = 0.0
     ALLOC_(l2_inv(G%isd:G%ied,G%jsd:G%jed)) 
     ALLOC_(rgauss(G%isd:G%ied,G%jsd:G%jed)) 
-    call get_param(param_file, "MOM", "SEED_STOCH_EOS", seed, &
+    call get_param(param_file, "MOM_stoch_eos", "SEED_STOCH_EOS", seed, &
                  "Specfied seed for random number sequence ", default=0)
     call random_2d_constructor(rn_CS, G%HI, Time, seed)
     call random_2d_norm(rn_CS, G%HI, rgauss)
@@ -98,13 +104,13 @@ subroutine MOM_stoch_eos_init(G,Time,param_file,stoch_eos_CS,restart_CS,diag)
 
     !register diagnostics
     stoch_eos_CS%id_tvar_sgs = register_diag_field('ocean_model', 'tvar_sgs', diag%axesTL, Time, &
-      'Parameterized SGS Temperature Variance ', 'None')	
-	if (stoch_eos_CS%use_stoch_eos) then
+      'Parameterized SGS Temperature Variance ', 'None')
+    if (stoch_eos_CS%use_stoch_eos) then
       stoch_eos_CS%id_stoch_eos = register_diag_field('ocean_model', 'stoch_eos', diag%axesT1, Time, &
         'random pattern for EOS', 'None')
       stoch_eos_CS%id_stoch_phi = register_diag_field('ocean_model', 'stoch_phi', diag%axesT1, Time, &
         'phi for EOS', 'None')
-	endif
+    endif
   endif
 
   end subroutine MOM_stoch_eos_init
@@ -153,13 +159,19 @@ subroutine MOM_stoch_eos_run(G,u,v,delt,Time,stoch_eos_CS,diag)
   end subroutine MOM_stoch_eos_run
 
 
-  subroutine MOM_calc_varT(G,Gv,h,tv,stoch_eos_CS)
+  subroutine MOM_calc_varT(G,GV,h,tv,stoch_eos_CS,dt)
   type(ocean_grid_type), intent(in)    :: G
   type(verticalGrid_type),                   intent(in)  :: GV  !< Vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)  :: h   !< Layer thickness [H ~> m]
   type(thermo_var_ptrs), intent(inout) :: tv   !< Thermodynamics structure
   type(MOM_stoch_eos_CS), intent(inout) :: stoch_eos_CS
+  real, intent(in)  :: dt    !< Time increment [T ~> s]
 ! locals
+  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: &
+    T, &          ! The temperature (or density) [degC], with the values in
+                  ! in massless layers filled vertically by diffusion.
+    S             ! The filled salinity [ppt], with the values in
+                  ! in massless layers filled vertically by diffusion.
   integer                                ::  i,j,k
   real :: Tl(5)              ! copy and T in local stencil [degC]
   real :: mn_T               ! mean of T in local stencil [degC]
@@ -171,6 +183,9 @@ subroutine MOM_calc_varT(G,Gv,h,tv,stoch_eos_CS)
   ! 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)
+
+  call vert_fill_TS(h, tv%T, tv%S, stoch_eos_CS%kappa_smooth*dt, T, S, G, GV, halo_here=1, larger_h_denom=.true.)
+
   do k=1,G%ke
      do j=G%jsc,G%jec
         do i=G%isc,G%iec
@@ -181,8 +196,8 @@ subroutine MOM_calc_varT(G,Gv,h,tv,stoch_eos_CS)
            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)
+           Tl(1) = T(i,j,k) ; Tl(2) = T(i-1,j,k) ; Tl(3) = T(i+1,j,k)
+           Tl(4) = T(i,j-1,k) ; Tl(5) = 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.

src/parameterizations/lateral/MOM_mixed_layer_restrat.F90

@@ -59,7 +59,9 @@ module MOM_mixed_layer_restrat
   logical :: debug = .false.       !< If true, calculate checksums of fields for debugging.
   type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
                                    !! timing of diagnostic output.
-
+  logical :: use_stanley_ml        !< If true, use the Stanley parameterization of SGS T variance
+                                   !! if false, MLE will calculate a MLD based on a density difference
+                                   !! based on the parameter MLE_DENSITY_DIFF.
   real, dimension(:,:), pointer :: &
          MLD_filtered => NULL(), &   !< Time-filtered MLD [H ~> m or kg m-2]
          MLD_filtered_slow => NULL() !< Slower time-filtered MLD [H ~> m or kg m-2]
@@ -177,6 +179,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
   real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer
                                        ! densities [R L2 T-2 ~> Pa].
   real, dimension(SZI_(G)) :: rhoAtK, rho1, d1, pRef_N2 ! Used for N2
+  real, dimension(SZI_(G)) :: covTS, varS !SGS TS covariance, S variance in Stanley param; currently 0
   real :: aFac, bFac ! Nondimensional ratios [nondim]
   real :: ddRho    ! A density difference [R ~> kg m-3]
   real :: hAtVel, zpa, zpb, dh, res_scaling_fac
@@ -196,6 +199,9 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
 
   is  = G%isc  ; ie  = G%iec  ; js  = G%jsc  ; je  = G%jec ; nz = GV%ke
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
+
+  covTS(:)=0.0 !!Functionality not implemented yet; in future, should be passed in tv
+  varS(:)=0.0
 
   if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// &
          "An equation of state must be used with this module.")
@@ -208,15 +214,25 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
     EOSdom(:) = EOS_domain(G%HI, halo=1)
     do j = js-1, je+1
       dK(:) = 0.5 * h(:,j,1) ! Depth of center of surface layer
-      call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom)
+      if (CS%use_stanley_ml) then
+        call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, tv%varT(:,j,1), covTS, varS, &
+          rhoSurf, tv%eqn_of_state, EOSdom)   
+      else
+        call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom)
+      endif
       deltaRhoAtK(:) = 0.
       MLD_fast(:,j) = 0.
       do k = 2, nz
         dKm1(:) = dK(:) ! Depth of center of layer K-1
         dK(:) = dK(:) + 0.5 * ( h(:,j,k) + h(:,j,k-1) ) ! Depth of center of layer K
         ! Mixed-layer depth, using sigma-0 (surface reference pressure)
         deltaRhoAtKm1(:) = deltaRhoAtK(:) ! Store value from previous iteration of K
-        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, tv%eqn_of_state, EOSdom)
+        if (CS%use_stanley_ml) then
+          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, tv%varT(:,j,k), covTS, varS, &
+            deltaRhoAtK, tv%eqn_of_state, EOSdom) 
+        else
+          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, tv%eqn_of_state, EOSdom)
+        endif
         do i = is-1,ie+1
           deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface
         enddo
@@ -303,6 +319,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
 !$OMP                               h_neglect,g_Rho0,I4dt,CS,uhml,uhtr,dt,vhml,vhtr,EOSdom,   &
 !$OMP                               utimescale_diag,vtimescale_diag,forces,dz_neglect, &
 !$OMP                               htot_slow,MLD_slow,Rml_av_slow,VarMix,I_LFront,    &
+!$OMP                               covTS, varS,                                       &
 !$OMP                               res_upscale, nz,MLD_fast,uDml_diag,vDml_diag)      &
 !$OMP                       private(rho_ml,h_vel,u_star,absf,mom_mixrate,timescale,    &
 !$OMP                               line_is_empty, keep_going,res_scaling_fac,         &
@@ -320,7 +337,12 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var
         h_avail(i,j,k) = max(I4dt*G%areaT(i,j)*(h(i,j,k)-GV%Angstrom_H),0.0)
       enddo
       if (keep_going) then
-        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), tv%eqn_of_state, EOSdom)
+        if (CS%use_stanley_ml) then
+          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, tv%varT(:,j,k), covTS, varS, &
+            rho_ml(:), tv%eqn_of_state, EOSdom)
+        else   
+          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), tv%eqn_of_state, EOSdom)
+        endif
         line_is_empty = .true.
         do i=is-1,ie+1
           if (htot_fast(i,j) < MLD_fast(i,j)) then
@@ -617,6 +639,10 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// &
          "Module must be initialized before it is used.")
   if ((nkml<2) .or. (CS%ml_restrat_coef<=0.0)) return
+
+  if (CS%use_stanley_ml) call MOM_error(FATAL, &
+       "MOM_mixedlayer_restrat: The Stanley parameterization is not"//&
+       "available with the BML.")  
 
   uDml(:)    = 0.0 ; vDml(:) = 0.0
   I4dt       = 0.25 / dt
@@ -843,6 +869,9 @@ 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. (2010)", units="nondim", default=0.0)
+  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.)
   ! We use GV%nkml to distinguish between the old and new implementation of MLE.
   ! The old implementation only works for the layer model with nkml>0.
   if (GV%nkml==0) then