mom_tidal_mixing Module Reference

Detailed Description

Interface to vertical tidal mixing schemes including CVMix tidal mixing.

Data Types
type	tidal_mixing_cs
	Control structure with parameters for the tidal mixing module. More...
type	tidal_mixing_diags
	Containers for tidal mixing diagnostics. More...

Functions/Subroutines
logical function, public	tidal_mixing_init (Time, G, GV, US, param_file, diag, CS)
	Initializes internal tidal dissipation scheme for diapycnal mixing. More...
subroutine, public	calculate_tidal_mixing (h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, CS, N2_lay, N2_int, Kd_lay, Kd_int, Kd_max, Kv)
	Depending on whether or not CVMix is active, calls the associated subroutine to compute internal tidal dissipation and to add the effect of internal-tide-driven mixing to the layer or interface diffusivities. More...
subroutine	calculate_cvmix_tidal (h, j, G, GV, US, CS, N2_int, Kd_lay, Kd_int, Kv)
	Calls the CVMix routines to compute tidal dissipation and to add the effect of internal-tide-driven mixing to the interface diffusivities. More...
subroutine	add_int_tide_diffusivity (h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, CS, N2_lay, Kd_lay, Kd_int, Kd_max)
	This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities. The mechanisms considered are (1) local dissipation of internal waves generated by the barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves. Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction, Froude-number-depending breaking, PSI, etc.). More...
subroutine, public	setup_tidal_diagnostics (G, CS)
	Sets up diagnostics arrays for tidal mixing. More...
subroutine, public	post_tidal_diagnostics (G, GV, h, CS)
	This subroutine offers up diagnostics of the tidal mixing. More...
subroutine, public	tidal_mixing_h_amp (h_amp, G, j, CS)
	This subroutine returns a zonal slice of the topographic roughness amplitudes. More...
subroutine	read_tidal_energy (G, US, tidal_energy_type, tidal_energy_file, CS)
	This subroutine read tidal energy inputs from a file. More...
subroutine	read_tidal_constituents (G, US, tidal_energy_file, CS)
	This subroutine reads tidal input energy from a file by constituent. More...
subroutine, public	tidal_mixing_end (CS)
	Clear pointers and deallocate memory. More...

Variables
character *(20), parameter	stlaurent_profile_string = "STLAURENT_02"
	Coded parmameters for specifying mixing schemes.
character *(20), parameter	polzin_profile_string = "POLZIN_09"
	Coded parmameters for specifying mixing schemes.
integer, parameter	stlaurent_02 = 1
	Coded parmameters for specifying mixing schemes.
integer, parameter	polzin_09 = 2
	Coded parmameters for specifying mixing schemes.
character *(20), parameter	simmons_scheme_string = "SIMMONS"
	Coded parmameters for specifying mixing schemes.
character *(20), parameter	schmittner_scheme_string = "SCHMITTNER"
	Coded parmameters for specifying mixing schemes.
integer, parameter	simmons = 1
	Coded parmameters for specifying mixing schemes.
integer, parameter	schmittner = 2
	Coded parmameters for specifying mixing schemes.

Function/Subroutine Documentation

add_int_tide_diffusivity()

subroutine mom_tidal_mixing::add_int_tide_diffusivity ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h, real, dimension( g %isd: g %ied), intent(in)  N2_bot, integer, intent(in)  j, real, dimension( g %isd: g %ied, g %ke), intent(in)  TKE_to_Kd, real, dimension( g %isd: g %ied, g %ke), intent(in)  max_TKE, type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, type(tidal_mixing_cs), pointer  CS, real, dimension( g %isd: g %ied, g %ke), intent(in)  N2_lay, real, dimension( g %isd: g %ied, g %ke), intent(inout), optional  Kd_lay, real, dimension( g %isd: g %ied, g %ke+1), intent(inout), optional  Kd_int, real, intent(in)  Kd_max  )

private

This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities. The mechanisms considered are (1) local dissipation of internal waves generated by the barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves. Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction, Froude-number-depending breaking, PSI, etc.).

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	n2_bot	The near-bottom squared buoyancy frequency frequency [T-2 ~> s-2].
[in]	n2_lay	The squared buoyancy frequency of the layers [T-2 ~> s-2].
[in]	j	The j-index to work on
[in]	tke_to_kd	The conversion rate between the TKE dissipated within a layer and the diapycnal diffusivity within that layer, usually (~Rho_0 / (G_Earth * dRho_lay)) [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
[in]	max_tke	The energy required to for a layer to entrain to its maximum realizable thickness [Z3 T-3 ~> m3 s-3]
	cs	The control structure for this module
[in,out]	kd_lay	The diapycnal diffusivity in layers [Z2 T-1 ~> m2 s-1].
[in,out]	kd_int	The diapycnal diffusivity at interfaces
[in]	kd_max	The maximum increment for diapycnal diffusivity due to TKE-based processes [Z2 T-1 ~> m2 s-1]. Set this to a negative value to have no limit.

Definition at line 962 of file MOM_tidal_mixing.F90.

  type(ocean_grid_type),            intent(in)    :: G      !< The ocean's grid structure
  type(verticalGrid_type),          intent(in)    :: GV     !< The ocean's 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 thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G)),         intent(in)    :: N2_bot !< The near-bottom squared buoyancy frequency
                                                            !! frequency [T-2 ~> s-2].
  real, dimension(SZI_(G),SZK_(G)), intent(in)    :: N2_lay !< The squared buoyancy frequency of the
                                                            !! layers [T-2 ~> s-2].
  integer,                          intent(in)    :: j      !< The j-index to work on
  real, dimension(SZI_(G),SZK_(G)), intent(in)    :: TKE_to_Kd !< The conversion rate between the TKE
                                                            !! dissipated within a layer and the
                                                            !! diapycnal diffusivity within that layer,
                                                            !! usually (~Rho_0 / (G_Earth * dRho_lay))
                                                            !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
  real, dimension(SZI_(G),SZK_(G)), intent(in)    :: max_TKE !< The energy required to for a layer to entrain
                                                            !! to its maximum realizable thickness [Z3 T-3 ~> m3 s-3]
  type(tidal_mixing_cs),            pointer       :: CS     !< The control structure for this module
  real, dimension(SZI_(G),SZK_(G)), &
                          optional, intent(inout) :: Kd_lay !< The diapycnal diffusivity in layers [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZK_(G)+1), &
                          optional, intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces
                                                            !! [Z2 T-1 ~> m2 s-1].
  real,                             intent(in)    :: Kd_max !< The maximum increment for diapycnal
                                                            !! diffusivity due to TKE-based processes
                                                            !! [Z2 T-1 ~> m2 s-1].
                                                            !! Set this to a negative value to have no limit.
 
  ! local
 
  real, dimension(SZI_(G)) :: &
    htot,             & ! total thickness above or below a layer, or the
                        ! integrated thickness in the BBL [Z ~> m].
    htot_wkb,         & ! WKB scaled distance from top to bottom [Z ~> m].
    tke_itidal_bot,   & ! internal tide TKE at ocean bottom [Z3 T-3 ~> m3 s-3]
    tke_niku_bot,     & ! lee-wave TKE at ocean bottom [Z3 T-3 ~> m3 s-3]
    tke_lowmode_bot,  & ! internal tide TKE at ocean bottom lost from all remote low modes [Z3 T-3 ~> m3 s-3] (BDM)
    inv_int,          & ! inverse of TKE decay for int tide over the depth of the ocean [nondim]
    inv_int_lee,      & ! inverse of TKE decay for lee waves over the depth of the ocean [nondim]
    inv_int_low,      & ! inverse of TKE decay for low modes over the depth of the ocean [nondim] (BDM)
    z0_polzin,        & ! TKE decay scale in Polzin formulation [Z ~> m].
    z0_polzin_scaled, & ! TKE decay scale in Polzin formulation [Z ~> m].
                        ! multiplied by N2_bot/N2_meanz to be coherent with the WKB scaled z
                        ! z*=int(N2/N2_bot) * N2_bot/N2_meanz = int(N2/N2_meanz)
                        ! z0_Polzin_scaled = z0_Polzin * N2_bot/N2_meanz
    n2_meanz,         & ! vertically averaged squared buoyancy frequency [T-2] for WKB scaling
    tke_itidal_rem,   & ! remaining internal tide TKE (from barotropic source) [Z3 T-3 ~> m3 s-3]
    tke_niku_rem,     & ! remaining lee-wave TKE [Z3 T-3 ~> m3 s-3]
    tke_lowmode_rem,  & ! remaining internal tide TKE (from propagating low mode source) [Z3 T-3 ~> m3 s-3] (BDM)
    tke_frac_top,     & ! fraction of bottom TKE that should appear at top of a layer [nondim]
    tke_frac_top_lee, & ! fraction of bottom TKE that should appear at top of a layer [nondim]
    tke_frac_top_lowmode, &
                        ! fraction of bottom TKE that should appear at top of a layer [nondim] (BDM)
    z_from_bot,       & ! distance from bottom [Z ~> m].
    z_from_bot_wkb      ! WKB scaled distance from bottom [Z ~> m].
 
  real :: I_rho0        ! Inverse of the Boussinesq reference density, i.e. 1 / RHO0 [R-1 ~> m3 kg-1]
  real :: Kd_add        ! diffusivity to add in a layer [Z2 T-1 ~> m2 s-1].
  real :: TKE_itide_lay ! internal tide TKE imparted to a layer (from barotropic) [Z3 T-3 ~> m3 s-3]
  real :: TKE_Niku_lay  ! lee-wave TKE imparted to a layer [Z3 T-3 ~> m3 s-3]
  real :: TKE_lowmode_lay ! internal tide TKE imparted to a layer (from low mode) [Z3 T-3 ~> m3 s-3] (BDM)
  real :: frac_used     ! fraction of TKE that can be used in a layer [nondim]
  real :: Izeta         ! inverse of TKE decay scale [Z-1 ~> m-1].
  real :: Izeta_lee     ! inverse of TKE decay scale for lee waves [Z-1 ~> m-1].
  real :: z0Ps_num      ! The numerator of the unlimited z0_Polzin_scaled [Z T-3 ~> m s-3].
  real :: z0Ps_denom    ! The denominator of the unlimited z0_Polzin_scaled [T-3 ~> s-3].
  real :: z0_psl        ! temporary variable [Z ~> m].
  real :: TKE_lowmode_tot ! TKE from all low modes [R Z3 T-3 ~> W m-2] (BDM)
 
  logical :: use_Polzin, use_Simmons
  character(len=160) :: mesg  ! The text of an error message
  integer :: i, k, is, ie, nz
  integer :: a, fr, m
  type(tidal_mixing_diags), pointer :: dd => null()
 
  is = g%isc ; ie = g%iec ; nz = g%ke
  dd => cs%dd
 
  if (.not.(cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation)) return
 
  do i=is,ie ; htot(i) = 0.0 ; inv_int(i) = 0.0 ; inv_int_lee(i) = 0.0 ; inv_int_low(i) = 0.0 ;enddo
  do k=1,nz ; do i=is,ie
    htot(i) = htot(i) + gv%H_to_Z*h(i,j,k)
  enddo ; enddo
 
  i_rho0 = 1.0 / (gv%Rho0)
 
  use_polzin = ((cs%Int_tide_dissipation .and. (cs%int_tide_profile == polzin_09)) .or. &
                (cs%lee_wave_dissipation .and. (cs%lee_wave_profile == polzin_09)) .or. &
                (cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == polzin_09)))
  use_simmons = ((cs%Int_tide_dissipation .and. (cs%int_tide_profile == stlaurent_02)) .or. &
                 (cs%lee_wave_dissipation .and. (cs%lee_wave_profile == stlaurent_02)) .or. &
                 (cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == stlaurent_02)))
 
  ! Calculate parameters for vertical structure of dissipation
  ! Simmons:
  if ( use_simmons ) then
    izeta = 1.0 / max(cs%Int_tide_decay_scale, gv%H_subroundoff*gv%H_to_Z)
    izeta_lee = 1.0 / max(cs%Int_tide_decay_scale*cs%Decay_scale_factor_lee, &
                          gv%H_subroundoff*gv%H_to_Z)
    do i=is,ie
      cs%Nb(i,j) = sqrt(n2_bot(i))
      if (associated(dd%N2_bot)) dd%N2_bot(i,j) = n2_bot(i)
      if ( cs%Int_tide_dissipation ) then
        if (izeta*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
          inv_int(i) = 1.0 / (1.0 - exp(-izeta*htot(i)))
        endif
      endif
      if ( cs%Lee_wave_dissipation ) then
        if (izeta_lee*htot(i) > 1.0e-14) then  ! L'Hospital's version of Adcroft's reciprocal rule.
          inv_int_lee(i) = 1.0 / (1.0 - exp(-izeta_lee*htot(i)))
        endif
      endif
      if ( cs%Lowmode_itidal_dissipation) then
        if (izeta*htot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
          inv_int_low(i) = 1.0 / (1.0 - exp(-izeta*htot(i)))
        endif
      endif
      z_from_bot(i) = gv%H_to_Z*h(i,j,nz)
    enddo
  endif ! Simmons
 
  ! Polzin:
  if ( use_polzin ) then
    ! WKB scaling of the vertical coordinate
    do i=is,ie ; n2_meanz(i) = 0.0 ; enddo
    do k=1,nz ; do i=is,ie
      n2_meanz(i) = n2_meanz(i) + n2_lay(i,k) * gv%H_to_Z * h(i,j,k)
    enddo ; enddo
    do i=is,ie
      n2_meanz(i) = n2_meanz(i) / (htot(i) + gv%H_subroundoff*gv%H_to_Z)
      if (associated(dd%N2_meanz))  dd%N2_meanz(i,j) = n2_meanz(i)
    enddo
 
    ! WKB scaled z*(z=H) z* at the surface using the modified Polzin WKB scaling
    do i=is,ie ; htot_wkb(i) = htot(i) ; enddo
!    do i=is,ie ; htot_WKB(i) = 0.0 ; enddo
!    do k=1,nz ; do i=is,ie
!      htot_WKB(i) = htot_WKB(i) + GV%H_to_Z*h(i,j,k) * N2_lay(i,k) / N2_meanz(i)
!    enddo ; enddo
    ! htot_WKB(i) = htot(i) ! Nearly equivalent and simpler
 
    do i=is,ie
      cs%Nb(i,j) = sqrt(n2_bot(i))
      if (cs%answers_2018) then
        if ((cs%tideamp(i,j) > 0.0) .and. &
            (cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j)**3 > 1.0e-14*us%T_to_s**3) ) then
          z0_polzin(i) = cs%Polzin_decay_scale_factor * cs%Nu_Polzin * &
                         cs%Nbotref_Polzin**2 * cs%tideamp(i,j) / &
                       ( cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j)**3 )
          if (z0_polzin(i) < cs%Polzin_min_decay_scale) &
            z0_polzin(i) = cs%Polzin_min_decay_scale
          if (n2_meanz(i) > 1.0e-14*us%T_to_s**2  ) then
            z0_polzin_scaled(i) = z0_polzin(i)*cs%Nb(i,j)**2 / n2_meanz(i)
          else
            z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
          endif
          if (z0_polzin_scaled(i) > (cs%Polzin_decay_scale_max_factor * htot(i)) ) &
            z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
        else
          z0_polzin(i) = cs%Polzin_decay_scale_max_factor * htot(i)
          z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
        endif
      else
        z0ps_num = (cs%Polzin_decay_scale_factor * cs%Nu_Polzin * cs%Nbotref_Polzin**2) * cs%tideamp(i,j)
        z0ps_denom = ( cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j) * n2_meanz(i) )
        if ((cs%tideamp(i,j) > 0.0) .and. &
            (z0ps_num < z0ps_denom * cs%Polzin_decay_scale_max_factor * htot(i))) then
          z0_polzin_scaled(i) = z0ps_num / z0ps_denom
 
          if (abs(n2_meanz(i) * z0_polzin_scaled(i)) < &
              cs%Nb(i,j)**2 * (cs%Polzin_decay_scale_max_factor * htot(i))) then
            z0_polzin(i) = z0_polzin_scaled(i) * (n2_meanz(i) / cs%Nb(i,j)**2)
          else
            z0_polzin(i) = cs%Polzin_decay_scale_max_factor * htot(i)
          endif
        else
          z0_polzin(i) = cs%Polzin_decay_scale_max_factor * htot(i)
          z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
        endif
      endif
 
      if (associated(dd%Polzin_decay_scale)) &
        dd%Polzin_decay_scale(i,j) = z0_polzin(i)
      if (associated(dd%Polzin_decay_scale_scaled)) &
        dd%Polzin_decay_scale_scaled(i,j) = z0_polzin_scaled(i)
      if (associated(dd%N2_bot)) dd%N2_bot(i,j) = cs%Nb(i,j)*cs%Nb(i,j)
 
      if (cs%answers_2018) then
        ! These expressions use dimensional constants to avoid NaN values.
        if ( cs%Int_tide_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
          if (htot_wkb(i) > 1.0e-14*us%m_to_Z) &
            inv_int(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1.0
        endif
        if ( cs%lee_wave_dissipation .and. (cs%lee_wave_profile == polzin_09) ) then
          if (htot_wkb(i) > 1.0e-14*us%m_to_Z) &
            inv_int_lee(i) = ( z0_polzin_scaled(i)*cs%Decay_scale_factor_lee / htot_wkb(i) ) + 1.0
        endif
        if ( cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
          if (htot_wkb(i) > 1.0e-14*us%m_to_Z) &
            inv_int_low(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1.0
        endif
      else
        ! These expressions give values of Inv_int < 10^14 using a variant of Adcroft's reciprocal rule.
        inv_int(i) = 0.0 ; inv_int_lee(i) = 0.0 ; inv_int_low(i) = 0.0
        if ( cs%Int_tide_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
          if (z0_polzin_scaled(i) < 1.0e14 * htot_wkb(i)) &
            inv_int(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1.0
        endif
        if ( cs%lee_wave_dissipation .and. (cs%lee_wave_profile == polzin_09) ) then
          if (z0_polzin_scaled(i) < 1.0e14 * htot_wkb(i)) &
            inv_int_lee(i) = ( z0_polzin_scaled(i)*cs%Decay_scale_factor_lee / htot_wkb(i) ) + 1.0
        endif
        if ( cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
          if (z0_polzin_scaled(i) < 1.0e14 * htot_wkb(i)) &
            inv_int_low(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1.0
        endif
      endif
 
      z_from_bot(i) = gv%H_to_Z*h(i,j,nz)
      ! Use the new formulation for WKB scaling.  N2 is referenced to its vertical mean.
      if (cs%answers_2018) then
        if (n2_meanz(i) > 1.0e-14*us%T_to_s**2 ) then
          z_from_bot_wkb(i) = gv%H_to_Z*h(i,j,nz) * n2_lay(i,nz) / n2_meanz(i)
        else ; z_from_bot_wkb(i) = 0 ; endif
      else
        if (gv%H_to_Z*h(i,j,nz) * n2_lay(i,nz) < n2_meanz(i) * (1.0e14 * htot_wkb(i))) then
          z_from_bot_wkb(i) = gv%H_to_Z*h(i,j,nz) * n2_lay(i,nz) / n2_meanz(i)
        else ; z_from_bot_wkb(i) = 0 ; endif
      endif
    enddo
  endif  ! Polzin
 
  ! Calculate/get dissipation values at bottom
  ! Both Polzin and Simmons:
  do i=is,ie
    ! Dissipation of locally trapped internal tide (non-propagating high modes)
    tke_itidal_bot(i) = min(cs%TKE_itidal(i,j)*cs%Nb(i,j), cs%TKE_itide_max)
    if (associated(dd%TKE_itidal_used)) &
      dd%TKE_itidal_used(i,j) = tke_itidal_bot(i)
    tke_itidal_bot(i) = (i_rho0 * cs%Mu_itides * cs%Gamma_itides) * tke_itidal_bot(i)
    ! Dissipation of locally trapped lee waves
    tke_niku_bot(i) = 0.0
    if (cs%Lee_wave_dissipation) then
      tke_niku_bot(i) = (i_rho0 * cs%Mu_itides * cs%Gamma_lee) * cs%TKE_Niku(i,j)
    endif
    ! Dissipation of propagating internal tide (baroclinic low modes; rays) (BDM)
    tke_lowmode_tot    = 0.0
    tke_lowmode_bot(i) = 0.0
    if (cs%Lowmode_itidal_dissipation) then
      ! get loss rate due to wave drag on low modes (already multiplied by q)
 
      ! TODO: uncomment the following call and fix it
      !call get_lowmode_loss(i,j,G,CS%int_tide_CSp,"WaveDrag",TKE_lowmode_tot)
      write (mesg,*) "========", __file__, __line__
      call mom_error(fatal,trim(mesg)//": this block not supported yet. (aa)")
 
      tke_lowmode_bot(i) = cs%Mu_itides * i_rho0 * tke_lowmode_tot
    endif
    ! Vertical energy flux at bottom
    tke_itidal_rem(i)  = inv_int(i)     * tke_itidal_bot(i)
    tke_niku_rem(i)    = inv_int_lee(i) * tke_niku_bot(i)
    tke_lowmode_rem(i) = inv_int_low(i) * tke_lowmode_bot(i)
 
    if (associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,nz) = tke_itidal_rem(i) !why is this here? BDM
  enddo
 
  ! Estimate the work that would be done by mixing in each layer.
  ! Simmons:
  if ( use_simmons ) then
    do k=nz-1,2,-1 ; do i=is,ie
      if (max_tke(i,k) <= 0.0) cycle
      z_from_bot(i) = z_from_bot(i) + gv%H_to_Z*h(i,j,k)
 
      ! Fraction of bottom flux predicted to reach top of this layer
      tke_frac_top(i)         = inv_int(i)     * exp(-izeta * z_from_bot(i))
      tke_frac_top_lee(i)     = inv_int_lee(i) * exp(-izeta_lee * z_from_bot(i))
      tke_frac_top_lowmode(i) = inv_int_low(i) * exp(-izeta * z_from_bot(i))
 
      ! Actual influx at bottom of layer minus predicted outflux at top of layer to give
      ! predicted power expended
      tke_itide_lay   = tke_itidal_rem(i)  - tke_itidal_bot(i) * tke_frac_top(i)
      tke_niku_lay    = tke_niku_rem(i)    - tke_niku_bot(i)   * tke_frac_top_lee(i)
      tke_lowmode_lay = tke_lowmode_rem(i) - tke_lowmode_bot(i)* tke_frac_top_lowmode(i)
 
      ! Actual power expended may be less than predicted if stratification is weak; adjust
      if (tke_itide_lay + tke_niku_lay + tke_lowmode_lay > (max_tke(i,k))) then
        frac_used = (max_tke(i,k)) / (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
        tke_itide_lay   = frac_used * tke_itide_lay
        tke_niku_lay    = frac_used * tke_niku_lay
        tke_lowmode_lay = frac_used * tke_lowmode_lay
      endif
 
      ! Calculate vertical flux available to bottom of layer above
      tke_itidal_rem(i)  = tke_itidal_rem(i)  - tke_itide_lay
      tke_niku_rem(i)    = tke_niku_rem(i)    - tke_niku_lay
      tke_lowmode_rem(i) = tke_lowmode_rem(i) - tke_lowmode_lay
 
      ! Convert power to diffusivity
      kd_add  = tke_to_kd(i,k) * (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
 
      if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
      if (present(kd_lay)) then
        kd_lay(i,k) = kd_lay(i,k) + kd_add
      endif
 
      if (present(kd_int)) then
        kd_int(i,k)   = kd_int(i,k)   + 0.5 * kd_add
        kd_int(i,k+1) = kd_int(i,k+1) + 0.5 * kd_add
      endif
 
      ! diagnostics
      if (associated(dd%Kd_itidal)) then
        ! If at layers, dd%Kd_itidal is just TKE_to_Kd(i,k) * TKE_itide_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_itide_lay
        if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
        if (k>1)  dd%Kd_itidal(i,j,k)   = dd%Kd_itidal(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_itidal(i,j,k+1) = dd%Kd_itidal(i,j,k+1) + 0.5*kd_add
      endif
      if (associated(dd%Kd_Itidal_work)) &
        dd%Kd_itidal_work(i,j,k) = gv%Rho0 * tke_itide_lay
      if (associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,k) = tke_itidal_rem(i)
 
      if (associated(dd%Kd_Niku)) then
        ! If at layers, dd%Kd_Niku(i,j,K) is just TKE_to_Kd(i,k) * TKE_Niku_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_niku_lay
        if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
        if (k>1) dd%Kd_Niku(i,j,k)    = dd%Kd_Niku(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_Niku(i,j,k+1) = dd%Kd_Niku(i,j,k+1) + 0.5*kd_add
      endif
!     if (associated(dd%Kd_Niku)) dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
      if (associated(dd%Kd_Niku_work)) &
        dd%Kd_Niku_work(i,j,k) = gv%Rho0 * tke_niku_lay
 
      if (associated(dd%Kd_lowmode)) then
        ! If at layers, dd%Kd_lowmode is just TKE_to_Kd(i,k) * TKE_lowmode_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_lowmode_lay
        if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
        if (k>1)  dd%Kd_lowmode(i,j,k)   = dd%Kd_lowmode(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_lowmode(i,j,k+1) = dd%Kd_lowmode(i,j,k+1) + 0.5*kd_add
      endif
      if (associated(dd%Kd_lowmode_work)) &
        dd%Kd_lowmode_work(i,j,k) = gv%Rho0 * tke_lowmode_lay
      if (associated(dd%Fl_lowmode)) dd%Fl_lowmode(i,j,k) = tke_lowmode_rem(i)
 
    enddo ; enddo
  endif ! Simmons
 
  ! Polzin:
  if ( use_polzin ) then
    do k=nz-1,2,-1 ; do i=is,ie
      if (max_tke(i,k) <= 0.0) cycle
      z_from_bot(i) = z_from_bot(i) + gv%H_to_Z*h(i,j,k)
      if (cs%answers_2018) then
        if (n2_meanz(i) > 1.0e-14*us%T_to_s**2 ) then
          z_from_bot_wkb(i) = z_from_bot_wkb(i) &
              + gv%H_to_Z * h(i,j,k) * n2_lay(i,k) / n2_meanz(i)
        else ; z_from_bot_wkb(i) = 0 ; endif
      else
        if (gv%H_to_Z*h(i,j,k) * n2_lay(i,k) < (1.0e14 * htot_wkb(i)) * n2_meanz(i)) then
          z_from_bot_wkb(i) = z_from_bot_wkb(i) + &
             gv%H_to_Z*h(i,j,k) * n2_lay(i,k) / n2_meanz(i)
        endif
      endif
 
      ! Fraction of bottom flux predicted to reach top of this layer
      tke_frac_top(i)     = ( inv_int(i) * z0_polzin_scaled(i) ) / &
                            ( z0_polzin_scaled(i) + z_from_bot_wkb(i) )
      z0_psl = z0_polzin_scaled(i)*cs%Decay_scale_factor_lee
      tke_frac_top_lee(i) = (inv_int_lee(i) * z0_psl) / (z0_psl + z_from_bot_wkb(i))
      tke_frac_top_lowmode(i) = ( inv_int_low(i) * z0_polzin_scaled(i) ) / &
                            ( z0_polzin_scaled(i) + z_from_bot_wkb(i) )
 
      ! Actual influx at bottom of layer minus predicted outflux at top of layer to give
      ! predicted power expended
      tke_itide_lay   = tke_itidal_rem(i)  - tke_itidal_bot(i) *tke_frac_top(i)
      tke_niku_lay    = tke_niku_rem(i)    - tke_niku_bot(i)   * tke_frac_top_lee(i)
      tke_lowmode_lay = tke_lowmode_rem(i) - tke_lowmode_bot(i)*tke_frac_top_lowmode(i)
 
      ! Actual power expended may be less than predicted if stratification is weak; adjust
      if (tke_itide_lay + tke_niku_lay + tke_lowmode_lay > (max_tke(i,k))) then
        frac_used = (max_tke(i,k)) / (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
        tke_itide_lay   = frac_used * tke_itide_lay
        tke_niku_lay    = frac_used * tke_niku_lay
        tke_lowmode_lay = frac_used * tke_lowmode_lay
      endif
 
      ! Calculate vertical flux available to bottom of layer above
      tke_itidal_rem(i)  = tke_itidal_rem(i)  - tke_itide_lay
      tke_niku_rem(i)    = tke_niku_rem(i)    - tke_niku_lay
      tke_lowmode_rem(i) = tke_lowmode_rem(i) - tke_lowmode_lay
 
      ! Convert power to diffusivity
      kd_add  = tke_to_kd(i,k) * (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
 
      if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
      if (present(kd_lay)) then
        kd_lay(i,k) = kd_lay(i,k) + kd_add
      endif
 
      if (present(kd_int)) then
        kd_int(i,k)   = kd_int(i,k)   + 0.5 * kd_add
        kd_int(i,k+1) = kd_int(i,k+1) + 0.5 * kd_add
      endif
 
      ! diagnostics
      if (associated(dd%Kd_itidal)) then
        ! If at layers, this is just dd%Kd_itidal(i,j,K) = TKE_to_Kd(i,k) * TKE_itide_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_itide_lay
        if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
        if (k>1)  dd%Kd_itidal(i,j,k)   = dd%Kd_itidal(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_itidal(i,j,k+1) = dd%Kd_itidal(i,j,k+1) + 0.5*kd_add
      endif
      if (associated(dd%Kd_Itidal_work)) &
        dd%Kd_itidal_work(i,j,k) = gv%Rho0 * tke_itide_lay
      if (associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,k) = tke_itidal_rem(i)
 
      if (associated(dd%Kd_Niku)) then
        ! If at layers, this is just dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_niku_lay
        if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
        if (k>1) dd%Kd_Niku(i,j,k)    = dd%Kd_Niku(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_Niku(i,j,k+1) = dd%Kd_Niku(i,j,k+1) + 0.5*kd_add
      endif
   !  if (associated(dd%Kd_Niku)) dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
      if (associated(dd%Kd_Niku_work)) dd%Kd_Niku_work(i,j,k) = gv%Rho0 * tke_niku_lay
 
      if (associated(dd%Kd_lowmode)) then
        ! If at layers, dd%Kd_lowmode is just TKE_to_Kd(i,k) * TKE_lowmode_lay
        ! The following sets the interface diagnostics.
        kd_add = tke_to_kd(i,k) * tke_lowmode_lay
        if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
        if (k>1)  dd%Kd_lowmode(i,j,k)   = dd%Kd_lowmode(i,j,k)   + 0.5*kd_add
        if (k<nz) dd%Kd_lowmode(i,j,k+1) = dd%Kd_lowmode(i,j,k+1) + 0.5*kd_add
      endif
      if (associated(dd%Kd_lowmode_work)) &
        dd%Kd_lowmode_work(i,j,k) = gv%Rho0 * tke_lowmode_lay
      if (associated(dd%Fl_lowmode)) dd%Fl_lowmode(i,j,k) = tke_lowmode_rem(i)
 
    enddo ; enddo
  endif ! Polzin
 

calculate_cvmix_tidal()

subroutine mom_tidal_mixing::calculate_cvmix_tidal ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h, integer, intent(in)  j, type(ocean_grid_type), intent(in)  G, type(verticalgrid_type), intent(in)  GV, type(unit_scale_type), intent(in)  US, type(tidal_mixing_cs), pointer  CS, real, dimension( g %isd: g %ied, g %ke+1), intent(in)  N2_int, real, dimension( g %isd: g %ied, g %ke), intent(inout), optional  Kd_lay, real, dimension( g %isd: g %ied, g %ke+1), intent(inout), optional  Kd_int, real, dimension(:,:,:), pointer  Kv  )

private

Calls the CVMix routines to compute tidal dissipation and to add the effect of internal-tide-driven mixing to the interface diffusivities.

Parameters

[in]	j	The j-index to work on
[in]	g	Grid structure.
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
	cs	This module's control structure.
[in]	n2_int	The squared buoyancy frequency at the interfaces [T-2 ~> s-2].
[in]	h	Layer thicknesses [H ~> m or kg m-2].
[in,out]	kd_lay	The diapycnal diffusivity in the layers [Z2 T-1 ~> m2 s-1].
[in,out]	kd_int	The diapycnal diffusivity at interfaces [Z2 T-1 ~> m2 s-1].
	kv	The "slow" vertical viscosity at each interface (not layer!) [Z2 T-1 ~> m2 s-1].

Definition at line 718 of file MOM_tidal_mixing.F90.

  integer,                 intent(in)    :: j     !< The j-index to work on
  type(ocean_grid_type),   intent(in)    :: G     !< Grid structure.
  type(verticalGrid_type), intent(in)    :: GV    !< ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US    !< A dimensional unit scaling type
  type(tidal_mixing_cs),   pointer       :: CS    !< This module's control structure.
  real, dimension(SZI_(G),SZK_(G)+1), intent(in) :: N2_int !< The squared buoyancy
                                                  !! frequency at the interfaces [T-2 ~> s-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: h     !< Layer thicknesses [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)), &
                 optional, intent(inout) :: Kd_lay!< The diapycnal diffusivity in the layers [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZK_(G)+1), &
                 optional, intent(inout) :: Kd_int!< The diapycnal diffusivity at interfaces [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].
  ! Local variables
  real, dimension(SZK_(G)+1) :: Kd_tidal    ! tidal diffusivity [m2 s-1]
  real, dimension(SZK_(G)+1) :: Kv_tidal    ! tidal viscosity [m2 s-1]
  real, dimension(SZK_(G)+1) :: vert_dep    ! vertical deposition
  real, dimension(SZK_(G)+1) :: iFaceHeight ! Height of interfaces [m]
  real, dimension(SZK_(G)+1) :: SchmittnerSocn
  real, dimension(SZK_(G))   :: cellHeight  ! Height of cell centers [m]
  real, dimension(SZK_(G))   :: tidal_qe_md ! Tidal dissipation energy interpolated from 3d input
                                            ! to model coordinates
  real, dimension(SZK_(G)+1) :: N2_int_i    ! De-scaled interface buoyancy frequency [s-2]
  real, dimension(SZK_(G))   :: Schmittner_coeff
  real, dimension(SZK_(G))   :: h_m         ! Cell thickness [m]
  real, allocatable, dimension(:,:) :: exp_hab_zetar
 
  integer :: i, k, is, ie
  real :: dh, hcorr, Simmons_coeff
  real, parameter :: rho_fw = 1000.0 ! fresh water density [kg/m^3]
                                     ! TODO: when coupled, get this from CESM (SHR_CONST_RHOFW)
  type(tidal_mixing_diags), pointer :: dd => null()
 
  is  = g%isc ; ie  = g%iec
  dd => cs%dd
 
  select case (cs%CVMix_tidal_scheme)
  case (simmons)
    do i=is,ie
 
      if (g%mask2dT(i,j)<1) cycle
 
      ifaceheight = 0.0 ! BBL is all relative to the surface
      hcorr = 0.0
      do k=1,g%ke
        ! cell center and cell bottom in meters (negative values in the ocean)
        dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment, rescaled to m for use by CVMix.
        dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
        hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
        dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
        cellheight(k)    = ifaceheight(k) - 0.5 * dh
        ifaceheight(k+1) = ifaceheight(k) - dh
      enddo
 
      call cvmix_compute_simmons_invariant( nlev                    = g%ke,                &
                                            energy_flux             = cs%tidal_qe_2d(i,j), &
                                            rho                     = rho_fw,              &
                                            simmonscoeff            = simmons_coeff,       &
                                            vertdep                 = vert_dep,            &
                                            zw                      = ifaceheight,         &
                                            zt                      = cellheight,          &
                                            cvmix_tidal_params_user = cs%CVMix_tidal_params)
 
      ! Since we pass tidal_qe_2d=(CS%Gamma_itides)*tidal_energy_flux_2d, and not tidal_energy_flux_2d in
      ! above subroutine call, we divide Simmons_coeff by CS%Gamma_itides as a corrective step:
      ! TODO: (CS%Gamma_itides)*tidal_energy_flux_2d is unnecessary, directly use tidal_energy_flux_2d
      simmons_coeff = simmons_coeff / cs%Gamma_itides
 
 
      ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable
      do k=1,g%ke+1
        n2_int_i(k) = us%s_to_T**2 * n2_int(i,k)
      enddo
 
      call cvmix_coeffs_tidal( mdiff_out               = kv_tidal,            &
                               tdiff_out               = kd_tidal,            &
                               nsqr                    = n2_int_i,            &
                               oceandepth              = -ifaceheight(g%ke+1),&
                               simmonscoeff            = simmons_coeff,       &
                               vert_dep                = vert_dep,            &
                               nlev                    = g%ke,                &
                               max_nlev                = g%ke,                &
                               cvmix_params            = cs%CVMix_glb_params, &
                               cvmix_tidal_params_user = cs%CVMix_tidal_params)
 
      ! Update diffusivity
      if (present(kd_lay)) then
        do k=1,g%ke
          kd_lay(i,k) = kd_lay(i,k) + 0.5 * us%m2_s_to_Z2_T * (kd_tidal(k) + kd_tidal(k+1))
        enddo
      endif
      if (present(kd_int)) then
        do k=1,g%ke+1
          kd_int(i,k) = kd_int(i,k) +  (us%m2_s_to_Z2_T * kd_tidal(k))
        enddo
      endif
      ! Update viscosity with the proper unit conversion.
      if (associated(kv)) then
        do k=1,g%ke+1
          kv(i,j,k) = kv(i,j,k) + us%m2_s_to_Z2_T * kv_tidal(k)  ! Rescale from m2 s-1 to Z2 T-1.
        enddo
      endif
 
      ! diagnostics
      if (associated(dd%Kd_itidal)) then
        dd%Kd_itidal(i,j,:) = us%m2_s_to_Z2_T*kd_tidal(:)
      endif
      if (associated(dd%N2_int)) then
        dd%N2_int(i,j,:) = n2_int(i,:)
      endif
      if (associated(dd%Simmons_coeff_2d)) then
        dd%Simmons_coeff_2d(i,j) = simmons_coeff
      endif
      if (associated(dd%vert_dep_3d)) then
        dd%vert_dep_3d(i,j,:) = vert_dep(:)
      endif
 
    enddo ! i=is,ie
 
  case (schmittner)
 
    ! TODO: correct exp_hab_zetar shapes in CVMix_compute_Schmittner_invariant
    ! and CVMix_compute_SchmittnerCoeff low subroutines
 
    allocate(exp_hab_zetar(g%ke+1,g%ke+1))
 
    do i=is,ie
 
      if (g%mask2dT(i,j)<1) cycle
 
      ifaceheight = 0.0 ! BBL is all relative to the surface
      hcorr = 0.0
      do k=1,g%ke
        h_m(k) = h(i,j,k)*gv%H_to_m  ! Rescale thicknesses to m for use by CVmix.
        ! cell center and cell bottom in meters (negative values in the ocean)
        dh = h_m(k) + hcorr ! Nominal thickness less the accumulated error (could temporarily make dh<0)
        hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
        dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
        cellheight(k)    = ifaceheight(k) - 0.5 * dh
        ifaceheight(k+1) = ifaceheight(k) - dh
      enddo
 
      schmittnersocn = 0.0 ! TODO: compute this
 
      ! form the time-invariant part of Schmittner coefficient term
      call cvmix_compute_schmittner_invariant(nlev                    = g%ke,           &
                                              vertdep                 = vert_dep,       &
                                              efficiency              = cs%Mu_itides,   &
                                              rho                     = rho_fw,         &
                                              exp_hab_zetar           = exp_hab_zetar,  &
                                              zw                      = ifaceheight,    &
                                              cvmix_tidal_params_user = cs%CVMix_tidal_params)
                  !TODO: in above call, there is no need to pass efficiency, since it gets
                  ! passed via CVMix_init_tidal and stored in CVMix_tidal_params. Change
                  ! CVMix API to prevent this redundancy.
 
      ! remap from input z coordinate to model coordinate:
      tidal_qe_md = 0.0
      call remapping_core_h(cs%remap_cs, size(cs%h_src), cs%h_src, cs%tidal_qe_3d_in(i,j,:), &
                            g%ke, h_m, tidal_qe_md)
 
      ! form the Schmittner coefficient that is based on 3D q*E, which is formed from
      ! summing q_i*TidalConstituent_i over the number of constituents.
      call cvmix_compute_schmittnercoeff( nlev                    = g%ke,               &
                                          energy_flux             = tidal_qe_md(:),     &
                                          rho                     = rho_fw,             &
                                          schmittnercoeff         = schmittner_coeff,   &
                                          exp_hab_zetar           = exp_hab_zetar,      &
                                          cvmix_tidal_params_user = cs%CVMix_tidal_params)
 
      ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable
      do k=1,g%ke+1
        n2_int_i(k) = us%s_to_T**2 * n2_int(i,k)
      enddo
 
      call cvmix_coeffs_tidal_schmittner( mdiff_out               = kv_tidal,             &
                                          tdiff_out               = kd_tidal,             &
                                          nsqr                    = n2_int_i,             &
                                          oceandepth              = -ifaceheight(g%ke+1), &
                                          vert_dep                = vert_dep,             &
                                          nlev                    = g%ke,                 &
                                          max_nlev                = g%ke,                 &
                                          schmittnercoeff         = schmittner_coeff,     &
                                          schmittnersouthernocean = schmittnersocn,       &
                                          cvmix_params            = cs%CVMix_glb_params,  &
                                          cvmix_tidal_params_user = cs%CVMix_tidal_params)
 
      ! Update diffusivity
      if (present(kd_lay)) then
        do k=1,g%ke
          kd_lay(i,k) = kd_lay(i,k) + 0.5 * us%m2_s_to_Z2_T * (kd_tidal(k) + kd_tidal(k+1))
        enddo
      endif
      if (present(kd_int)) then
        do k=1,g%ke+1
          kd_int(i,k) = kd_int(i,k) +  (us%m2_s_to_Z2_T * kd_tidal(k))
        enddo
      endif
 
      ! Update viscosity
      if (associated(kv)) then
        do k=1,g%ke+1
          kv(i,j,k) = kv(i,j,k) + us%m2_s_to_Z2_T * kv_tidal(k)   ! Rescale from m2 s-1 to Z2 T-1.
        enddo
      endif
 
      ! diagnostics
      if (associated(dd%Kd_itidal)) then
        dd%Kd_itidal(i,j,:) = us%m2_s_to_Z2_T*kd_tidal(:)
      endif
      if (associated(dd%N2_int)) then
        dd%N2_int(i,j,:) = n2_int(i,:)
      endif
      if (associated(dd%Schmittner_coeff_3d)) then
        dd%Schmittner_coeff_3d(i,j,:) = schmittner_coeff(:)
      endif
      if (associated(dd%tidal_qe_md)) then
        dd%tidal_qe_md(i,j,:) = tidal_qe_md(:)
      endif
      if (associated(dd%vert_dep_3d)) then
        dd%vert_dep_3d(i,j,:) = vert_dep(:)
      endif
    enddo ! i=is,ie
 
    deallocate(exp_hab_zetar)
 
  case default
     call mom_error(fatal, "tidal_mixing_init: Unrecognized setting "// &
         "#define CVMIX_TIDAL_SCHEME found in input file.")
  end select
 

calculate_tidal_mixing()

subroutine, public mom_tidal_mixing::calculate_tidal_mixing	(	real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		real, dimension(szi_(g)), intent(in)	N2_bot,
		integer, intent(in)	j,
		real, dimension(szi_(g),szk_(g)), intent(in)	TKE_to_Kd,
		real, dimension(szi_(g),szk_(g)), intent(in)	max_TKE,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(tidal_mixing_cs), pointer	CS,
		real, dimension(szi_(g),szk_(g)), intent(in)	N2_lay,
		real, dimension(szi_(g),szk_(g)+1), intent(in)	N2_int,
		real, dimension(szi_(g),szk_(g)), intent(inout), optional	Kd_lay,
		real, dimension(szi_(g),szk_(g)+1), intent(inout), optional	Kd_int,
		real, intent(in)	Kd_max,
		real, dimension(:,:,:), pointer	Kv
	)

Depending on whether or not CVMix is active, calls the associated subroutine to compute internal tidal dissipation and to add the effect of internal-tide-driven mixing to the layer or interface diffusivities.

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	n2_bot	The near-bottom squared buoyancy frequency [T-2 ~> s-2].
[in]	n2_lay	The squared buoyancy frequency of the layers [T-2 ~> s-2].
[in]	n2_int	The squared buoyancy frequency at the interfaces [T-2 ~> s-2].
[in]	j	The j-index to work on
[in]	tke_to_kd	The conversion rate between the TKE dissipated within a layer and the diapycnal diffusivity within that layer, usually (~Rho_0 / (G_Earth * dRho_lay)) [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
[in]	max_tke	The energy required to for a layer to entrain to its maximum realizable thickness [Z3 T-3 ~> m3 s-3]
	cs	The control structure for this module
[in,out]	kd_lay	The diapycnal diffusivity in layers [Z2 T-1 ~> m2 s-1].
[in,out]	kd_int	The diapycnal diffusivity at interfaces,
[in]	kd_max	The maximum increment for diapycnal diffusivity due to TKE-based processes, [Z2 T-1 ~> m2 s-1]. Set this to a negative value to have no limit.
	kv	The "slow" vertical viscosity at each interface (not layer!) [Z2 T-1 ~> m2 s-1].

Definition at line 672 of file MOM_tidal_mixing.F90.

  type(ocean_grid_type),            intent(in)    :: G      !< The ocean's grid structure
  type(verticalGrid_type),          intent(in)    :: GV     !< The ocean's 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 thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G)),         intent(in)    :: N2_bot !< The near-bottom squared buoyancy
                                                            !! frequency [T-2 ~> s-2].
  real, dimension(SZI_(G),SZK_(G)), intent(in)    :: N2_lay !< The squared buoyancy frequency of the
                                                            !! layers [T-2 ~> s-2].
  real, dimension(SZI_(G),SZK_(G)+1), intent(in)  :: N2_int !< The squared buoyancy frequency at the
                                                            !! interfaces [T-2 ~> s-2].
  integer,                          intent(in)    :: j      !< The j-index to work on
  real, dimension(SZI_(G),SZK_(G)), intent(in)    :: TKE_to_Kd !< The conversion rate between the TKE
                                                            !! dissipated within a layer and the
                                                            !! diapycnal diffusivity within that layer,
                                                            !! usually (~Rho_0 / (G_Earth * dRho_lay))
                                                            !! [Z2 T-1 / Z3 T-3 = T2 Z-1 ~> s2 m-1]
  real, dimension(SZI_(G),SZK_(G)), intent(in)    :: max_TKE !< The energy required to for a layer to entrain
                                                            !! to its maximum realizable thickness [Z3 T-3 ~> m3 s-3]
  type(tidal_mixing_cs),            pointer       :: CS     !< The control structure for this module
  real, dimension(SZI_(G),SZK_(G)), &
                          optional, intent(inout) :: Kd_lay !< The diapycnal diffusivity in layers [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZK_(G)+1), &
                          optional, intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces,
                                                            !! [Z2 T-1 ~> m2 s-1].
  real,                             intent(in)    :: Kd_max !< The maximum increment for diapycnal
                                                            !! diffusivity due to TKE-based processes,
                                                            !! [Z2 T-1 ~> m2 s-1].
                                                            !! Set this to a negative value to have no limit.
  real, dimension(:,:,:),           pointer       :: Kv     !< The "slow" vertical viscosity at each interface
                                                            !! (not layer!) [Z2 T-1 ~> m2 s-1].
 
  if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. cs%Lowmode_itidal_dissipation) then
    if (cs%use_CVMix_tidal) then
      call calculate_cvmix_tidal(h, j, g, gv, us, cs, n2_int, kd_lay, kd_int, kv)
    else
      call add_int_tide_diffusivity(h, n2_bot, j, tke_to_kd, max_tke, &
                                    g, gv, us, cs, n2_lay, kd_lay, kd_int, kd_max)
    endif
  endif

post_tidal_diagnostics()

subroutine, public mom_tidal_mixing::post_tidal_diagnostics	(	type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		type(tidal_mixing_cs), pointer	CS
	)

This subroutine offers up diagnostics of the tidal mixing.

Parameters

[in]	g	The ocean's grid structure
[in]	gv	The ocean's vertical grid structure.
[in]	h	Layer thicknesses [H ~> m or kg m-2].
	cs	The control structure for this module

Definition at line 1498 of file MOM_tidal_mixing.F90.

  type(ocean_grid_type),    intent(in)   :: G   !< The ocean's grid structure
  type(verticalGrid_type),  intent(in)   :: GV  !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                            intent(in)   :: h   !< Layer thicknesses [H ~> m or kg m-2].
  type(tidal_mixing_cs),    pointer      :: CS  !< The control structure for this module
 
  ! local
  type(tidal_mixing_diags), pointer :: dd => null()
 
  dd => cs%dd
 
  if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. cs%Lowmode_itidal_dissipation) then
    if (cs%id_TKE_itidal  > 0) call post_data(cs%id_TKE_itidal,  dd%TKE_itidal_used, cs%diag)
    if (cs%id_TKE_leewave > 0) call post_data(cs%id_TKE_leewave, cs%TKE_Niku,        cs%diag)
    if (cs%id_Nb          > 0) call post_data(cs%id_Nb,      cs%Nb,      cs%diag)
    if (cs%id_N2_bot      > 0) call post_data(cs%id_N2_bot,  dd%N2_bot,  cs%diag)
    if (cs%id_N2_meanz    > 0) call post_data(cs%id_N2_meanz,dd%N2_meanz,cs%diag)
 
    if (cs%id_Fl_itidal > 0) call post_data(cs%id_Fl_itidal, dd%Fl_itidal, cs%diag)
    if (cs%id_Kd_itidal > 0) call post_data(cs%id_Kd_itidal, dd%Kd_itidal, cs%diag)
    if (cs%id_Kd_Niku   > 0) call post_data(cs%id_Kd_Niku,   dd%Kd_Niku,   cs%diag)
    if (cs%id_Kd_lowmode> 0) call post_data(cs%id_Kd_lowmode, dd%Kd_lowmode, cs%diag)
    if (cs%id_Fl_lowmode> 0) call post_data(cs%id_Fl_lowmode, dd%Fl_lowmode, cs%diag)
 
    if (cs%id_N2_int> 0)        call post_data(cs%id_N2_int, dd%N2_int, cs%diag)
    if (cs%id_vert_dep> 0)      call post_data(cs%id_vert_dep, dd%vert_dep_3d, cs%diag)
    if (cs%id_Simmons_coeff> 0) call post_data(cs%id_Simmons_coeff, dd%Simmons_coeff_2d, cs%diag)
    if (cs%id_Schmittner_coeff> 0) call post_data(cs%id_Schmittner_coeff, dd%Schmittner_coeff_3d, cs%diag)
    if (cs%id_tidal_qe_md> 0) call post_data(cs%id_tidal_qe_md, dd%tidal_qe_md, cs%diag)
 
    if (cs%id_Kd_Itidal_Work > 0) &
      call post_data(cs%id_Kd_Itidal_Work, dd%Kd_Itidal_Work, cs%diag)
    if (cs%id_Kd_Niku_Work > 0) call post_data(cs%id_Kd_Niku_Work, dd%Kd_Niku_Work, cs%diag)
    if (cs%id_Kd_Lowmode_Work > 0) &
      call post_data(cs%id_Kd_Lowmode_Work, dd%Kd_Lowmode_Work, cs%diag)
 
    if (cs%id_Polzin_decay_scale > 0 ) &
      call post_data(cs%id_Polzin_decay_scale, dd%Polzin_decay_scale, cs%diag)
    if (cs%id_Polzin_decay_scale_scaled > 0 ) &
      call post_data(cs%id_Polzin_decay_scale_scaled, dd%Polzin_decay_scale_scaled, cs%diag)
  endif
 
  if (associated(dd%Kd_itidal)) deallocate(dd%Kd_itidal)
  if (associated(dd%Kd_lowmode)) deallocate(dd%Kd_lowmode)
  if (associated(dd%Fl_itidal)) deallocate(dd%Fl_itidal)
  if (associated(dd%Fl_lowmode)) deallocate(dd%Fl_lowmode)
  if (associated(dd%Polzin_decay_scale)) deallocate(dd%Polzin_decay_scale)
  if (associated(dd%Polzin_decay_scale_scaled)) deallocate(dd%Polzin_decay_scale_scaled)
  if (associated(dd%N2_bot)) deallocate(dd%N2_bot)
  if (associated(dd%N2_meanz)) deallocate(dd%N2_meanz)
  if (associated(dd%Kd_Niku)) deallocate(dd%Kd_Niku)
  if (associated(dd%Kd_Niku_work)) deallocate(dd%Kd_Niku_work)
  if (associated(dd%Kd_Itidal_Work))  deallocate(dd%Kd_Itidal_Work)
  if (associated(dd%Kd_Lowmode_Work)) deallocate(dd%Kd_Lowmode_Work)
  if (associated(dd%TKE_itidal_used)) deallocate(dd%TKE_itidal_used)
  if (associated(dd%N2_int)) deallocate(dd%N2_int)
  if (associated(dd%vert_dep_3d)) deallocate(dd%vert_dep_3d)
  if (associated(dd%Simmons_coeff_2d)) deallocate(dd%Simmons_coeff_2d)
  if (associated(dd%Schmittner_coeff_3d)) deallocate(dd%Schmittner_coeff_3d)
  if (associated(dd%tidal_qe_md)) deallocate(dd%tidal_qe_md)
 

read_tidal_constituents()

subroutine mom_tidal_mixing::read_tidal_constituents ( type(ocean_grid_type), intent(in)  G, type(unit_scale_type), intent(in)  US, character(len=200), intent(in)  tidal_energy_file, type(tidal_mixing_cs), pointer  CS  )

private

This subroutine reads tidal input energy from a file by constituent.

Parameters

[in]	g	The ocean's grid structure
[in]	us	A dimensional unit scaling type
[in]	tidal_energy_file	The file from which to read tidal energy inputs
	cs	The control structure for this module

Definition at line 1612 of file MOM_tidal_mixing.F90.

  type(ocean_grid_type), intent(in) :: G    !< The ocean's grid structure
  type(unit_scale_type), intent(in) :: US   !< A dimensional unit scaling type
  character(len=200),    intent(in) :: tidal_energy_file !< The file from which to read tidal energy inputs
  type(tidal_mixing_cs), pointer    :: CS   !< The control structure for this module
 
  ! local variables
  real, parameter :: C1_3 = 1.0/3.0
  real, dimension(SZI_(G),SZJ_(G)) :: &
    tidal_qk1, &  ! qk1 coefficient used in Schmittner & Egbert
    tidal_qo1     ! qo1 coefficient used in Schmittner & Egbert
  real, allocatable, dimension(:) :: &
    z_t, &        ! depth from surface to midpoint of input layer [Z]
    z_w           ! depth from surface to top of input layer [Z]
  real, allocatable, dimension(:,:,:) :: &
    tc_m2, &      ! input lunar semidiurnal tidal energy flux [W/m^2]
    tc_s2, &      ! input solar semidiurnal tidal energy flux [W/m^2]
    tc_k1, &      ! input lunar diurnal tidal energy flux [W/m^2]
    tc_o1         ! input lunar diurnal tidal energy flux [W/m^2]
  integer, dimension(4) :: nz_in
  integer               :: k, is, ie, js, je, isd, ied, jsd, jed, i, j
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
 
  ! get number of input levels:
  call field_size(tidal_energy_file, 'z_t', nz_in)
 
  ! allocate local variables
  allocate(z_t(nz_in(1)), z_w(nz_in(1)) )
  allocate(tc_m2(isd:ied,jsd:jed,nz_in(1)), &
           tc_s2(isd:ied,jsd:jed,nz_in(1)), &
           tc_k1(isd:ied,jsd:jed,nz_in(1)), &
           tc_o1(isd:ied,jsd:jed,nz_in(1)) )
 
  ! allocate CS variables associated with 3d tidal energy dissipation
  if (.not. allocated(cs%tidal_qe_3d_in)) allocate(cs%tidal_qe_3d_in(isd:ied,jsd:jed,nz_in(1)))
  if (.not. allocated(cs%h_src))          allocate(cs%h_src(nz_in(1)))
 
  ! read in tidal constituents
  call mom_read_data(tidal_energy_file, 'M2', tc_m2, g%domain)
  call mom_read_data(tidal_energy_file, 'S2', tc_s2, g%domain)
  call mom_read_data(tidal_energy_file, 'K1', tc_k1, g%domain)
  call mom_read_data(tidal_energy_file, 'O1', tc_o1, g%domain)
  ! Note the hard-coded assumption that z_t and z_w in the file are in centimeters.
  call mom_read_data(tidal_energy_file, 'z_t', z_t, scale=100.0*us%m_to_Z)
  call mom_read_data(tidal_energy_file, 'z_w', z_w, scale=100.0*us%m_to_Z)
 
  do j=js,je ; do i=is,ie
    if (abs(g%geoLatT(i,j)) < 30.0) then
      tidal_qk1(i,j) = c1_3
      tidal_qo1(i,j) = c1_3
    else
      tidal_qk1(i,j) = 1.0
      tidal_qo1(i,j) = 1.0
    endif
  enddo ; enddo
 
  cs%tidal_qe_3d_in(:,:,:) = 0.0
  do k=1,nz_in(1)
    ! Store the input cell thickness in m for use with CVmix.
    cs%h_src(k) = us%Z_to_m*(z_t(k)-z_w(k))*2.0
    ! form tidal_qe_3d_in from weighted tidal constituents
    do j=js,je ; do i=is,ie
      if ((z_t(k) <= g%bathyT(i,j)) .and. (z_w(k) > cs%tidal_diss_lim_tc)) &
        cs%tidal_qe_3d_in(i,j,k) = c1_3*tc_m2(i,j,k) + c1_3*tc_s2(i,j,k) + &
                tidal_qk1(i,j)*tc_k1(i,j,k) + tidal_qo1(i,j)*tc_o1(i,j,k)
    enddo ; enddo
  enddo
 
  !open(unit=1905,file="out_1905.txt",access="APPEND")
  !do j=G%jsd,G%jed
  !  do i=isd,ied
  !    if ( i+G%idg_offset .eq. 90 .and. j+G%jdg_offset .eq. 126) then
  !      write(1905,*) "-------------------------------------------"
  !      do k=50,nz_in(1)
  !          write(1905,*) i,j,k
  !          write(1905,*) CS%tidal_qe_3d_in(i,j,k), tc_m2(i,j,k)
  !          write(1905,*) z_t(k), G%bathyT(i,j), z_w(k),CS%tidal_diss_lim_tc
  !      end do
  !    endif
  !  enddo
  !enddo
  !close(1905)
 
  ! test if qE is positive
  if (any(cs%tidal_qe_3d_in<0.0)) then
    call mom_error(fatal, "read_tidal_constituents: Negative tidal_qe_3d_in terms.")
  endif
 
  !! collapse 3D q*E to 2D q*E
  !CS%tidal_qe_2d(:,:) = 0.0
  !do k=1,nz_in(1) ; do j=js,je ; do i=is,ie
  !  if (z_t(k) <= G%bathyT(i,j)) &
  !    CS%tidal_qe_2d(i,j) = CS%tidal_qe_2d(i,j) + CS%tidal_qe_3d_in(i,j,k)
  !enddo ; enddo ; enddo
 
  ! initialize input remapping:
  call initialize_remapping(cs%remap_cs, remapping_scheme="PLM", &
                            boundary_extrapolation=.false., check_remapping=cs%debug, &
                            answers_2018=cs%remap_answers_2018)
 
  deallocate(tc_m2)
  deallocate(tc_s2)
  deallocate(tc_k1)
  deallocate(tc_o1)
  deallocate(z_t)
  deallocate(z_w)
 

read_tidal_energy()

subroutine mom_tidal_mixing::read_tidal_energy ( type(ocean_grid_type), intent(in)  G, type(unit_scale_type), intent(in)  US, character(len=20), intent(in)  tidal_energy_type, character(len=200), intent(in)  tidal_energy_file, type(tidal_mixing_cs), pointer  CS  )

private

This subroutine read tidal energy inputs from a file.

Parameters

[in]	g	The ocean's grid structure
[in]	us	A dimensional unit scaling type
[in]	tidal_energy_type	The type of tidal energy inputs to read
[in]	tidal_energy_file	The file from which to read tidalinputs
	cs	The control structure for this module

Definition at line 1582 of file MOM_tidal_mixing.F90.

  type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure
  type(unit_scale_type),   intent(in) :: US   !< A dimensional unit scaling type
  character(len=20),       intent(in) :: tidal_energy_type !< The type of tidal energy inputs to read
  character(len=200),      intent(in) :: tidal_energy_file !< The file from which to read tidalinputs
  type(tidal_mixing_cs),   pointer    :: CS   !< The control structure for this module
  ! local
  integer :: i, j, isd, ied, jsd, jed, nz
  real, allocatable, dimension(:,:) :: tidal_energy_flux_2d ! input tidal energy flux at T-grid points [W m-2]
 
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
 
  select case (uppercase(tidal_energy_type(1:4)))
  case ('JAYN') ! Jayne 2009
    if (.not. allocated(cs%tidal_qe_2d)) allocate(cs%tidal_qe_2d(isd:ied,jsd:jed))
    allocate(tidal_energy_flux_2d(isd:ied,jsd:jed))
    call mom_read_data(tidal_energy_file,'wave_dissipation',tidal_energy_flux_2d, g%domain)
    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      cs%tidal_qe_2d(i,j) = cs%Gamma_itides * tidal_energy_flux_2d(i,j)
    enddo ; enddo
    deallocate(tidal_energy_flux_2d)
  case ('ER03') ! Egbert & Ray 2003
    call read_tidal_constituents(g, us, tidal_energy_file, cs)
  case default
    call mom_error(fatal, "read_tidal_energy: Unknown tidal energy file type.")
  end select
 

setup_tidal_diagnostics()

subroutine, public mom_tidal_mixing::setup_tidal_diagnostics	(	type(ocean_grid_type), intent(in)	G,
		type(tidal_mixing_cs), pointer	CS
	)

Sets up diagnostics arrays for tidal mixing.

Parameters

[in]	g	The ocean's grid structure
	cs	The control structure for this module

Definition at line 1413 of file MOM_tidal_mixing.F90.

  type(ocean_grid_type), intent(in) :: G  !< The ocean's grid structure
  type(tidal_mixing_cs), pointer    :: CS !< The control structure for this module
 
  ! local
  integer :: isd, ied, jsd, jed, nz
  type(tidal_mixing_diags), pointer :: dd => null()
 
  isd = g%isd; ied = g%ied; jsd = g%jsd; jed = g%jed; nz = g%ke
  dd => cs%dd
 
  if ((cs%id_Kd_itidal > 0) .or. (cs%id_Kd_Itidal_work > 0)) then
    allocate(dd%Kd_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Kd_itidal(:,:,:) = 0.0
  endif
  if ((cs%id_Kd_lowmode > 0) .or. (cs%id_Kd_lowmode_work > 0)) then
    allocate(dd%Kd_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Kd_lowmode(:,:,:) = 0.0
  endif
  if ( (cs%id_Fl_itidal > 0) ) then
    allocate(dd%Fl_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Fl_itidal(:,:,:) = 0.0
  endif
  if ( (cs%id_Fl_lowmode > 0) ) then
    allocate(dd%Fl_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Fl_lowmode(:,:,:) = 0.0
  endif
  if ( (cs%id_Polzin_decay_scale > 0) ) then
    allocate(dd%Polzin_decay_scale(isd:ied,jsd:jed))
    dd%Polzin_decay_scale(:,:) = 0.0
  endif
  if ( (cs%id_N2_bot > 0) ) then
    allocate(dd%N2_bot(isd:ied,jsd:jed)) ; dd%N2_bot(:,:) = 0.0
  endif
  if ( (cs%id_N2_meanz > 0) ) then
    allocate(dd%N2_meanz(isd:ied,jsd:jed)) ; dd%N2_meanz(:,:) = 0.0
  endif
  if ( (cs%id_Polzin_decay_scale_scaled > 0) ) then
    allocate(dd%Polzin_decay_scale_scaled(isd:ied,jsd:jed))
    dd%Polzin_decay_scale_scaled(:,:) = 0.0
  endif
  if ((cs%id_Kd_Niku > 0) .or. (cs%id_Kd_Niku_work > 0)) then
    allocate(dd%Kd_Niku(isd:ied,jsd:jed,nz+1)) ; dd%Kd_Niku(:,:,:) = 0.0
  endif
  if (cs%id_Kd_Niku_work > 0) then
    allocate(dd%Kd_Niku_work(isd:ied,jsd:jed,nz)) ; dd%Kd_Niku_work(:,:,:) = 0.0
  endif
  if (cs%id_Kd_Itidal_work > 0) then
    allocate(dd%Kd_Itidal_work(isd:ied,jsd:jed,nz))
    dd%Kd_Itidal_work(:,:,:) = 0.0
  endif
  if (cs%id_Kd_Lowmode_Work > 0) then
    allocate(dd%Kd_Lowmode_Work(isd:ied,jsd:jed,nz))
    dd%Kd_Lowmode_Work(:,:,:) = 0.0
  endif
  if (cs%id_TKE_itidal > 0) then
    allocate(dd%TKE_Itidal_used(isd:ied,jsd:jed)) ; dd%TKE_Itidal_used(:,:) = 0.
  endif
  ! additional diags for CVMix
  if (cs%id_N2_int > 0) then
    allocate(dd%N2_int(isd:ied,jsd:jed,nz+1)) ; dd%N2_int(:,:,:) = 0.0
  endif
  if (cs%id_Simmons_coeff > 0) then
    if (cs%CVMix_tidal_scheme .ne. simmons) then
      call mom_error(fatal, "setup_tidal_diagnostics: Simmons_coeff diagnostics is available "//&
                            "only when CVMix_tidal_scheme is Simmons")
    endif
    allocate(dd%Simmons_coeff_2d(isd:ied,jsd:jed)) ; dd%Simmons_coeff_2d(:,:) = 0.0
  endif
  if (cs%id_vert_dep > 0) then
    allocate(dd%vert_dep_3d(isd:ied,jsd:jed,nz+1)) ; dd%vert_dep_3d(:,:,:) = 0.0
  endif
  if (cs%id_Schmittner_coeff > 0) then
    if (cs%CVMix_tidal_scheme .ne. schmittner) then
      call mom_error(fatal, "setup_tidal_diagnostics: Schmittner_coeff diagnostics is available "//&
                            "only when CVMix_tidal_scheme is Schmittner.")
    endif
    allocate(dd%Schmittner_coeff_3d(isd:ied,jsd:jed,nz)) ; dd%Schmittner_coeff_3d(:,:,:) = 0.0
  endif
  if (cs%id_tidal_qe_md > 0) then
    if (cs%CVMix_tidal_scheme .ne. schmittner) then
      call mom_error(fatal, "setup_tidal_diagnostics: tidal_qe_md diagnostics is available "//&
                            "only when CVMix_tidal_scheme is Schmittner.")
    endif
    allocate(dd%tidal_qe_md(isd:ied,jsd:jed,nz)) ; dd%tidal_qe_md(:,:,:) = 0.0
  endif

tidal_mixing_end()

subroutine, public mom_tidal_mixing::tidal_mixing_end

(

type(tidal_mixing_cs), pointer

CS

)

Clear pointers and deallocate memory.

Parameters

cs	This module's control structure, which will be deallocated in this routine.

Definition at line 1724 of file MOM_tidal_mixing.F90.

  type(tidal_mixing_cs), pointer :: CS !< This module's control structure, which
                                       !! will be deallocated in this routine.
 
  if (.not.associated(cs)) return
 
  !TODO deallocate all the dynamically allocated members here ...
  if (allocated(cs%tidal_qe_2d))    deallocate(cs%tidal_qe_2d)
  if (allocated(cs%tidal_qe_3d_in)) deallocate(cs%tidal_qe_3d_in)
  if (allocated(cs%h_src))          deallocate(cs%h_src)
  deallocate(cs%dd)
  deallocate(cs)
 

tidal_mixing_h_amp()

subroutine, public mom_tidal_mixing::tidal_mixing_h_amp	(	real, dimension( g %isd: g %ied), intent(out)	h_amp,
		type(ocean_grid_type), intent(in)	G,
		integer, intent(in)	j,
		type(tidal_mixing_cs), pointer	CS
	)

This subroutine returns a zonal slice of the topographic roughness amplitudes.

Parameters

[in]	g	The ocean's grid structure
[out]	h_amp	The topographic roughness amplitude [Z ~> m]
[in]	j	j-index of the row to work on
	cs	The control structure for this module

Definition at line 1563 of file MOM_tidal_mixing.F90.

  type(ocean_grid_type),    intent(in)  :: G     !< The ocean's grid structure
  real, dimension(SZI_(G)), intent(out) :: h_amp !< The topographic roughness amplitude [Z ~> m]
  integer,                  intent(in)  :: j     !< j-index of the row to work on
  type(tidal_mixing_cs),    pointer     :: CS    !< The control structure for this module
 
  integer :: i
 
  h_amp(:) = 0.0
  if ( cs%Int_tide_dissipation .and. .not. cs%use_CVMix_tidal ) then
    do i=g%isc,g%iec
      h_amp(i) = sqrt(cs%h2(i,j))
    enddo
  endif
 

tidal_mixing_init()

logical function, public mom_tidal_mixing::tidal_mixing_init	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(tidal_mixing_cs), pointer	CS
	)

Initializes internal tidal dissipation scheme for diapycnal mixing.

Parameters

[in]	time	The current time.
[in]	g	Grid structure.
[in]	gv	Vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	param_file	Run-time parameter file handle
[in,out]	diag	Diagnostics control structure.
	cs	This module's control structure.

TODO: add units

Definition at line 213 of file MOM_tidal_mixing.F90.

  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(tidal_mixing_cs),    pointer       :: CS         !< This module's control structure.
 
  ! Local variables
  logical :: read_tideamp
  logical :: default_2018_answers
  character(len=20)  :: tmpstr, int_tide_profile_str
  character(len=20)  :: CVMix_tidal_scheme_str, tidal_energy_type
  character(len=200) :: filename, h2_file, Niku_TKE_input_file
  character(len=200) :: tidal_energy_file, tideamp_file
  real :: utide, hamp, prandtl_tidal, max_frac_rough
  real :: Niku_scale ! local variable for scaling the Nikurashin TKE flux data
  integer :: i, j, is, ie, js, je
  integer :: isd, ied, jsd, jed
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "MOM_tidal_mixing"     !< This module's name.
 
  if (associated(cs)) then
    call mom_error(warning, "tidal_mixing_init called when control structure "// &
                            "is already associated.")
    return
  endif
  allocate(cs)
  allocate(cs%dd)
 
  cs%debug = cs%debug.and.is_root_pe()
 
  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
 
  cs%diag => diag
 
  ! Read parameters
  call get_param(param_file, mdl, "USE_CVMix_TIDAL", cs%use_CVMix_tidal, &
                 default=.false., do_not_log=.true.)
  call get_param(param_file, mdl, "INT_TIDE_DISSIPATION", cs%int_tide_dissipation, &
                 default=cs%use_CVMix_tidal, do_not_log=.true.)
  call log_version(param_file, mdl, version, &
                 "Vertical Tidal Mixing Parameterization", &
                 all_default=.not.(cs%use_CVMix_tidal .or. cs%int_tide_dissipation))
  call get_param(param_file, mdl, "USE_CVMix_TIDAL", cs%use_CVMix_tidal, &
                 "If true, turns on tidal mixing via CVMix", &
                 default=.false.)
 
  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".",do_not_log=.true.)
  cs%inputdir = slasher(cs%inputdir)
  call get_param(param_file, mdl, "INT_TIDE_DISSIPATION", cs%int_tide_dissipation, &
                 "If true, use an internal tidal dissipation scheme to "//&
                 "drive diapycnal mixing, along the lines of St. Laurent "//&
                 "et al. (2002) and Simmons et al. (2004).", default=cs%use_CVMix_tidal)
 
  ! return if tidal mixing is inactive
  tidal_mixing_init = cs%int_tide_dissipation
  if (.not. tidal_mixing_init) return
 
  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDAL_MIXING_2018_ANSWERS", cs%answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)
  call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", cs%remap_answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)
 
  if (cs%int_tide_dissipation) then
 
    ! Read in CVMix tidal scheme if CVMix tidal mixing is on
    if (cs%use_CVMix_tidal) then
      call get_param(param_file, mdl, "CVMIX_TIDAL_SCHEME", cvmix_tidal_scheme_str, &
                 "CVMIX_TIDAL_SCHEME selects the CVMix tidal mixing "//&
                 "scheme with INT_TIDE_DISSIPATION. Valid values are:\n"//&
                 "\t SIMMONS - Use the Simmons et al (2004) tidal \n"//&
                 "\t                mixing scheme.\n"//&
                 "\t SCHMITTNER - Use the Schmittner et al (2014) tidal \n"//&
                 "\t                mixing scheme.", &
                 default=simmons_scheme_string)
      cvmix_tidal_scheme_str = uppercase(cvmix_tidal_scheme_str)
 
      select case (cvmix_tidal_scheme_str)
        case (simmons_scheme_string)    ; cs%CVMix_tidal_scheme = simmons
        case (schmittner_scheme_string) ; cs%CVMix_tidal_scheme = schmittner
        case default
        call mom_error(fatal, "tidal_mixing_init: Unrecognized setting "// &
            "#define CVMIX_TIDAL_SCHEME "//trim(cvmix_tidal_scheme_str)//" found in input file.")
      end select
    endif ! CS%use_CVMix_tidal
 
    ! Read in vertical profile of tidal energy dissipation
    if ( cs%CVMix_tidal_scheme.eq.schmittner .or. .not. cs%use_CVMix_tidal) then
      call get_param(param_file, mdl, "INT_TIDE_PROFILE", int_tide_profile_str, &
                   "INT_TIDE_PROFILE selects the vertical profile of energy "//&
                   "dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//&
                   "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//&
                   "\t                decay profile.\n"//&
                   "\t POLZIN_09 - Use the Polzin WKB-stretched algebraic \n"//&
                   "\t                decay profile.", &
                   default=stlaurent_profile_string)
      int_tide_profile_str = uppercase(int_tide_profile_str)
 
      select case (int_tide_profile_str)
        case (stlaurent_profile_string)   ; cs%int_tide_profile = stlaurent_02
        case (polzin_profile_string)      ; cs%int_tide_profile = polzin_09
        case default
          call mom_error(fatal, "tidal_mixing_init: Unrecognized setting "// &
              "#define INT_TIDE_PROFILE "//trim(int_tide_profile_str)//" found in input file.")
      end select
    endif
 
  elseif (cs%use_CVMix_tidal) then
        call mom_error(fatal, "tidal_mixing_init: Cannot set INT_TIDE_DISSIPATION to False "// &
            "when USE_CVMix_TIDAL is set to True.")
  endif
 
  call get_param(param_file, mdl, "LEE_WAVE_DISSIPATION", cs%Lee_wave_dissipation, &
                 "If true, use an lee wave driven dissipation scheme to "//&
                 "drive diapycnal mixing, along the lines of Nikurashin "//&
                 "(2010) and using the St. Laurent et al. (2002) "//&
                 "and Simmons et al. (2004) vertical profile", default=.false.)
  if (cs%lee_wave_dissipation) then
    if (cs%use_CVMix_tidal) then
        call mom_error(fatal, "tidal_mixing_init: Lee wave driven dissipation scheme cannot "// &
            "be used when CVMix tidal mixing scheme is active.")
    endif
    call get_param(param_file, mdl, "LEE_WAVE_PROFILE", tmpstr, &
                 "LEE_WAVE_PROFILE selects the vertical profile of energy "//&
                 "dissipation with LEE_WAVE_DISSIPATION. Valid values are:\n"//&
                 "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//&
                 "\t                decay profile.\n"//&
                 "\t POLZIN_09 - Use the Polzin WKB-stretched algebraic \n"//&
                 "\t                decay profile.", &
                 default=stlaurent_profile_string)
    tmpstr = uppercase(tmpstr)
    select case (tmpstr)
      case (stlaurent_profile_string) ; cs%lee_wave_profile = stlaurent_02
      case (polzin_profile_string) ; cs%lee_wave_profile = polzin_09
      case default
        call mom_error(fatal, "tidal_mixing_init: Unrecognized setting "// &
            "#define LEE_WAVE_PROFILE "//trim(tmpstr)//" found in input file.")
    end select
  endif
 
  call get_param(param_file, mdl, "INT_TIDE_LOWMODE_DISSIPATION", cs%Lowmode_itidal_dissipation, &
                 "If true, consider mixing due to breaking low modes that "//&
                 "have been remotely generated; as with itidal drag on the "//&
                 "barotropic tide, use an internal tidal dissipation scheme to "//&
                 "drive diapycnal mixing, along the lines of St. Laurent "//&
                 "et al. (2002) and Simmons et al. (2004).", default=.false.)
 
  if ((cs%Int_tide_dissipation .and. (cs%int_tide_profile == polzin_09)) .or. &
      (cs%lee_wave_dissipation .and. (cs%lee_wave_profile == polzin_09))) then
    if (cs%use_CVMix_tidal) then
        call mom_error(fatal, "tidal_mixing_init: Polzin scheme cannot "// &
            "be used when CVMix tidal mixing scheme is active.")
    endif
    call get_param(param_file, mdl, "NU_POLZIN", cs%Nu_Polzin, &
                 "When the Polzin decay profile is used, this is a "//&
                 "non-dimensional constant in the expression for the "//&
                 "vertical scale of decay for the tidal energy dissipation.", &
                 units="nondim", default=0.0697)
    call get_param(param_file, mdl, "NBOTREF_POLZIN", cs%Nbotref_Polzin, &
                 "When the Polzin decay profile is used, this is the "//&
                 "reference value of the buoyancy frequency at the ocean "//&
                 "bottom in the Polzin formulation for the vertical "//&
                 "scale of decay for the tidal energy dissipation.", &
                 units="s-1", default=9.61e-4, scale=us%T_to_s)
    call get_param(param_file, mdl, "POLZIN_DECAY_SCALE_FACTOR", &
                 cs%Polzin_decay_scale_factor, &
                 "When the Polzin decay profile is used, this is a "//&
                 "scale factor for the vertical scale of decay of the tidal "//&
                 "energy dissipation.", default=1.0, units="nondim")
    call get_param(param_file, mdl, "POLZIN_SCALE_MAX_FACTOR", &
                 cs%Polzin_decay_scale_max_factor, &
                 "When the Polzin decay profile is used, this is a factor "//&
                 "to limit the vertical scale of decay of the tidal "//&
                 "energy dissipation to POLZIN_DECAY_SCALE_MAX_FACTOR "//&
                 "times the depth of the ocean.", units="nondim", default=1.0)
    call get_param(param_file, mdl, "POLZIN_MIN_DECAY_SCALE", cs%Polzin_min_decay_scale, &
                 "When the Polzin decay profile is used, this is the "//&
                 "minimum vertical decay scale for the vertical profile\n"//&
                 "of internal tide dissipation with the Polzin (2009) formulation", &
                 units="m", default=0.0, scale=us%m_to_Z)
  endif
 
  if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation) then
    call get_param(param_file, mdl, "INT_TIDE_DECAY_SCALE", cs%Int_tide_decay_scale, &
                 "The decay scale away from the bottom for tidal TKE with "//&
                 "the new coding when INT_TIDE_DISSIPATION is used.", &
                 !units="m", default=0.0)
                 units="m", default=500.0, scale=us%m_to_Z)  ! TODO: confirm this new default
    call get_param(param_file, mdl, "MU_ITIDES", cs%Mu_itides, &
                 "A dimensionless turbulent mixing efficiency used with "//&
                 "INT_TIDE_DISSIPATION, often 0.2.", units="nondim", default=0.2)
    call get_param(param_file, mdl, "GAMMA_ITIDES", cs%Gamma_itides, &
                 "The fraction of the internal tidal energy that is "//&
                 "dissipated locally with INT_TIDE_DISSIPATION. "//&
                 "THIS NAME COULD BE BETTER.", &
                 units="nondim", default=0.3333)
    call get_param(param_file, mdl, "MIN_ZBOT_ITIDES", cs%min_zbot_itides, &
                 "Turn off internal tidal dissipation when the total "//&
                 "ocean depth is less than this value.", units="m", default=0.0, scale=us%m_to_Z)
  endif
 
  if ( (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation) .and. &
        .not. cs%use_CVMix_tidal) then
 
    call safe_alloc_ptr(cs%Nb,isd,ied,jsd,jed)
    call safe_alloc_ptr(cs%h2,isd,ied,jsd,jed)
    call safe_alloc_ptr(cs%TKE_itidal,isd,ied,jsd,jed)
    call safe_alloc_ptr(cs%mask_itidal,isd,ied,jsd,jed) ; cs%mask_itidal(:,:) = 1.0
 
    call get_param(param_file, mdl, "KAPPA_ITIDES", cs%kappa_itides, &
                 "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
                 "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
                 units="m-1", default=8.e-4*atan(1.0), scale=us%Z_to_m)
 
    call get_param(param_file, mdl, "UTIDE", cs%utide, &
                 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
                 units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
    call safe_alloc_ptr(cs%tideamp,is,ie,js,je) ; cs%tideamp(:,:) = cs%utide
 
    call get_param(param_file, mdl, "KAPPA_H2_FACTOR", cs%kappa_h2_factor, &
                 "A scaling factor for the roughness amplitude with "//&
                 "INT_TIDE_DISSIPATION.",  units="nondim", default=1.0)
    call get_param(param_file, mdl, "TKE_ITIDE_MAX", cs%TKE_itide_max, &
                 "The maximum internal tide energy source available to mix "//&
                 "above the bottom boundary layer with INT_TIDE_DISSIPATION.", &
                 units="W m-2", default=1.0e3, scale=us%W_m2_to_RZ3_T3)
 
    call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, &
                 "If true, read a file (given by TIDEAMP_FILE) containing "//&
                 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
    if (read_tideamp) then
      if (cs%use_CVMix_tidal) then
          call mom_error(fatal, "tidal_mixing_init: Tidal amplitude files are "// &
              "not compatible with CVMix tidal mixing. ")
      endif
      call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
                 "The path to the file containing the spatially varying "//&
                 "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc")
      filename = trim(cs%inputdir) // trim(tideamp_file)
      call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename)
      call mom_read_data(filename, 'tideamp', cs%tideamp, g%domain, timelevel=1, scale=us%m_to_Z*us%T_to_s)
    endif
 
    call get_param(param_file, mdl, "H2_FILE", h2_file, &
                 "The path to the file containing the sub-grid-scale "//&
                 "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
                 fail_if_missing=(.not.cs%use_CVMix_tidal))
    filename = trim(cs%inputdir) // trim(h2_file)
    call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename)
    call mom_read_data(filename, 'h2', cs%h2, g%domain, timelevel=1, scale=us%m_to_Z**2)
 
    call get_param(param_file, mdl, "FRACTIONAL_ROUGHNESS_MAX", max_frac_rough, &
                 "The maximum topographic roughness amplitude as a fraction of the mean depth, "//&
                 "or a negative value for no limitations on roughness.", &
                 units="nondim", default=0.1)
 
    do j=js,je ; do i=is,ie
      if (g%bathyT(i,j) < cs%min_zbot_itides) cs%mask_itidal(i,j) = 0.0
      cs%tideamp(i,j) = cs%tideamp(i,j) * cs%mask_itidal(i,j) * g%mask2dT(i,j)
 
      ! Restrict rms topo to a fraction (often 10 percent) of the column depth.
      if (cs%answers_2018 .and. (max_frac_rough >= 0.0)) then
        hamp = min(max_frac_rough*g%bathyT(i,j), sqrt(cs%h2(i,j)))
        cs%h2(i,j) = hamp*hamp
      else
        if (max_frac_rough >= 0.0) &
          cs%h2(i,j) = min((max_frac_rough*g%bathyT(i,j))**2, cs%h2(i,j))
      endif
 
      utide = cs%tideamp(i,j)
      ! Compute the fixed part of internal tidal forcing.
      ! The units here are [R Z3 T-2 ~> J m-2 = kg s-2] here.
      cs%TKE_itidal(i,j) = 0.5 * cs%kappa_h2_factor * gv%Rho0 * &
           cs%kappa_itides * cs%h2(i,j) * utide*utide
    enddo ; enddo
 
  endif
 
  if (cs%Lee_wave_dissipation) then
 
    call get_param(param_file, mdl, "NIKURASHIN_TKE_INPUT_FILE",niku_tke_input_file, &
                 "The path to the file containing the TKE input from lee "//&
                 "wave driven mixing. Used with LEE_WAVE_DISSIPATION.", &
                 fail_if_missing=.true.)
    call get_param(param_file, mdl, "NIKURASHIN_SCALE",niku_scale, &
                 "A non-dimensional factor by which to scale the lee-wave "//&
                 "driven TKE input. Used with LEE_WAVE_DISSIPATION.", &
                 units="nondim", default=1.0)
 
    filename = trim(cs%inputdir) // trim(niku_tke_input_file)
    call log_param(param_file, mdl, "INPUTDIR/NIKURASHIN_TKE_INPUT_FILE", &
                   filename)
    call safe_alloc_ptr(cs%TKE_Niku,is,ie,js,je) ; cs%TKE_Niku(:,:) = 0.0
    call mom_read_data(filename, 'TKE_input', cs%TKE_Niku, g%domain, timelevel=1, &  ! ??? timelevel -aja
                       scale=niku_scale*us%W_m2_to_RZ3_T3)
 
    call get_param(param_file, mdl, "GAMMA_NIKURASHIN",cs%Gamma_lee, &
                 "The fraction of the lee wave energy that is dissipated "//&
                 "locally with LEE_WAVE_DISSIPATION.", units="nondim", &
                 default=0.3333)
    call get_param(param_file, mdl, "DECAY_SCALE_FACTOR_LEE",cs%Decay_scale_factor_lee, &
                 "Scaling for the vertical decay scaleof the local "//&
                 "dissipation of lee waves dissipation.", units="nondim", &
                 default=1.0)
  else
    cs%Decay_scale_factor_lee = -9.e99 ! This should never be used if CS%Lee_wave_dissipation = False
  endif
 
  ! Configure CVMix
  if (cs%use_CVMix_tidal) then
 
    ! Read in CVMix params
    !call openParameterBlock(param_file,'CVMix_TIDAL')
    call get_param(param_file, mdl, "TIDAL_MAX_COEF", cs%tidal_max_coef, &
                   "largest acceptable value for tidal diffusivity", &
                   units="m^2/s", default=50e-4) ! the default is 50e-4 in CVMix, 100e-4 in POP.
    call get_param(param_file, mdl, "TIDAL_DISS_LIM_TC", cs%tidal_diss_lim_tc, &
                   "Min allowable depth for dissipation for tidal-energy-constituent data. "//&
                   "No dissipation contribution is applied above TIDAL_DISS_LIM_TC.", &
                   units="m", default=0.0, scale=us%m_to_Z)
    call get_param(param_file, mdl, "TIDAL_ENERGY_FILE",tidal_energy_file, &
                 "The path to the file containing tidal energy "//&
                 "dissipation. Used with CVMix tidal mixing schemes.", &
                 fail_if_missing=.true.)
    call get_param(param_file, mdl, 'MIN_THICKNESS', cs%min_thickness, default=0.001, &
                   do_not_log=.true.)
    call get_param(param_file, mdl, "PRANDTL_TIDAL", prandtl_tidal, &
                   "Prandtl number used by CVMix tidal mixing schemes "//&
                   "to convert vertical diffusivities into viscosities.", &
                    units="nondim", default=1.0, &
                   do_not_log=.true.)
    call cvmix_put(cs%CVMix_glb_params,'Prandtl',prandtl_tidal)
 
    tidal_energy_file = trim(cs%inputdir) // trim(tidal_energy_file)
    call get_param(param_file, mdl, "TIDAL_ENERGY_TYPE",tidal_energy_type, &
                 "The type of input tidal energy flux dataset. Valid values are"//&
                   "\t Jayne\n"//&
                   "\t ER03 \n",&
                 fail_if_missing=.true.)
    ! Check whether tidal energy input format and CVMix tidal mixing scheme are consistent
    if ( .not. ( &
          (uppercase(tidal_energy_type(1:4)).eq.'JAYN' .and. cs%CVMix_tidal_scheme.eq.simmons).or. &
          (uppercase(tidal_energy_type(1:4)).eq.'ER03' .and. cs%CVMix_tidal_scheme.eq.schmittner) ) )then
        call mom_error(fatal, "tidal_mixing_init: Tidal energy file type ("//&
                      trim(tidal_energy_type)//") is incompatible with CVMix tidal "//&
                      " mixing scheme: "//trim(cvmix_tidal_scheme_str) )
    endif
    cvmix_tidal_scheme_str = lowercase(cvmix_tidal_scheme_str)
 
    ! Set up CVMix
    call cvmix_init_tidal(cvmix_tidal_params_user = cs%CVMix_tidal_params,    &
                          mix_scheme              = cvmix_tidal_scheme_str,   &
                          efficiency              = cs%Mu_itides,             &
                          vertical_decay_scale    = cs%int_tide_decay_scale*us%Z_to_m,  &
                          max_coefficient         = cs%tidal_max_coef,        &
                          local_mixing_frac       = cs%Gamma_itides,          &
                          depth_cutoff            = cs%min_zbot_itides*us%Z_to_m)
 
    call read_tidal_energy(g, us, tidal_energy_type, tidal_energy_file, cs)
 
    !call closeParameterBlock(param_file)
 
  endif ! CVMix on
 
  ! Register Diagnostics fields
 
  if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. &
      cs%Lowmode_itidal_dissipation) then
 
    cs%id_Kd_itidal = register_diag_field('ocean_model','Kd_itides',diag%axesTi,time, &
         'Internal Tide Driven Diffusivity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
 
    if (cs%use_CVMix_tidal) then
      cs%id_N2_int = register_diag_field('ocean_model','N2_int',diag%axesTi,time, &
          'Bouyancy frequency squared, at interfaces', 's-2', conversion=us%s_to_T**2)
      !> TODO: add units
      cs%id_Simmons_coeff = register_diag_field('ocean_model','Simmons_coeff',diag%axesT1,time, &
           'time-invariant portion of the tidal mixing coefficient using the Simmons', '')
      cs%id_Schmittner_coeff = register_diag_field('ocean_model','Schmittner_coeff',diag%axesTL,time, &
           'time-invariant portion of the tidal mixing coefficient using the Schmittner', '')
      cs%id_tidal_qe_md = register_diag_field('ocean_model','tidal_qe_md',diag%axesTL,time, &
           'input tidal energy dissipated locally interpolated to model vertical coordinates', '')
      cs%id_vert_dep = register_diag_field('ocean_model','vert_dep',diag%axesTi,time, &
           'vertical deposition function needed for Simmons et al tidal  mixing', '')
 
    else
      cs%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal',diag%axesT1,time, &
          'Internal Tide Driven Turbulent Kinetic Energy', &
          'W m-2', conversion=us%RZ3_T3_to_W_m2)
      cs%id_Nb = register_diag_field('ocean_model','Nb',diag%axesT1,time, &
           'Bottom Buoyancy Frequency', 's-1', conversion=us%s_to_T)
 
      cs%id_Kd_lowmode = register_diag_field('ocean_model','Kd_lowmode',diag%axesTi,time, &
           'Internal Tide Driven Diffusivity (from propagating low modes)', &
           'm2 s-1', conversion=us%Z2_T_to_m2_s)
 
      cs%id_Fl_itidal = register_diag_field('ocean_model','Fl_itides',diag%axesTi,time, &
          'Vertical flux of tidal turbulent dissipation', &
          'm3 s-3', conversion=(us%Z_to_m**3*us%s_to_T**3))
 
      cs%id_Fl_lowmode = register_diag_field('ocean_model','Fl_lowmode',diag%axesTi,time, &
           'Vertical flux of tidal turbulent dissipation (from propagating low modes)', &
           'm3 s-3', conversion=(us%Z_to_m**3*us%s_to_T**3))
 
      cs%id_Polzin_decay_scale = register_diag_field('ocean_model','Polzin_decay_scale',diag%axesT1,time, &
           'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', &
           'm', conversion=us%Z_to_m)
 
      cs%id_Polzin_decay_scale_scaled = register_diag_field('ocean_model', &
           'Polzin_decay_scale_scaled', diag%axesT1, time, &
           'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme, '// &
           'scaled by N2_bot/N2_meanz', 'm', conversion=us%Z_to_m)
 
      cs%id_N2_bot = register_diag_field('ocean_model','N2_b',diag%axesT1,time, &
           'Bottom Buoyancy frequency squared', 's-2', conversion=us%s_to_T**2)
 
      cs%id_N2_meanz = register_diag_field('ocean_model','N2_meanz', diag%axesT1, time, &
           'Buoyancy frequency squared averaged over the water column', 's-2', conversion=us%s_to_T**2)
 
      cs%id_Kd_Itidal_Work = register_diag_field('ocean_model','Kd_Itidal_Work',diag%axesTL,time, &
           'Work done by Internal Tide Diapycnal Mixing', &
           'W m-2', conversion=us%RZ3_T3_to_W_m2)
 
      cs%id_Kd_Niku_Work = register_diag_field('ocean_model','Kd_Nikurashin_Work',diag%axesTL,time, &
           'Work done by Nikurashin Lee Wave Drag Scheme', &
           'W m-2', conversion=us%RZ3_T3_to_W_m2)
 
      cs%id_Kd_Lowmode_Work = register_diag_field('ocean_model','Kd_Lowmode_Work',diag%axesTL,time, &
           'Work done by Internal Tide Diapycnal Mixing (low modes)', &
           'W m-2', conversion=us%RZ3_T3_to_W_m2)
 
      if (cs%Lee_wave_dissipation) then
        cs%id_TKE_leewave = register_diag_field('ocean_model','TKE_leewave',diag%axesT1,time, &
            'Lee wave Driven Turbulent Kinetic Energy', &
            'W m-2', conversion=us%RZ3_T3_to_W_m2)
        cs%id_Kd_Niku = register_diag_field('ocean_model','Kd_Nikurashin',diag%axesTi,time, &
            'Lee Wave Driven Diffusivity', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
      endif
    endif ! S%use_CVMix_tidal
  endif