MOM6
mom_internal_tides Module Reference

Detailed Description

Subroutines that use the ray-tracing equations to propagate the internal tide energy density.

Author
Benjamin Mater & Robert Hallberg, 2015

Data Types

type  int_tide_cs
 This control structure has parameters for the MOM_internal_tides module. More...
 
type  loop_bounds_type
 A structure with the active energy loop bounds. More...
 

Functions/Subroutines

subroutine, public propagate_int_tide (h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, G, GV, US, CS)
 Calls subroutines in this file that are needed to refract, propagate, and dissipate energy density of the internal tide. More...
 
subroutine sum_en (G, CS, En, label)
 Checks for energy conservation on computational domain. More...
 
subroutine itidal_lowmode_loss (G, US, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss, dt, full_halos)
 Calculates the energy lost from the propagating internal tide due to scattering over small-scale roughness along the lines of Jayne & St. Laurent (2001). More...
 
subroutine, public get_lowmode_loss (i, j, G, CS, mechanism, TKE_loss_sum)
 This subroutine extracts the energy lost from the propagating internal which has been summed across all angles, frequencies, and modes for a given mechanism and location. More...
 
subroutine refract (En, cn, freq, dt, G, US, NAngle, use_PPMang)
 Implements refraction on the internal waves at a single frequency. More...
 
subroutine ppm_angular_advect (En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang)
 This subroutine calculates the 1-d flux for advection in angular space using a monotonic piecewise parabolic scheme. This needs to be called from within i and j spatial loops. More...
 
subroutine propagate (En, cn, freq, dt, G, US, CS, NAngle)
 Propagates internal waves at a single frequency. More...
 
subroutine propagate_corner_spread (En, energized_wedge, NAngle, speed, dt, G, CS, LB)
 This subroutine does first-order corner advection. It was written with the hopes of smoothing out the garden sprinkler effect, but is too numerically diffusive to be of much use as of yet. It is not yet compatible with reflection schemes (BDM). More...
 
subroutine propagate_x (En, speed_x, Cgx_av, dCgx, dt, G, US, Nangle, CS, LB)
 Propagates the internal wave energy in the logical x-direction. More...
 
subroutine propagate_y (En, speed_y, Cgy_av, dCgy, dt, G, US, Nangle, CS, LB)
 Propagates the internal wave energy in the logical y-direction. More...
 
subroutine zonal_flux_en (u, h, hL, hR, uh, dt, G, US, j, ish, ieh, vol_CFL)
 Evaluates the zonal mass or volume fluxes in a layer. More...
 
subroutine merid_flux_en (v, h, hL, hR, vh, dt, G, US, J, ish, ieh, vol_CFL)
 Evaluates the meridional mass or volume fluxes in a layer. More...
 
subroutine reflect (En, NAngle, CS, G, LB)
 Reflection of the internal waves at a single frequency. More...
 
subroutine teleport (En, NAngle, CS, G, LB)
 Moves energy across lines of partial reflection to prevent reflection of energy that is supposed to get across. More...
 
subroutine correct_halo_rotation (En, test, G, NAngle)
 Rotates points in the halos where required to accommodate changes in grid orientation, such as at the tripolar fold. More...
 
subroutine ppm_reconstruction_x (h_in, h_l, h_r, G, LB, simple_2nd)
 Calculates left/right edge values for PPM reconstruction in x-direction. More...
 
subroutine ppm_reconstruction_y (h_in, h_l, h_r, G, LB, simple_2nd)
 Calculates left/right edge valus for PPM reconstruction in y-direction. More...
 
subroutine ppm_limit_pos (h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)
 Limits the left/right edge values of the PPM reconstruction to give a reconstruction that is positive-definite. Here this is reinterpreted as giving a constant thickness if the mean thickness is less than h_min, with a minimum of h_min otherwise. More...
 
subroutine, public internal_tides_init (Time, G, GV, US, param_file, diag, CS)
 This subroutine initializes the internal tides module. More...
 
subroutine, public internal_tides_end (CS)
 This subroutine deallocates the memory associated with the internal tides control structure. More...
 

Function/Subroutine Documentation

◆ correct_halo_rotation()

subroutine mom_internal_tides::correct_halo_rotation ( real, dimension(:,:,:,:,:), intent(inout)  En,
real, dimension(szi_(g),szj_(g),2), intent(in)  test,
type(ocean_grid_type), intent(in)  G,
integer, intent(in)  NAngle 
)
private

Rotates points in the halos where required to accommodate changes in grid orientation, such as at the tripolar fold.

Parameters
[in]gThe ocean's grid structure
[in,out]enThe internal gravity wave energy density as a function of space, angular orientation, frequency, and vertical mode [R Z3 T-2 ~> J m-2].
[in]testAn x-unit vector that has been passed through
[in]nangleThe number of wave orientations in the discretized wave energy spectrum.

Definition at line 1810 of file MOM_internal_tides.F90.

1810  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1811  real, dimension(:,:,:,:,:), intent(inout) :: En !< The internal gravity wave energy density as a
1812  !! function of space, angular orientation, frequency,
1813  !! and vertical mode [R Z3 T-2 ~> J m-2].
1814  real, dimension(SZI_(G),SZJ_(G),2), &
1815  intent(in) :: test !< An x-unit vector that has been passed through
1816  !! the halo updates, to enable the rotation of the
1817  !! wave energies in the halo region to be corrected.
1818  integer, intent(in) :: NAngle !< The number of wave orientations in the
1819  !! discretized wave energy spectrum.
1820  ! Local variables
1821  real, dimension(G%isd:G%ied,NAngle) :: En2d
1822  integer, dimension(G%isd:G%ied) :: a_shift
1823  integer :: i_first, i_last, a_new
1824  integer :: a, i, j, isd, ied, jsd, jed, m, fr
1825  character(len=160) :: mesg ! The text of an error message
1826  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1827 
1828  do j=jsd,jed
1829  i_first = ied+1 ; i_last = isd-1
1830  do i=isd,ied
1831  a_shift(i) = 0
1832  if (test(i,j,1) /= 1.0) then
1833  if (i<i_first) i_first = i
1834  if (i>i_last) i_last = i
1835 
1836  if (test(i,j,1) == -1.0) then ; a_shift(i) = nangle/2
1837  elseif (test(i,j,2) == 1.0) then ; a_shift(i) = -nangle/4
1838  elseif (test(i,j,2) == -1.0) then ; a_shift(i) = nangle/4
1839  else
1840  write(mesg,'("Unrecognized rotation test vector ",2ES9.2," at ",F7.2," E, ",&
1841  &F7.2," N; i,j=",2i4)') &
1842  test(i,j,1), test(i,j,2), g%GeoLonT(i,j), g%GeoLatT(i,j), i, j
1843  call mom_error(fatal, mesg)
1844  endif
1845  endif
1846  enddo
1847 
1848  if (i_first <= i_last) then
1849  ! At least one point in this row needs to be rotated.
1850  do m=1,size(en,5) ; do fr=1,size(en,4)
1851  do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
1852  a_new = a + a_shift(i)
1853  if (a_new < 1) a_new = a_new + nangle
1854  if (a_new > nangle) a_new = a_new - nangle
1855  en2d(i,a_new) = en(i,j,a,fr,m)
1856  endif ; enddo ; enddo
1857  do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
1858  en(i,j,a,fr,m) = en2d(i,a)
1859  endif ; enddo ; enddo
1860  enddo ; enddo
1861  endif
1862  enddo

References mom_error_handler::mom_error().

Referenced by propagate_int_tide().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_lowmode_loss()

subroutine, public mom_internal_tides::get_lowmode_loss ( integer, intent(in)  i,
integer, intent(in)  j,
type(ocean_grid_type), intent(in)  G,
type(int_tide_cs), pointer  CS,
character(len=*), intent(in)  mechanism,
real, intent(out)  TKE_loss_sum 
)

This subroutine extracts the energy lost from the propagating internal which has been summed across all angles, frequencies, and modes for a given mechanism and location.

It can be called from another module to get values from this module's (private) CS.

Parameters
[in]iThe i-index of the value to be reported.
[in]jThe j-index of the value to be reported.
[in]gThe ocean's grid structure
csThe control structure returned by a previous call to int_tide_init.
[in]mechanismThe named mechanism of loss to return
[out]tke_loss_sumTotal energy loss rate due to specified mechanism [R Z3 T-3 ~> W m-2].

Definition at line 728 of file MOM_internal_tides.F90.

728  integer, intent(in) :: i !< The i-index of the value to be reported.
729  integer, intent(in) :: j !< The j-index of the value to be reported.
730  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
731  type(int_tide_CS), pointer :: CS !< The control structure returned by a
732  !! previous call to int_tide_init.
733  character(len=*), intent(in) :: mechanism !< The named mechanism of loss to return
734  real, intent(out) :: TKE_loss_sum !< Total energy loss rate due to specified
735  !! mechanism [R Z3 T-3 ~> W m-2].
736 
737  if (mechanism == 'LeakDrag') tke_loss_sum = cs%tot_leak_loss(i,j) ! not used for mixing yet
738  if (mechanism == 'QuadDrag') tke_loss_sum = cs%tot_quad_loss(i,j) ! not used for mixing yet
739  if (mechanism == 'WaveDrag') tke_loss_sum = cs%tot_itidal_loss(i,j) ! currently used for mixing
740  if (mechanism == 'Froude') tke_loss_sum = cs%tot_Froude_loss(i,j) ! not used for mixing yet
741 

◆ internal_tides_end()

subroutine, public mom_internal_tides::internal_tides_end ( type(int_tide_cs), pointer  CS)

This subroutine deallocates the memory associated with the internal tides control structure.

Parameters
csA pointer to the control structure returned by a previous call to internal_tides_init, it will be deallocated here.

Definition at line 2532 of file MOM_internal_tides.F90.

2532  type(int_tide_CS), pointer :: CS !< A pointer to the control structure returned by a previous
2533  !! call to internal_tides_init, it will be deallocated here.
2534 
2535  if (associated(cs)) then
2536  if (associated(cs%En)) deallocate(cs%En)
2537  if (allocated(cs%frequency)) deallocate(cs%frequency)
2538  if (allocated(cs%id_En_mode)) deallocate(cs%id_En_mode)
2539  if (allocated(cs%id_Ub_mode)) deallocate(cs%id_Ub_mode)
2540  if (allocated(cs%id_cp_mode)) deallocate(cs%id_cp_mode)
2541  deallocate(cs)
2542  endif
2543  cs => null()

◆ internal_tides_init()

subroutine, public mom_internal_tides::internal_tides_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(inout)  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(in), target  diag,
type(int_tide_cs), pointer  CS 
)

This subroutine initializes the internal tides module.

Parameters
[in]timeThe current model time.
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]param_fileA structure to parse for run-time parameters.
[in]diagA structure that is used to regulate diagnostic output.
csA pointer that is set to point to the control structure for this module.

Definition at line 2100 of file MOM_internal_tides.F90.

2100  type(time_type), target, intent(in) :: Time !< The current model time.
2101  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2102  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2103  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2104  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2105  !! parameters.
2106  type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
2107  !! diagnostic output.
2108  type(int_tide_CS),pointer :: CS !< A pointer that is set to point to the control
2109  !! structure for this module.
2110  ! Local variables
2111  real :: Angle_size ! size of wedges, rad
2112  real, allocatable :: angles(:) ! orientations of wedge centers, rad
2113  real, allocatable, dimension(:,:) :: h2 ! topographic roughness scale, m^2
2114  real :: kappa_itides, kappa_h2_factor
2115  ! characteristic topographic wave number
2116  ! and a scaling factor
2117  real, allocatable :: ridge_temp(:,:)
2118  ! array for temporary storage of flags
2119  ! of cells with double-reflecting ridges
2120  logical :: use_int_tides, use_temperature
2121  real :: period_1 ! The period of the gravest modeled mode [T ~> s]
2122  integer :: num_angle, num_freq, num_mode, m, fr
2123  integer :: isd, ied, jsd, jed, a, id_ang, i, j
2124  type(axes_grp) :: axes_ang
2125  ! This include declares and sets the variable "version".
2126 #include "version_variable.h"
2127  character(len=40) :: mdl = "MOM_internal_tides" ! This module's name.
2128  character(len=16), dimension(8) :: freq_name
2129  character(len=40) :: var_name
2130  character(len=160) :: var_descript
2131  character(len=200) :: filename
2132  character(len=200) :: refl_angle_file, land_mask_file
2133  character(len=200) :: refl_pref_file, refl_dbl_file
2134  character(len=200) :: dy_Cu_file, dx_Cv_file
2135  character(len=200) :: h2_file
2136 
2137  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2138 
2139  if (associated(cs)) then
2140  call mom_error(warning, "internal_tides_init called "//&
2141  "with an associated control structure.")
2142  return
2143  else
2144  allocate(cs)
2145  endif
2146 
2147  use_int_tides = .false.
2148  call read_param(param_file, "INTERNAL_TIDES", use_int_tides)
2149  cs%do_int_tides = use_int_tides
2150  if (.not.use_int_tides) return
2151 
2152  use_temperature = .true.
2153  call read_param(param_file, "ENABLE_THERMODYNAMICS", use_temperature)
2154  if (.not.use_temperature) call mom_error(fatal, &
2155  "register_int_tide_restarts: internal_tides only works with "//&
2156  "ENABLE_THERMODYNAMICS defined.")
2157 
2158  ! Set number of frequencies, angles, and modes to consider
2159  num_freq = 1 ; num_angle = 24 ; num_mode = 1
2160  call read_param(param_file, "INTERNAL_TIDE_FREQS", num_freq)
2161  call read_param(param_file, "INTERNAL_TIDE_ANGLES", num_angle)
2162  call read_param(param_file, "INTERNAL_TIDE_MODES", num_mode)
2163  if (.not.((num_freq > 0) .and. (num_angle > 0) .and. (num_mode > 0))) return
2164  cs%nFreq = num_freq ; cs%nAngle = num_angle ; cs%nMode = num_mode
2165 
2166  ! Allocate energy density array
2167  allocate(cs%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode))
2168  cs%En(:,:,:,:,:) = 0.0
2169 
2170  ! Allocate phase speed array
2171  allocate(cs%cp(isd:ied, jsd:jed, num_freq, num_mode))
2172  cs%cp(:,:,:,:) = 0.0
2173 
2174  ! Allocate and populate frequency array (each a multiple of first for now)
2175  allocate(cs%frequency(num_freq))
2176  call get_param(param_file, mdl, "FIRST_MODE_PERIOD", period_1, units="s", scale=us%s_to_T)
2177  do fr=1,num_freq
2178  cs%frequency(fr) = (8.0*atan(1.0) * (real(fr)) / period_1) ! ADDED BDM
2179  enddo
2180 
2181  ! Read all relevant parameters and write them to the model log.
2182 
2183  cs%Time => time ! direct a pointer to the current model time target
2184 
2185  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
2186  cs%inputdir = slasher(cs%inputdir)
2187 
2188  call log_version(param_file, mdl, version, "")
2189 
2190  call get_param(param_file, mdl, "INTERNAL_TIDE_FREQS", num_freq, &
2191  "The number of distinct internal tide frequency bands "//&
2192  "that will be calculated.", default=1)
2193  call get_param(param_file, mdl, "INTERNAL_TIDE_MODES", num_mode, &
2194  "The number of distinct internal tide modes "//&
2195  "that will be calculated.", default=1)
2196  call get_param(param_file, mdl, "INTERNAL_TIDE_ANGLES", num_angle, &
2197  "The number of angular resolution bands for the internal "//&
2198  "tide calculations.", default=24)
2199 
2200  if (use_int_tides) then
2201  if ((num_freq <= 0) .and. (num_mode <= 0) .and. (num_angle <= 0)) then
2202  call mom_error(warning, "Internal tides were enabled, but the number "//&
2203  "of requested frequencies, modes and angles were not all positive.")
2204  return
2205  endif
2206  else
2207  if ((num_freq > 0) .and. (num_mode > 0) .and. (num_angle > 0)) then
2208  call mom_error(warning, "Internal tides were not enabled, even though "//&
2209  "the number of requested frequencies, modes and angles were all "//&
2210  "positive.")
2211  return
2212  endif
2213  endif
2214 
2215  if (cs%NFreq /= num_freq) call mom_error(fatal, "Internal_tides_init: "//&
2216  "Inconsistent number of frequencies.")
2217  if (cs%NAngle /= num_angle) call mom_error(fatal, "Internal_tides_init: "//&
2218  "Inconsistent number of angles.")
2219  if (cs%NMode /= num_mode) call mom_error(fatal, "Internal_tides_init: "//&
2220  "Inconsistent number of modes.")
2221  if (4*(num_angle/4) /= num_angle) call mom_error(fatal, &
2222  "Internal_tides_init: INTERNAL_TIDE_ANGLES must be a multiple of 4.")
2223 
2224  cs%diag => diag
2225 
2226  call get_param(param_file, mdl, "INTERNAL_TIDE_DECAY_RATE", cs%decay_rate, &
2227  "The rate at which internal tide energy is lost to the "//&
2228  "interior ocean internal wave field.", &
2229  units="s-1", default=0.0, scale=us%T_to_s)
2230  call get_param(param_file, mdl, "INTERNAL_TIDE_VOLUME_BASED_CFL", cs%vol_CFL, &
2231  "If true, use the ratio of the open face lengths to the "//&
2232  "tracer cell areas when estimating CFL numbers in the "//&
2233  "internal tide code.", default=.false.)
2234  call get_param(param_file, mdl, "INTERNAL_TIDE_CORNER_ADVECT", cs%corner_adv, &
2235  "If true, internal tide ray-tracing advection uses a "//&
2236  "corner-advection scheme rather than PPM.", default=.false.)
2237  call get_param(param_file, mdl, "INTERNAL_TIDE_SIMPLE_2ND_PPM", cs%simple_2nd, &
2238  "If true, CONTINUITY_PPM uses a simple 2nd order "//&
2239  "(arithmetic mean) interpolation of the edge values. "//&
2240  "This may give better PV conservation properties. While "//&
2241  "it formally reduces the accuracy of the continuity "//&
2242  "solver itself in the strongly advective limit, it does "//&
2243  "not reduce the overall order of accuracy of the dynamic "//&
2244  "core.", default=.false.)
2245  call get_param(param_file, mdl, "INTERNAL_TIDE_UPWIND_1ST", cs%upwind_1st, &
2246  "If true, the internal tide ray-tracing advection uses "//&
2247  "1st-order upwind advection. This scheme is highly "//&
2248  "continuity solver. This scheme is highly "//&
2249  "diffusive but may be useful for debugging.", default=.false.)
2250  call get_param(param_file, mdl, "INTERNAL_TIDE_BACKGROUND_DRAG", &
2251  cs%apply_background_drag, "If true, the internal tide "//&
2252  "ray-tracing advection uses a background drag term as a sink.",&
2253  default=.false.)
2254  call get_param(param_file, mdl, "INTERNAL_TIDE_QUAD_DRAG", cs%apply_bottom_drag, &
2255  "If true, the internal tide ray-tracing advection uses "//&
2256  "a quadratic bottom drag term as a sink.", default=.false.)
2257  call get_param(param_file, mdl, "INTERNAL_TIDE_WAVE_DRAG", cs%apply_wave_drag, &
2258  "If true, apply scattering due to small-scale roughness as a sink.", &
2259  default=.false.)
2260  call get_param(param_file, mdl, "INTERNAL_TIDE_FROUDE_DRAG", cs%apply_Froude_drag, &
2261  "If true, apply wave breaking as a sink.", &
2262  default=.false.)
2263  call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
2264  "CDRAG is the drag coefficient relating the magnitude of "//&
2265  "the velocity field to the bottom stress.", units="nondim", &
2266  default=0.003)
2267  call get_param(param_file, mdl, "INTERNAL_TIDE_ENERGIZED_ANGLE", cs%energized_angle, &
2268  "If positive, only one angular band of the internal tides "//&
2269  "gets all of the energy. (This is for debugging.)", default=-1)
2270  call get_param(param_file, mdl, "USE_PPM_ANGULAR", cs%use_PPMang, &
2271  "If true, use PPM for advection of energy in angular space.", &
2272  default=.false.)
2273  call get_param(param_file, mdl, "GAMMA_ITIDES", cs%q_itides, &
2274  "The fraction of the internal tidal energy that is "//&
2275  "dissipated locally with INT_TIDE_DISSIPATION. "//&
2276  "THIS NAME COULD BE BETTER.", &
2277  units="nondim", default=0.3333)
2278  call get_param(param_file, mdl, "KAPPA_ITIDES", kappa_itides, &
2279  "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
2280  "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
2281  units="m-1", default=8.e-4*atan(1.0), scale=us%L_to_m)
2282  call get_param(param_file, mdl, "KAPPA_H2_FACTOR", kappa_h2_factor, &
2283  "A scaling factor for the roughness amplitude with n"//&
2284  "INT_TIDE_DISSIPATION.", units="nondim", default=1.0)
2285 
2286  ! Allocate various arrays needed for loss rates
2287  allocate(h2(isd:ied,jsd:jed)) ; h2(:,:) = 0.0
2288  allocate(cs%TKE_itidal_loss_fixed(isd:ied,jsd:jed))
2289  cs%TKE_itidal_loss_fixed = 0.0
2290  allocate(cs%TKE_leak_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2291  cs%TKE_leak_loss(:,:,:,:,:) = 0.0
2292  allocate(cs%TKE_quad_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2293  cs%TKE_quad_loss(:,:,:,:,:) = 0.0
2294  allocate(cs%TKE_itidal_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2295  cs%TKE_itidal_loss(:,:,:,:,:) = 0.0
2296  allocate(cs%TKE_Froude_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2297  cs%TKE_Froude_loss(:,:,:,:,:) = 0.0
2298  allocate(cs%tot_leak_loss(isd:ied,jsd:jed)) ; cs%tot_leak_loss(:,:) = 0.0
2299  allocate(cs%tot_quad_loss(isd:ied,jsd:jed) ) ; cs%tot_quad_loss(:,:) = 0.0
2300  allocate(cs%tot_itidal_loss(isd:ied,jsd:jed)) ; cs%tot_itidal_loss(:,:) = 0.0
2301  allocate(cs%tot_Froude_loss(isd:ied,jsd:jed)) ; cs%tot_Froude_loss(:,:) = 0.0
2302 
2303  ! Compute the fixed part of the bottom drag loss from baroclinic modes
2304  call get_param(param_file, mdl, "H2_FILE", h2_file, &
2305  "The path to the file containing the sub-grid-scale "//&
2306  "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
2307  fail_if_missing=.true.)
2308  filename = trim(cs%inputdir) // trim(h2_file)
2309  call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename)
2310  call mom_read_data(filename, 'h2', h2, g%domain, timelevel=1, scale=us%m_to_Z)
2311  do j=g%jsc,g%jec ; do i=g%isc,g%iec
2312  ! Restrict rms topo to 10 percent of column depth.
2313  h2(i,j) = min(0.01*(g%bathyT(i,j))**2, h2(i,j))
2314  ! Compute the fixed part; units are [R L-2 Z3 ~> kg m-2] here
2315  ! will be multiplied by N and the squared near-bottom velocity to get into [R Z3 T-3 ~> W m-2]
2316  cs%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*gv%Rho0 * us%L_to_Z*kappa_itides * h2(i,j)
2317  enddo ; enddo
2318 
2319  deallocate(h2)
2320 
2321  ! Read in prescribed coast/ridge/shelf angles from file
2322  call get_param(param_file, mdl, "REFL_ANGLE_FILE", refl_angle_file, &
2323  "The path to the file containing the local angle of "//&
2324  "the coastline/ridge/shelf with respect to the equator.", &
2325  fail_if_missing=.false.)
2326  filename = trim(cs%inputdir) // trim(refl_angle_file)
2327  call log_param(param_file, mdl, "INPUTDIR/REFL_ANGLE_FILE", filename)
2328  allocate(cs%refl_angle(isd:ied,jsd:jed)) ; cs%refl_angle(:,:) = cs%nullangle
2329  call mom_read_data(filename, 'refl_angle', cs%refl_angle, &
2330  g%domain, timelevel=1)
2331  ! replace NANs with null value
2332  do j=g%jsc,g%jec ; do i=g%isc,g%iec
2333  if (is_nan(cs%refl_angle(i,j))) cs%refl_angle(i,j) = cs%nullangle
2334  enddo ; enddo
2335  call pass_var(cs%refl_angle,g%domain)
2336 
2337  ! Read in prescribed partial reflection coefficients from file
2338  call get_param(param_file, mdl, "REFL_PREF_FILE", refl_pref_file, &
2339  "The path to the file containing the reflection coefficients.", &
2340  fail_if_missing=.false.)
2341  filename = trim(cs%inputdir) // trim(refl_pref_file)
2342  call log_param(param_file, mdl, "INPUTDIR/REFL_PREF_FILE", filename)
2343  allocate(cs%refl_pref(isd:ied,jsd:jed)) ; cs%refl_pref(:,:) = 1.0
2344  call mom_read_data(filename, 'refl_pref', cs%refl_pref, g%domain, timelevel=1)
2345  !CS%refl_pref = CS%refl_pref*1 ! adjust partial reflection if desired
2346  call pass_var(cs%refl_pref,g%domain)
2347 
2348  ! Tag reflection cells with partial reflection (done here for speed)
2349  allocate(cs%refl_pref_logical(isd:ied,jsd:jed)) ; cs%refl_pref_logical(:,:) = .false.
2350  do j=jsd,jed
2351  do i=isd,ied
2352  ! flag cells with partial reflection
2353  if (cs%refl_angle(i,j) /= cs%nullangle .and. &
2354  cs%refl_pref(i,j) < 1.0 .and. cs%refl_pref(i,j) > 0.0) then
2355  cs%refl_pref_logical(i,j) = .true.
2356  endif
2357  enddo
2358  enddo
2359 
2360  ! Read in double-reflective (ridge) tags from file
2361  call get_param(param_file, mdl, "REFL_DBL_FILE", refl_dbl_file, &
2362  "The path to the file containing the double-reflective ridge tags.", &
2363  fail_if_missing=.false.)
2364  filename = trim(cs%inputdir) // trim(refl_dbl_file)
2365  call log_param(param_file, mdl, "INPUTDIR/REFL_DBL_FILE", filename)
2366  allocate(ridge_temp(isd:ied,jsd:jed)) ; ridge_temp(:,:) = 0.0
2367  call mom_read_data(filename, 'refl_dbl', ridge_temp, g%domain, timelevel=1)
2368  call pass_var(ridge_temp,g%domain)
2369  allocate(cs%refl_dbl(isd:ied,jsd:jed)) ; cs%refl_dbl(:,:) = .false.
2370  do i=isd,ied; do j=jsd,jed
2371  if (ridge_temp(i,j) == 1) then; cs%refl_dbl(i,j) = .true.
2372  else ; cs%refl_dbl(i,j) = .false. ; endif
2373  enddo ; enddo
2374 
2375  ! Read in prescribed land mask from file (if overwriting -BDM).
2376  ! This should be done in MOM_initialize_topography subroutine
2377  ! defined in MOM_fixed_initialization.F90 (BDM)
2378  !call get_param(param_file, mdl, "LAND_MASK_FILE", land_mask_file, &
2379  ! "The path to the file containing the land mask.", &
2380  ! fail_if_missing=.false.)
2381  !filename = trim(CS%inputdir) // trim(land_mask_file)
2382  !call log_param(param_file, mdl, "INPUTDIR/LAND_MASK_FILE", filename)
2383  !G%mask2dCu(:,:) = 1 ; G%mask2dCv(:,:) = 1 ; G%mask2dT(:,:) = 1
2384  !call MOM_read_data(filename, 'land_mask', G%mask2dCu, G%domain, timelevel=1)
2385  !call MOM_read_data(filename, 'land_mask', G%mask2dCv, G%domain, timelevel=1)
2386  !call MOM_read_data(filename, 'land_mask', G%mask2dT, G%domain, timelevel=1)
2387  !call pass_vector(G%mask2dCu, G%mask2dCv, G%domain, To_All+Scalar_Pair, CGRID_NE)
2388  !call pass_var(G%mask2dT,G%domain)
2389 
2390  ! Read in prescribed partial east face blockages from file (if overwriting -BDM)
2391  !call get_param(param_file, mdl, "dy_Cu_FILE", dy_Cu_file, &
2392  ! "The path to the file containing the east face blockages.", &
2393  ! fail_if_missing=.false.)
2394  !filename = trim(CS%inputdir) // trim(dy_Cu_file)
2395  !call log_param(param_file, mdl, "INPUTDIR/dy_Cu_FILE", filename)
2396  !G%dy_Cu(:,:) = 0.0
2397  !call MOM_read_data(filename, 'dy_Cu', G%dy_Cu, G%domain, timelevel=1, scale=US%m_to_L)
2398 
2399  ! Read in prescribed partial north face blockages from file (if overwriting -BDM)
2400  !call get_param(param_file, mdl, "dx_Cv_FILE", dx_Cv_file, &
2401  ! "The path to the file containing the north face blockages.", &
2402  ! fail_if_missing=.false.)
2403  !filename = trim(CS%inputdir) // trim(dx_Cv_file)
2404  !call log_param(param_file, mdl, "INPUTDIR/dx_Cv_FILE", filename)
2405  !G%dx_Cv(:,:) = 0.0
2406  !call MOM_read_data(filename, 'dx_Cv', G%dx_Cv, G%domain, timelevel=1, scale=US%m_to_L)
2407  !call pass_vector(G%dy_Cu, G%dx_Cv, G%domain, To_All+Scalar_Pair, CGRID_NE)
2408 
2409  ! Register maps of reflection parameters
2410  cs%id_refl_ang = register_diag_field('ocean_model', 'refl_angle', diag%axesT1, &
2411  time, 'Local angle of coastline/ridge/shelf with respect to equator', 'rad')
2412  cs%id_refl_pref = register_diag_field('ocean_model', 'refl_pref', diag%axesT1, &
2413  time, 'Partial reflection coefficients', '')
2414  cs%id_dx_Cv = register_diag_field('ocean_model', 'dx_Cv', diag%axesT1, &
2415  time, 'North face unblocked width', 'm', conversion=us%L_to_m)
2416  cs%id_dy_Cu = register_diag_field('ocean_model', 'dy_Cu', diag%axesT1, &
2417  time, 'East face unblocked width', 'm', conversion=us%L_to_m)
2418  cs%id_land_mask = register_diag_field('ocean_model', 'land_mask', diag%axesT1, &
2419  time, 'Land mask', 'logical') ! used if overriding (BDM)
2420  ! Output reflection parameters as diags here (not needed every timestep)
2421  if (cs%id_refl_ang > 0) call post_data(cs%id_refl_ang, cs%refl_angle, cs%diag)
2422  if (cs%id_refl_pref > 0) call post_data(cs%id_refl_pref, cs%refl_pref, cs%diag)
2423  if (cs%id_dx_Cv > 0) call post_data(cs%id_dx_Cv, g%dx_Cv, cs%diag)
2424  if (cs%id_dy_Cu > 0) call post_data(cs%id_dy_Cu, g%dy_Cu, cs%diag)
2425  if (cs%id_land_mask > 0) call post_data(cs%id_land_mask, g%mask2dT, cs%diag)
2426 
2427  ! Register 2-D energy density (summed over angles, freq, modes)
2428  cs%id_tot_En = register_diag_field('ocean_model', 'ITide_tot_En', diag%axesT1, &
2429  time, 'Internal tide total energy density', &
2430  'J m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**2)
2431  ! Register 2-D drag scale used for quadratic bottom drag
2432  cs%id_itide_drag = register_diag_field('ocean_model', 'ITide_drag', diag%axesT1, &
2433  time, 'Interior and bottom drag internal tide decay timescale', 's-1', conversion=us%s_to_T)
2434  !Register 2-D energy input into internal tides
2435  cs%id_TKE_itidal_input = register_diag_field('ocean_model', 'TKE_itidal_input', diag%axesT1, &
2436  time, 'Conversion from barotropic to baroclinic tide, '//&
2437  'a fraction of which goes into rays', &
2438  'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2439  ! Register 2-D energy losses (summed over angles, freq, modes)
2440  cs%id_tot_leak_loss = register_diag_field('ocean_model', 'ITide_tot_leak_loss', diag%axesT1, &
2441  time, 'Internal tide energy loss to background drag', &
2442  'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2443  cs%id_tot_quad_loss = register_diag_field('ocean_model', 'ITide_tot_quad_loss', diag%axesT1, &
2444  time, 'Internal tide energy loss to bottom drag', &
2445  'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2446  cs%id_tot_itidal_loss = register_diag_field('ocean_model', 'ITide_tot_itidal_loss', diag%axesT1, &
2447  time, 'Internal tide energy loss to wave drag', &
2448  'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2449  cs%id_tot_Froude_loss = register_diag_field('ocean_model', 'ITide_tot_Froude_loss', diag%axesT1, &
2450  time, 'Internal tide energy loss to wave breaking', &
2451  'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2452  cs%id_tot_allprocesses_loss = register_diag_field('ocean_model', 'ITide_tot_allprocesses_loss', diag%axesT1, &
2453  time, 'Internal tide energy loss summed over all processes', &
2454  'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2455 
2456  allocate(cs%id_En_mode(cs%nFreq,cs%nMode)) ; cs%id_En_mode(:,:) = -1
2457  allocate(cs%id_En_ang_mode(cs%nFreq,cs%nMode)) ; cs%id_En_ang_mode(:,:) = -1
2458  allocate(cs%id_itidal_loss_mode(cs%nFreq,cs%nMode)) ; cs%id_itidal_loss_mode(:,:) = -1
2459  allocate(cs%id_allprocesses_loss_mode(cs%nFreq,cs%nMode)) ; cs%id_allprocesses_loss_mode(:,:) = -1
2460  allocate(cs%id_itidal_loss_ang_mode(cs%nFreq,cs%nMode)) ; cs%id_itidal_loss_ang_mode(:,:) = -1
2461  allocate(cs%id_Ub_mode(cs%nFreq,cs%nMode)) ; cs%id_Ub_mode(:,:) = -1
2462  allocate(cs%id_cp_mode(cs%nFreq,cs%nMode)) ; cs%id_cp_mode(:,:) = -1
2463 
2464  allocate(angles(cs%NAngle)) ; angles(:) = 0.0
2465  angle_size = (8.0*atan(1.0)) / (real(num_angle))
2466  do a=1,num_angle ; angles(a) = (real(a) - 1) * angle_size ; enddo
2467 
2468  id_ang = diag_axis_init("angle", angles, "Radians", "N", "Angular Orienation of Fluxes")
2469  call define_axes_group(diag, (/ diag%axesT1%handles(1), diag%axesT1%handles(2), id_ang /), axes_ang)
2470 
2471  do fr=1,cs%nFreq ; write(freq_name(fr), '("freq",i1)') fr ; enddo
2472  do m=1,cs%nMode ; do fr=1,cs%nFreq
2473  ! Register 2-D energy density (summed over angles) for each freq and mode
2474  write(var_name, '("Itide_En_freq",i1,"_mode",i1)') fr, m
2475  write(var_descript, '("Internal tide energy density in frequency ",i1," mode ",i1)') fr, m
2476  cs%id_En_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2477  diag%axesT1, time, var_descript, 'J m-2', conversion=us%R_to_kg_m3*us%Z_to_m**2*us%s_to_T**3)
2478  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2479 
2480  ! Register 3-D (i,j,a) energy density for each freq and mode
2481  write(var_name, '("Itide_En_ang_freq",i1,"_mode",i1)') fr, m
2482  write(var_descript, '("Internal tide angular energy density in frequency ",i1," mode ",i1)') fr, m
2483  cs%id_En_ang_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2484  axes_ang, time, var_descript, 'J m-2 band-1', conversion=us%R_to_kg_m3*us%Z_to_m**2*us%s_to_T**3)
2485  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2486 
2487  ! Register 2-D energy loss (summed over angles) for each freq and mode
2488  ! wave-drag only
2489  write(var_name, '("Itide_wavedrag_loss_freq",i1,"_mode",i1)') fr, m
2490  write(var_descript, '("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
2491  cs%id_itidal_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2492  diag%axesT1, time, var_descript, 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2493  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2494  ! all loss processes
2495  write(var_name, '("Itide_allprocesses_loss_freq",i1,"_mode",i1)') fr, m
2496  write(var_descript, '("Internal tide energy loss due to all processes from frequency ",i1," mode ",i1)') fr, m
2497  cs%id_allprocesses_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2498  diag%axesT1, time, var_descript, 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2499  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2500 
2501  ! Register 3-D (i,j,a) energy loss for each freq and mode
2502  ! wave-drag only
2503  write(var_name, '("Itide_wavedrag_loss_ang_freq",i1,"_mode",i1)') fr, m
2504  write(var_descript, '("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
2505  cs%id_itidal_loss_ang_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2506  axes_ang, time, var_descript, 'W m-2 band-1', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
2507  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2508 
2509  ! Register 2-D period-averaged near-bottom horizonal velocity for each freq and mode
2510  write(var_name, '("Itide_Ub_freq",i1,"_mode",i1)') fr, m
2511  write(var_descript, '("Near-bottom horizonal velocity for frequency ",i1," mode ",i1)') fr, m
2512  cs%id_Ub_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2513  diag%axesT1, time, var_descript, 'm s-1', conversion=us%L_T_to_m_s)
2514  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2515 
2516  ! Register 2-D horizonal phase velocity for each freq and mode
2517  write(var_name, '("Itide_cp_freq",i1,"_mode",i1)') fr, m
2518  write(var_descript, '("Horizonal phase velocity for frequency ",i1," mode ",i1)') fr, m
2519  cs%id_cp_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2520  diag%axesT1, time, var_descript, 'm s-1', conversion=us%L_T_to_m_s)
2521  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2522 
2523  enddo ; enddo
2524 
2525  ! Initialize wave_structure (not sure if this should be here - BDM)
2526  call wave_structure_init(time, g, param_file, diag, cs%wave_structure_CSp)
2527 

References mom_diag_mediator::define_axes_group(), mom_error_handler::mom_error(), mom_error_handler::mom_mesg(), and mom_wave_structure::wave_structure_init().

Referenced by mom_diabatic_driver::diabatic_driver_init().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ itidal_lowmode_loss()

subroutine mom_internal_tides::itidal_lowmode_loss ( type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(int_tide_cs), pointer  CS,
real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(in)  Nb,
real, dimension(g%isd:g%ied,g%jsd:g%jed,cs%nfreq,cs%nmode), intent(inout)  Ub,
real, dimension(g%isd:g%ied,g%jsd:g%jed,cs%nangle,cs%nfreq,cs%nmode), intent(inout)  En,
real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(in)  TKE_loss_fixed,
real, dimension(g%isd:g%ied,g%jsd:g%jed,cs%nangle,cs%nfreq,cs%nmode), intent(out)  TKE_loss,
real, intent(in)  dt,
logical, intent(in), optional  full_halos 
)
private

Calculates the energy lost from the propagating internal tide due to scattering over small-scale roughness along the lines of Jayne & St. Laurent (2001).

Parameters
[in]gThe ocean's grid structure.
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to int_tide_init.
[in]nbNear-bottom stratification [T-1 ~> s-1].
[in,out]ubRMS (over one period) near-bottom horizontal
[in]tke_loss_fixedFixed part of energy loss [R L-2 Z3 ~> kg m-2]
[in,out]enEnergy density of the internal waves [R Z3 T-2 ~> J m-2].
[out]tke_lossEnergy loss rate [R Z3 T-3 ~> W m-2]
[in]dtTime increment [T ~> s].
[in]full_halosIf true, do the calculation over the entirecomputational domain.

Definition at line 636 of file MOM_internal_tides.F90.

636  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
637  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
638  type(int_tide_CS), pointer :: CS !< The control structure returned by a
639  !! previous call to int_tide_init.
640  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
641  intent(in) :: Nb !< Near-bottom stratification [T-1 ~> s-1].
642  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%nFreq,CS%nMode), &
643  intent(inout) :: Ub !< RMS (over one period) near-bottom horizontal
644  !! mode velocity [L T-1 ~> m s-1].
645  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
646  intent(in) :: TKE_loss_fixed !< Fixed part of energy loss [R L-2 Z3 ~> kg m-2]
647  !! (rho*kappa*h^2).
648  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
649  intent(inout) :: En !< Energy density of the internal waves [R Z3 T-2 ~> J m-2].
650  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
651  intent(out) :: TKE_loss !< Energy loss rate [R Z3 T-3 ~> W m-2]
652  !! (q*rho*kappa*h^2*N*U^2).
653  real, intent(in) :: dt !< Time increment [T ~> s].
654  logical,optional, intent(in) :: full_halos !< If true, do the calculation over the
655  !! entirecomputational domain.
656  ! Local variables
657  integer :: j,i,m,fr,a, is, ie, js, je
658  real :: En_tot ! energy for a given mode, frequency, and point summed over angles [R Z3 T-2 ~> J m-2]
659  real :: TKE_loss_tot ! dissipation for a given mode, frequency, and point summed over angles [R Z3 T-3 ~> W m-2]
660  real :: TKE_sum_check ! temporary for check summing
661  real :: frac_per_sector ! fraction of energy in each wedge
662  real :: q_itides ! fraction of energy actually lost to mixing (remainder, 1-q, is
663  ! assumed to stay in propagating mode for now - BDM)
664  real :: loss_rate ! approximate loss rate for implicit calc [T-1 ~> s-1]
665  real :: En_negl ! negilibly small number to prevent division by zero
666 
667  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
668 
669  q_itides = cs%q_itides
670  en_negl = 1e-30*us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**2
671 
672  if (present(full_halos)) then ; if (full_halos) then
673  is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
674  endif ; endif
675 
676  do j=js,je ; do i=is,ie ; do m=1,cs%nMode ; do fr=1,cs%nFreq
677 
678  ! Sum energy across angles
679  en_tot = 0.0
680  do a=1,cs%nAngle
681  en_tot = en_tot + en(i,j,a,fr,m)
682  enddo
683 
684  ! Calculate TKE loss rate; units of [R Z3 T-3 ~> W m-2] here.
685  tke_loss_tot = q_itides * tke_loss_fixed(i,j) * nb(i,j) * ub(i,j,fr,m)**2
686 
687  ! Update energy remaining (this is a pseudo implicit calc)
688  ! (E(t+1)-E(t))/dt = -TKE_loss(E(t+1)/E(t)), which goes to zero as E(t+1) goes to zero
689  if (en_tot > 0.0) then
690  do a=1,cs%nAngle
691  frac_per_sector = en(i,j,a,fr,m)/en_tot
692  tke_loss(i,j,a,fr,m) = frac_per_sector*tke_loss_tot ! Wm-2
693  loss_rate = tke_loss(i,j,a,fr,m) / (en(i,j,a,fr,m) + en_negl) ! [T-1 ~> s-1]
694  en(i,j,a,fr,m) = en(i,j,a,fr,m) / (1.0 + dt*loss_rate)
695  enddo
696  else
697  ! no loss if no energy
698  tke_loss(i,j,:,fr,m) = 0.0
699  endif
700 
701  ! Update energy remaining (this is the old explicit calc)
702  !if (En_tot > 0.0) then
703  ! do a=1,CS%nAngle
704  ! frac_per_sector = En(i,j,a,fr,m)/En_tot
705  ! TKE_loss(i,j,a,fr,m) = frac_per_sector*TKE_loss_tot
706  ! if (TKE_loss(i,j,a,fr,m)*dt <= En(i,j,a,fr,m))then
707  ! En(i,j,a,fr,m) = En(i,j,a,fr,m) - TKE_loss(i,j,a,fr,m)*dt
708  ! else
709  ! call MOM_error(WARNING, "itidal_lowmode_loss: energy loss greater than avalable, "// &
710  ! " setting En to zero.", all_print=.true.)
711  ! En(i,j,a,fr,m) = 0.0
712  ! endif
713  ! enddo
714  !else
715  ! ! no loss if no energy
716  ! TKE_loss(i,j,:,fr,m) = 0.0
717  !endif
718 
719  enddo ; enddo ; enddo ; enddo
720 

Referenced by propagate_int_tide().

Here is the caller graph for this function:

◆ merid_flux_en()

subroutine mom_internal_tides::merid_flux_en ( real, dimension( g %isd: g %ied), intent(in)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  hL,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  hR,
real, dimension( g %isd: g %ied), intent(inout)  vh,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
integer, intent(in)  J,
integer, intent(in)  ish,
integer, intent(in)  ieh,
logical, intent(in)  vol_CFL 
)
private

Evaluates the meridional mass or volume fluxes in a layer.

Parameters
[in]gThe ocean's grid structure.
[in]vThe meridional velocity [L T-1 ~> m s-1].
[in]hEnergy density used to calculate the fluxes [R Z3 T-2 ~> J m-2].
[in]hlLeft- Energy densities in the reconstruction [R Z3 T-2 ~> J m-2].
[in]hrRight- Energy densities in the reconstruction [R Z3 T-2 ~> J m-2].
[in,out]vhThe meridional energy transport [R Z3 L2 T-3 ~> J s-1].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
[in]jThe j-index to work on.
[in]ishThe start i-index range to work on.
[in]iehThe end i-index range to work on.
[in]vol_cflIf true, rescale the ratio of face areas to the cell areas when estimating the CFL number.

Definition at line 1550 of file MOM_internal_tides.F90.

1550  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1551  real, dimension(SZI_(G)), intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
1552  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h !< Energy density used to calculate the
1553  !! fluxes [R Z3 T-2 ~> J m-2].
1554  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hL !< Left- Energy densities in the
1555  !! reconstruction [R Z3 T-2 ~> J m-2].
1556  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hR !< Right- Energy densities in the
1557  !! reconstruction [R Z3 T-2 ~> J m-2].
1558  real, dimension(SZI_(G)), intent(inout) :: vh !< The meridional energy transport [R Z3 L2 T-3 ~> J s-1].
1559  real, intent(in) :: dt !< Time increment [T ~> s].
1560  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1561  integer, intent(in) :: J !< The j-index to work on.
1562  integer, intent(in) :: ish !< The start i-index range to work on.
1563  integer, intent(in) :: ieh !< The end i-index range to work on.
1564  logical, intent(in) :: vol_CFL !< If true, rescale the ratio of face
1565  !! areas to the cell areas when estimating
1566  !! the CFL number.
1567  ! Local variables
1568  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim].
1569  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1570  ! with the same units as h_in.
1571  integer :: i
1572 
1573  do i=ish,ieh
1574  if (v(i) > 0.0) then
1575  if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1576  else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif
1577  curv_3 = hl(i,j) + hr(i,j) - 2.0*h(i,j)
1578  vh(i) = g%dx_Cv(i,j) * v(i) * ( hr(i,j) + cfl * &
1579  (0.5*(hl(i,j) - hr(i,j)) + curv_3*(cfl - 1.5)) )
1580  elseif (v(i) < 0.0) then
1581  if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1582  else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif
1583  curv_3 = hl(i,j+1) + hr(i,j+1) - 2.0*h(i,j+1)
1584  vh(i) = g%dx_Cv(i,j) * v(i) * ( hl(i,j+1) + cfl * &
1585  (0.5*(hr(i,j+1)-hl(i,j+1)) + curv_3*(cfl - 1.5)) )
1586  else
1587  vh(i) = 0.0
1588  endif
1589  enddo

Referenced by propagate_y().

Here is the caller graph for this function:

◆ ppm_angular_advect()

subroutine mom_internal_tides::ppm_angular_advect ( real, dimension(1-halo_ang:nangle+halo_ang), intent(in)  En2d,
real, dimension(1-halo_ang:nangle+halo_ang), intent(in)  CFL_ang,
real, dimension(0:nangle), intent(out)  Flux_En,
integer, intent(in)  NAngle,
real, intent(in)  dt,
integer, intent(in)  halo_ang 
)
private

This subroutine calculates the 1-d flux for advection in angular space using a monotonic piecewise parabolic scheme. This needs to be called from within i and j spatial loops.

Parameters
[in]nangleThe number of wave orientations in the discretized wave energy spectrum.
[in]dtTime increment [T ~> s].
[in]halo_angThe halo size in angular space
[in]en2dThe internal gravity wave energy density as a
[in]cfl_angThe CFL number of the energy advection across angles
[out]flux_enThe time integrated internal wave energy flux across angles [R Z3 T-2 ~> J m-2].

Definition at line 876 of file MOM_internal_tides.F90.

876  integer, intent(in) :: NAngle !< The number of wave orientations in the
877  !! discretized wave energy spectrum.
878  real, intent(in) :: dt !< Time increment [T ~> s].
879  integer, intent(in) :: halo_ang !< The halo size in angular space
880  real, dimension(1-halo_ang:NAngle+halo_ang), &
881  intent(in) :: En2d !< The internal gravity wave energy density as a
882  !! function of angular resolution [R Z3 T-2 ~> J m-2].
883  real, dimension(1-halo_ang:NAngle+halo_ang), &
884  intent(in) :: CFL_ang !< The CFL number of the energy advection across angles
885  real, dimension(0:NAngle), intent(out) :: Flux_En !< The time integrated internal wave energy flux
886  !! across angles [R Z3 T-2 ~> J m-2].
887  ! Local variables
888  real :: flux
889  real :: u_ang
890  real :: Angle_size
891  real :: I_Angle_size
892  real :: I_dt
893  integer :: a
894  real :: aR, aL, dMx, dMn, Ep, Ec, Em, dA, mA, a6
895 
896  i_dt = 1 / dt
897  angle_size = (8.0*atan(1.0)) / (real(nangle))
898  i_angle_size = 1 / angle_size
899  flux_en(:) = 0
900 
901  do a=0,nangle
902  u_ang = cfl_ang(a)*angle_size*i_dt
903  if (u_ang >= 0.0) then
904  ! Implementation of PPM-H3
905  ep = en2d(a+1)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
906  ec = en2d(a) *i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
907  em = en2d(a-1)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
908  al = ( 5.*ec + ( 2.*em - ep ) )/6. ! H3 estimate
909  al = max( min(ec,em), al) ; al = min( max(ec,em), al) ! Bound
910  ar = ( 5.*ec + ( 2.*ep - em ) )/6. ! H3 estimate
911  ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar) ! Bound
912  da = ar - al ; ma = 0.5*( ar + al )
913  if ((ep-ec)*(ec-em) <= 0.) then
914  al = ec ; ar = ec ! PCM for local extremum
915  elseif ( da*(ec-ma) > (da*da)/6. ) then
916  al = 3.*ec - 2.*ar !?
917  elseif ( da*(ec-ma) < - (da*da)/6. ) then
918  ar = 3.*ec - 2.*al !?
919  endif
920  a6 = 6.*ec - 3. * (ar + al) ! Curvature
921  ! CALCULATE FLUX RATE (Jm-2/s)
922  flux = u_ang*( ar + 0.5 * cfl_ang(a) * ( ( al - ar ) + a6 * ( 1. - 2./3. * cfl_ang(a) ) ) )
923  !flux = u_ang*( aR - 0.5 * CFL_ang(A) * ( ( aR - aL ) - a6 * ( 1. - 2./3. * CFL_ang(A) ) ) )
924  ! CALCULATE AMOUNT FLUXED (Jm-2)
925  flux_en(a) = dt * flux
926  !Flux_En(A) = (dt * I_Angle_size) * flux
927  else
928  ! Implementation of PPM-H3
929  ep = en2d(a+2)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
930  ec = en2d(a+1)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
931  em = en2d(a) *i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
932  al = ( 5.*ec + ( 2.*em - ep ) )/6. ! H3 estimate
933  al = max( min(ec,em), al) ; al = min( max(ec,em), al) ! Bound
934  ar = ( 5.*ec + ( 2.*ep - em ) )/6. ! H3 estimate
935  ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar) ! Bound
936  da = ar - al ; ma = 0.5*( ar + al )
937  if ((ep-ec)*(ec-em) <= 0.) then
938  al = ec ; ar = ec ! PCM for local extremum
939  elseif ( da*(ec-ma) > (da*da)/6. ) then
940  al = 3.*ec - 2.*ar
941  elseif ( da*(ec-ma) < - (da*da)/6. ) then
942  ar = 3.*ec - 2.*al
943  endif
944  a6 = 6.*ec - 3. * (ar + al) ! Curvature
945  ! CALCULATE FLUX RATE (Jm-2/s)
946  flux = u_ang*( ar + 0.5 * cfl_ang(a) * ( ( al - ar ) + a6 * ( 1. - 2./3. * cfl_ang(a) ) ) )
947  !flux = u_ang*( aL + 0.5 * CFL_ang(A) * ( ( aR - aL ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) )
948  ! CALCULATE AMOUNT FLUXED (Jm-2)
949  flux_en(a) = dt * flux
950  !Flux_En(A) = (dt * I_Angle_size) * flux
951  endif
952  enddo

Referenced by refract().

Here is the caller graph for this function:

◆ ppm_limit_pos()

subroutine mom_internal_tides::ppm_limit_pos ( real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_in,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout)  h_L,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout)  h_R,
real, intent(in)  h_min,
type(ocean_grid_type), intent(in)  G,
integer, intent(in)  iis,
integer, intent(in)  iie,
integer, intent(in)  jis,
integer, intent(in)  jie 
)
private

Limits the left/right edge values of the PPM reconstruction to give a reconstruction that is positive-definite. Here this is reinterpreted as giving a constant thickness if the mean thickness is less than h_min, with a minimum of h_min otherwise.

Parameters
[in]gThe ocean's grid structure.
[in]h_inThickness of layer (2D).
[in,out]h_lLeft edge value (2D).
[in,out]h_rRight edge value (2D).
[in]h_minThe minimum thickness that can be obtained by a concave parabolic fit.
[in]iisStart i-index for computations
[in]iieEnd i-index for computations
[in]jisStart j-index for computations
[in]jieEnd j-index for computations

Definition at line 2020 of file MOM_internal_tides.F90.

2020  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2021  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Thickness of layer (2D).
2022  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left edge value (2D).
2023  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right edge value (2D).
2024  real, intent(in) :: h_min !< The minimum thickness that can be
2025  !! obtained by a concave parabolic fit.
2026  integer, intent(in) :: iis !< Start i-index for computations
2027  integer, intent(in) :: iie !< End i-index for computations
2028  integer, intent(in) :: jis !< Start j-index for computations
2029  integer, intent(in) :: jie !< End j-index for computations
2030  ! Local variables
2031  real :: curv, dh, scale
2032  character(len=256) :: mesg ! The text of an error message
2033  integer :: i,j
2034 
2035  do j=jis,jie ; do i=iis,iie
2036  ! This limiter prevents undershooting minima within the domain with
2037  ! values less than h_min.
2038  curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
2039  if (curv > 0.0) then ! Only minima are limited.
2040  dh = h_r(i,j) - h_l(i,j)
2041  if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
2042  if (h_in(i,j) <= h_min) then
2043  h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
2044  elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then
2045  ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
2046  ! be limited in this case. 0 < scale < 1.
2047  scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
2048  h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
2049  h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
2050  endif
2051  endif
2052  endif
2053  enddo ; enddo

Referenced by ppm_reconstruction_x(), and ppm_reconstruction_y().

Here is the caller graph for this function:

◆ ppm_reconstruction_x()

subroutine mom_internal_tides::ppm_reconstruction_x ( real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)  h_in,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  h_l,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)  h_r,
type(ocean_grid_type), intent(in)  G,
type(loop_bounds_type), intent(in)  LB,
logical, intent(in), optional  simple_2nd 
)
private

Calculates left/right edge values for PPM reconstruction in x-direction.

Parameters
[in]gThe ocean's grid structure.
[in]h_inEnergy density in a sector (2D).
[out]h_lLeft edge value of reconstruction (2D).
[out]h_rRight edge value of reconstruction (2D).
[in]lbA structure with the active loop bounds.
[in]simple_2ndIf true, use the arithmetic mean energy densities as default edge values for a simple 2nd order scheme.

Definition at line 1867 of file MOM_internal_tides.F90.

1867  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1868  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Energy density in a sector (2D).
1869  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D).
1870  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D).
1871  type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.
1872  logical, optional, intent(in) :: simple_2nd !< If true, use the arithmetic mean
1873  !! energy densities as default edge values
1874  !! for a simple 2nd order scheme.
1875  ! Local variables
1876  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
1877  real, parameter :: oneSixth = 1./6.
1878  real :: h_ip1, h_im1
1879  real :: dMx, dMn
1880  logical :: use_CW84, use_2nd
1881  character(len=256) :: mesg ! The text of an error message
1882  integer :: i, j, isl, iel, jsl, jel, stencil
1883 
1884  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
1885  isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh
1886 
1887  ! This is the stencil of the reconstruction, not the scheme overall.
1888  stencil = 2 ; if (use_2nd) stencil = 1
1889 
1890  if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then
1891  write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
1892  & "x-halo that needs to be increased by ",i2,".")') &
1893  stencil + max(g%isd-isl,iel-g%ied)
1894  call mom_error(fatal,mesg)
1895  endif
1896  if ((jsl < g%jsd) .or. (jel > g%jed)) then
1897  write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
1898  & "y-halo that needs to be increased by ",i2,".")') &
1899  max(g%jsd-jsl,jel-g%jed)
1900  call mom_error(fatal,mesg)
1901  endif
1902 
1903  if (use_2nd) then
1904  do j=jsl,jel ; do i=isl,iel
1905  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1906  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1907  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
1908  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
1909  enddo ; enddo
1910  else
1911  do j=jsl,jel ; do i=isl-1,iel+1
1912  if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then
1913  slp(i,j) = 0.0
1914  else
1915  ! This uses a simple 2nd order slope.
1916  slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
1917  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
1918  dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
1919  dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
1920  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
1921  ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
1922  endif
1923  enddo ; enddo
1924 
1925  do j=jsl,jel ; do i=isl,iel
1926  ! Neighboring values should take into account any boundaries. The 3
1927  ! following sets of expressions are equivalent.
1928  ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
1929  ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
1930  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1931  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1932  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
1933  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
1934  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
1935  enddo ; enddo
1936  endif
1937 
1938  call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)

References mom_error_handler::mom_error(), and ppm_limit_pos().

Referenced by propagate_x().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ppm_reconstruction_y()

subroutine mom_internal_tides::ppm_reconstruction_y ( real, dimension(szi_(g),szj_(g)), intent(in)  h_in,
real, dimension(szi_(g),szj_(g)), intent(out)  h_l,
real, dimension(szi_(g),szj_(g)), intent(out)  h_r,
type(ocean_grid_type), intent(in)  G,
type(loop_bounds_type), intent(in)  LB,
logical, intent(in), optional  simple_2nd 
)
private

Calculates left/right edge valus for PPM reconstruction in y-direction.

Parameters
[in]gThe ocean's grid structure.
[in]h_inEnergy density in a sector (2D).
[out]h_lLeft edge value of reconstruction (2D).
[out]h_rRight edge value of reconstruction (2D).
[in]lbA structure with the active loop bounds.
[in]simple_2ndIf true, use the arithmetic mean energy densities as default edge values for a simple 2nd order scheme.

Definition at line 1943 of file MOM_internal_tides.F90.

1943  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1944  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Energy density in a sector (2D).
1945  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D).
1946  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D).
1947  type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.
1948  logical, optional, intent(in) :: simple_2nd !< If true, use the arithmetic mean
1949  !! energy densities as default edge values
1950  !! for a simple 2nd order scheme.
1951  ! Local variables
1952  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
1953  real, parameter :: oneSixth = 1./6.
1954  real :: h_jp1, h_jm1
1955  real :: dMx, dMn
1956  logical :: use_2nd
1957  character(len=256) :: mesg ! The text of an error message
1958  integer :: i, j, isl, iel, jsl, jel, stencil
1959 
1960  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
1961  isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1
1962 
1963  ! This is the stencil of the reconstruction, not the scheme overall.
1964  stencil = 2 ; if (use_2nd) stencil = 1
1965 
1966  if ((isl < g%isd) .or. (iel > g%ied)) then
1967  write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
1968  & "x-halo that needs to be increased by ",i2,".")') &
1969  max(g%isd-isl,iel-g%ied)
1970  call mom_error(fatal,mesg)
1971  endif
1972  if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then
1973  write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
1974  & "y-halo that needs to be increased by ",i2,".")') &
1975  stencil + max(g%jsd-jsl,jel-g%jed)
1976  call mom_error(fatal,mesg)
1977  endif
1978 
1979  if (use_2nd) then
1980  do j=jsl,jel ; do i=isl,iel
1981  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
1982  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
1983  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
1984  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
1985  enddo ; enddo
1986  else
1987  do j=jsl-1,jel+1 ; do i=isl,iel
1988  if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then
1989  slp(i,j) = 0.0
1990  else
1991  ! This uses a simple 2nd order slope.
1992  slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
1993  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
1994  dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
1995  dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
1996  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
1997  ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
1998  endif
1999  enddo ; enddo
2000 
2001  do j=jsl,jel ; do i=isl,iel
2002  ! Neighboring values should take into account any boundaries. The 3
2003  ! following sets of expressions are equivalent.
2004  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2005  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2006  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
2007  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
2008  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
2009  enddo ; enddo
2010  endif
2011 
2012  call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)

References mom_error_handler::mom_error(), and ppm_limit_pos().

Referenced by propagate_y().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ propagate()

subroutine mom_internal_tides::propagate ( real, dimension(g%isd:g%ied,g%jsd:g%jed,nangle), intent(inout)  En,
real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(in)  cn,
real, intent(in)  freq,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(int_tide_cs), pointer  CS,
integer, intent(in)  NAngle 
)
private

Propagates internal waves at a single frequency.

Parameters
[in,out]gThe ocean's grid structure.
[in]nangleThe number of wave orientations in the discretized wave energy spectrum.
[in,out]enThe internal gravity wave energy density as a
[in]cnBaroclinic mode speed [L T-1 ~> m s-1].
[in]freqWave frequency [T-1 ~> s-1].
[in]dtTime step [T ~> s].
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to int_tide_init.

Definition at line 957 of file MOM_internal_tides.F90.

957  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
958  integer, intent(in) :: NAngle !< The number of wave orientations in the
959  !! discretized wave energy spectrum.
960  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
961  intent(inout) :: En !< The internal gravity wave energy density as a
962  !! function of space and angular resolution,
963  !! [R Z3 T-2 ~> J m-2].
964  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
965  intent(in) :: cn !< Baroclinic mode speed [L T-1 ~> m s-1].
966  real, intent(in) :: freq !< Wave frequency [T-1 ~> s-1].
967  real, intent(in) :: dt !< Time step [T ~> s].
968  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
969  type(int_tide_CS), pointer :: CS !< The control structure returned by a
970  !! previous call to int_tide_init.
971  ! Local variables
972  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: &
973  speed ! The magnitude of the group velocity at the q points for corner adv [L T-1 ~> m s-1].
974  integer, parameter :: stencil = 2
975  real, dimension(SZIB_(G),SZJ_(G)) :: &
976  speed_x ! The magnitude of the group velocity at the Cu points [L T-1 ~> m s-1].
977  real, dimension(SZI_(G),SZJB_(G)) :: &
978  speed_y ! The magnitude of the group velocity at the Cv points [L T-1 ~> m s-1].
979  real, dimension(0:NAngle) :: &
980  cos_angle, sin_angle
981  real, dimension(NAngle) :: &
982  Cgx_av, Cgy_av, dCgx, dCgy
983  real :: f2 ! The squared Coriolis parameter [T-2 ~> s-2].
984  real :: Angle_size, I_Angle_size, angle
985  real :: Ifreq ! The inverse of the frequency [T ~> s]
986  real :: freq2 ! The frequency squared [T-2 ~> s-2]
987  type(loop_bounds_type) :: LB
988  integer :: is, ie, js, je, asd, aed, na
989  integer :: ish, ieh, jsh, jeh
990  integer :: i, j, a
991 
992  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na = size(en,3)
993  asd = 1-stencil ; aed = nangle+stencil
994 
995  ifreq = 1.0 / freq
996  freq2 = freq**2
997 
998  ! Define loop bounds: Need extensions on j-loop so propagate_y
999  ! (done after propagate_x) will have updated values in the halo
1000  ! for correct PPM reconstruction. Use if no teleporting and
1001  ! no pass_var between propagate_x and propagate_y.
1002  !jsh = js-3 ; jeh = je+3 ; ish = is ; ieh = ie
1003 
1004  ! Define loop bounds: Need 1-pt extensions on loops because
1005  ! teleporting eats up a halo point. Use if teleporting.
1006  ! Also requires pass_var before propagate_y.
1007  jsh = js-1 ; jeh = je+1 ; ish = is-1 ; ieh = ie+1
1008 
1009  angle_size = (8.0*atan(1.0)) / real(nangle)
1010  i_angle_size = 1.0 / angle_size
1011 
1012  if (cs%corner_adv) then
1013  ! IMPLEMENT CORNER ADVECTION IN HORIZONTAL--------------------
1014  ! FIND AVERAGE GROUP VELOCITY (SPEED) AT CELL CORNERS
1015  ! NOTE: THIS HAS NOT BE ADAPTED FOR REFLECTION YET (BDM)!!
1016  ! Fix indexing here later
1017  speed(:,:) = 0
1018  do j=jsh-1,jeh ; do i=ish-1,ieh
1019  f2 = g%CoriolisBu(i,j)**2
1020  speed(i,j) = 0.25*(cn(i,j) + cn(i+1,j) + cn(i+1,j+1) + cn(i,j+1)) * &
1021  sqrt(max(freq2 - f2, 0.0)) * ifreq
1022  enddo ; enddo
1023  do a=1,na
1024  ! Apply the propagation WITH CORNER ADVECTION/FINITE VOLUME APPROACH.
1025  lb%jsh = js ; lb%jeh = je ; lb%ish = is ; lb%ieh = ie
1026  call propagate_corner_spread(en(:,:,a), a, nangle, speed, dt, g, cs, lb)
1027  enddo ! a-loop
1028  else
1029  ! IMPLEMENT PPM ADVECTION IN HORIZONTAL-----------------------
1030  ! These could be in the control structure, as they do not vary.
1031  do a=0,na
1032  ! These are the angles at the cell edges...
1033  angle = (real(a) - 0.5) * angle_size
1034  cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
1035  enddo
1036 
1037  do a=1,na
1038  cgx_av(a) = (sin_angle(a) - sin_angle(a-1)) * i_angle_size
1039  cgy_av(a) = -(cos_angle(a) - cos_angle(a-1)) * i_angle_size
1040  dcgx(a) = sqrt(0.5 + 0.5*(sin_angle(a)*cos_angle(a) - &
1041  sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
1042  cgx_av(a)**2)
1043  dcgy(a) = sqrt(0.5 - 0.5*(sin_angle(a)*cos_angle(a) - &
1044  sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
1045  cgy_av(a)**2)
1046  enddo
1047 
1048  do j=jsh,jeh ; do i=ish-1,ieh
1049  f2 = 0.5 * (g%CoriolisBu(i,j)**2 + g%CoriolisBu(i,j-1)**2)
1050  speed_x(i,j) = 0.5*(cn(i,j) + cn(i+1,j)) * g%mask2dCu(i,j) * &
1051  sqrt(max(freq2 - f2, 0.0)) * ifreq
1052  enddo ; enddo
1053  do j=jsh-1,jeh ; do i=ish,ieh
1054  f2 = 0.5 * (g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j)**2)
1055  speed_y(i,j) = 0.5*(cn(i,j) + cn(i,j+1)) * g%mask2dCv(i,j) * &
1056  sqrt(max(freq2 - f2, 0.0)) * ifreq
1057  enddo ; enddo
1058 
1059  ! Apply propagation in x-direction (reflection included)
1060  lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
1061  call propagate_x(en(:,:,:), speed_x, cgx_av(:), dcgx(:), dt, g, us, cs%nAngle, cs, lb)
1062 
1063  ! Check for energy conservation on computational domain (for debugging)
1064  !call sum_En(G,CS,En(:,:,:),'post-propagate_x')
1065 
1066  ! Update halos
1067  call pass_var(en(:,:,:),g%domain)
1068 
1069  ! Apply propagation in y-direction (reflection included)
1070  ! LB%jsh = js ; LB%jeh = je ; LB%ish = is ; LB%ieh = ie ! Use if no teleport
1071  lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
1072  call propagate_y(en(:,:,:), speed_y, cgy_av(:), dcgy(:), dt, g, us, cs%nAngle, cs, lb)
1073 
1074  ! Check for energy conservation on computational domain (for debugging)
1075  !call sum_En(G,CS,En(:,:,:),'post-propagate_y')
1076  endif
1077 

References propagate_corner_spread(), propagate_x(), and propagate_y().

Referenced by propagate_int_tide().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ propagate_corner_spread()

subroutine mom_internal_tides::propagate_corner_spread ( real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(inout)  En,
integer, intent(in)  energized_wedge,
integer, intent(in)  NAngle,
real, dimension(g%isdb:g%iedb,g%jsd:g%jed), intent(in)  speed,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(int_tide_cs), pointer  CS,
type(loop_bounds_type), intent(in)  LB 
)
private

This subroutine does first-order corner advection. It was written with the hopes of smoothing out the garden sprinkler effect, but is too numerically diffusive to be of much use as of yet. It is not yet compatible with reflection schemes (BDM).

Parameters
[in]gThe ocean's grid structure.
[in,out]enThe energy density integrated over an angular
[in]speedThe magnitude of the group velocity at the cell
[in]energized_wedgeIndex of current ray direction.
[in]nangleThe number of wave orientations in the discretized wave energy spectrum.
[in]dtTime increment [T ~> s].
csThe control structure returned by a previous call to continuity_PPM_init.
[in]lbA structure with the active energy loop bounds.

Definition at line 1084 of file MOM_internal_tides.F90.

1084  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1085  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
1086  intent(inout) :: En !< The energy density integrated over an angular
1087  !! band [R Z3 T-2 ~> J m-2], intent in/out.
1088  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed), &
1089  intent(in) :: speed !< The magnitude of the group velocity at the cell
1090  !! corner points [L T-1 ~> m s-1].
1091  integer, intent(in) :: energized_wedge !< Index of current ray direction.
1092  integer, intent(in) :: NAngle !< The number of wave orientations in the
1093  !! discretized wave energy spectrum.
1094  real, intent(in) :: dt !< Time increment [T ~> s].
1095  type(int_tide_CS), pointer :: CS !< The control structure returned by a previous
1096  !! call to continuity_PPM_init.
1097  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1098  ! Local variables
1099  integer :: i, j, k, ish, ieh, jsh, jeh, m
1100  real :: TwoPi, Angle_size
1101  real :: energized_angle ! angle through center of current wedge
1102  real :: theta ! angle at edge of wedge
1103  real :: Nsubrays ! number of sub-rays for averaging
1104  ! count includes the two rays that bound the current wedge,
1105  ! i.e. those at -dtheta/2 and +dtheta/2 from energized angle
1106  real :: I_Nsubwedges ! inverse of number of sub-wedges
1107  real :: cos_thetaDT, sin_thetaDT ! cos(theta)*dt, sin(theta)*dt
1108  real :: xNE,xNW,xSW,xSE,yNE,yNW,ySW,ySE ! corner point coordinates of advected fluid parcel
1109  real :: CFL_xNE,CFL_xNW,CFL_xSW,CFL_xSE,CFL_yNE,CFL_yNW,CFL_ySW,CFL_ySE,CFL_max
1110  real :: xN,xS,xE,xW,yN,yS,yE,yW ! intersection point coordinates of parcel edges and grid
1111  real :: xCrn,yCrn ! grid point contained within advected fluid parcel
1112  real :: xg,yg ! grid point of interest
1113  real :: slopeN,slopeW,slopeS,slopeE, bN,bW,bS,bE ! parameters defining parcel sides
1114  real :: aNE,aN,aNW,aW,aSW,aS,aSE,aE,aC ! sub-areas of advected parcel
1115  real :: a_total ! total area of advected parcel
1116  real :: a1,a2,a3,a4 ! areas used in calculating polygon areas (sub-areas) of advected parcel
1117  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: x,y ! coordinates of cell corners
1118  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: Idx,Idy ! inverse of dx,dy at cell corners
1119  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: dx,dy ! dx,dy at cell corners
1120  real, dimension(2) :: E_new ! Energy in cell after advection for subray [R Z3 T-2 ~> J m-2]; set size
1121  ! here to define Nsubrays - this should be made an input option later!
1122 
1123  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1124  twopi = (8.0*atan(1.0))
1125  nsubrays = real(size(e_new))
1126  i_nsubwedges = 1./(nsubrays - 1)
1127 
1128  angle_size = twopi / real(nangle)
1129  energized_angle = angle_size * real(energized_wedge - 1) ! for a=1 aligned with x-axis
1130  !energized_angle = Angle_size * real(energized_wedge - 1) + 2.0*Angle_size !
1131  !energized_angle = Angle_size * real(energized_wedge - 1) + 0.5*Angle_size !
1132  do j=jsh-1,jeh ; do i=ish-1,ieh
1133  ! This will only work for a Cartesian grid for which G%geoLonBu is in the same units has dx.
1134  ! This needs to be extensively revised to work for a general grid.
1135  x(i,j) = g%US%m_to_L*g%geoLonBu(i,j)
1136  y(i,j) = g%US%m_to_L*g%geoLatBu(i,j)
1137  idx(i,j) = g%IdxBu(i,j) ; dx(i,j) = g%dxBu(i,j)
1138  idy(i,j) = g%IdyBu(i,j) ; dy(i,j) = g%dyBu(i,j)
1139  enddo ; enddo
1140 
1141  do j=jsh,jeh ; do i=ish,ieh
1142  do m=1,int(nsubrays)
1143  theta = energized_angle - 0.5*angle_size + real(m - 1)*angle_size*i_nsubwedges
1144  if (theta < 0.0) then
1145  theta = theta + twopi
1146  elseif (theta > twopi) then
1147  theta = theta - twopi
1148  endif
1149  cos_thetadt = cos(theta)*dt
1150  sin_thetadt = sin(theta)*dt
1151 
1152  ! corner point coordinates of advected fluid parcel ----------
1153  xg = x(i,j); yg = y(i,j)
1154  xne = xg - speed(i,j)*cos_thetadt
1155  yne = yg - speed(i,j)*sin_thetadt
1156  cfl_xne = (xg-xne)*idx(i,j)
1157  cfl_yne = (yg-yne)*idy(i,j)
1158 
1159  xg = x(i-1,j); yg = y(i-1,j)
1160  xnw = xg - speed(i-1,j)*cos_thetadt
1161  ynw = yg - speed(i-1,j)*sin_thetadt
1162  cfl_xnw = (xg-xnw)*idx(i-1,j)
1163  cfl_ynw = (yg-ynw)*idy(i-1,j)
1164 
1165  xg = x(i-1,j-1); yg = y(i-1,j-1)
1166  xsw = xg - speed(i-1,j-1)*cos_thetadt
1167  ysw = yg - speed(i-1,j-1)*sin_thetadt
1168  cfl_xsw = (xg-xsw)*idx(i-1,j-1)
1169  cfl_ysw = (yg-ysw)*idy(i-1,j-1)
1170 
1171  xg = x(i,j-1); yg = y(i,j-1)
1172  xse = xg - speed(i,j-1)*cos_thetadt
1173  yse = yg - speed(i,j-1)*sin_thetadt
1174  cfl_xse = (xg-xse)*idx(i,j-1)
1175  cfl_yse = (yg-yse)*idy(i,j-1)
1176 
1177  cfl_max = max(abs(cfl_xne),abs(cfl_xnw),abs(cfl_xsw), &
1178  abs(cfl_xse),abs(cfl_yne),abs(cfl_ynw), &
1179  abs(cfl_ysw),abs(cfl_yse))
1180  if (cfl_max > 1.0) then
1181  call mom_error(warning, "propagate_corner_spread: CFL exceeds 1.", .true.)
1182  endif
1183 
1184  ! intersection point coordinates of parcel edges and cell edges ---
1185  if (0.0 <= theta .and. theta < 0.25*twopi) then
1186  xn = x(i-1,j-1)
1187  yw = y(i-1,j-1)
1188  elseif (0.25*twopi <= theta .and. theta < 0.5*twopi) then
1189  xn = x(i,j-1)
1190  yw = y(i,j-1)
1191  elseif (0.5*twopi <= theta .and. theta < 0.75*twopi) then
1192  xn = x(i,j)
1193  yw = y(i,j)
1194  elseif (0.75*twopi <= theta .and. theta <= 1.00*twopi) then
1195  xn = x(i-1,j)
1196  yw = y(i-1,j)
1197  endif
1198  xs = xn
1199  ye = yw
1200 
1201  ! north intersection
1202  slopen = (yne - ynw)/(xne - xnw)
1203  bn = -slopen*xne + yne
1204  yn = slopen*xn + bn
1205  ! west intersection
1206  if (xnw == xsw) then
1207  xw = xnw
1208  else
1209  slopew = (ynw - ysw)/(xnw - xsw)
1210  bw = -slopew*xnw + ynw
1211  xw = (yw - bw)/slopew
1212  endif
1213  ! south intersection
1214  slopes = (ysw - yse)/(xsw - xse)
1215  bs = -slopes*xsw + ysw
1216  ys = slopes*xs + bs
1217  ! east intersection
1218  if (xne == xse) then
1219  xe = xne
1220  else
1221  slopee = (yse - yne)/(xse - xne)
1222  be = -slopee*xse + yse
1223  xe = (ye - be)/slopee
1224  endif
1225 
1226  ! areas --------------------------------------------
1227  ane = 0.0; an = 0.0; anw = 0.0; ! initialize areas
1228  aw = 0.0; asw = 0.0; as = 0.0; ! initialize areas
1229  ase = 0.0; ae = 0.0; ac = 0.0; ! initialize areas
1230  if (0.0 <= theta .and. theta < 0.25*twopi) then
1231  xcrn = x(i-1,j-1); ycrn = y(i-1,j-1)
1232  ! west area
1233  a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1234  a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1235  a3 = (yw - ynw)*(0.5*(xw + xnw))
1236  a4 = (ynw - yn)*(0.5*(xnw + xn))
1237  aw = a1 + a2 + a3 + a4
1238  ! southwest area
1239  a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1240  a2 = (ys - ysw)*(0.5*(xs + xsw))
1241  a3 = (ysw - yw)*(0.5*(xsw + xw))
1242  a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1243  asw = a1 + a2 + a3 + a4
1244  ! south area
1245  a1 = (ye - yse)*(0.5*(xe + xse))
1246  a2 = (yse - ys)*(0.5*(xse + xs))
1247  a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1248  a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1249  as = a1 + a2 + a3 + a4
1250  ! area within cell
1251  a1 = (yne - ye)*(0.5*(xne + xe))
1252  a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1253  a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1254  a4 = (yn - yne)*(0.5*(xn + xne))
1255  ac = a1 + a2 + a3 + a4
1256  elseif (0.25*twopi <= theta .and. theta < 0.5*twopi) then
1257  xcrn = x(i,j-1); ycrn = y(i,j-1)
1258  ! south area
1259  a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1260  a2 = (ys - ysw)*(0.5*(xs + xsw))
1261  a3 = (ysw - yw)*(0.5*(xsw + xw))
1262  a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1263  as = a1 + a2 + a3 + a4
1264  ! southeast area
1265  a1 = (ye - yse)*(0.5*(xe + xse))
1266  a2 = (yse - ys)*(0.5*(xse + xs))
1267  a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1268  a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1269  ase = a1 + a2 + a3 + a4
1270  ! east area
1271  a1 = (yne - ye)*(0.5*(xne + xe))
1272  a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1273  a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1274  a4 = (yn - yne)*(0.5*(xn + xne))
1275  ae = a1 + a2 + a3 + a4
1276  ! area within cell
1277  a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1278  a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1279  a3 = (yw - ynw)*(0.5*(xw + xnw))
1280  a4 = (ynw - yn)*(0.5*(xnw + xn))
1281  ac = a1 + a2 + a3 + a4
1282  elseif (0.5*twopi <= theta .and. theta < 0.75*twopi) then
1283  xcrn = x(i,j); ycrn = y(i,j)
1284  ! east area
1285  a1 = (ye - yse)*(0.5*(xe + xse))
1286  a2 = (yse - ys)*(0.5*(xse + xs))
1287  a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1288  a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1289  ae = a1 + a2 + a3 + a4
1290  ! northeast area
1291  a1 = (yne - ye)*(0.5*(xne + xe))
1292  a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1293  a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1294  a4 = (yn - yne)*(0.5*(xn + xne))
1295  ane = a1 + a2 + a3 + a4
1296  ! north area
1297  a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1298  a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1299  a3 = (yw - ynw)*(0.5*(xw + xnw))
1300  a4 = (ynw - yn)*(0.5*(xnw + xn))
1301  an = a1 + a2 + a3 + a4
1302  ! area within cell
1303  a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1304  a2 = (ys - ysw)*(0.5*(xs + xsw))
1305  a3 = (ysw - yw)*(0.5*(xsw + xw))
1306  a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1307  ac = a1 + a2 + a3 + a4
1308  elseif (0.75*twopi <= theta .and. theta <= 1.00*twopi) then
1309  xcrn = x(i-1,j); ycrn = y(i-1,j)
1310  ! north area
1311  a1 = (yne - ye)*(0.5*(xne + xe))
1312  a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1313  a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1314  a4 = (yn - yne)*(0.5*(xn + xne))
1315  an = a1 + a2 + a3 + a4
1316  ! northwest area
1317  a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1318  a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1319  a3 = (yw - ynw)*(0.5*(xw + xnw))
1320  a4 = (ynw - yn)*(0.5*(xnw + xn))
1321  anw = a1 + a2 + a3 + a4
1322  ! west area
1323  a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1324  a2 = (ys - ysw)*(0.5*(xs + xsw))
1325  a3 = (ysw - yw)*(0.5*(xsw + xw))
1326  a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1327  aw = a1 + a2 + a3 + a4
1328  ! area within cell
1329  a1 = (ye - yse)*(0.5*(xe + xse))
1330  a2 = (yse - ys)*(0.5*(xse + xs))
1331  a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1332  a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1333  ac = a1 + a2 + a3 + a4
1334  endif
1335 
1336  ! energy weighting ----------------------------------------
1337  a_total = ane + an + anw + aw + asw + as + ase + ae + ac
1338  e_new(m) = (ane*en(i+1,j+1) + an*en(i,j+1) + anw*en(i-1,j+1) + &
1339  aw*en(i-1,j) + asw*en(i-1,j-1) + as*en(i,j-1) + &
1340  ase*en(i+1,j-1) + ae*en(i+1,j) + ac*en(i,j)) / (dx(i,j)*dy(i,j))
1341  enddo ! m-loop
1342  ! update energy in cell
1343  en(i,j) = sum(e_new)/nsubrays
1344  enddo ; enddo

References mom_error_handler::mom_error().

Referenced by propagate().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ propagate_int_tide()

subroutine, public mom_internal_tides::propagate_int_tide ( real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
real, dimension(szi_(g),szj_(g),cs%nmode), intent(in)  cn,
real, dimension(szi_(g),szj_(g)), intent(in)  TKE_itidal_input,
real, dimension(szi_(g),szj_(g)), intent(in)  vel_btTide,
real, dimension(szi_(g),szj_(g)), intent(in)  Nb,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(int_tide_cs), pointer  CS 
)

Calls subroutines in this file that are needed to refract, propagate, and dissipate energy density of the internal tide.

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]tvPointer to thermodynamic variables (needed for wave structure).
[in]tke_itidal_inputThe energy input to the internal waves [R Z3 T-3 ~> W m-2].
[in]vel_bttideBarotropic velocity read from file [L T-1 ~> m s-1].
[in]nbNear-bottom buoyancy frequency [T-1 ~> s-1].
[in]dtLength of time over which to advance the internal tides [T ~> s].
csThe control structure returned by a previous call to int_tide_init.
[in]cnThe internal wave speeds of each

Definition at line 154 of file MOM_internal_tides.F90.

154  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
155  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
156  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
157  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
158  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
159  type(thermo_var_ptrs), intent(in) :: tv !< Pointer to thermodynamic variables
160  !! (needed for wave structure).
161  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: TKE_itidal_input !< The energy input to the
162  !! internal waves [R Z3 T-3 ~> W m-2].
163  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: vel_btTide !< Barotropic velocity read
164  !! from file [L T-1 ~> m s-1].
165  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Nb !< Near-bottom buoyancy frequency [T-1 ~> s-1].
166  real, intent(in) :: dt !< Length of time over which to advance
167  !! the internal tides [T ~> s].
168  type(int_tide_CS), pointer :: CS !< The control structure returned by a
169  !! previous call to int_tide_init.
170  real, dimension(SZI_(G),SZJ_(G),CS%nMode), &
171  intent(in) :: cn !< The internal wave speeds of each
172  !! mode [L T-1 ~> m s-1].
173  ! Local variables
174  real, dimension(SZI_(G),SZJ_(G),2) :: &
175  test
176  real, dimension(SZI_(G),SZJ_(G),CS%nFreq,CS%nMode) :: &
177  tot_En_mode, & ! energy summed over angles only [R Z3 T-2 ~> J m-2]
178  Ub, & ! near-bottom horizontal velocity of wave (modal) [L T-1 ~> m s-1]
179  Umax ! Maximum horizontal velocity of wave (modal) [L T-1 ~> m s-1]
180  real, dimension(SZI_(G),SZJB_(G)) :: &
181  flux_heat_y, &
182  flux_prec_y
183  real, dimension(SZI_(G),SZJ_(G)) :: &
184  tot_En, & ! energy summed over angles, modes, frequencies [R Z3 T-2 ~> J m-2]
185  tot_leak_loss, tot_quad_loss, tot_itidal_loss, tot_Froude_loss, tot_allprocesses_loss, &
186  ! energy loss rates summed over angle, freq, and mode [R Z3 T-3 ~> W m-2]
187  drag_scale, & ! bottom drag scale [T-1 ~> s-1]
188  itidal_loss_mode, allprocesses_loss_mode
189  ! energy loss rates for a given mode and frequency (summed over angles) [R Z3 T-3 ~> W m-2]
190  real :: frac_per_sector, f2, Kmag2
191  real :: I_D_here ! The inverse of the local depth [Z-1 ~> m-1]
192  real :: I_rho0 ! The inverse fo the Boussinesq density [R-1 ~> m3 kg-1]
193  real :: freq2 ! The frequency squared [T-2 ~> s-2]
194  real :: c_phase ! The phase speed [L T-1 ~> m s-1]
195  real :: loss_rate ! An energy loss rate [T-1 ~> s-1]
196  real :: Fr2_max
197  real :: cn_subRO ! A tiny wave speed to prevent division by zero [L T-1 ~> m s-1]
198  real :: En_new, En_check ! Energies for debugging [R Z3 T-2 ~> J m-2]
199  real :: En_initial, Delta_E_check ! Energies for debugging [R Z3 T-2 ~> J m-2]
200  real :: TKE_Froude_loss_check, TKE_Froude_loss_tot ! Energy losses for debugging [R Z3 T-3 ~> W m-2]
201  character(len=160) :: mesg ! The text of an error message
202  integer :: a, m, fr, i, j, is, ie, js, je, isd, ied, jsd, jed, nAngle, nzm
203  integer :: id_g, jd_g ! global (decomp-invar) indices (for debugging)
204  type(group_pass_type), save :: pass_test, pass_En
205 
206  if (.not.associated(cs)) return
207  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
208  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nangle = cs%NAngle
209  i_rho0 = 1.0 / gv%Rho0
210  cn_subro = 1e-100*us%m_s_to_L_T ! The hard-coded value here might need to increase.
211 
212  ! Set the wave speeds for the modes, using cg(n) ~ cg(1)/n.**********************
213  ! This is wrong, of course, but it works reasonably in some cases.
214  ! Uncomment if wave_speed is not used to calculate the true values (BDM).
215  !do m=1,CS%nMode ; do j=jsd,jed ; do i=isd,ied
216  ! cn(i,j,m) = cn(i,j,1) / real(m)
217  !enddo ; enddo ; enddo
218 
219  ! Add the forcing.***************************************************************
220  if (cs%energized_angle <= 0) then
221  frac_per_sector = 1.0 / real(cs%nAngle * cs%nMode * cs%nFreq)
222  do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
223  f2 = 0.25*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
224  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
225  if (cs%frequency(fr)**2 > f2) &
226  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + dt*frac_per_sector*(1.0-cs%q_itides) * &
227  tke_itidal_input(i,j)
228  enddo ; enddo ; enddo ; enddo ; enddo
229  elseif (cs%energized_angle <= cs%nAngle) then
230  frac_per_sector = 1.0 / real(cs%nMode * cs%nFreq)
231  a = cs%energized_angle
232  do m=1,cs%nMode ; do fr=1,cs%nFreq ; do j=js,je ; do i=is,ie
233  f2 = 0.25*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
234  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
235  if (cs%frequency(fr)**2 > f2) &
236  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + dt*frac_per_sector*(1.0-cs%q_itides) * &
237  tke_itidal_input(i,j)
238  enddo ; enddo ; enddo ; enddo
239  else
240  call mom_error(warning, "Internal tide energy is being put into a angular "//&
241  "band that does not exist.")
242  endif
243 
244  ! Pass a test vector to check for grid rotation in the halo updates.
245  do j=jsd,jed ; do i=isd,ied ; test(i,j,1) = 1.0 ; test(i,j,2) = 0.0 ; enddo ; enddo
246  do m=1,cs%nMode ; do fr=1,cs%nFreq
247  call create_group_pass(pass_en, cs%En(:,:,:,fr,m), g%domain)
248  enddo ; enddo
249  call create_group_pass(pass_test, test(:,:,1), test(:,:,2), g%domain, stagger=agrid)
250  call start_group_pass(pass_test, g%domain)
251 
252  ! Apply half the refraction.
253  do m=1,cs%nMode ; do fr=1,cs%nFreq
254  call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, &
255  g, us, cs%nAngle, cs%use_PPMang)
256  enddo ; enddo
257 
258  ! Check for En<0 - for debugging, delete later
259  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
260  do j=js,je ; do i=is,ie
261  if (cs%En(i,j,a,fr,m)<0.0) then
262  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
263  write(mesg,*) 'After first refraction: En<0.0 at ig=', id_g, ', jg=', jd_g, &
264  'En=',cs%En(i,j,a,fr,m)
265  call mom_error(warning, "propagate_int_tide: "//trim(mesg))
266  cs%En(i,j,a,fr,m) = 0.0
267 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
268  endif
269  enddo ; enddo
270  enddo ; enddo ; enddo
271 
272  call do_group_pass(pass_en, g%domain)
273 
274  call complete_group_pass(pass_test, g%domain)
275 
276  ! Rotate points in the halos as necessary.
277  call correct_halo_rotation(cs%En, test, g, cs%nAngle)
278 
279  ! Propagate the waves.
280  do m=1,cs%NMode ; do fr=1,cs%Nfreq
281  call propagate(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), dt, &
282  g, us, cs, cs%NAngle)
283  enddo ; enddo
284 
285  ! Check for En<0 - for debugging, delete later
286  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
287  do j=js,je ; do i=is,ie
288  if (cs%En(i,j,a,fr,m)<0.0) then
289  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
290  if (abs(cs%En(i,j,a,fr,m))>1.0) then ! only print if large
291  write(mesg,*) 'After propagation: En<0.0 at ig=', id_g, ', jg=', jd_g, &
292  'En=', cs%En(i,j,a,fr,m)
293  call mom_error(warning, "propagate_int_tide: "//trim(mesg))
294 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
295  endif
296  cs%En(i,j,a,fr,m) = 0.0
297  endif
298  enddo ; enddo
299  enddo ; enddo ; enddo
300 
301  ! Apply the other half of the refraction.
302  do m=1,cs%NMode ; do fr=1,cs%Nfreq
303  call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, &
304  g, us, cs%NAngle, cs%use_PPMang)
305  enddo ; enddo
306 
307  ! Check for En<0 - for debugging, delete later
308  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
309  do j=js,je ; do i=is,ie
310  if (cs%En(i,j,a,fr,m)<0.0) then
311  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
312  write(mesg,*) 'After second refraction: En<0.0 at ig=', id_g, ', jg=', jd_g, &
313  'En=',cs%En(i,j,a,fr,m)
314  call mom_error(warning, "propagate_int_tide: "//trim(mesg))
315  cs%En(i,j,a,fr,m) = 0.0
316 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
317  endif
318  enddo ; enddo
319  enddo ; enddo ; enddo
320 
321  ! Apply various dissipation mechanisms.
322  if (cs%apply_background_drag .or. cs%apply_bottom_drag &
323  .or. cs%apply_wave_drag .or. cs%apply_Froude_drag &
324  .or. (cs%id_tot_En > 0)) then
325  tot_en(:,:) = 0.0
326  tot_en_mode(:,:,:,:) = 0.0
327  do m=1,cs%NMode ; do fr=1,cs%Nfreq
328  do j=jsd,jed ; do i=isd,ied ; do a=1,cs%nAngle
329  tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
330  tot_en_mode(i,j,fr,m) = tot_en_mode(i,j,fr,m) + cs%En(i,j,a,fr,m)
331  enddo ; enddo ; enddo
332  enddo ; enddo
333  endif
334 
335  ! Extract the energy for mixing due to misc. processes (background leakage)------
336  if (cs%apply_background_drag) then
337  do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
338  ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale
339  ! to each En component (technically not correct; fix later)
340  cs%TKE_leak_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * cs%decay_rate ! loss rate [R Z3 T-3 ~> W m-2]
341  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt * cs%decay_rate) ! implicit update
342  enddo ; enddo ; enddo ; enddo ; enddo
343  endif
344  ! Check for En<0 - for debugging, delete later
345  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
346  do j=js,je ; do i=is,ie
347  if (cs%En(i,j,a,fr,m)<0.0) then
348  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
349  write(mesg,*) 'After leak loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
350  'En=',cs%En(i,j,a,fr,m)
351  call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
352  cs%En(i,j,a,fr,m) = 0.0
353 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
354  endif
355  enddo ; enddo
356  enddo ; enddo ; enddo
357 
358  ! Extract the energy for mixing due to bottom drag-------------------------------
359  if (cs%apply_bottom_drag) then
360  do j=jsd,jed ; do i=isd,ied
361  ! Note the 1 m dimensional scale here. Should this be a parameter?
362  i_d_here = 1.0 / (max(g%bathyT(i,j), 1.0*us%m_to_Z))
363  drag_scale(i,j) = cs%cdrag * sqrt(max(0.0, us%L_to_Z**2*vel_bttide(i,j)**2 + &
364  tot_en(i,j) * i_rho0 * i_d_here)) * i_d_here
365  enddo ; enddo
366  do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
367  ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale
368  ! to each En component (technically not correct; fix later)
369  cs%TKE_quad_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * drag_scale(i,j) ! loss rate
370  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt * drag_scale(i,j)) ! implicit update
371  enddo ; enddo ; enddo ; enddo ; enddo
372  endif
373  ! Check for En<0 - for debugging, delete later
374  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
375  do j=js,je ; do i=is,ie
376  if (cs%En(i,j,a,fr,m)<0.0) then
377  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
378  write(mesg,*) 'After bottom loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
379  'En=',cs%En(i,j,a,fr,m)
380  call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
381  cs%En(i,j,a,fr,m) = 0.0
382 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
383  !stop
384  endif
385  enddo ; enddo
386  enddo ; enddo ; enddo
387 
388  ! Extract the energy for mixing due to scattering (wave-drag)--------------------
389  ! still need to allow a portion of the extracted energy to go to higher modes.
390  ! First, find velocity profiles
391  if (cs%apply_wave_drag .or. cs%apply_Froude_drag) then
392  do m=1,cs%NMode ; do fr=1,cs%Nfreq
393  ! Calculate modal structure for given mode and frequency
394  call wave_structure(h, tv, g, gv, us, cn(:,:,m), m, cs%frequency(fr), &
395  cs%wave_structure_CSp, tot_en_mode(:,:,fr,m), full_halos=.true.)
396  ! Pick out near-bottom and max horizontal baroclinic velocity values at each point
397  do j=jsd,jed ; do i=isd,ied
398  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
399  nzm = cs%wave_structure_CSp%num_intfaces(i,j)
400  ub(i,j,fr,m) = cs%wave_structure_CSp%Uavg_profile(i,j,nzm)
401  umax(i,j,fr,m) = maxval(cs%wave_structure_CSp%Uavg_profile(i,j,1:nzm))
402  enddo ; enddo ! i-loop, j-loop
403  enddo ; enddo ! fr-loop, m-loop
404  endif ! apply_wave or _Froude_drag (Ub or Umax needed)
405  ! Finally, apply loss
406  if (cs%apply_wave_drag) then
407  ! Calculate loss rate and apply loss over the time step
408  call itidal_lowmode_loss(g, us, cs, nb, ub, cs%En, cs%TKE_itidal_loss_fixed, &
409  cs%TKE_itidal_loss, dt, full_halos=.false.)
410  endif
411  ! Check for En<0 - for debugging, delete later
412  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
413  do j=js,je ; do i=is,ie
414  if (cs%En(i,j,a,fr,m)<0.0) then
415  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
416  write(mesg,*) 'After wave drag loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
417  'En=',cs%En(i,j,a,fr,m)
418  call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
419  cs%En(i,j,a,fr,m) = 0.0
420 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
421  endif
422  enddo ; enddo
423  enddo ; enddo ; enddo
424 
425  ! Extract the energy for mixing due to wave breaking-----------------------------
426  if (cs%apply_Froude_drag) then
427  ! Pick out maximum baroclinic velocity values; calculate Fr=max(u)/cg
428  do m=1,cs%NMode ; do fr=1,cs%Nfreq
429  freq2 = cs%frequency(fr)**2
430  do j=jsd,jed ; do i=isd,ied
431  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
432  ! Calculate horizontal phase velocity magnitudes
433  f2 = 0.25*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
434  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
435  kmag2 = (freq2 - f2) / (cn(i,j,m)**2 + cn_subro**2)
436  c_phase = 0.0
437  if (kmag2 > 0.0) then
438  c_phase = sqrt(freq2/kmag2)
439  nzm = cs%wave_structure_CSp%num_intfaces(i,j)
440  fr2_max = (umax(i,j,fr,m) / c_phase)**2
441  ! Dissipate energy if Fr>1; done here with an arbitrary time scale
442  if (fr2_max > 1.0) then
443  en_initial = sum(cs%En(i,j,:,fr,m)) ! for debugging
444  ! Calculate effective decay rate [s-1] if breaking occurs over a time step
445  loss_rate = (1.0 - fr2_max) / (fr2_max * dt)
446  do a=1,cs%nAngle
447  ! Determine effective dissipation rate (Wm-2)
448  cs%TKE_Froude_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * abs(loss_rate)
449  ! Update energy
450  en_new = cs%En(i,j,a,fr,m)/fr2_max ! for debugging
451  en_check = cs%En(i,j,a,fr,m) - cs%TKE_Froude_loss(i,j,a,fr,m)*dt ! for debugging
452  ! Re-scale (reduce) energy due to breaking
453  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m)/fr2_max
454  ! Check (for debugging only)
455  if (abs(en_new - en_check) > 1e-10*us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**2) then
456  call mom_error(warning, "MOM_internal_tides: something is wrong with Fr-breaking.", &
457  all_print=.true.)
458  write(mesg,*) "En_new=", en_new , "En_check=", en_check
459  call mom_error(warning, "MOM_internal_tides: "//trim(mesg), all_print=.true.)
460  endif
461  enddo
462  ! Check (for debugging)
463  delta_e_check = en_initial - sum(cs%En(i,j,:,fr,m))
464  tke_froude_loss_check = abs(delta_e_check)/dt
465  tke_froude_loss_tot = sum(cs%TKE_Froude_loss(i,j,:,fr,m))
466  if (abs(tke_froude_loss_check - tke_froude_loss_tot) > 1e-10) then
467  call mom_error(warning, "MOM_internal_tides: something is wrong with Fr energy update.", &
468  all_print=.true.)
469  write(mesg,*) "TKE_Froude_loss_check=", tke_froude_loss_check, &
470  "TKE_Froude_loss_tot=", tke_froude_loss_tot
471  call mom_error(warning, "MOM_internal_tides: "//trim(mesg), all_print=.true.)
472  endif
473  endif ! Fr2>1
474  endif ! Kmag2>0
475  cs%cp(i,j,fr,m) = c_phase
476  enddo ; enddo
477  enddo ; enddo
478  endif
479  ! Check for En<0 - for debugging, delete later
480  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
481  do j=js,je ; do i=is,ie
482  if (cs%En(i,j,a,fr,m)<0.0) then
483  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
484  write(mesg,*) 'After Froude loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
485  'En=',cs%En(i,j,a,fr,m)
486  call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
487  cs%En(i,j,a,fr,m) = 0.0
488 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
489  !stop
490  endif
491  enddo ; enddo
492  enddo ; enddo ; enddo
493 
494  ! Check for energy conservation on computational domain.*************************
495  do m=1,cs%NMode ; do fr=1,cs%Nfreq
496  call sum_en(g,cs,cs%En(:,:,:,fr,m),'prop_int_tide')
497  enddo ; enddo
498 
499  ! Output diagnostics.************************************************************
500  if (query_averaging_enabled(cs%diag)) then
501  ! Output two-dimensional diagnostistics
502  if (cs%id_tot_En > 0) call post_data(cs%id_tot_En, tot_en, cs%diag)
503  if (cs%id_itide_drag > 0) call post_data(cs%id_itide_drag, drag_scale, cs%diag)
504  if (cs%id_TKE_itidal_input > 0) call post_data(cs%id_TKE_itidal_input, &
505  tke_itidal_input, cs%diag)
506 
507  ! Output 2-D energy density (summed over angles) for each freq and mode
508  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_En_mode(fr,m) > 0) then
509  tot_en(:,:) = 0.0
510  do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
511  tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
512  enddo ; enddo ; enddo
513  call post_data(cs%id_En_mode(fr,m), tot_en, cs%diag)
514  endif ; enddo ; enddo
515 
516  ! Output 3-D (i,j,a) energy density for each freq and mode
517  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_En_ang_mode(fr,m) > 0) then
518  call post_data(cs%id_En_ang_mode(fr,m), cs%En(:,:,:,fr,m) , cs%diag)
519  endif ; enddo ; enddo
520 
521  ! Output 2-D energy loss (summed over angles, freq, modes)
522  tot_leak_loss(:,:) = 0.0
523  tot_quad_loss(:,:) = 0.0
524  tot_itidal_loss(:,:) = 0.0
525  tot_froude_loss(:,:) = 0.0
526  tot_allprocesses_loss(:,:) = 0.0
527  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
528  tot_leak_loss(i,j) = tot_leak_loss(i,j) + cs%TKE_leak_loss(i,j,a,fr,m)
529  tot_quad_loss(i,j) = tot_quad_loss(i,j) + cs%TKE_quad_loss(i,j,a,fr,m)
530  tot_itidal_loss(i,j) = tot_itidal_loss(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
531  tot_froude_loss(i,j) = tot_froude_loss(i,j) + cs%TKE_Froude_loss(i,j,a,fr,m)
532  enddo ; enddo ; enddo ; enddo ; enddo
533  do j=js,je ; do i=is,ie
534  tot_allprocesses_loss(i,j) = tot_leak_loss(i,j) + tot_quad_loss(i,j) + &
535  tot_itidal_loss(i,j) + tot_froude_loss(i,j)
536  enddo ; enddo
537  cs%tot_leak_loss = tot_leak_loss
538  cs%tot_quad_loss = tot_quad_loss
539  cs%tot_itidal_loss = tot_itidal_loss
540  cs%tot_Froude_loss = tot_froude_loss
541  cs%tot_allprocesses_loss = tot_allprocesses_loss
542  if (cs%id_tot_leak_loss > 0) then
543  call post_data(cs%id_tot_leak_loss, tot_leak_loss, cs%diag)
544  endif
545  if (cs%id_tot_quad_loss > 0) then
546  call post_data(cs%id_tot_quad_loss, tot_quad_loss, cs%diag)
547  endif
548  if (cs%id_tot_itidal_loss > 0) then
549  call post_data(cs%id_tot_itidal_loss, tot_itidal_loss, cs%diag)
550  endif
551  if (cs%id_tot_Froude_loss > 0) then
552  call post_data(cs%id_tot_Froude_loss, tot_froude_loss, cs%diag)
553  endif
554  if (cs%id_tot_allprocesses_loss > 0) then
555  call post_data(cs%id_tot_allprocesses_loss, tot_allprocesses_loss, cs%diag)
556  endif
557 
558  ! Output 2-D energy loss (summed over angles) for each freq and mode
559  do m=1,cs%NMode ; do fr=1,cs%Nfreq
560  if (cs%id_itidal_loss_mode(fr,m) > 0 .or. cs%id_allprocesses_loss_mode(fr,m) > 0) then
561  itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well)
562  allprocesses_loss_mode(:,:) = 0.0 ! all processes summed together
563  do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
564  itidal_loss_mode(i,j) = itidal_loss_mode(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
565  allprocesses_loss_mode(i,j) = allprocesses_loss_mode(i,j) + &
566  cs%TKE_leak_loss(i,j,a,fr,m) + cs%TKE_quad_loss(i,j,a,fr,m) + &
567  cs%TKE_itidal_loss(i,j,a,fr,m) + cs%TKE_Froude_loss(i,j,a,fr,m)
568  enddo ; enddo ; enddo
569  call post_data(cs%id_itidal_loss_mode(fr,m), itidal_loss_mode, cs%diag)
570  call post_data(cs%id_allprocesses_loss_mode(fr,m), allprocesses_loss_mode, cs%diag)
571  endif ; enddo ; enddo
572 
573  ! Output 3-D (i,j,a) energy loss for each freq and mode
574  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_itidal_loss_ang_mode(fr,m) > 0) then
575  call post_data(cs%id_itidal_loss_ang_mode(fr,m), cs%TKE_itidal_loss(:,:,:,fr,m) , cs%diag)
576  endif ; enddo ; enddo
577 
578  ! Output 2-D period-averaged horizontal near-bottom mode velocity for each freq and mode
579  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_Ub_mode(fr,m) > 0) then
580  call post_data(cs%id_Ub_mode(fr,m), ub(:,:,fr,m), cs%diag)
581  endif ; enddo ; enddo
582 
583  ! Output 2-D horizontal phase velocity for each freq and mode
584  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_cp_mode(fr,m) > 0) then
585  call post_data(cs%id_cp_mode(fr,m), cs%cp(:,:,fr,m), cs%diag)
586  endif ; enddo ; enddo
587 
588  endif
589 

References mom_domains::complete_group_pass(), correct_halo_rotation(), itidal_lowmode_loss(), mom_error_handler::mom_error(), propagate(), refract(), mom_domains::start_group_pass(), sum_en(), and mom_wave_structure::wave_structure().

Here is the call graph for this function:

◆ propagate_x()

subroutine mom_internal_tides::propagate_x ( real, dimension(g%isd:g%ied,g%jsd:g%jed,nangle), intent(inout)  En,
real, dimension(g%isdb:g%iedb,g%jsd:g%jed), intent(in)  speed_x,
real, dimension(nangle), intent(in)  Cgx_av,
real, dimension(nangle), intent(in)  dCgx,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
integer, intent(in)  Nangle,
type(int_tide_cs), pointer  CS,
type(loop_bounds_type), intent(in)  LB 
)
private

Propagates the internal wave energy in the logical x-direction.

Parameters
[in]gThe ocean's grid structure.
[in]nangleThe number of wave orientations in the discretized wave energy spectrum.
[in,out]enThe energy density integrated over an angular
[in]speed_xThe magnitude of the group velocity at the
[in]cgx_avThe average x-projection in each angular band.
[in]dcgxThe difference in x-projections between the edges of each angular band.
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to continuity_PPM_init.
[in]lbA structure with the active energy loop bounds.

Definition at line 1349 of file MOM_internal_tides.F90.

1349  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1350  integer, intent(in) :: NAngle !< The number of wave orientations in the
1351  !! discretized wave energy spectrum.
1352  real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
1353  intent(inout) :: En !< The energy density integrated over an angular
1354  !! band [R Z3 T-2 ~> J m-2], intent in/out.
1355  real, dimension(G%IsdB:G%IedB,G%jsd:G%jed), &
1356  intent(in) :: speed_x !< The magnitude of the group velocity at the
1357  !! Cu points [L T-1 ~> m s-1].
1358  real, dimension(Nangle), intent(in) :: Cgx_av !< The average x-projection in each angular band.
1359  real, dimension(Nangle), intent(in) :: dCgx !< The difference in x-projections between the
1360  !! edges of each angular band.
1361  real, intent(in) :: dt !< Time increment [T ~> s].
1362  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1363  type(int_tide_CS), pointer :: CS !< The control structure returned by a previous call
1364  !! to continuity_PPM_init.
1365  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1366  ! Local variables
1367  real, dimension(SZI_(G),SZJ_(G)) :: &
1368  EnL, EnR ! Left and right face energy densities [R Z3 T-2 ~> J m-2].
1369  real, dimension(SZIB_(G),SZJ_(G)) :: &
1370  flux_x ! The internal wave energy flux [J T-1 ~> J s-1].
1371  real, dimension(SZIB_(G)) :: &
1372  cg_p, cg_m, flux1, flux2
1373  !real, dimension(SZI_(G),SZJB_(G),Nangle) :: En_m, En_p
1374  real, dimension(SZI_(G),SZJB_(G),Nangle) :: &
1375  Fdt_m, Fdt_p! Left and right energy fluxes [J]
1376  integer :: i, j, k, ish, ieh, jsh, jeh, a
1377 
1378  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1379  do a=1,nangle
1380  ! This sets EnL and EnR.
1381  if (cs%upwind_1st) then
1382  do j=jsh,jeh ; do i=ish-1,ieh+1
1383  enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
1384  enddo ; enddo
1385  else
1386  call ppm_reconstruction_x(en(:,:,a), enl, enr, g, lb, simple_2nd=cs%simple_2nd)
1387  endif
1388 
1389  do j=jsh,jeh
1390  ! This is done once with single speed (GARDEN SPRINKLER EFFECT POSSIBLE)
1391  do i=ish-1,ieh
1392  cg_p(i) = speed_x(i,j) * (cgx_av(a))
1393  enddo
1394  call zonal_flux_en(cg_p, en(:,j,a), enl(:,j), enr(:,j), flux1, &
1395  dt, g, us, j, ish, ieh, cs%vol_CFL)
1396  do i=ish-1,ieh ; flux_x(i,j) = flux1(i); enddo
1397  enddo
1398 
1399  do j=jsh,jeh ; do i=ish,ieh
1400  fdt_m(i,j,a) = dt*flux_x(i-1,j) ! left face influx (J)
1401  fdt_p(i,j,a) = -dt*flux_x(i,j) ! right face influx (J)
1402  enddo ; enddo
1403 
1404  enddo ! a-loop
1405 
1406  ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected
1407  ! and will eventually propagate out of cell. (Thid code only reflects if En > 0)
1408  call reflect(fdt_m(:,:,:), nangle, cs, g, lb)
1409  call teleport(fdt_m(:,:,:), nangle, cs, g, lb)
1410  call reflect(fdt_p(:,:,:), nangle, cs, g, lb)
1411  call teleport(fdt_p(:,:,:), nangle, cs, g, lb)
1412 
1413  ! Update reflected energy (Jm-2)
1414  do j=jsh,jeh ; do i=ish,ieh
1415  ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) & ! for debugging
1416  ! call MOM_error(FATAL, "propagate_x: OutFlux>Available")
1417  en(i,j,:) = en(i,j,:) + g%IareaT(i,j)*(fdt_m(i,j,:) + fdt_p(i,j,:))
1418  enddo ; enddo
1419 

References ppm_reconstruction_x(), reflect(), teleport(), and zonal_flux_en().

Referenced by propagate().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ propagate_y()

subroutine mom_internal_tides::propagate_y ( real, dimension(g%isd:g%ied,g%jsd:g%jed,nangle), intent(inout)  En,
real, dimension(g%isd:g%ied,g%jsdb:g%jedb), intent(in)  speed_y,
real, dimension(nangle), intent(in)  Cgy_av,
real, dimension(nangle), intent(in)  dCgy,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
integer, intent(in)  Nangle,
type(int_tide_cs), pointer  CS,
type(loop_bounds_type), intent(in)  LB 
)
private

Propagates the internal wave energy in the logical y-direction.

Parameters
[in]gThe ocean's grid structure.
[in]nangleThe number of wave orientations in the discretized wave energy spectrum.
[in,out]enThe energy density integrated over an angular
[in]speed_yThe magnitude of the group velocity at the
[in]cgy_avThe average y-projection in each angular band.
[in]dcgyThe difference in y-projections between the edges of each angular band.
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
csThe control structure returned by a previous call to continuity_PPM_init.
[in]lbA structure with the active energy loop bounds.

Definition at line 1424 of file MOM_internal_tides.F90.

1424  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1425  integer, intent(in) :: NAngle !< The number of wave orientations in the
1426  !! discretized wave energy spectrum.
1427  real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
1428  intent(inout) :: En !< The energy density integrated over an angular
1429  !! band [R Z3 T-2 ~> J m-2], intent in/out.
1430  real, dimension(G%isd:G%ied,G%JsdB:G%JedB), &
1431  intent(in) :: speed_y !< The magnitude of the group velocity at the
1432  !! Cv points [L T-1 ~> m s-1].
1433  real, dimension(Nangle), intent(in) :: Cgy_av !< The average y-projection in each angular band.
1434  real, dimension(Nangle), intent(in) :: dCgy !< The difference in y-projections between the
1435  !! edges of each angular band.
1436  real, intent(in) :: dt !< Time increment [T ~> s].
1437  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1438  type(int_tide_CS), pointer :: CS !< The control structure returned by a previous call
1439  !! to continuity_PPM_init.
1440  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1441  ! Local variables
1442  real, dimension(SZI_(G),SZJ_(G)) :: &
1443  EnL, EnR ! South and north face energy densities [R Z3 T-2 ~> J m-2].
1444  real, dimension(SZI_(G),SZJB_(G)) :: &
1445  flux_y ! The internal wave energy flux [J T-1 ~> J s-1].
1446  real, dimension(SZI_(G)) :: &
1447  cg_p, cg_m, flux1, flux2
1448  !real, dimension(SZI_(G),SZJB_(G),Nangle) :: En_m, En_p
1449  real, dimension(SZI_(G),SZJB_(G),Nangle) :: &
1450  Fdt_m, Fdt_p! South and north energy fluxes [J]
1451  character(len=160) :: mesg ! The text of an error message
1452  integer :: i, j, k, ish, ieh, jsh, jeh, a
1453 
1454  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1455  do a=1,nangle
1456  ! This sets EnL and EnR.
1457  if (cs%upwind_1st) then
1458  do j=jsh-1,jeh+1 ; do i=ish,ieh
1459  enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
1460  enddo ; enddo
1461  else
1462  call ppm_reconstruction_y(en(:,:,a), enl, enr, g, lb, simple_2nd=cs%simple_2nd)
1463  endif
1464 
1465  do j=jsh-1,jeh
1466  ! This is done once with single speed (GARDEN SPRINKLER EFFECT POSSIBLE)
1467  do i=ish,ieh
1468  cg_p(i) = speed_y(i,j) * (cgy_av(a))
1469  enddo
1470  call merid_flux_en(cg_p, en(:,:,a), enl(:,:), enr(:,:), flux1, &
1471  dt, g, us, j, ish, ieh, cs%vol_CFL)
1472  do i=ish,ieh ; flux_y(i,j) = flux1(i); enddo
1473  enddo
1474 
1475  do j=jsh,jeh ; do i=ish,ieh
1476  fdt_m(i,j,a) = dt*flux_y(i,j-1) ! south face influx (J)
1477  fdt_p(i,j,a) = -dt*flux_y(i,j) ! north face influx (J)
1478  !if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) then ! for debugging
1479  ! call MOM_error(WARNING, "propagate_y: OutFlux>Available prior to reflection", .true.)
1480  ! write(mesg,*) "flux_y_south=",flux_y(i,J-1),"flux_y_north=",flux_y(i,J),"En=",En(i,j,a), &
1481  ! "cn_south=", speed_y(i,J-1) * (Cgy_av(a)), "cn_north=", speed_y(i,J) * (Cgy_av(a))
1482  ! call MOM_error(WARNING, mesg, .true.)
1483  !endif
1484  enddo ; enddo
1485 
1486  enddo ! a-loop
1487 
1488  ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected
1489  ! and will eventually propagate out of cell. (Thid code only reflects if En > 0)
1490  call reflect(fdt_m(:,:,:), nangle, cs, g, lb)
1491  call teleport(fdt_m(:,:,:), nangle, cs, g, lb)
1492  call reflect(fdt_p(:,:,:), nangle, cs, g, lb)
1493  call teleport(fdt_p(:,:,:), nangle, cs, g, lb)
1494 
1495  ! Update reflected energy (Jm-2)
1496  do a=1,nangle ; do j=jsh,jeh ; do i=ish,ieh
1497  ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) & ! for debugging
1498  ! call MOM_error(FATAL, "propagate_y: OutFlux>Available", .true.)
1499  en(i,j,a) = en(i,j,a) + g%IareaT(i,j)*(fdt_m(i,j,a) + fdt_p(i,j,a))
1500  enddo ; enddo ; enddo
1501 

References merid_flux_en(), ppm_reconstruction_y(), reflect(), and teleport().

Referenced by propagate().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ reflect()

subroutine mom_internal_tides::reflect ( real, dimension(g%isd:g%ied,g%jsd:g%jed,nangle), intent(inout)  En,
integer, intent(in)  NAngle,
type(int_tide_cs), pointer  CS,
type(ocean_grid_type), intent(in)  G,
type(loop_bounds_type), intent(in)  LB 
)
private

Reflection of the internal waves at a single frequency.

Parameters
[in]gThe ocean's grid structure
[in]nangleThe number of wave orientations in the discretized wave energy spectrum.
[in,out]enThe internal gravity wave energy density as a
csThe control structure returned by a previous call to int_tide_init.
[in]lbA structure with the active energy loop bounds.

Definition at line 1594 of file MOM_internal_tides.F90.

1594  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1595  integer, intent(in) :: NAngle !< The number of wave orientations in the
1596  !! discretized wave energy spectrum.
1597  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
1598  intent(inout) :: En !< The internal gravity wave energy density as a
1599  !! function of space and angular resolution
1600  !! [R Z3 T-2 ~> J m-2].
1601  type(int_tide_CS), pointer :: CS !< The control structure returned by a
1602  !! previous call to int_tide_init.
1603  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1604  ! Local variables
1605  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
1606  ! angle of boudary wrt equator
1607  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl
1608  ! fraction of wave energy reflected
1609  ! values should collocate with angle_c
1610  logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
1611  ! tags of cells with double reflection
1612 
1613  real :: TwoPi ! 2*pi
1614  real :: Angle_size ! size of beam wedge (rad)
1615  real :: angle_wall ! angle of coast/ridge/shelf wrt equator
1616  real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator
1617  real :: angle_r ! angle of reflected ray wrt equator
1618  real, dimension(1:Nangle) :: En_reflected
1619  integer :: i, j, a, a_r, na
1620  !integer :: isd, ied, jsd, jed ! start and end local indices on data domain
1621  ! ! (values include halos)
1622  integer :: isc, iec, jsc, jec ! start and end local indices on PE
1623  ! (values exclude halos)
1624  integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain
1625  ! leaving out outdated halo points (march in)
1626  integer :: id_g, jd_g ! global (decomp-invar) indices
1627 
1628  !isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1629  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1630  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1631 
1632  twopi = 8.0*atan(1.0)
1633  angle_size = twopi / (real(nangle))
1634 
1635  do a=1,nangle
1636  ! These are the angles at the cell centers
1637  ! (should do this elsewhere since doesn't change with time)
1638  angle_i(a) = angle_size * real(a - 1) ! for a=1 aligned with x-axis
1639  enddo
1640 
1641  angle_c = cs%refl_angle
1642  part_refl = cs%refl_pref
1643  ridge = cs%refl_dbl
1644  en_reflected(:) = 0.0
1645 
1646  !do j=jsc-1,jec+1
1647  do j=jsh,jeh
1648  !do i=isc-1,iec+1
1649  do i=ish,ieh
1650  ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset
1651  ! redistribute energy in angular space if ray will hit boundary
1652  ! i.e., if energy is in a reflecting cell
1653  if (angle_c(i,j) /= cs%nullangle) then
1654  do a=1,nangle
1655  if (en(i,j,a) > 0.0) then
1656  ! if ray is incident, keep specified boundary angle
1657  if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then
1658  angle_wall = angle_c(i,j)
1659  ! if ray is not incident but in ridge cell, use complementary angle
1660  elseif (ridge(i,j)) then
1661  angle_wall = angle_c(i,j) + 0.5*twopi
1662  if (angle_wall > twopi) then
1663  angle_wall = angle_wall - twopi*floor(abs(angle_wall)/twopi)
1664  elseif (angle_wall < 0.0) then
1665  angle_wall = angle_wall + twopi*ceiling(abs(angle_wall)/twopi)
1666  endif
1667  ! if ray is not incident and not in a ridge cell, keep specified angle
1668  else
1669  angle_wall = angle_c(i,j)
1670  endif
1671  ! do reflection
1672  if (sin(angle_i(a) - angle_wall) >= 0.0) then
1673  angle_r = 2.0*angle_wall - angle_i(a)
1674  if (angle_r > twopi) then
1675  angle_r = angle_r - twopi*floor(abs(angle_r)/twopi)
1676  elseif (angle_r < 0.0) then
1677  angle_r = angle_r + twopi*ceiling(abs(angle_r)/twopi)
1678  endif
1679  a_r = nint(angle_r/angle_size) + 1
1680  do while (a_r > nangle) ; a_r = a_r - nangle ; enddo
1681  if (a /= a_r) then
1682  en_reflected(a_r) = part_refl(i,j)*en(i,j,a)
1683  en(i,j,a) = (1.0-part_refl(i,j))*en(i,j,a)
1684  endif
1685  endif
1686  endif
1687  enddo ! a-loop
1688  en(i,j,:) = en(i,j,:) + en_reflected(:)
1689  en_reflected(:) = 0.0
1690  endif
1691  enddo ! i-loop
1692  enddo ! j-loop
1693 
1694  ! Check to make sure no energy gets onto land (only run for debugging)
1695  ! do a=1,NAngle ; do j=jsc,jec ; do i=isc,iec
1696  ! if (En(i,j,a) > 0.001 .and. G%mask2dT(i,j) == 0) then
1697  ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset
1698  ! write (mesg,*) 'En=', En(i,j,a), 'a=', a, 'ig_g=',id_g, 'jg_g=',jd_g
1699  ! call MOM_error(FATAL, "reflect: Energy detected out of bounds: "//trim(mesg), .true.)
1700  ! endif
1701  ! enddo ; enddo ; enddo
1702 

Referenced by propagate_x(), and propagate_y().

Here is the caller graph for this function:

◆ refract()

subroutine mom_internal_tides::refract ( real, dimension(g%isd:g%ied,g%jsd:g%jed,nangle), intent(inout)  En,
real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(in)  cn,
real, intent(in)  freq,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
integer, intent(in)  NAngle,
logical, intent(in)  use_PPMang 
)
private

Implements refraction on the internal waves at a single frequency.

Parameters
[in]gThe ocean's grid structure.
[in]nangleThe number of wave orientations in the discretized wave energy spectrum.
[in,out]enThe internal gravity wave energy density as a
[in]cnBaroclinic mode speed [L T-1 ~> m s-1].
[in]freqWave frequency [T-1 ~> s-1].
[in]dtTime step [T ~> s].
[in]usA dimensional unit scaling type
[in]use_ppmangIf true, use PPM for advection rather than upwind.

Definition at line 746 of file MOM_internal_tides.F90.

746  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
747  integer, intent(in) :: NAngle !< The number of wave orientations in the
748  !! discretized wave energy spectrum.
749  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
750  intent(inout) :: En !< The internal gravity wave energy density as a
751  !! function of space and angular resolution,
752  !! [R Z3 T-2 ~> J m-2].
753  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
754  intent(in) :: cn !< Baroclinic mode speed [L T-1 ~> m s-1].
755  real, intent(in) :: freq !< Wave frequency [T-1 ~> s-1].
756  real, intent(in) :: dt !< Time step [T ~> s].
757  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
758  logical, intent(in) :: use_PPMang !< If true, use PPM for advection rather
759  !! than upwind.
760  ! Local variables
761  integer, parameter :: stencil = 2
762  real, dimension(SZI_(G),1-stencil:NAngle+stencil) :: &
763  En2d
764  real, dimension(1-stencil:NAngle+stencil) :: &
765  cos_angle, sin_angle
766  real, dimension(SZI_(G)) :: &
767  Dk_Dt_Kmag, Dl_Dt_Kmag
768  real, dimension(SZI_(G),0:nAngle) :: &
769  Flux_E
770  real, dimension(SZI_(G),SZJ_(G),1-stencil:NAngle+stencil) :: &
771  CFL_ang
772  real :: f2 ! The squared Coriolis parameter [T-2 ~> s-2].
773  real :: favg ! The average Coriolis parameter at a point [T-1 ~> s-1].
774  real :: df_dy, df_dx ! The x- and y- gradients of the Coriolis parameter [T-1 L-1 ~> s-1 m-1].
775  real :: dlnCn_dx ! The x-gradient of the wave speed divided by itself [L-1 ~> m-1].
776  real :: dlnCn_dy ! The y-gradient of the wave speed divided by itself [L-1 ~> m-1].
777  real :: Angle_size, dt_Angle_size, angle
778  real :: Ifreq, Kmag2, I_Kmag
779  real :: cn_subRO ! A tiny wave speed to prevent division by zero [L T-1 ~> m s-1]
780  integer :: is, ie, js, je, asd, aed, na
781  integer :: i, j, a
782 
783  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na = size(en,3)
784  asd = 1-stencil ; aed = nangle+stencil
785 
786  ifreq = 1.0 / freq
787  cn_subro = 1e-100*us%m_s_to_L_T ! The hard-coded value here might need to increase.
788  angle_size = (8.0*atan(1.0)) / (real(nangle))
789  dt_angle_size = dt / angle_size
790 
791  do a=asd,aed
792  angle = (real(a) - 0.5) * angle_size
793  cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
794  enddo
795 
796  !### There should also be refraction due to cn.grad(grid_orientation).
797  cfl_ang(:,:,:) = 0.0
798  do j=js,je
799  ! Copy En into angle space with halos.
800  do a=1,na ; do i=is,ie
801  en2d(i,a) = en(i,j,a)
802  enddo ; enddo
803  do a=asd,0 ; do i=is,ie
804  en2d(i,a) = en2d(i,a+nangle)
805  en2d(i,nangle+stencil+a) = en2d(i,stencil+a)
806  enddo ; enddo
807 
808  ! Do the refraction.
809  do i=is,ie
810  f2 = 0.25* ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
811  (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
812  favg = 0.25*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
813  (g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j)))
814  df_dx = 0.5*((g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1)) - &
815  (g%CoriolisBu(i-1,j) + g%CoriolisBu(i-1,j-1))) * g%IdxT(i,j)
816  dlncn_dx = 0.5*( g%IdxCu(i,j) * (cn(i+1,j) - cn(i,j)) / &
817  (0.5*(cn(i+1,j) + cn(i,j)) + cn_subro) + &
818  g%IdxCu(i-1,j) * (cn(i,j) - cn(i-1,j)) / &
819  (0.5*(cn(i,j) + cn(i-1,j)) + cn_subro) )
820  df_dy = 0.5*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) - &
821  (g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j-1))) * g%IdyT(i,j)
822  dlncn_dy = 0.5*( g%IdyCv(i,j) * (cn(i,j+1) - cn(i,j)) / &
823  (0.5*(cn(i,j+1) + cn(i,j)) + cn_subro) + &
824  g%IdyCv(i,j-1) * (cn(i,j) - cn(i,j-1)) / &
825  (0.5*(cn(i,j) + cn(i,j-1)) + cn_subro) )
826  kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cn_subro**2)
827  if (kmag2 > 0.0) then
828  i_kmag = 1.0 / sqrt(kmag2)
829  dk_dt_kmag(i) = -ifreq * (favg*df_dx + (freq**2 - f2) * dlncn_dx) * i_kmag
830  dl_dt_kmag(i) = -ifreq * (favg*df_dy + (freq**2 - f2) * dlncn_dy) * i_kmag
831  else
832  dk_dt_kmag(i) = 0.0
833  dl_dt_kmag(i) = 0.0
834  endif
835  enddo
836 
837  ! Determine the energy fluxes in angular orientation space.
838  do a=asd,aed ; do i=is,ie
839  cfl_ang(i,j,a) = (cos_angle(a) * dl_dt_kmag(i) - sin_angle(a) * dk_dt_kmag(i)) * dt_angle_size
840  if (abs(cfl_ang(i,j,a)) > 1.0) then
841  call mom_error(warning, "refract: CFL exceeds 1.", .true.)
842  if (cfl_ang(i,j,a) > 0.0) then ; cfl_ang(i,j,a) = 1.0 ; else ; cfl_ang(i,j,a) = -1.0 ; endif
843  endif
844  enddo ; enddo
845 
846  ! Advect in angular space
847  if (.not.use_ppmang) then
848  ! Use simple upwind
849  do a=0,na ; do i=is,ie
850  if (cfl_ang(i,j,a) > 0.0) then
851  flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a)
852  else
853  flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a+1)
854  endif
855  enddo ; enddo
856  else
857  ! Use PPM
858  do i=is,ie
859  call ppm_angular_advect(en2d(i,:),cfl_ang(i,j,:),flux_e(i,:),nangle,dt,stencil)
860  enddo
861  endif
862 
863  ! Update and copy back to En.
864  do a=1,na ; do i=is,ie
865  !if (En2d(i,a)+(Flux_E(i,A-1)-Flux_E(i,A)) < 0.0) then ! for debugging
866  ! call MOM_error(FATAL, "refract: OutFlux>Available")
867  !endif
868  en(i,j,a) = en2d(i,a) + (flux_e(i,a-1) - flux_e(i,a))
869  enddo ; enddo
870  enddo ! j-loop

References mom_error_handler::mom_error(), and ppm_angular_advect().

Referenced by propagate_int_tide().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ sum_en()

subroutine mom_internal_tides::sum_en ( type(ocean_grid_type), intent(in)  G,
type(int_tide_cs), pointer  CS,
real, dimension(g%isd:g%ied,g%jsd:g%jed,cs%nangle), intent(in)  En,
character(len=*), intent(in)  label 
)
private

Checks for energy conservation on computational domain.

Parameters
[in]gThe ocean's grid structure.
csThe control structure returned by a previous call to int_tide_init.
[in]enThe energy density of the internal tides [R Z3 T-2 ~> J m-2].
[in]labelA label to use in error messages

Definition at line 594 of file MOM_internal_tides.F90.

594  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
595  type(int_tide_CS), pointer :: CS !< The control structure returned by a
596  !! previous call to int_tide_init.
597  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle), &
598  intent(in) :: En !< The energy density of the internal tides [R Z3 T-2 ~> J m-2].
599  character(len=*), intent(in) :: label !< A label to use in error messages
600  ! Local variables
601  real :: En_sum ! The total energy [R Z3 T-2 ~> J m-2]
602  real :: tmpForSumming
603  integer :: m,fr,a
604  ! real :: En_sum_diff, En_sum_pdiff
605  ! character(len=160) :: mesg ! The text of an error message
606  ! real :: days
607 
608  en_sum = 0.0
609  tmpforsumming = 0.0
610  do a=1,cs%nAngle
611  tmpforsumming = global_area_mean(en(:,:,a),g)*g%areaT_global
612  en_sum = en_sum + tmpforsumming
613  enddo
614  cs%En_sum = en_sum
615  !En_sum_diff = En_sum - CS%En_sum
616  !if (CS%En_sum /= 0.0) then
617  ! En_sum_pdiff= (En_sum_diff/CS%En_sum)*100.0
618  !else
619  ! En_sum_pdiff= 0.0
620  !endif
621  !! Print to screen
622  !if (is_root_pe()) then
623  ! days = time_type_to_real(CS%Time) / 86400.0
624  ! write(mesg,*) trim(label)//': days =', days, ', En_sum=', En_sum, &
625  ! ', En_sum_diff=', En_sum_diff, ', Percent change=', En_sum_pdiff, '%'
626  ! call MOM_mesg(mesg)
627  !if (is_root_pe() .and. (abs(En_sum_pdiff) > 1.0)) &
628  ! call MOM_error(FATAL, "Run stopped due to excessive internal tide energy change.")
629  !endif
630 

References mom_spatial_means::global_area_mean().

Referenced by propagate_int_tide().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ teleport()

subroutine mom_internal_tides::teleport ( real, dimension(g%isd:g%ied,g%jsd:g%jed,nangle), intent(inout)  En,
integer, intent(in)  NAngle,
type(int_tide_cs), pointer  CS,
type(ocean_grid_type), intent(in)  G,
type(loop_bounds_type), intent(in)  LB 
)
private

Moves energy across lines of partial reflection to prevent reflection of energy that is supposed to get across.

Parameters
[in]gThe ocean's grid structure.
[in]nangleThe number of wave orientations in the discretized wave energy spectrum.
[in,out]enThe internal gravity wave energy density as a
csThe control structure returned by a previous call to int_tide_init.
[in]lbA structure with the active energy loop bounds.

Definition at line 1708 of file MOM_internal_tides.F90.

1708  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1709  integer, intent(in) :: NAngle !< The number of wave orientations in the
1710  !! discretized wave energy spectrum.
1711  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
1712  intent(inout) :: En !< The internal gravity wave energy density as a
1713  !! function of space and angular resolution
1714  !! [R Z3 T-2 ~> J m-2].
1715  type(int_tide_CS), pointer :: CS !< The control structure returned by a
1716  !! previous call to int_tide_init.
1717  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1718  ! Local variables
1719  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
1720  ! angle of boudary wrt equator
1721  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl
1722  ! fraction of wave energy reflected
1723  ! values should collocate with angle_c
1724  logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: pref_cell
1725  ! flag for partial reflection
1726  logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
1727  ! tags of cells with double reflection
1728  real :: TwoPi ! size of beam wedge (rad)
1729  real :: Angle_size ! size of beam wedge (rad)
1730  real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator
1731  real, dimension(1:NAngle) :: cos_angle, sin_angle
1732  real :: En_tele ! energy to be "teleported" [R Z3 T-2 ~> J m-2]
1733  character(len=160) :: mesg ! The text of an error message
1734  integer :: i, j, a
1735  !integer :: isd, ied, jsd, jed ! start and end local indices on data domain
1736  ! ! (values include halos)
1737  !integer :: isc, iec, jsc, jec ! start and end local indices on PE
1738  ! ! (values exclude halos)
1739  integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain
1740  ! leaving out outdated halo points (march in)
1741  integer :: id_g, jd_g ! global (decomp-invar) indices
1742  integer :: jos, ios ! offsets
1743  real :: cos_normal, sin_normal, angle_wall
1744  ! cos/sin of cross-ridge normal, ridge angle
1745 
1746  !isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1747  !isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
1748  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1749 
1750  twopi = 8.0*atan(1.0)
1751  angle_size = twopi / (real(nangle))
1752 
1753  do a=1,nangle
1754  ! These are the angles at the cell centers
1755  ! (should do this elsewhere since doesn't change with time)
1756  angle_i(a) = angle_size * real(a - 1) ! for a=1 aligned with x-axis
1757  cos_angle(a) = cos(angle_i(a)) ; sin_angle(a) = sin(angle_i(a))
1758  enddo
1759 
1760  angle_c = cs%refl_angle
1761  part_refl = cs%refl_pref
1762  pref_cell = cs%refl_pref_logical
1763  ridge = cs%refl_dbl
1764 
1765  do j=jsh,jeh
1766  do i=ish,ieh
1767  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
1768  if (pref_cell(i,j)) then
1769  do a=1,nangle
1770  if (en(i,j,a) > 0) then
1771  ! if ray is incident, keep specified boundary angle
1772  if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then
1773  angle_wall = angle_c(i,j)
1774  ! if ray is not incident but in ridge cell, use complementary angle
1775  elseif (ridge(i,j)) then
1776  angle_wall = angle_c(i,j) + 0.5*twopi
1777  ! if ray is not incident and not in a ridge cell, keep specified angle
1778  else
1779  angle_wall = angle_c(i,j)
1780  endif
1781  ! teleport if incident
1782  if (sin(angle_i(a) - angle_wall) >= 0.0) then
1783  en_tele = en(i,j,a)
1784  cos_normal = cos(angle_wall + 0.25*twopi)
1785  sin_normal = sin(angle_wall + 0.25*twopi)
1786  ! find preferred zonal offset based on shelf/ridge angle
1787  ios = int(sign(1.,cos_normal))
1788  ! find preferred meridional offset based on shelf/ridge angle
1789  jos = int(sign(1.,sin_normal))
1790  ! find receptive ocean cell in direction of offset
1791  if (.not. pref_cell(i+ios,j+jos)) then
1792  en(i,j,a) = en(i,j,a) - en_tele
1793  en(i+ios,j+jos,a) = en(i+ios,j+jos,a) + en_tele
1794  else
1795  write(mesg,*) 'idg=',id_g,'jd_g=',jd_g,'a=',a
1796  call mom_error(fatal, "teleport: no receptive ocean cell at "//trim(mesg), .true.)
1797  endif
1798  endif ! incidence check
1799  endif ! energy check
1800  enddo ! a-loop
1801  endif ! pref check
1802  enddo ! i-loop
1803  enddo ! j-loop
1804 

References mom_error_handler::mom_error().

Referenced by propagate_x(), and propagate_y().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ zonal_flux_en()

subroutine mom_internal_tides::zonal_flux_en ( real, dimension(szib_(g)), intent(in)  u,
real, dimension(szi_(g)), intent(in)  h,
real, dimension(szi_(g)), intent(in)  hL,
real, dimension(szi_(g)), intent(in)  hR,
real, dimension(szib_(g)), intent(inout)  uh,
real, intent(in)  dt,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
integer, intent(in)  j,
integer, intent(in)  ish,
integer, intent(in)  ieh,
logical, intent(in)  vol_CFL 
)
private

Evaluates the zonal mass or volume fluxes in a layer.

Parameters
[in]gThe ocean's grid structure.
[in]uThe zonal velocity [L T-1 ~> m s-1].
[in]hEnergy density used to calculate the fluxes [R Z3 T-2 ~> J m-2].
[in]hlLeft- Energy densities in the reconstruction [R Z3 T-2 ~> J m-2].
[in]hrRight- Energy densities in the reconstruction [R Z3 T-2 ~> J m-2].
[in,out]uhThe zonal energy transport [R Z3 L2 T-3 ~> J s-1].
[in]dtTime increment [T ~> s].
[in]usA dimensional unit scaling type
[in]jThe j-index to work on.
[in]ishThe start i-index range to work on.
[in]iehThe end i-index range to work on.
[in]vol_cflIf true, rescale the ratio of face areas to the cell areas when estimating the CFL number.

Definition at line 1506 of file MOM_internal_tides.F90.

1506  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1507  real, dimension(SZIB_(G)), intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
1508  real, dimension(SZI_(G)), intent(in) :: h !< Energy density used to calculate the fluxes
1509  !! [R Z3 T-2 ~> J m-2].
1510  real, dimension(SZI_(G)), intent(in) :: hL !< Left- Energy densities in the reconstruction
1511  !! [R Z3 T-2 ~> J m-2].
1512  real, dimension(SZI_(G)), intent(in) :: hR !< Right- Energy densities in the reconstruction
1513  !! [R Z3 T-2 ~> J m-2].
1514  real, dimension(SZIB_(G)), intent(inout) :: uh !< The zonal energy transport [R Z3 L2 T-3 ~> J s-1].
1515  real, intent(in) :: dt !< Time increment [T ~> s].
1516  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1517  integer, intent(in) :: j !< The j-index to work on.
1518  integer, intent(in) :: ish !< The start i-index range to work on.
1519  integer, intent(in) :: ieh !< The end i-index range to work on.
1520  logical, intent(in) :: vol_CFL !< If true, rescale the ratio of face areas to
1521  !! the cell areas when estimating the CFL number.
1522  ! Local variables
1523  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim].
1524  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1525  ! with the same units as h_in.
1526  integer :: i
1527 
1528  do i=ish-1,ieh
1529  ! Set new values of uh and duhdu.
1530  if (u(i) > 0.0) then
1531  if (vol_cfl) then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
1532  else ; cfl = u(i) * dt * g%IdxT(i,j) ; endif
1533  curv_3 = (hl(i) + hr(i)) - 2.0*h(i)
1534  uh(i) = g%dy_Cu(i,j) * u(i) * &
1535  (hr(i) + cfl * (0.5*(hl(i) - hr(i)) + curv_3*(cfl - 1.5)))
1536  elseif (u(i) < 0.0) then
1537  if (vol_cfl) then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
1538  else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ; endif
1539  curv_3 = (hl(i+1) + hr(i+1)) - 2.0*h(i+1)
1540  uh(i) = g%dy_Cu(i,j) * u(i) * &
1541  (hl(i+1) + cfl * (0.5*(hr(i+1)-hl(i+1)) + curv_3*(cfl - 1.5)))
1542  else
1543  uh(i) = 0.0
1544  endif
1545  enddo

Referenced by propagate_x().

Here is the caller graph for this function: