MOM_bkgnd_mixing.F90

!> Interface to background mixing schemes, including the Bryan and Lewis (1979)
!! which is applied via CVMix.
 
module mom_bkgnd_mixing
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_debugging,       only : hchksum
use mom_diag_mediator,   only : diag_ctrl, time_type, register_diag_field
use mom_diag_mediator,   only : post_data
use mom_eos,             only : calculate_density, calculate_density_derivs
use mom_error_handler,   only : mom_error, fatal, warning, note
use mom_file_parser,     only : get_param, log_version, param_file_type
use mom_file_parser,     only : openparameterblock, closeparameterblock
use mom_forcing_type,    only : forcing
use mom_grid,            only : ocean_grid_type
use mom_unit_scaling,    only : unit_scale_type
use mom_verticalgrid,    only : verticalgrid_type
use mom_variables,       only : thermo_var_ptrs,  vertvisc_type
use mom_intrinsic_functions, only : invcosh
use cvmix_background,    only : cvmix_init_bkgnd, cvmix_coeffs_bkgnd
 
implicit none ; private
 
#include <MOM_memory.h>
 
public bkgnd_mixing_init
public bkgnd_mixing_end
public calculate_bkgnd_mixing
public sfc_bkgnd_mixing
 
! 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 this module.
type, public :: bkgnd_mixing_cs  ! TODO: private
 
  ! Parameters
  real    :: bryan_lewis_c1         !< The vertical diffusivity values for  Bryan-Lewis profile
                                    !! at |z|=D [m2 s-1]
  real    :: bryan_lewis_c2         !< The amplitude of variation in diffusivity for the
                                    !! Bryan-Lewis diffusivity profile [m2 s-1]
  real    :: bryan_lewis_c3         !< The inverse length scale for transition region in the
                                    !! Bryan-Lewis diffusivity profile [m-1]
  real    :: bryan_lewis_c4         !< The depth where diffusivity is Bryan_Lewis_bl1 in the
                                    !! Bryan-Lewis profile [m]
  real    :: bckgrnd_vdc1           !< Background diffusivity (Ledwell) when
                                    !! horiz_varying_background=.true. [Z2 T-1 ~> m2 s-1]
  real    :: bckgrnd_vdc_eq         !< Equatorial diffusivity (Gregg) when
                                    !! horiz_varying_background=.true. [Z2 T-1 ~> m2 s-1]
  real    :: bckgrnd_vdc_psim       !< Max. PSI induced diffusivity (MacKinnon) when
                                    !! horiz_varying_background=.true. [Z2 T-1 ~> m2 s-1]
  real    :: bckgrnd_vdc_banda      !< Banda Sea diffusivity (Gordon) when
                                    !! horiz_varying_background=.true. [Z2 T-1 ~> m2 s-1]
  real    :: kd_min                 !< minimum diapycnal diffusivity [Z2 T-1 ~> m2 s-1]
  real    :: kd                     !< interior diapycnal diffusivity [Z2 T-1 ~> m2 s-1]
  real    :: omega                  !< The Earth's rotation rate [T-1 ~> s-1].
  real    :: n0_2omega              !< ratio of the typical Buoyancy frequency to
                                    !! twice the Earth's rotation period, used with the
                                    !! Henyey scaling from the mixing
  real    :: prandtl_bkgnd          !< Turbulent Prandtl number used to convert
                                    !! vertical background diffusivity into viscosity
  real    :: kd_tanh_lat_scale      !< A nondimensional scaling for the range of
                                    !! diffusivities with Kd_tanh_lat_fn. Valid values
                                    !! are in the range of -2 to 2; 0.4 reproduces CM2M.
  real    :: kdml                   !< mixed layer diapycnal diffusivity [Z2 T-1 ~> m2 s-1]
                                    !! when bulkmixedlayer==.false.
  real    :: hmix                   !< mixed layer thickness [Z ~> m] when bulkmixedlayer==.false.
  logical :: kd_tanh_lat_fn         !< If true, use the tanh dependence of Kd_sfc on
                                    !! latitude, like GFDL CM2.1/CM2M.  There is no
                                    !! physical justification for this form, and it can
                                    !! not be used with Henyey_IGW_background.
  logical :: bryan_lewis_diffusivity!< If true, background vertical diffusivity
                                    !! uses Bryan-Lewis (1979) like tanh profile.
  logical :: horiz_varying_background !< If true, apply vertically uniform, latitude-dependent
                                    !! background diffusivity, as described in Danabasoglu et al., 2012
  logical :: henyey_igw_background  !< If true, use a simplified variant of the
             !! Henyey et al, JGR (1986) latitudinal scaling for the background diapycnal diffusivity,
             !! which gives a marked decrease in the diffusivity near the equator.  The simplification
             !! here is to assume that the in-situ stratification is the same as the reference stratificaiton.
  logical :: henyey_igw_background_new !< same as Henyey_IGW_background
             !! but incorporate the effect of stratification on TKE dissipation,
             !! e = f/f_0 * acosh(N/f) / acosh(N_0/f_0) * e_0
             !! where e is the TKE dissipation, and N_0 and f_0
             !! are the reference buoyancy frequency and inertial frequencies respectively.
             !! e_0 is the reference dissipation at (N_0,f_0). In the previous version, N=N_0.
             !! Additionally, the squared inverse relationship between  diapycnal diffusivities
             !! and stratification is included:
             !!
             !! kd = e/N^2
             !!
             !! where kd is the diapycnal diffusivity. This approach assumes that work done
             !! against gravity is uniformly distributed throughout the column. Whereas, kd=kd_0*e,
             !! as in the original version, concentrates buoyancy work in regions of strong stratification.
  logical :: bulkmixedlayer !< If true, a refined bulk mixed layer scheme is used
  logical :: debug !< If true, turn on debugging in this module
  ! Daignostic handles and pointers
  type(diag_ctrl), pointer :: diag => null() !< A structure that regulates diagnostic output
  integer :: id_kd_bkgnd = -1 !< Diagnotic IDs
  integer :: id_kv_bkgnd = -1 !< Diagnostic IDs
 
  real, allocatable, dimension(:,:)   ::  kd_sfc !< surface value of the diffusivity [Z2 T-1 ~> m2 s-1]
  ! Diagnostics arrays
  real, allocatable, dimension(:,:,:) :: kd_bkgnd !< Background diffusivity [Z2 T-1 ~> m2 s-1]
  real, allocatable, dimension(:,:,:) :: kv_bkgnd !< Background viscosity  [Z2 s-1 ~> m2 s-1]
 
  character(len=40)  :: bkgnd_scheme_str = "none" !< Background scheme identifier
 
end type bkgnd_mixing_cs
 
character(len=40)  :: mdl = "MOM_bkgnd_mixing" !< This module's name.
 
contains
 
!> Initialize the background mixing routine.
subroutine bkgnd_mixing_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(bkgnd_mixing_cs),    pointer      :: cs         !< This module's control structure.
 
  ! Local variables
  real :: kv                    ! The interior vertical viscosity [Z2 T-1 ~> m2 s-1] - read to set prandtl
                                ! number unless it is provided as a parameter
  real :: prandtl_bkgnd_comp    ! Kv/CS%Kd. Gets compared with user-specified prandtl_bkgnd.
 
! This include declares and sets the variable "version".
#include "version_variable.h"
 
  if (associated(cs)) then
    call mom_error(warning, "bkgnd_mixing_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)
 
  ! Read parameters
  call log_version(param_file, mdl, version, &
    "Adding static vertical background mixing coefficients")
 
  call get_param(param_file, mdl, "KD", cs%Kd, &
                 "The background diapycnal diffusivity of density in the "//&
                 "interior. Zero or the molecular value, ~1e-7 m2 s-1, "//&
                 "may be used.", units="m2 s-1", scale=us%m2_s_to_Z2_T, fail_if_missing=.true.)
 
  call get_param(param_file, mdl, "KV", kv, &
                 "The background kinematic viscosity in the interior. "//&
                 "The molecular value, ~1e-6 m2 s-1, may be used.", &
                 units="m2 s-1", scale=us%m2_s_to_Z2_T, fail_if_missing=.true.)
 
  call get_param(param_file, mdl, "KD_MIN", cs%Kd_min, &
                 "The minimum diapycnal diffusivity.", &
                 units="m2 s-1", default=0.01*cs%Kd*us%Z2_T_to_m2_s, scale=us%m2_s_to_Z2_T)
 
  ! The following is needed to set one of the choices of vertical background mixing
 
  ! BULKMIXEDLAYER is not always defined (e.g., CM2G63L), so the following by pass
  ! the need to include BULKMIXEDLAYER in MOM_input
  cs%bulkmixedlayer = (gv%nkml > 0)
  if (cs%bulkmixedlayer) then
    ! Check that Kdml is not set when using bulk mixed layer
    call get_param(param_file, mdl, "KDML", cs%Kdml, default=-1.)
    if (cs%Kdml>0.) call mom_error(fatal, &
                 "bkgnd_mixing_init: KDML cannot be set when using"// &
                 "bulk mixed layer.")
    cs%Kdml = cs%Kd ! This is not used with a bulk mixed layer, but also
                    ! cannot be a NaN.
  else
    call get_param(param_file, mdl, "KDML", cs%Kdml, &
                 "If BULKMIXEDLAYER is false, KDML is the elevated "//&
                 "diapycnal diffusivity in the topmost HMIX of fluid. "//&
                 "KDML is only used if BULKMIXEDLAYER is false.", &
                 units="m2 s-1", default=cs%Kd*us%Z2_T_to_m2_s, scale=us%m2_s_to_Z2_T)
    call get_param(param_file, mdl, "HMIX_FIXED", cs%Hmix, &
                 "The prescribed depth over which the near-surface "//&
                 "viscosity and diffusivity are elevated when the bulk "//&
                 "mixed layer is not used.", units="m", scale=us%m_to_Z, fail_if_missing=.true.)
  endif
 
  call get_param(param_file, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)
 
!  call openParameterBlock(param_file,'MOM_BACKGROUND_MIXING')
 
  call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", cs%Bryan_Lewis_diffusivity, &
                 "If true, use a Bryan & Lewis (JGR 1979) like tanh "//&
                 "profile of background diapycnal diffusivity with depth. "//&
                 "This is done via CVMix.", default=.false.)
 
  if (cs%Bryan_Lewis_diffusivity) then
    call check_bkgnd_scheme(cs, "BRYAN_LEWIS_DIFFUSIVITY")
 
    call get_param(param_file, mdl, "BRYAN_LEWIS_C1", cs%Bryan_Lewis_c1, &
                   "The vertical diffusivity values for Bryan-Lewis profile at |z|=D.", &
                   units="m2 s-1", fail_if_missing=.true.)
 
    call get_param(param_file, mdl, "BRYAN_LEWIS_C2", cs%Bryan_Lewis_c2, &
                   "The amplitude of variation in diffusivity for the Bryan-Lewis profile", &
                   units="m2 s-1", fail_if_missing=.true.)
 
    call get_param(param_file, mdl, "BRYAN_LEWIS_C3", cs%Bryan_Lewis_c3, &
                   "The inverse length scale for transition region in the Bryan-Lewis profile", &
                   units="m-1", fail_if_missing=.true.)
 
    call get_param(param_file, mdl, "BRYAN_LEWIS_C4", cs%Bryan_Lewis_c4, &
                   "The depth where diffusivity is BRYAN_LEWIS_C1 in the Bryan-Lewis profile",&
                   units="m", fail_if_missing=.true.)
 
  endif ! CS%Bryan_Lewis_diffusivity
 
  call get_param(param_file, mdl, "HORIZ_VARYING_BACKGROUND", cs%horiz_varying_background, &
                 "If true, apply vertically uniform, latitude-dependent background "//&
                 "diffusivity, as described in Danabasoglu et al., 2012", &
                 default=.false.)
 
  if (cs%horiz_varying_background) then
    call check_bkgnd_scheme(cs, "HORIZ_VARYING_BACKGROUND")
 
    call get_param(param_file, mdl, "BCKGRND_VDC1", cs%bckgrnd_vdc1, &
                   "Background diffusivity (Ledwell) when HORIZ_VARYING_BACKGROUND=True", &
                   units="m2 s-1",default = 0.16e-04, scale=us%m2_s_to_Z2_T)
 
    call get_param(param_file, mdl, "BCKGRND_VDC_EQ", cs%bckgrnd_vdc_eq, &
                   "Equatorial diffusivity (Gregg) when HORIZ_VARYING_BACKGROUND=True", &
                   units="m2 s-1",default = 0.01e-04, scale=us%m2_s_to_Z2_T)
 
    call get_param(param_file, mdl, "BCKGRND_VDC_PSIM", cs%bckgrnd_vdc_psim, &
                   "Max. PSI induced diffusivity (MacKinnon) when HORIZ_VARYING_BACKGROUND=True", &
                   units="m2 s-1",default = 0.13e-4, scale=us%m2_s_to_Z2_T)
 
    call get_param(param_file, mdl, "BCKGRND_VDC_BAN", cs%bckgrnd_vdc_Banda, &
                   "Banda Sea diffusivity (Gordon) when HORIZ_VARYING_BACKGROUND=True", &
                   units="m2 s-1",default = 1.0e-4, scale=us%m2_s_to_Z2_T)
  endif
 
  call get_param(param_file, mdl, "PRANDTL_BKGND", cs%prandtl_bkgnd, &
                 "Turbulent Prandtl number used to convert vertical "//&
                 "background diffusivities into viscosities.", &
                 units="nondim", default=1.0)
 
  if (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background) then
 
    prandtl_bkgnd_comp = cs%prandtl_bkgnd
    if (cs%Kd /= 0.0) prandtl_bkgnd_comp = kv/cs%Kd
 
    if ( abs(cs%prandtl_bkgnd - prandtl_bkgnd_comp)>1.e-14) then
      call mom_error(fatal,"set_diffusivity_init: The provided KD, KV,"//&
                           "and PRANDTL_BKGND values are incompatible. The following "//&
                           "must hold: KD*PRANDTL_BKGND==KV")
    endif
 
  endif
 
  call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND", cs%Henyey_IGW_background, &
                 "If true, use a latitude-dependent scaling for the near "//&
                 "surface background diffusivity, as described in "//&
                 "Harrison & Hallberg, JPO 2008.", default=.false.)
  if (cs%Henyey_IGW_background) call check_bkgnd_scheme(cs, "HENYEY_IGW_BACKGROUND")
 
 
  call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND_NEW", cs%Henyey_IGW_background_new, &
                 "If true, use a better latitude-dependent scaling for the "//&
                 "background diffusivity, as described in "//&
                 "Harrison & Hallberg, JPO 2008.", default=.false.)
  if (cs%Henyey_IGW_background_new) call check_bkgnd_scheme(cs, "HENYEY_IGW_BACKGROUND_NEW")
 
  if (cs%Kd>0.0 .and. (trim(cs%bkgnd_scheme_str)=="BRYAN_LEWIS_DIFFUSIVITY" .or.&
                          trim(cs%bkgnd_scheme_str)=="HORIZ_VARYING_BACKGROUND" )) then
    call mom_error(warning, "set_diffusivity_init: a nonzero constant background "//&
         "diffusivity (KD) is specified along with "//trim(cs%bkgnd_scheme_str))
  endif
 
  if (cs%Henyey_IGW_background) then
    call get_param(param_file, mdl, "HENYEY_N0_2OMEGA", cs%N0_2Omega, &
                  "The ratio of the typical Buoyancy frequency to twice "//&
                  "the Earth's rotation period, used with the Henyey "//&
                  "scaling from the mixing.", units="nondim", default=20.0)
    call get_param(param_file, mdl, "OMEGA", cs%omega, &
                 "The rotation rate of the earth.", units="s-1", &
                 default=7.2921e-5, scale=us%T_to_s)
  endif
 
  call get_param(param_file, mdl, "KD_TANH_LAT_FN", &
                  cs%Kd_tanh_lat_fn, &
                 "If true, use a tanh dependence of Kd_sfc on latitude, "//&
                 "like CM2.1/CM2M.  There is no physical justification "//&
                 "for this form, and it can not be used with "//&
                 "HENYEY_IGW_BACKGROUND.", default=.false.)
 
  if (cs%Kd_tanh_lat_fn) &
    call get_param(param_file, mdl, "KD_TANH_LAT_SCALE", cs%Kd_tanh_lat_scale, &
                 "A nondimensional scaling for the range ofdiffusivities "//&
                 "with KD_TANH_LAT_FN. Valid values are in the range of "//&
                 "-2 to 2; 0.4 reproduces CM2M.", units="nondim", default=0.0)
 
  if (cs%Henyey_IGW_background .and. cs%Kd_tanh_lat_fn) call mom_error(fatal, &
    "MOM_bkgnd_mixing: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.")
 
!  call closeParameterBlock(param_file)
 
  ! allocate arrays and set them to zero
  allocate(cs%Kd_bkgnd(szi_(g), szj_(g), szk_(g)+1)); cs%kd_bkgnd(:,:,:) = 0.
  allocate(cs%kv_bkgnd(szi_(g), szj_(g), szk_(g)+1)); cs%kv_bkgnd(:,:,:) = 0.
  allocate(cs%Kd_sfc(szi_(g), szj_(g))); cs%Kd_sfc(:,:) = 0.
 
  ! Register diagnostics
  cs%diag => diag
  cs%id_kd_bkgnd = register_diag_field('ocean_model', 'Kd_bkgnd', diag%axesTi, time, &
      'Background diffusivity added by MOM_bkgnd_mixing module', 'm2/s', conversion=us%Z2_T_to_m2_s)
  cs%id_kv_bkgnd = register_diag_field('ocean_model', 'Kv_bkgnd', diag%axesTi, time, &
      'Background viscosity added by MOM_bkgnd_mixing module', 'm2/s', conversion=us%Z2_T_to_m2_s)
 
end subroutine bkgnd_mixing_init
 
!> Get surface vertical background diffusivities/viscosities.
subroutine sfc_bkgnd_mixing(G, US, CS)
 
  type(ocean_grid_type),          intent(in)    :: g  !< Grid structure.
  type(unit_scale_type),          intent(in)    :: us !< A dimensional unit scaling type
  type(bkgnd_mixing_cs), pointer, intent(inout) :: cs !< The control structure returned by
                                                      !! a previous call to bkgnd_mixing_init.
  ! local variables
  real :: i_x30  !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg)))
  real :: deg_to_rad !< factor converting degrees to radians, pi/180.
  real :: abs_sin    !< absolute value of sine of latitude [nondim]
  real :: epsilon
  integer :: i, j, k, is, ie, js, je
 
  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec
 
  ! set some parameters
  deg_to_rad = atan(1.0)/45.0 ! = PI/180
  epsilon = 1.e-10
 
 
  if (.not. (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background)) then
!$OMP parallel do default(none) shared(is,ie,js,je,CS)
    do j=js,je ; do i=is,ie
      cs%Kd_sfc(i,j) = cs%Kd
    enddo ; enddo
  endif
 
  if (cs%Henyey_IGW_background) then
    i_x30 = 2.0 / invcosh(cs%N0_2Omega*2.0) ! This is evaluated at 30 deg.
!$OMP parallel do default(none) &
!$OMP shared(is,ie,js,je,CS,G,deg_to_rad,epsilon,I_x30) &
!$OMP private(abs_sin)
    do j=js,je ; do i=is,ie
      abs_sin = abs(sin(g%geoLatT(i,j)*deg_to_rad))
      cs%Kd_sfc(i,j) = max(cs%Kd_min, cs%Kd_sfc(i,j) * &
           ((abs_sin * invcosh(cs%N0_2Omega/max(epsilon,abs_sin))) * i_x30) )
    enddo ; enddo
  elseif (cs%Kd_tanh_lat_fn) then
!$OMP parallel do default(none) shared(is,ie,js,je,CS,G)
    do j=js,je ; do i=is,ie
      !   The transition latitude and latitude range are hard-scaled here, since
      ! this is not really intended for wide-spread use, but rather for
      ! comparison with CM2M / CM2.1 settings.
      cs%Kd_sfc(i,j) = max(cs%Kd_min, cs%Kd_sfc(i,j) * (1.0 + &
          cs%Kd_tanh_lat_scale * 0.5*tanh((abs(g%geoLatT(i,j)) - 35.0)/5.0) ))
    enddo ; enddo
  endif
 
  if (cs%debug) call hchksum(cs%Kd_sfc,"After sfc_bkgnd_mixing: Kd_sfc",g%HI,haloshift=0, scale=us%Z2_T_to_m2_s)
 
end subroutine sfc_bkgnd_mixing
 
 
!> Calculates the vertical background diffusivities/viscosities
subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kv, j, 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)    :: h   !< Layer thickness [H ~> m or kg m-2].
  type(thermo_var_ptrs),                    intent(in)    :: tv  !< Thermodynamics structure.
  real, dimension(SZI_(G),SZK_(G)),         intent(in)    :: n2_lay !< squared buoyancy frequency associated
                                                                 !! with layers [T-2 ~> s-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: kd_lay !< Diapycnal diffusivity of each layer
                                                                 !! [Z2 T-1 ~> m2 s-1].
  real, dimension(:,:,:),                   pointer       :: kv  !< The "slow" vertical viscosity at each interface
                                                                 !! (not layer!) [Z2 T-1 ~> m2 s-1]
  integer,                                  intent(in)    :: j   !< Meridional grid index
  type(bkgnd_mixing_cs),                    pointer       :: cs  !< The control structure returned by
                                                                 !! a previous call to bkgnd_mixing_init.
 
  ! local variables
  real, dimension(SZK_(G)+1) :: depth_int  !< distance from surface of the interfaces [m]
  real, dimension(SZK_(G)+1) :: kd_col     !< Diffusivities at the interfaces [m2 s-1]
  real, dimension(SZK_(G)+1) :: kv_col     !< Viscosities at the interfaces [m2 s-1]
  real, dimension(SZI_(G)) :: depth        !< distance from surface of an interface [Z ~> m]
  real :: depth_c    !< depth of the center of a layer [Z ~> m]
  real :: i_hmix     !< inverse of fixed mixed layer thickness [Z-1 ~> m-1]
  real :: i_2omega   !< 1/(2 Omega) [T ~> s]
  real :: n_2omega   !  The ratio of the stratification to the Earth's rotation rate [nondim]
  real :: n02_n2     !  The ratio a reference stratification to the actual stratification [nondim]
  real :: i_x30      !< 2/acos(2) = 1/(sin(30 deg) * acosh(1/sin(30 deg)))
  real :: deg_to_rad !< factor converting degrees to radians, pi/180.
  real :: abs_sin    !< absolute value of sine of latitude [nondim]
  real :: epsilon    ! The minimum value of the sine of latitude [nondim]
  real :: bckgrnd_vdc_psin !< PSI diffusivity in northern hemisphere [Z2 T-1 ~> m2 s-1]
  real :: bckgrnd_vdc_psis !< PSI diffusivity in southern hemisphere [Z2 T-1 ~> m2 s-1]
  integer :: i, k, is, ie, js, je, nz
 
  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec ; nz = g%ke
 
  ! set some parameters
  deg_to_rad = atan(1.0)/45.0 ! = PI/180
  epsilon = 1.e-10
 
  ! Set up the background diffusivity.
  if (cs%Bryan_Lewis_diffusivity) then
 
    do i=is,ie
      depth_int(1) = 0.0
      do k=2,nz+1
        depth_int(k) = depth_int(k-1) + gv%H_to_m*h(i,j,k-1)
      enddo
 
      call cvmix_init_bkgnd(max_nlev=nz, &
                            zw = depth_int(:), &  !< interface depths relative to the surface in m, must be positive.
                            bl1 = cs%Bryan_Lewis_c1, &
                            bl2 = cs%Bryan_Lewis_c2, &
                            bl3 = cs%Bryan_Lewis_c3, &
                            bl4 = cs%Bryan_Lewis_c4, &
                            prandtl = cs%prandtl_bkgnd)
 
      kd_col(:) = 0.0 ; kv_col(:) = 0.0  ! Is this line necessary?
      call cvmix_coeffs_bkgnd(mdiff_out=kv_col, tdiff_out=kd_col, nlev=nz, max_nlev=nz)
 
      ! Update Kd and Kv.
      do k=1,nz+1
        cs%Kv_bkgnd(i,j,k) = us%m2_s_to_Z2_T*kv_col(k)
        cs%Kd_bkgnd(i,j,k) = us%m2_s_to_Z2_T*kd_col(k)
      enddo
      do k=1,nz
        kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * us%m2_s_to_Z2_T * (kd_col(k) + kd_col(k+1))
      enddo
    enddo ! i loop
 
  elseif ((.not. cs%Bryan_Lewis_diffusivity) .and. (.not.cs%bulkmixedlayer) .and. &
          (.not. cs%horiz_varying_background) .and. (cs%Kd /= cs%Kdml)) then
    i_hmix = 1.0 / cs%Hmix
    do i=is,ie ; depth(i) = 0.0 ; enddo
    do k=1,nz ; do i=is,ie
      depth_c = depth(i) + 0.5*gv%H_to_Z*h(i,j,k)
      if (depth_c <= cs%Hmix) then ; cs%Kd_bkgnd(i,j,k) = cs%Kdml
      elseif (depth_c >= 2.0*cs%Hmix) then ; cs%Kd_bkgnd(i,j,k) = cs%Kd_sfc(i,j)
      else
        kd_lay(i,j,k) = ((cs%Kd_sfc(i,j) - cs%Kdml) * i_hmix) * depth_c + &
                        (2.0*cs%Kdml - cs%Kd_sfc(i,j))
      endif
 
      depth(i) = depth(i) + gv%H_to_Z*h(i,j,k)
    enddo ; enddo
 
  elseif (cs%horiz_varying_background) then
    !### Note that there are lots of hard-coded parameters (mostly latitudes and longitudes) here.
    do i=is,ie
      bckgrnd_vdc_psis = cs%bckgrnd_vdc_psim * exp(-(0.4*(g%geoLatT(i,j)+28.9))**2)
      bckgrnd_vdc_psin = cs%bckgrnd_vdc_psim * exp(-(0.4*(g%geoLatT(i,j)-28.9))**2)
      !### Add parentheses.
      cs%Kd_bkgnd(i,j,:) = cs%bckgrnd_vdc_eq + bckgrnd_vdc_psin + bckgrnd_vdc_psis
 
      if (g%geoLatT(i,j) < -10.0) then
        cs%Kd_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) + cs%bckgrnd_vdc1
      elseif (g%geoLatT(i,j) <= 10.0) then
        cs%Kd_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) + cs%bckgrnd_vdc1 * (g%geoLatT(i,j)/10.0)**2
      else
        cs%Kd_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) + cs%bckgrnd_vdc1
      endif
 
      ! North Banda Sea
      if ( (g%geoLatT(i,j) < -1.0)  .and. (g%geoLatT(i,j) > -4.0) .and. &
           ( mod(g%geoLonT(i,j)+360.0,360.0) > 103.0) .and. &
           ( mod(g%geoLonT(i,j)+360.0,360.0) < 134.0) ) then
        cs%Kd_bkgnd(i,j,:) = cs%bckgrnd_vdc_Banda
      endif
 
      ! Middle Banda Sea
      if ( (g%geoLatT(i,j) <= -4.0) .and. (g%geoLatT(i,j) > -7.0) .and. &
           ( mod(g%geoLonT(i,j)+360.0,360.0) > 106.0) .and. &
           ( mod(g%geoLonT(i,j)+360.0,360.0) < 140.0) ) then
        cs%Kd_bkgnd(i,j,:) = cs%bckgrnd_vdc_Banda
      endif
 
      ! South Banda Sea
      if ( (g%geoLatT(i,j) <= -7.0) .and. (g%geoLatT(i,j) > -8.3) .and. &
           ( mod(g%geoLonT(i,j)+360.0,360.0) > 111.0) .and. &
           ( mod(g%geoLonT(i,j)+360.0,360.0) < 142.0) ) then
        cs%Kd_bkgnd(i,j,:) = cs%bckgrnd_vdc_Banda
      endif
 
      ! Compute kv_bkgnd
      cs%kv_bkgnd(i,j,:) = cs%Kd_bkgnd(i,j,:) * cs%prandtl_bkgnd
 
      ! Update Kd (uniform profile; no interpolation needed)
      kd_lay(i,j,:) = cs%Kd_bkgnd(i,j,1)
 
    enddo
 
  elseif (cs%Henyey_IGW_background_new) then
    i_x30 = 2.0 / invcosh(cs%N0_2Omega*2.0) ! This is evaluated at 30 deg.
    i_2omega = 0.5 / cs%omega
    do k=1,nz ; do i=is,ie
      abs_sin = max(epsilon, abs(sin(g%geoLatT(i,j)*deg_to_rad)))
      n_2omega = max(abs_sin, sqrt(n2_lay(i,k))*i_2omega)
      n02_n2 = (cs%N0_2Omega/n_2omega)**2
      kd_lay(i,j,k) = max(cs%Kd_min, cs%Kd_sfc(i,j) * &
           ((abs_sin * invcosh(n_2omega/abs_sin)) * i_x30)*n02_n2)
    enddo ; enddo
 
  else
    do k=1,nz ; do i=is,ie
      kd_lay(i,j,k) = cs%Kd_sfc(i,j)
    enddo ; enddo
  endif
 
  ! Update CS%kd_bkgnd and CS%kv_bkgnd for diagnostic purposes
  if (.not. (cs%Bryan_Lewis_diffusivity .or. cs%horiz_varying_background)) then
    do i=is,ie
      cs%kd_bkgnd(i,j,1) = 0.0; cs%kv_bkgnd(i,j,1) = 0.0
      cs%kd_bkgnd(i,j,nz+1) = 0.0; cs%kv_bkgnd(i,j,nz+1) = 0.0
      do k=2,nz
        cs%Kd_bkgnd(i,j,k) = 0.5*(kd_lay(i,j,k-1) + kd_lay(i,j,k))
        cs%Kv_bkgnd(i,j,k) = cs%Kd_bkgnd(i,j,k) * cs%prandtl_bkgnd
      enddo
    enddo
  endif
 
  ! Update Kv
  if (associated(kv)) then
    do k=1,nz+1 ; do i=is,ie
      kv(i,j,k) = kv(i,j,k) + cs%Kv_bkgnd(i,j,k)
    enddo ; enddo
  endif
 
  ! TODO: In both CS%Bryan_Lewis_diffusivity and CS%horiz_varying_background, KV and KD at surface
  ! and bottom interfaces are set to be nonzero. Make sure this is not problematic.
 
end subroutine calculate_bkgnd_mixing
 
!> Reads the parameter "USE_CVMix_BACKGROUND" 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_bkgnd_is_used(param_file)
  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
  call get_param(param_file, mdl, "USE_CVMix_BACKGROUND", cvmix_bkgnd_is_used, &
                 default=.false., do_not_log = .true.)
 
end function cvmix_bkgnd_is_used
 
!> Sets CS%bkgnd_scheme_str to check whether multiple background diffusivity schemes are activated.
!! The string is also for error/log messages.
subroutine check_bkgnd_scheme(CS,str)
  type(bkgnd_mixing_cs), pointer :: CS  !< Control structure
  character(len=*), intent(in)   :: str !< Background scheme identifier deducted from MOM_input
                                        !! parameters
 
  if (trim(cs%bkgnd_scheme_str)=="none") then
    cs%bkgnd_scheme_str = str
  else
     call mom_error(fatal, "set_diffusivity_init: Cannot activate both "//trim(str)//" and "//&
          trim(cs%bkgnd_scheme_str)//".")
  endif
 
end subroutine
 
!> Clear pointers and dealocate memory
subroutine bkgnd_mixing_end(CS)
  type(bkgnd_mixing_cs), pointer :: cs !< Control structure for this module that
                                       !! will be deallocated in this subroutine
 
  if (.not. associated(cs)) return
 
  deallocate(cs%kd_bkgnd)
  deallocate(cs%kv_bkgnd)
  deallocate(cs)
 
end subroutine bkgnd_mixing_end
 
 
end module mom_bkgnd_mixing