MOM_cvmix_shear.F90

Go to the documentation of this file.

!> Interface to CVMix interior shear schemes
module mom_cvmix_shear

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

!---------------------------------------------------
! module MOM_cvmix_shear
! Author: Brandon Reichl
! Date: Aug 31, 2016
! Purpose: Interface to CVMix interior shear schemes
! Further information to be added at a later time.
!---------------------------------------------------

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_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, eos_type
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

!> Control structure including parameters for CVMix interior shear schemes.
type, public :: cvmix_shear_cs
  logical :: use_lmd94, use_pp81            !< Flags for various schemes
  real    :: ri_zero                        !< LMD94 critical Richardson number
  real    :: nu_zero                        !< LMD94 maximum interior diffusivity
  real    :: kpp_exp                        !<
  real, allocatable, dimension(:,:,:) :: n2 !< Squared Brunt-Vaisala frequency (1/s2)
  real, allocatable, dimension(:,:,:) :: s2 !< Squared shear frequency (1/s2)
  character(10) :: mix_scheme               !< Mixing scheme name (string)
end type cvmix_shear_cs

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

contains

!> Subroutine for calculating (internal) diffusivity
subroutine calculate_cvmix_shear(u_H, v_H, h, tv, KH,  &
                                 KM, G, GV, CS )
  type(ocean_grid_type),                      intent(in)  :: G !< Grid structure.
  type(verticalgrid_type),                    intent(in)  :: GV !< Vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: u_H !< Initial zonal velocity on T points, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: v_H !< Initial meridional velocity on T points, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h !< Layer thickness, in 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) :: KH !< The vertical viscosity at each interface
                                                                !! (not layer!) in m2 s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: KM !< The vertical viscosity at each interface
                                                                !! (not layer!) in 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
  real :: pref, DU, DV, DRHO, DZ, N2, S2
  real, dimension(2*(G%ke)) :: pres_1d, temp_1d, salt_1d, rho_1d
  real, dimension(G%ke+1) ::  Ri_Grad !< Gradient Richardson number

  ! some constants
  gorho = gv%g_Earth / gv%Rho0

  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.
      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%H_to_Pa * h(i,j,k)

      enddo ! k-loop finishes

      ! compute in-situ density
      call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, 1, 2*g%ke, 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_m)
        n2 = drho/dz
        s2 = (du*du+dv*dv)/(dz*dz)
        ri_grad(k) = max(0.,n2)/max(s2,1.e-16)
      enddo

      ! Call to CVMix wrapper for computing interior mixing coefficients.
      call  cvmix_coeffs_shear(mdiff_out=km(i,j,:), &
                                   tdiff_out=kh(i,j,:), &
                                   rich=ri_grad, &
                                   nlev=g%ke,    &
                                   max_nlev=g%ke)
    enddo
  enddo

end subroutine calculate_cvmix_shear


!> Initialized the cvmix internal shear mixing routine.
!! \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, 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(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 log_version(param_file, mdl, version, &
    "Parameterization of shear-driven turbulence via CVMix (various options)")
  call get_param(param_file, mdl, "USE_LMD94", cs%use_LMD94, &
                 "If true, use the Large-McWilliams-Doney (JGR 1994) \n"//&
                 "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) \n"//&
                 "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).gt.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.7)
  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 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)
  ! Allocation and initialization
  allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) );cs%N2(:,:,:) = 0.
  allocate( cs%S2( szi_(g), szj_(g), szk_(g)+1 ) );cs%S2(:,:,:) = 0.

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

end module mom_cvmix_shear