MOM_internal_tide_input.F90

Go to the documentation of this file.

module mom_int_tide_input
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, January 2013                                   *
!*                                                                     *
!*    This file contains the subroutines that sets the energy input    *
!*  to the internal tides.                                             *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v                                        *
!*    j    x ^ x ^ x   At >:  u                                        *
!*    j    > o > o >   At o:  h, buoy, ustar, T, S, Kd, ea, eb, etc.   *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
use mom_diag_mediator, only : diag_ctrl, time_type
use mom_diag_mediator, only : safe_alloc_ptr, post_data, register_diag_field
use mom_diag_to_z, only : diag_to_z_cs, register_zint_diag, calc_zint_diags
use mom_debugging, only : hchksum
use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_forcing_type, only : forcing
use mom_grid, only : ocean_grid_type
use mom_io, only : slasher, vardesc
use mom_thickness_diffuse, only : vert_fill_ts
use mom_variables, only : thermo_var_ptrs, vertvisc_type, p3d
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs

use fms_mod, only : read_data

implicit none ; private

#include <MOM_memory.h>

public set_int_tide_input, int_tide_input_init, int_tide_input_end

type, public :: int_tide_input_cs ; private
  logical :: debug           ! If true, write verbose checksums for debugging.
  type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                             ! timing of diagnostic output.
  real :: tke_itide_max ! Maximum Internal tide conversion (W m-2)
                        ! available to mix above the BBL

  real, allocatable, dimension(:,:) :: &
    tke_itidal_coef     ! The time-invariant field that enters the TKE_itidal
                        ! input calculation, in J m-2.

  integer :: id_tke_itidal = -1, id_nb = -1, id_n2_bot = -1
  character(len=200) :: inputdir
end type int_tide_input_cs

type, public :: int_tide_input_type
  real, allocatable, dimension(:,:) :: &
    tke_itidal_input, & ! The internal tide TKE input at the bottom of
                        ! the ocean, in W m-2.
    h2, &               ! The squared topographic roughness height, in m2.
    tideamp, &          ! The amplitude of the tidal velocities, in m s-1.
    nb                  ! The bottom stratification, in s-1.
end type int_tide_input_type

contains

subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, CS)
  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(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< The zonal velocity, in m s-1
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< The meridional velocity, in m s-1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                     intent(in)    :: tv
  type(forcing),                             intent(in)    :: fluxes
  type(int_tide_input_type),                 intent(inout) :: itide
  real,                                      intent(in)    :: dt
  type(int_tide_input_cs),                   pointer       :: CS

! Arguments: u - Zonal velocity, in m s-1.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m or kg m-2.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      fluxes - A structure of surface fluxes that may be used.
!  (inout)   itide - A structure containing fields related to the internal
!                    tide sources.
!  (in)      dt - The time increment in s.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - This module's control structure.

  real, dimension(SZI_(G),SZJ_(G)) :: &
    N2_bot        ! The bottom squared buoyancy frequency, in s-2.

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    T_f, S_f      ! The temperature and salinity in C and PSU with the values in
                  ! the massless layers filled vertically by diffusion.
  logical :: use_EOS    ! If true, density is calculated from T & S using an
                        ! equation of state.
  integer :: i, j, k, is, ie, js, je, nz
  integer :: isd, ied, jsd, jed

  real :: kappa_fill  ! diffusivity used to fill massless layers
  real :: dt_fill     ! timestep used to fill massless layers

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  if (.not.associated(cs)) call mom_error(fatal,"set_diffusivity: "//&
         "Module must be initialized before it is used.")

  kappa_fill = 1.e-3 ! m2 s-1
  dt_fill = 7200.

  use_eos = associated(tv%eqn_of_state)

  ! Smooth the properties through massless layers.
  if (use_eos) then
    call vert_fill_ts(h, tv%T, tv%S, kappa_fill, dt_fill, t_f, s_f, g, gv)
  endif

  call find_n2_bottom(h, tv, t_f, s_f, itide%h2, fluxes, g, gv, n2_bot)

!$OMP parallel do default(none) shared(is,ie,js,je,G,itide,N2_bot,CS)
  do j=js,je ; do i=is,ie
    itide%Nb(i,j) = g%mask2dT(i,j) * sqrt(n2_bot(i,j))
    itide%TKE_itidal_input(i,j) = min(cs%TKE_itidal_coef(i,j)*itide%Nb(i,j),cs%TKE_itide_max)
  enddo ; enddo

  if (cs%debug) then
    call hchksum(n2_bot,"N2_bot",g%HI,haloshift=0)
    call hchksum(itide%TKE_itidal_input,"TKE_itidal_input",g%HI,haloshift=0)
  endif

  if (cs%id_TKE_itidal > 0) call post_data(cs%id_TKE_itidal, itide%TKE_itidal_input, cs%diag)
  if (cs%id_Nb > 0) call post_data(cs%id_Nb, itide%Nb, cs%diag)
  if (cs%id_N2_bot > 0 ) call post_data(cs%id_N2_bot,n2_bot,cs%diag)

end subroutine set_int_tide_input

subroutine find_n2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, N2_bot)
  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, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                    intent(in)   :: tv
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)   :: T_f, S_f
  real, dimension(SZI_(G),SZJ_(G)),         intent(in)   :: h2
  type(forcing),                            intent(in)   :: fluxes
  type(int_tide_input_cs),                  pointer      :: CS
  real, dimension(SZI_(G),SZJ_(G)),         intent(out)  :: N2_bot

  real, dimension(SZI_(G),SZK_(G)+1) :: &
    dRho_int      ! The unfiltered density differences across interfaces.
  real, dimension(SZI_(G)) :: &
    pres, &       ! The pressure at each interface, in Pa.
    Temp_int, &   ! The temperature at each interface, in degC.
    Salin_int, &  ! The salinity at each interface, in PSU.
    drho_bot, &
    h_amp, &
    hb, &
    z_from_bot, &
    dRho_dT, &    ! The partial derivatives of density with temperature and
    dRho_dS       ! salinity, in kg m-3 degC-1 and kg m-3 PSU-1.

  real :: dz_int  ! The thickness associated with an interface, in m.
  real :: G_Rho0  ! The gravitation acceleration divided by the Boussinesq
                  ! density, in m4 s-2 kg-1.
  logical :: do_i(szi_(g)), do_any
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  g_rho0 = gv%g_Earth / gv%Rho0

  ! Find the (limited) density jump across each interface.
  do i=is,ie
    drho_int(i,1) = 0.0 ; drho_int(i,nz+1) = 0.0
  enddo
!$OMP parallel do default(none) shared(is,ie,js,je,nz,tv,fluxes,G,GV,h,T_f,S_f, &
!$OMP                                  h2,N2_bot,G_Rho0) &
!$OMP                          private(pres,Temp_Int,Salin_Int,dRho_dT,dRho_dS, &
!$OMP                                  hb,dRho_bot,z_from_bot,do_i,h_amp,       &
!$OMP                                  do_any,dz_int) &
!$OMP                     firstprivate(dRho_int)
  do j=js,je
    if (associated(tv%eqn_of_state)) then
      if (associated(fluxes%p_surf)) then
        do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo
      else
        do i=is,ie ; pres(i) = 0.0 ; enddo
      endif
      do k=2,nz
        do i=is,ie
          pres(i) = pres(i) + gv%H_to_Pa*h(i,j,k-1)
          temp_int(i) = 0.5 * (t_f(i,j,k) + t_f(i,j,k-1))
          salin_int(i) = 0.5 * (s_f(i,j,k) + s_f(i,j,k-1))
        enddo
        call calculate_density_derivs(temp_int, salin_int, pres, &
                 drho_dt(:), drho_ds(:), is, ie-is+1, tv%eqn_of_state)
        do i=is,ie
          drho_int(i,k) = max(drho_dt(i)*(t_f(i,j,k) - t_f(i,j,k-1)) + &
                              drho_ds(i)*(s_f(i,j,k) - s_f(i,j,k-1)), 0.0)
        enddo
      enddo
    else
      do k=2,nz ; do i=is,ie
        drho_int(i,k) = gv%Rlay(k) - gv%Rlay(k-1)
      enddo ; enddo
    endif

    ! Find the bottom boundary layer stratification.
    do i=is,ie
      hb(i) = 0.0 ; drho_bot(i) = 0.0
      z_from_bot(i) = 0.5*gv%H_to_m*h(i,j,nz)
      do_i(i) = (g%mask2dT(i,j) > 0.5)
      h_amp(i) = sqrt(h2(i,j))
    enddo

    do k=nz,2,-1
      do_any = .false.
      do i=is,ie ; if (do_i(i)) then
        dz_int = 0.5*gv%H_to_m*(h(i,j,k) + h(i,j,k-1))
        z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above

        hb(i) = hb(i) + dz_int
        drho_bot(i) = drho_bot(i) + drho_int(i,k)

        if (z_from_bot(i) > h_amp(i)) then
          if (k>2) then
            ! Always include at least one full layer.
            hb(i) = hb(i) + 0.5*gv%H_to_m*(h(i,j,k-1) + h(i,j,k-2))
            drho_bot(i) = drho_bot(i) + drho_int(i,k-1)
          endif
          do_i(i) = .false.
        else
          do_any = .true.
        endif
      endif ; enddo
      if (.not.do_any) exit
    enddo

    do i=is,ie
      if (hb(i) > 0.0) then
        n2_bot(i,j) = (g_rho0 * drho_bot(i)) / hb(i)
      else ;  n2_bot(i,j) = 0.0 ; endif
    enddo
  enddo

end subroutine find_n2_bottom

subroutine int_tide_input_init(Time, G, GV, param_file, diag, CS, itide)
  type(time_type),          intent(in)    :: Time
  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(param_file_type),    intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(diag_ctrl), target,  intent(inout) :: diag
  type(int_tide_input_cs),   pointer      :: CS
  type(int_tide_input_type), pointer      :: itide
! Arguments: Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module
!  (in)      diag_to_Z_CSp - A pointer to the Z-diagnostics control structure.
  type(vardesc) :: vd
  logical :: read_tideamp
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_int_tide_input"  ! This module's name.
  character(len=20)  :: tmpstr
  character(len=200) :: filename, tideamp_file, h2_file

  real :: mask_itidal
  real :: utide              ! constant tidal amplitude (m s-1) to be used if
                             ! tidal amplitude file is not present.
  real :: kappa_h2_factor    ! factor for the product of wavenumber * rms sgs height.
  real :: kappa_itides       ! topographic wavenumber and non-dimensional scaling
  real :: min_zbot_itides    ! Minimum ocean depth for internal tide conversion,
                             ! in m.
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed

  if (associated(cs)) then
    call mom_error(warning, "int_tide_input_init called with an associated "// &
                            "control structure.")
    return
  endif
  if (associated(itide)) then
    call mom_error(warning, "int_tide_input_init called with an associated "// &
                            "internal tide input type.")
    return
  endif
  allocate(cs)
  allocate(itide)

  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 all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")

  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
  cs%inputdir = slasher(cs%inputdir)

  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)

  call get_param(param_file, mdl, "MIN_ZBOT_ITIDES", min_zbot_itides, &
               "Turn off internal tidal dissipation when the total \n"//&
               "ocean depth is less than this value.", units="m", default=0.0)

  call get_param(param_file, mdl, "UTIDE", utide, &
               "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
               units="m s-1", default=0.0)

  allocate(itide%Nb(isd:ied,jsd:jed))  ; itide%Nb(:,:) = 0.0
  allocate(itide%h2(isd:ied,jsd:jed))  ; itide%h2(:,:) = 0.0
  allocate(itide%TKE_itidal_input(isd:ied,jsd:jed)) ; itide%TKE_itidal_input(:,:) = 0.0
  allocate(itide%tideamp(isd:ied,jsd:jed)) ; itide%tideamp(:,:) = utide
  allocate(cs%TKE_itidal_coef(isd:ied,jsd:jed)) ; cs%TKE_itidal_coef(:,:) = 0.0

  call get_param(param_file, mdl, "KAPPA_ITIDES", kappa_itides, &
               "A topographic wavenumber used with INT_TIDE_DISSIPATION. \n"//&
               "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
               units="m-1", default=8.e-4*atan(1.0))

  call get_param(param_file, mdl, "KAPPA_H2_FACTOR", kappa_h2_factor, &
               "A scaling factor for the roughness amplitude with n"//&
               "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 availble to mix \n"//&
               "above the bottom boundary layer with INT_TIDE_DISSIPATION.", &
               units="W m-2",  default=1.0e3)

  call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, &
               "If true, read a file (given by TIDEAMP_FILE) containing \n"//&
               "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
  if (read_tideamp) then
    call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
               "The path to the file containing the spatially varying \n"//&
               "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 read_data(filename, 'tideamp', itide%tideamp, &
                   domain=g%domain%mpp_domain, timelevel=1)
  endif

  call get_param(param_file, mdl, "H2_FILE", h2_file, &
               "The path to the file containing the sub-grid-scale \n"//&
               "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
               fail_if_missing=.true.)
  filename = trim(cs%inputdir) // trim(h2_file)
  call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename)
  call read_data(filename, 'h2', itide%h2, domain=g%domain%mpp_domain, &
                 timelevel=1)

  do j=js,je ; do i=is,ie
    mask_itidal = 1.0
    if (g%bathyT(i,j) < min_zbot_itides) mask_itidal = 0.0

    itide%tideamp(i,j) = itide%tideamp(i,j) * mask_itidal * g%mask2dT(i,j)

    ! Restrict rms topo to 10 percent of column depth.
    itide%h2(i,j) = min(0.01*g%bathyT(i,j)**2, itide%h2(i,j))

    ! Compute the fixed part of internal tidal forcing; units are [J m-2] here.
    cs%TKE_itidal_coef(i,j) = 0.5*kappa_h2_factor*gv%Rho0*&
         kappa_itides * itide%h2(i,j) * itide%tideamp(i,j)**2
  enddo; enddo


  cs%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal_itide',diag%axesT1,time, &
      'Internal Tide Driven Turbulent Kinetic Energy', 'Watt meter-2')

  cs%id_Nb = register_diag_field('ocean_model','Nb_itide',diag%axesT1,time, &
       'Bottom Buoyancy Frequency', 'sec-1')

  cs%id_N2_bot = register_diag_field('ocean_model','N2_b_itide',diag%axesT1,time, &
       'Bottom Buoyancy frequency squared', 's-2')

end subroutine int_tide_input_init

subroutine int_tide_input_end(CS)
  type(int_tide_input_cs), pointer :: CS

  if (associated(cs)) deallocate(cs)

end subroutine int_tide_input_end

end module mom_int_tide_input