MOM_CVMix_shear.F90

!> Interface to CVMix interior shear schemes
module mom_cvmix_shear

! This file is part of MOM6. See LICENSE.md for the license.

!> \author Brandon Reichl

use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : diag_ctrl, time_type
use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_grid, only : ocean_grid_type
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density
use cvmix_shear, only : cvmix_init_shear, cvmix_coeffs_shear
use mom_kappa_shear, only : kappa_shear_is_used
implicit none ; private

#include <MOM_memory.h>

public calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_is_used, cvmix_shear_end

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

!> Control structure including parameters for CVMix interior shear schemes.
type, public :: cvmix_shear_cs ! TODO: private
  logical :: use_lmd94                      !< Flags to use the LMD94 scheme
  logical :: use_pp81                       !< Flags to use Pacanowski and Philander (JPO 1981)
  logical :: smooth_ri                      !< If true, smooth Ri using a 1-2-1 filter
  real    :: ri_zero                        !< LMD94 critical Richardson number
  real    :: nu_zero                        !< LMD94 maximum interior diffusivity
  real    :: kpp_exp                        !< Exponent of unitless factor of diff.
                                            !! for KPP internal shear mixing scheme.
  real, allocatable, dimension(:,:,:) :: n2 !< Squared Brunt-Vaisala frequency [T-2 ~> s-2]
  real, allocatable, dimension(:,:,:) :: s2 !< Squared shear frequency [T-2 ~> s-2]
  real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number
  real, allocatable, dimension(:,:,:) :: ri_grad_smooth !< Gradient Richardson number
                                                        !! after smoothing
  character(10) :: mix_scheme               !< Mixing scheme name (string)

  type(diag_ctrl), pointer :: diag => null() !< Pointer to the diagnostics control structure
  !>@{ Diagnostic handles
  integer :: id_n2 = -1, id_s2 = -1, id_ri_grad = -1, id_kv = -1, id_kd = -1
  integer :: id_ri_grad_smooth = -1
  !>@}

end type cvmix_shear_cs

character(len=40)  :: mdl = "MOM_CVMix_shear"  !< This module's name.

contains

!> Subroutine for calculating (internal) vertical diffusivities/viscosities
subroutine calculate_cvmix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS )
  type(ocean_grid_type),                      intent(in)  :: G   !< Grid structure.
  type(verticalgrid_type),                    intent(in)  :: GV  !< Vertical grid structure.
  type(unit_scale_type),                      intent(in)  :: US  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: u_H !< Initial zonal velocity on T points [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: v_H !< Initial meridional velocity on T
                                                                 !! points [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h   !< Layer thickness [H ~> m or kg m-2].
  type(thermo_var_ptrs),                      intent(in)  :: tv  !< Thermodynamics structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd  !< The vertical diffusivity at each interface
                                                                 !! (not layer!) [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv  !< The vertical viscosity at each interface
                                                                 !! (not layer!) [Z2 T-1 ~> m2 s-1].
  type(cvmix_shear_cs),                       pointer     :: CS  !< The control structure returned by a previous
                                                                 !! call to CVMix_shear_init.
  ! Local variables
  integer :: i, j, k, kk, km1
  real :: GoRho  ! Gravitational acceleration divided by density [Z T-2 R-1 ~> m4 s-2 kg-2]
  real :: pref   ! Interface pressures [R L2 T-2 ~> Pa]
  real :: DU, DV ! Velocity differences [L T-1 ~> m s-1]
  real :: DZ     ! Grid spacing around an interface [Z ~> m]
  real :: N2     ! Buoyancy frequency at an interface [T-2 ~> s-2]
  real :: S2     ! Shear squared at an interface [T-2 ~> s-2]
  real :: dummy  ! A dummy variable [nondim]
  real :: dRho   ! Buoyancy differences [Z T-2 ~> m s-2]
  real, dimension(2*(G%ke)) :: pres_1d ! A column of interface pressures [R L2 T-2 ~> Pa]
  real, dimension(2*(G%ke)) :: temp_1d ! A column of temperatures [degC]
  real, dimension(2*(G%ke)) :: salt_1d ! A column of salinities [ppt]
  real, dimension(2*(G%ke)) :: rho_1d  ! A column of densities at interface pressures [R ~> kg m-3]
  real, dimension(G%ke+1) :: Ri_Grad !< Gradient Richardson number [nondim]
  real, dimension(G%ke+1) :: Kvisc   !< Vertical viscosity at interfaces [m2 s-1]
  real, dimension(G%ke+1) :: Kdiff   !< Diapycnal diffusivity at interfaces [m2 s-1]
  real :: epsln  !< Threshold to identify vanished layers [H ~> m or kg m-2]

  ! some constants
  gorho = us%L_to_Z**2 * gv%g_Earth / gv%Rho0
  epsln = 1.e-10 * gv%m_to_H

  do j = g%jsc, g%jec
    do i = g%isc, g%iec

      ! skip calling for land points
      if (g%mask2dT(i,j)==0.) cycle

      ! Richardson number computed for each cell in a column.
      pref = 0. ; if (associated(tv%p_surf)) pref = tv%p_surf(i,j)
      ri_grad(:)=1.e8 !Initialize w/ large Richardson value
      do k=1,g%ke
        ! pressure, temp, and saln for EOS
        ! kk+1 = k fields
        ! kk+2 = km1 fields
        km1  = max(1, k-1)
        kk   = 2*(k-1)
        pres_1d(kk+1) = pref
        pres_1d(kk+2) = pref
        temp_1d(kk+1) = tv%T(i,j,k)
        temp_1d(kk+2) = tv%T(i,j,km1)
        salt_1d(kk+1) = tv%S(i,j,k)
        salt_1d(kk+2) = tv%S(i,j,km1)

        ! pRef is pressure at interface between k and km1.
        ! iterate pRef for next pass through k-loop.
        pref = pref + (gv%g_Earth * gv%H_to_RZ) * h(i,j,k)

      enddo ! k-loop finishes

      ! compute in-situ density [R ~> kg m-3]
      call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, tv%eqn_of_state)

      ! N2 (can be negative) on interface
      do k = 1, g%ke
        km1 = max(1, k-1)
        kk = 2*(k-1)
        du = u_h(i,j,k) - u_h(i,j,km1)
        dv = v_h(i,j,k) - v_h(i,j,km1)
        drho = gorho * (rho_1d(kk+1) - rho_1d(kk+2))
        dz = (0.5*(h(i,j,km1) + h(i,j,k))+gv%H_subroundoff)*gv%H_to_Z
        n2 = drho / dz
        s2 = us%L_to_Z**2*(du*du+dv*dv)/(dz*dz)
        ri_grad(k) = max(0., n2) / max(s2, 1.e-10*us%T_to_s**2)

        ! fill 3d arrays, if user asks for diagsnostics
        if (cs%id_N2 > 0) cs%N2(i,j,k) = n2
        if (cs%id_S2 > 0) cs%S2(i,j,k) = s2

      enddo

      ri_grad(g%ke+1) = ri_grad(g%ke)

      if (cs%id_ri_grad > 0) cs%ri_grad(i,j,:) = ri_grad(:)

      if (cs%smooth_ri) then
        ! 1) fill Ri_grad in vanished layers with adjacent value
        do k = 2, g%ke
          if (h(i,j,k) <= epsln) ri_grad(k) = ri_grad(k-1)
        enddo

        ri_grad(g%ke+1) = ri_grad(g%ke)

        ! 2) vertically smooth Ri with 1-2-1 filter
        dummy =  0.25 * ri_grad(2)
        ri_grad(g%ke+1) = ri_grad(g%ke)
        do k = 3, g%ke
          ri_grad(k) = dummy + 0.5 * ri_grad(k) + 0.25 * ri_grad(k+1)
          dummy = 0.25 * ri_grad(k)
        enddo

        if (cs%id_ri_grad_smooth > 0) cs%ri_grad_smooth(i,j,:) = ri_grad(:)
      endif

      do k=1,g%ke+1
        kvisc(k) = us%Z2_T_to_m2_s * kv(i,j,k)
        kdiff(k) = us%Z2_T_to_m2_s * kd(i,j,k)
      enddo

      ! Call to CVMix wrapper for computing interior mixing coefficients.
      call  cvmix_coeffs_shear(mdiff_out=kvisc(:), &
                                   tdiff_out=kdiff(:), &
                                   rich=ri_grad(:), &
                                   nlev=g%ke,    &
                                   max_nlev=g%ke)
      do k=1,g%ke+1
        kv(i,j,k) = us%m2_s_to_Z2_T * kvisc(k)
        kd(i,j,k) = us%m2_s_to_Z2_T * kdiff(k)
      enddo
    enddo
  enddo

  ! write diagnostics
  if (cs%id_kd > 0) call post_data(cs%id_kd, kd, cs%diag)
  if (cs%id_kv > 0) call post_data(cs%id_kv, kv, cs%diag)
  if (cs%id_N2 > 0) call post_data(cs%id_N2, cs%N2, cs%diag)
  if (cs%id_S2 > 0) call post_data(cs%id_S2, cs%S2, cs%diag)
  if (cs%id_ri_grad > 0) call post_data(cs%id_ri_grad, cs%ri_grad, cs%diag)
  if (cs%id_ri_grad_smooth > 0) call post_data(cs%id_ri_grad_smooth ,cs%ri_grad_smooth, cs%diag)

end subroutine calculate_cvmix_shear


!> Initialized the CVMix internal shear mixing routine.
!! \todo Does this note require emphasis?
!! \note *This is where we test to make sure multiple internal shear
!!       mixing routines (including JHL) are not enabled at the same time.*
!! (returns) CVMix_shear_init - True if module is to be used, False otherwise
logical function cvmix_shear_init(Time, G, GV, US, param_file, diag, CS)
  type(time_type),         intent(in)    :: Time !< The current time.
  type(ocean_grid_type),   intent(in)    :: G  !< Grid structure.
  type(verticalgrid_type), intent(in)    :: GV !< Vertical grid structure.
  type(unit_scale_type),   intent(in)    :: US !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< Run-time parameter file handle
  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure.
  type(cvmix_shear_cs),    pointer       :: CS !< This module's control structure.
  ! Local variables
  integer :: NumberTrue=0
  logical :: use_JHL
! This include declares and sets the variable "version".
#include "version_variable.h"

  if (associated(cs)) then
    call mom_error(warning, "CVMix_shear_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)

! Set default, read and log parameters
  call get_param(param_file, mdl, "USE_LMD94", cs%use_LMD94, default=.false., do_not_log=.true.)
  call get_param(param_file, mdl, "USE_PP81", cs%use_PP81, default=.false., do_not_log=.true.)
  call log_version(param_file, mdl, version, &
           "Parameterization of shear-driven turbulence via CVMix (various options)", &
            all_default=.not.(cs%use_PP81.or.cs%use_LMD94))
  call get_param(param_file, mdl, "USE_LMD94", cs%use_LMD94, &
                 "If true, use the Large-McWilliams-Doney (JGR 1994) "//&
                 "shear mixing parameterization.", default=.false.)
  if (cs%use_LMD94) then
     numbertrue=numbertrue + 1
     cs%Mix_Scheme='KPP'
  endif
  call get_param(param_file, mdl, "USE_PP81", cs%use_PP81, &
                 "If true, use the Pacanowski and Philander (JPO 1981) "//&
                 "shear mixing parameterization.", default=.false.)
  if (cs%use_PP81) then
     numbertrue = numbertrue + 1
     cs%Mix_Scheme='PP'
  endif
  use_jhl=kappa_shear_is_used(param_file)
  if (use_jhl) numbertrue = numbertrue + 1
  ! After testing for interior schemes, make sure only 0 or 1 are enabled.
  ! Otherwise, warn user and kill job.
  if ((numbertrue) > 1) then
     call mom_error(fatal, 'MOM_CVMix_shear_init: '// &
           'Multiple shear driven internal mixing schemes selected,'//&
           ' please disable all but one scheme to proceed.')
  endif
  cvmix_shear_init=(cs%use_PP81.or.cs%use_LMD94)

! Forego remainder of initialization if not using this scheme
  if (.not. cvmix_shear_init) return
  call get_param(param_file, mdl, "NU_ZERO", cs%Nu_Zero, &
                 "Leading coefficient in KPP shear mixing.", &
                 units="nondim", default=5.e-3)
  call get_param(param_file, mdl, "RI_ZERO", cs%Ri_Zero, &
                 "Critical Richardson for KPP shear mixing, "// &
                 "NOTE this the internal mixing and this is "// &
                 "not for setting the boundary layer depth." &
                 ,units="nondim", default=0.8)
  call get_param(param_file, mdl, "KPP_EXP", cs%KPP_exp, &
                 "Exponent of unitless factor of diffusivities, "// &
                 "for KPP internal shear mixing scheme." &
                 ,units="nondim", default=3.0)
  call get_param(param_file, mdl, "SMOOTH_RI", cs%smooth_ri, &
                 "If true, vertically smooth the Richardson "// &
                 "number by applying a 1-2-1 filter once.", &
                 default = .false.)
  call cvmix_init_shear(mix_scheme=cs%Mix_Scheme, &
                        kpp_nu_zero=cs%Nu_Zero,   &
                        kpp_ri_zero=cs%Ri_zero,   &
                        kpp_exp=cs%KPP_exp)

  ! Register diagnostics; allocation and initialization
  cs%diag => diag

  cs%id_N2 = register_diag_field('ocean_model', 'N2_shear', diag%axesTi, time, &
      'Square of Brunt-Vaisala frequency used by MOM_CVMix_shear module', '1/s2', conversion=us%s_to_T**2)
  if (cs%id_N2 > 0) then
    allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) ) ; cs%N2(:,:,:) = 0.
  endif

  cs%id_S2 = register_diag_field('ocean_model', 'S2_shear', diag%axesTi, time, &
      'Square of vertical shear used by MOM_CVMix_shear module','1/s2', conversion=us%s_to_T**2)
  if (cs%id_S2 > 0) then
    allocate( cs%S2( szi_(g), szj_(g), szk_(g)+1 ) ) ; cs%S2(:,:,:) = 0.
  endif

  cs%id_ri_grad = register_diag_field('ocean_model', 'ri_grad_shear', diag%axesTi, time, &
      'Gradient Richarson number used by MOM_CVMix_shear module','nondim')
  if (cs%id_ri_grad > 0) then !Initialize w/ large Richardson value
    allocate( cs%ri_grad( szi_(g), szj_(g), szk_(g)+1 )) ; cs%ri_grad(:,:,:) = 1.e8
  endif

  cs%id_ri_grad_smooth = register_diag_field('ocean_model', 'ri_grad_shear_smooth', &
       diag%axesTi, time, &
      'Smoothed gradient Richarson number used by MOM_CVMix_shear module','nondim')
  if (cs%id_ri_grad_smooth > 0) then !Initialize w/ large Richardson value
    allocate( cs%ri_grad_smooth( szi_(g), szj_(g), szk_(g)+1 )) ; cs%ri_grad_smooth(:,:,:) = 1.e8
  endif

  cs%id_kd = register_diag_field('ocean_model', 'kd_shear_CVMix', diag%axesTi, time, &
      'Vertical diffusivity added by MOM_CVMix_shear module', 'm2/s', conversion=us%Z2_T_to_m2_s)
  cs%id_kv = register_diag_field('ocean_model', 'kv_shear_CVMix', diag%axesTi, time, &
      'Vertical viscosity added by MOM_CVMix_shear module', 'm2/s', conversion=us%Z2_T_to_m2_s)

end function cvmix_shear_init

!> Reads the parameters "LMD94" and "PP81" and returns state.
!!   This function allows other modules to know whether this parameterization will
!! be used without needing to duplicate the log entry.
logical function cvmix_shear_is_used(param_file)
  type(param_file_type), intent(in) :: param_file !< Run-time parameter files handle.
  ! Local variables
  logical :: LMD94, PP81
  call get_param(param_file, mdl, "USE_LMD94", lmd94, &
       default=.false., do_not_log = .true.)
  call get_param(param_file, mdl, "Use_PP81", pp81, &
       default=.false., do_not_log = .true.)
  cvmix_shear_is_used = (lmd94 .or. pp81)
end function cvmix_shear_is_used

!> Clear pointers and dealocate memory
subroutine cvmix_shear_end(CS)
  type(cvmix_shear_cs), pointer :: CS !< Control structure for this module that
                                      !! will be deallocated in this subroutine

  if (.not. associated(cs)) return

  if (cs%id_N2 > 0) deallocate(cs%N2)
  if (cs%id_S2 > 0) deallocate(cs%S2)
  if (cs%id_ri_grad > 0) deallocate(cs%ri_grad)
  deallocate(cs)

end subroutine cvmix_shear_end

end module mom_cvmix_shear