MOM_wave_interface.F90

!> Interface for surface waves
module mom_wave_interface

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

use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc
use mom_diag_mediator, only : diag_ctrl
use mom_domains,       only : pass_var, pass_vector, agrid
use mom_domains,       only : to_south, to_west, to_all
use mom_error_handler, only : mom_error, fatal, warning
use mom_file_parser,   only : get_param, log_version, param_file_type
use mom_grid,          only : ocean_grid_type
use mom_safe_alloc,    only : safe_alloc_ptr
use mom_time_manager,  only : time_type, operator(+), operator(/)
use mom_unit_scaling,  only : unit_scale_type
use mom_variables,     only : thermo_var_ptrs, surface
use mom_verticalgrid,  only : verticalgrid_type
use data_override_mod, only : data_override_init, data_override

implicit none ; private

#include <MOM_memory.h>

public mom_wave_interface_init ! Public interface to fully initialize the wave routines.
public mom_wave_interface_init_lite ! Public interface to quick initialize this module.
public update_surface_waves ! Public interface to update wave information at the
                            ! coupler/driver level.
public update_stokes_drift ! Public interface to update the Stokes drift profiles
                           ! called in step_mom.
public get_langmuir_number ! Public interface to compute Langmuir number called from
                           ! ePBL or KPP routines.
public stokesmixing ! NOT READY - Public interface to add down-Stokes gradient
                    ! momentum mixing (e.g. the approach of Harcourt 2013/2015)
public coriolisstokes ! NOT READY - Public interface to add Coriolis-Stokes acceleration
                      ! of the mean currents, needed for comparison with LES.  It is
                      ! presently advised against implementing in non-1d settings without
                      ! serious consideration of the full 3d wave-averaged Navier-Stokes
                      ! CL2 effects.
public waves_end ! public interface to deallocate and free wave related memory.

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

!> Container for all surface wave related parameters
type, public :: wave_parameters_cs ; private

  !Main surface wave options
  logical, public :: usewaves         !< Flag to enable surface gravity wave feature
  logical, public :: lagrangianmixing !< This feature is in development and not ready
                                      !! True if Stokes drift is present and mixing
                                      !! should be applied to Lagrangian current
                                      !! (mean current + Stokes drift).
                                      !! See Reichl et al., 2016 KPP-LT approach
  logical, public :: stokesmixing     !< This feature is in development and not ready.
                                      !! True if vertical mixing of momentum
                                      !! should be applied directly to Stokes current
                                      !! (with separate mixing parameter for Eulerian
                                      !! mixing contribution).
                                      !! See Harcourt 2013, 2015 Second-Moment approach
  logical, public :: coriolisstokes   !< This feature is in development and not ready.
                                      ! True if Coriolis-Stokes acceleration should be applied.
  integer, public :: stklevelmode=1   !< Sets if Stokes drift is defined at mid-points
                                      !! or layer averaged.  Set to 0 if mid-point and set to
                                      !! 1 if average value of Stokes drift over level.
                                      !! If advecting with Stokes transport, 1 is the correct
                                      !! approach.

  ! Surface Wave Dependent 1d/2d/3d vars
  real, allocatable, dimension(:), public :: &
       wavenum_cen        !< Wavenumber bands for read/coupled [m-1]
  real, allocatable, dimension(:), public :: &
       freq_cen           !< Frequency bands for read/coupled [s-1]
  real, allocatable, dimension(:), public :: &
       prescribedsurfstkx !< Surface Stokes drift if prescribed [m s-1]
  real, allocatable, dimension(:), public :: &
       prescribedsurfstky !< Surface Stokes drift if prescribed [m s-1]
  real, allocatable, dimension(:,:,:), public :: &
       us_x               !< 3d zonal Stokes drift profile [m s-1]
                          !! Horizontal -> U points
                          !! Vertical -> Mid-points
  real, allocatable, dimension(:,:,:), public :: &
       us_y               !< 3d meridional Stokes drift profile [m s-1]
                          !! Horizontal -> V points
                          !! Vertical -> Mid-points
  real, allocatable, dimension(:,:), public :: &
       la_sl,&            !< SL Langmuir number (directionality factored later)
                          !! Horizontal -> H points
       la_turb            !< Aligned Turbulent Langmuir number
                          !! Horizontal -> H points
  real, allocatable, dimension(:,:), public :: &
       us0_x              !< Surface Stokes Drift (zonal, m/s)
                          !! Horizontal -> U points
  real, allocatable, dimension(:,:), public :: &
       us0_y              !< Surface Stokes Drift (meridional, m/s)
                          !! Horizontal -> V points
  real, allocatable, dimension(:,:,:), public :: &
       stkx0              !< Stokes Drift spectrum (zonal, m/s)
                          !! Horizontal -> U points
                          !! 3rd dimension -> Freq/Wavenumber
  real, allocatable, dimension(:,:,:), public :: &
       stky0              !< Stokes Drift spectrum (meridional, m/s)
                          !! Horizontal -> V points
                          !! 3rd dimension -> Freq/Wavenumber
  real, allocatable, dimension(:,:,:), public :: &
       kvs                !< Viscosity for Stokes Drift shear [Z2 T-1 ~> m2 s-1]

  ! Pointers to auxiliary fields
  type(time_type), pointer, public :: time !< A pointer to the ocean model's clock.
  type(diag_ctrl), pointer, public :: diag !< A structure that is used to regulate the
                                           !! timing of diagnostic output.

  !> An arbitrary lower-bound on the Langmuir number.  Run-time parameter.
  !! Langmuir number is sqrt(u_star/u_stokes). When both are small
  !! but u_star is orders of magnitude smaller the Langmuir number could
  !! have unintended consequences.  Since both are small it can be safely capped
  !! to avoid such consequences.
  real :: la_min = 0.05

  !>@{ Diagnostic handles
  integer, public :: id_surfacestokes_x = -1 , id_surfacestokes_y = -1
  integer, public :: id_3dstokes_x = -1 , id_3dstokes_y = -1
  integer, public :: id_la_turb = -1
  !>@}

end type wave_parameters_cs

! Options not needed outside of this module

integer :: wavemethod=-99 !< Options for including wave information
                          !! Valid (tested) choices are:
                          !!   0 - Test Profile
                          !!   1 - Surface Stokes Drift Bands
                          !!   2 - DHH85
                          !!   3 - LF17
                          !! -99 - No waves computed, but empirical Langmuir number used.
                          !! \todo Module variable! Move into a control structure.

! Options if WaveMethod is Surface Stokes Drift Bands (1)
integer, public :: numbands =0 !< Number of wavenumber/frequency partitions to receive
                               !! This needs to match the number of bands provided
                               !! via either coupling or file.
                               !! \todo Module variable! Move into a control structure.
integer, public :: partitionmode !< Method for partition mode (meant to check input)
                                 !! 0 - wavenumbers
                                 !! 1 - frequencies
                                 !! \todo Module variable! Move into a control structure.
integer :: datasource !< Integer that specifies where the Model Looks for Data
                      !! Valid choices are:
                      !! 1 - FMS DataOverride Routine
                      !! 2 - Reserved For Coupler
                      !! 3 - User input (fixed values, useful for 1d testing)
                      !! \todo Module variable! Move into a control structure.

! Options if using FMS DataOverride Routine
character(len=40)  :: surfbandfilename !< Filename if using DataOverride
                                       !! \todo Module variable! Move into a control structure.
logical :: dataoverrideisinitialized !< Flag for DataOverride Initialization
                                     !! \todo Module variable! Move into a control structure.

! Options for computing Langmuir number
real :: la_frachbl         !< Fraction of OSBL for averaging Langmuir number
                           !! \todo Module variable! Move into a control structure.
logical :: la_misalignment = .false. !< Flag to use misalignment in Langmuir number
                                     !! \todo Module variable! Move into a control structure.

! This include declares and sets the variable "version".
#include "version_variable.h"

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

!>@{ Undocumented parameters.
!! \todo These module variables need to be documented as static/private variables or moved
!! into a control structure.
! Switches needed in import_stokes_drift
integer, parameter :: testprof = 0, surfbands = 1, &
                      dhh85 = 2, lf17 = 3, null_wavemethod=-99, &
                      dataovr = 1, coupler = 2, input = 3

! Options For Test Prof
Real    :: tp_stkx0, tp_stky0, tp_wvl
logical :: waveagepeakfreq ! Flag to use W
logical :: staticwaves, dhh85_is_set
real    :: waveage, wavewind
real    :: pi
!>@}

contains

!> Initializes parameters related to MOM_wave_interface
subroutine mom_wave_interface_init(time, G, GV, US, param_file, CS, diag )
  type(time_type), target, intent(in)    :: Time       !< Model time
  type(ocean_grid_type),   intent(inout) :: G          !< Grid structure
  type(verticalgrid_type), intent(in)    :: GV         !< Vertical grid structure
  type(unit_scale_type),   intent(in)    :: US         !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< Input parameter structure
  type(wave_parameters_cs), pointer      :: CS         !< Wave parameter control structure
  type(diag_ctrl), target, intent(inout) :: diag       !< Diagnostic Pointer
  ! Local variables
  ! I/O
  character*(13) :: TMPSTRING1,TMPSTRING2
  character*(5), parameter  :: NULL_STRING      = "EMPTY"
  character*(12), parameter :: TESTPROF_STRING  = "TEST_PROFILE"
  character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS"
  character*(5), parameter  :: DHH85_STRING     = "DHH85"
  character*(4), parameter  :: LF17_STRING      = "LF17"
  character*(12), parameter :: DATAOVR_STRING   = "DATAOVERRIDE"
  character*(7), parameter  :: COUPLER_STRING   = "COUPLER"
  character*(5), parameter  :: INPUT_STRING     = "INPUT"

  ! Dummy Check
  if (associated(cs)) then
     call mom_error(fatal, "wave_interface_init called with an associated"//&
                             "control structure.")
     return
  endif

  pi=4.0*atan(1.0)

  ! Allocate CS and set pointers
  allocate(cs)

  cs%diag => diag
  cs%Time => time

  ! Add any initializations needed here
  dataoverrideisinitialized = .false.

  ! The only way to get here is with UseWaves enabled.
  cs%UseWaves=.true.

  call log_version(param_file, mdl, version)

  ! Wave modified physics
  !  Presently these are all in research mode
  call get_param(param_file, mdl, "LAGRANGIAN_MIXING", cs%LagrangianMixing, &
       "Flag to use Lagrangian Mixing of momentum", units="", &
       default=.false.)
  if (cs%LagrangianMixing) then
    ! Force Code Intervention
    call mom_error(fatal,"Should you be enabling Lagrangian Mixing? Code not ready.")
  endif
  call get_param(param_file, mdl, "STOKES_MIXING", cs%StokesMixing, &
       "Flag to use Stokes Mixing of momentum", units="", &
       default=.false.)
  if (cs%StokesMixing) then
    ! Force Code Intervention
    call mom_error(fatal,"Should you be enabling Stokes Mixing? Code not ready.")
  endif
  call get_param(param_file, mdl, "CORIOLIS_STOKES", cs%CoriolisStokes, &
       "Flag to use Coriolis Stokes acceleration", units="", &
       default=.false.)
  if (cs%CoriolisStokes) then
    ! Force Code Intervention
    call mom_error(fatal,"Should you be enabling Coriolis-Stokes? Code not ready.")
  endif

  ! Get Wave Method and write to integer WaveMethod
  call get_param(param_file,mdl,"WAVE_METHOD",tmpstring1,             &
       "Choice of wave method, valid options include: \n"//           &
       " TEST_PROFILE  - Prescribed from surface Stokes drift \n"//   &
       "                 and a decay wavelength.\n"//                 &
       " SURFACE_BANDS - Computed from multiple surface values \n"//  &
       "                 and decay wavelengths.\n"//                  &
       " DHH85         - Uses Donelan et al. 1985 empirical \n"//     &
       "                 wave spectrum with prescribed values. \n"//  &
       " LF17          - Infers Stokes drift profile from wind \n"//  &
       "                 speed following Li and Fox-Kemper 2017.\n",  &
       units='', default=null_string)
  select case (trim(tmpstring1))
  case (null_string)! No Waves
    call mom_error(fatal, "wave_interface_init called with no specified "//&
                           "WAVE_METHOD.")
  case (testprof_string)! Test Profile
    wavemethod = testprof
    call get_param(param_file,mdl,"TP_STKX_SURF",tp_stkx0,&
         'Surface Stokes (x) for test profile',&
         units='m/s',default=0.1)
    call get_param(param_file,mdl,"TP_STKY_SURF",tp_stky0,&
         'Surface Stokes (y) for test profile',&
         units='m/s',default=0.0)
    call get_param(param_file,mdl,"TP_WVL",tp_wvl,&
         units='m', default=50.0, scale=us%m_to_Z)
  case (surfbands_string)! Surface Stokes Drift Bands
    wavemethod = surfbands
    call get_param(param_file, mdl, "SURFBAND_SOURCE",tmpstring2,       &
       "Choice of SURFACE_BANDS data mode, valid options include: \n"// &
       " DATAOVERRIDE  - Read from NetCDF using FMS DataOverride. \n"// &
       " COUPLER       - Look for variables from coupler pass \n"//     &
       " INPUT         - Testing with fixed values.",                   &
       units='', default=null_string)
    select case (trim(tmpstring2))
    case (null_string)! Default
      call mom_error(fatal, "wave_interface_init called with SURFACE_BANDS"//&
                           " but no SURFBAND_SOURCE.")
    case (dataovr_string)! Using Data Override
      datasource = dataovr
      call get_param(param_file, mdl, "SURFBAND_FILENAME", surfbandfilename, &
           "Filename of surface Stokes drift input band data.", default="StkSpec.nc")
    case (coupler_string)! Reserved for coupling
      datasource = coupler
    case (input_string)! A method to input the Stokes band (globally uniform)
      datasource = input
      call get_param(param_file,mdl,"SURFBAND_NB",numbands,                 &
         "Prescribe number of wavenumber bands for Stokes drift. "//      &
         "Make sure this is consistnet w/ WAVENUMBERS, STOKES_X, and "// &
         "STOKES_Y, there are no safety checks in the code.",              &
         units='', default=1)
      allocate( cs%WaveNum_Cen(1:numbands) )
      cs%WaveNum_Cen(:) = 0.0
      allocate( cs%PrescribedSurfStkX(1:numbands))
      cs%PrescribedSurfStkX(:) = 0.0
      allocate( cs%PrescribedSurfStkY(1:numbands))
      cs%PrescribedSurfStkY(:) = 0.0
      allocate( cs%STKx0(g%isdB:g%iedB,g%jsd:g%jed,1:numbands))
      cs%STKx0(:,:,:) = 0.0
      allocate( cs%STKy0(g%isd:g%ied,g%jsdB:g%jedB,1:numbands))
      cs%STKy0(:,:,:) = 0.0
      partitionmode=0
      call get_param(param_file,mdl,"SURFBAND_WAVENUMBERS",cs%WaveNum_Cen,      &
           "Central wavenumbers for surface Stokes drift bands.",units='rad/m', &
           default=0.12566)
      call get_param(param_file,mdl,"SURFBAND_STOKES_X",cs%PrescribedSurfStkX,      &
           "X-direction surface Stokes drift for bands.",units='m/s', &
           default=0.15)
      call get_param(param_file,mdl,"SURFBAND_STOKES_Y",cs%PrescribedSurfStkY,      &
           "Y-direction surface Stokes drift for bands.",units='m/s', &
           default=0.0)
   case default! No method provided
     call mom_error(fatal,'Check WAVE_METHOD.')
   end select

  case (dhh85_string)!Donelan et al., 1985 spectrum
    wavemethod = dhh85
    call mom_error(warning,"DHH85 only ever set-up for uniform cases w/"//&
                           " Stokes drift in x-direction.")
    call get_param(param_file,mdl,"DHH85_AGE_FP",waveagepeakfreq,   &
         "Choose true to use waveage in peak frequency.", &
         units='', default=.false.)
    call get_param(param_file,mdl,"DHH85_AGE",waveage,   &
         "Wave Age for DHH85 spectrum.", &
         units='', default=1.2)
    call get_param(param_file,mdl,"DHH85_WIND",wavewind,   &
         "Wind speed for DHH85 spectrum.", &
         units='', default=10.0)
    call get_param(param_file,mdl,"STATIC_DHH85",staticwaves,   &
         "Flag to disable updating DHH85 Stokes drift.", &
          default=.false.)
  case (lf17_string)!Li and Fox-Kemper 17 wind-sea Langmuir number
    wavemethod = lf17
   case default
     call mom_error(fatal,'Check WAVE_METHOD.')
  end select

  ! Langmuir number Options
  call get_param(param_file, mdl, "LA_DEPTH_RATIO", la_frachbl,              &
         "The depth (normalized by BLD) to average Stokes drift over in "//&
         "Langmuir number calculation, where La = sqrt(ust/Stokes).",       &
         units="nondim",default=0.04)
  call get_param(param_file, mdl, "LA_MISALIGNMENT", la_misalignment,    &
         "Flag (logical) if using misalignment bt shear and waves in LA",&
         default=.false.)
  call get_param(param_file, mdl, "MIN_LANGMUIR", cs%La_min,    &
         "A minimum value for all Langmuir numbers that is not physical, "//&
         "but is likely only encountered when the wind is very small and "//&
         "therefore its effects should be mostly benign.",units="nondim",&
         default=0.05)

  ! Allocate and initialize
  ! a. Stokes driftProfiles
  allocate(cs%Us_x(g%isdB:g%IedB,g%jsd:g%jed,g%ke))
  cs%Us_x(:,:,:) = 0.0
  allocate(cs%Us_y(g%isd:g%Ied,g%jsdB:g%jedB,g%ke))
  cs%Us_y(:,:,:) = 0.0
  ! b. Surface Values
  allocate(cs%US0_x(g%isdB:g%iedB,g%jsd:g%jed))
  cs%US0_x(:,:) = 0.0
  allocate(cs%US0_y(g%isd:g%ied,g%jsdB:g%jedB))
  cs%US0_y(:,:) = 0.0
  ! c. Langmuir number
  allocate(cs%La_SL(g%isc:g%iec,g%jsc:g%jec))
  allocate(cs%La_turb(g%isc:g%iec,g%jsc:g%jec))
  cs%La_SL(:,:) = 0.0
  cs%La_turb (:,:) = 0.0
  ! d. Viscosity for Stokes drift
  if (cs%StokesMixing) then
    allocate(cs%KvS(g%isd:g%Ied,g%jsd:g%jed,g%ke))
    cs%KvS(:,:,:) = 0.0
  endif

  ! Initialize Wave related outputs
  cs%id_surfacestokes_y = register_diag_field('ocean_model','surface_stokes_y', &
       cs%diag%axesCu1,time,'Surface Stokes drift (y)','m s-1')
  cs%id_surfacestokes_x = register_diag_field('ocean_model','surface_stokes_x', &
       cs%diag%axesCv1,time,'Surface Stokes drift (x)','m s-1')
  cs%id_3dstokes_y = register_diag_field('ocean_model','3d_stokes_y', &
       cs%diag%axesCvL,time,'3d Stokes drift (y)','m s-1')
  cs%id_3dstokes_x = register_diag_field('ocean_model','3d_stokes_x', &
       cs%diag%axesCuL,time,'3d Stokes drift (y)','m s-1')
  cs%id_La_turb = register_diag_field('ocean_model','La_turbulent',&
       cs%diag%axesT1,time,'Surface (turbulent) Langmuir number','nondim')

  return
end subroutine mom_wave_interface_init

!> A 'lite' init subroutine to initialize a few inputs needed if using wave information
!! with the wind-speed dependent Stokes drift formulation of LF17
subroutine mom_wave_interface_init_lite(param_file)
  type(param_file_type), intent(in) :: param_file !< Input parameter structure
  character*(5), parameter  :: NULL_STRING      = "EMPTY"
  character*(4), parameter  :: LF17_STRING      = "LF17"
  character*(13) :: TMPSTRING1
  logical :: StatisticalWaves

  ! Langmuir number Options
  call get_param(param_file, mdl, "LA_DEPTH_RATIO", la_frachbl,              &
       "The depth (normalized by BLD) to average Stokes drift over in "//&
       "Langmuir number calculation, where La = sqrt(ust/Stokes).",       &
       units="nondim",default=0.04)

  ! Check if using LA_LI2016
  call get_param(param_file,mdl,"USE_LA_LI2016",statisticalwaves,     &
                 do_not_log=.true.,default=.false.)
  if (statisticalwaves) then
    wavemethod = lf17
    pi=4.0*atan(1.0)
  else
    wavemethod = null_wavemethod
  end if

  return
end subroutine mom_wave_interface_init_lite

!> Subroutine that handles updating of surface wave/Stokes drift related properties
subroutine update_surface_waves(G, GV, US, Day, dt, CS)
  type(wave_parameters_cs), pointer    :: CS  !< Wave parameter Control structure
  type(ocean_grid_type), intent(inout) :: G   !< Grid structure
  type(verticalgrid_type), intent(in)  :: GV  !< Vertical grid structure
  type(unit_scale_type),   intent(in)  :: US   !< A dimensional unit scaling type
  type(time_type),         intent(in)  :: Day !< Current model time
  type(time_type),         intent(in)  :: dt  !< Timestep as a time-type
  ! Local variables
  integer :: ii, jj, kk, b
  type(time_type) :: Day_Center

  ! Computing central time of time step
  day_center = day + dt/2

  if (wavemethod == testprof) then
    ! Do nothing
  elseif (wavemethod==surfbands) then
    if (datasource==dataovr) then
      call surface_bands_by_data_override(day_center, g, gv, us, cs)
    elseif (datasource==coupler) then
      ! Reserve for coupler hooks
    elseif (datasource==input) then
      do b=1,numbands
        do ii=g%isdB,g%iedB
          do jj=g%jsd,g%jed
            cs%STKx0(ii,jj,b) = cs%PrescribedSurfStkX(b)
          enddo
        enddo
        do ii=g%isd,g%ied
          do jj=g%jsdB, g%jedB
            cs%STKY0(ii,jj,b) = cs%PrescribedSurfStkY(b)
          enddo
        enddo
      enddo
    endif
  endif

  return
end subroutine update_surface_waves

!> Constructs the Stokes Drift profile on the model grid based on
!! desired coupling options
subroutine update_stokes_drift(G, GV, US, CS, h, ustar)
  type(wave_parameters_cs),  pointer     :: CS    !< Wave parameter Control structure
  type(ocean_grid_type),   intent(inout) :: G     !< Grid structure
  type(verticalgrid_type), intent(in)    :: GV    !< Vertical grid structure
  type(unit_scale_type),   intent(in)    :: US    !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
       intent(in)    :: h     !< Thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G)), &
       intent(in)    :: ustar !< Wind friction velocity [Z T-1 ~> m s-1].
  ! Local Variables
  real    :: Top, MidPoint, Bottom, one_cm
  real    :: DecayScale
  real    :: CMN_FAC, WN, UStokes
  real    :: La
  integer :: ii, jj, kk, b, iim1, jjm1

  one_cm = 0.01*us%m_to_Z

  ! 1. If Test Profile Option is chosen
  !    Computing mid-point value from surface value and decay wavelength
  if (wavemethod==testprof) then
    decayscale = 4.*pi / tp_wvl !4pi
    do ii = g%isdB,g%iedB
      do jj = g%jsd,g%jed
        iim1 = max(1,ii-1)
        bottom = 0.0
        midpoint = 0.0
        do kk = 1,g%ke
          top = bottom
          midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(iim1,jj,kk))
          bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(iim1,jj,kk))
          cs%Us_x(ii,jj,kk) = tp_stkx0*exp(midpoint*decayscale)
        enddo
      enddo
    enddo
    do ii = g%isd,g%ied
      do jj = g%jsdB,g%jedB
        jjm1 = max(1,jj-1)
        bottom = 0.0
        midpoint = 0.0
        do kk = 1,g%ke
          top = bottom
          midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(ii,jjm1,kk))
          bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(ii,jjm1,kk))
          cs%Us_y(ii,jj,kk) = tp_stky0*exp(midpoint*decayscale)
        enddo
      enddo
    enddo
  ! 2. If Surface Bands is chosen
  !    In wavenumber mode compute integral for layer averaged Stokes drift.
  !    In frequency mode compuate value at midpoint.
  elseif (wavemethod==surfbands) then
    cs%Us_x(:,:,:) = 0.0
    cs%Us_y(:,:,:) = 0.0
    cs%Us0_x(:,:) = 0.0
    cs%Us0_y(:,:) = 0.0
    ! Computing X direction Stokes drift
    do ii = g%isdB,g%iedB
      do jj = g%jsd,g%jed
        ! 1. First compute the surface Stokes drift
        !    by integrating over the partitionas.
        do b = 1,numbands
          if (partitionmode==0) then
            ! In wavenumber we are averaging over (small) level
            cmn_fac = (1.0-exp(-one_cm*2.*cs%WaveNum_Cen(b))) / &
                      (one_cm*2.*cs%WaveNum_Cen(b))
          elseif (partitionmode==1) then
             ! In frequency we are not averaging over level and taking top
            cmn_fac = 1.0
          endif
          cs%US0_x(ii,jj) = cs%US0_x(ii,jj) + cs%STKx0(ii,jj,b)*cmn_fac
        enddo
        ! 2. Second compute the level averaged Stokes drift
        bottom = 0.0
        do kk = 1,g%ke
          top = bottom
          iim1 = max(ii-1,1)
          midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(iim1,jj,kk))
          bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(iim1,jj,kk))
          do b = 1,numbands
            if (partitionmode==0) then
              ! In wavenumber we are averaging over level
              cmn_fac = (exp(top*2.*cs%WaveNum_Cen(b))-exp(bottom*2.*cs%WaveNum_Cen(b)))&
                        / ((top-bottom)*(2.*cs%WaveNum_Cen(b)))
            elseif (partitionmode==1) then
              if (cs%StkLevelMode==0) then
                ! Take the value at the midpoint
                cmn_fac = exp(midpoint*2.*(2.*pi*cs%Freq_Cen(b)*us%T_to_s)**2/(us%L_to_Z**2*gv%g_Earth))
              elseif (cs%StkLevelMode==1) then
                ! Use a numerical integration and then
                ! divide by layer thickness
                wn = (2.*pi*cs%Freq_Cen(b)*us%T_to_s)**2 / (us%L_to_Z**2*gv%g_Earth) !bgr bug-fix missing g
                cmn_fac = (exp(2.*wn*top)-exp(2.*wn*bottom)) / (2.*wn*(top-bottom))
              endif
            endif
            cs%US_x(ii,jj,kk) = cs%US_x(ii,jj,kk) + cs%STKx0(ii,jj,b)*cmn_fac
          enddo
        enddo
      enddo
    enddo
    ! Computing Y direction Stokes drift
    do ii = g%isd,g%ied
      do jj = g%jsdB,g%jedB
        ! Compute the surface values.
        do b = 1,numbands
          if (partitionmode==0) then
            ! In wavenumber we are averaging over (small) level
            cmn_fac = (1.0-exp(-one_cm*2.*cs%WaveNum_Cen(b))) / &
                      (one_cm*2.*cs%WaveNum_Cen(b))
          elseif (partitionmode==1) then
            ! In frequency we are not averaging over level and taking top
            cmn_fac = 1.0
          endif
          cs%US0_y(ii,jj) = cs%US0_y(ii,jj) + cs%STKy0(ii,jj,b)*cmn_fac
        enddo
        ! Compute the level averages.
        bottom = 0.0
        do kk = 1,g%ke
          top = bottom
          jjm1 = max(jj-1,1)
          midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(ii,jjm1,kk))
          bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(ii,jjm1,kk))
          do b = 1,numbands
            if (partitionmode==0) then
              ! In wavenumber we are averaging over level
              cmn_fac = (exp(top*2.*cs%WaveNum_Cen(b)) - &
                         exp(bottom*2.*cs%WaveNum_Cen(b))) / &
                        ((top-bottom)*(2.*cs%WaveNum_Cen(b)))
            elseif (partitionmode==1) then
              if (cs%StkLevelMode==0) then
                ! Take the value at the midpoint
                cmn_fac = exp(midpoint*2.*(2.*pi*cs%Freq_Cen(b)*us%T_to_s)**2/(us%L_to_Z**2*gv%g_Earth))
              elseif (cs%StkLevelMode==1) then
                ! Use a numerical integration and then
                ! divide by layer thickness
                wn = (2.*pi*cs%Freq_Cen(b)*us%T_to_s)**2 / (us%L_to_Z**2*gv%g_Earth)
                cmn_fac = (exp(2.*wn*top)-exp(2.*wn*bottom)) / (2.*wn*(top-bottom))
              endif
            endif
            cs%US_y(ii,jj,kk) = cs%US_y(ii,jj,kk) + cs%STKy0(ii,jj,b)*cmn_fac
          enddo
        enddo
      enddo
    enddo
  elseif (wavemethod==dhh85) then
    if (.not.(staticwaves .and. dhh85_is_set)) then
      do ii = g%isdB,g%iedB
        do jj = g%jsd,g%jed
          bottom = 0.0
          do kk = 1,g%ke
            top = bottom
            iim1 = max(ii-1,1)
            midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(iim1,jj,kk))
            bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(iim1,jj,kk))
            !bgr note that this is using a u-point ii on h-point ustar
            !    this code has only been previous used for uniform
            !    grid cases.  This needs fixed if DHH85 is used for non
            !    uniform cases.
            call dhh85_mid(gv, us, midpoint, ustokes)
            ! Putting into x-direction (no option for direction
            cs%US_x(ii,jj,kk) = ustokes
          enddo
        enddo
      enddo
      do ii = g%isd,g%ied
        do jj = g%jsdB,g%jedB
          bottom = 0.0
          do kk=1, g%ke
            top = bottom
            jjm1 = max(jj-1,1)
            midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(ii,jjm1,kk))
            bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(ii,jjm1,kk))
            !bgr note that this is using a v-point jj on h-point ustar
            !    this code has only been previous used for uniform
            !    grid cases.  This needs fixed if DHH85 is used for non
            !    uniform cases.
            ! call DHH85_mid(GV, US, Midpoint, UStokes)
            ! Putting into x-direction, so setting y direction to 0
            cs%US_y(ii,jj,kk) = 0.0
            ! For rotational symmetry there should be the option for this to become = UStokes
            !    bgr - see note above, but this is true
            !          if this is used for anything
            !          other than simple LES comparison
          enddo
        enddo
      enddo
      dhh85_is_set = .true.
    endif
  else! Keep this else, fallback to 0 Stokes drift
    do kk= 1,g%ke
      do ii = g%isdB,g%iedB
        do jj = g%jsd,g%jed
          cs%Us_x(ii,jj,kk) = 0.
        enddo
      enddo
      do ii = g%isd,g%ied
        do jj = g%jsdB,g%jedB
          cs%Us_y(ii,jj,kk) = 0.
        enddo
      enddo
    enddo
  endif

  ! Turbulent Langmuir number is computed here and available to use anywhere.
  ! SL Langmuir number requires mixing layer depth, and therefore is computed
  ! in the routine it is needed by (e.g. KPP or ePBL).
  do ii = g%isc,g%iec
    do jj = g%jsc, g%jec
      top = h(ii,jj,1)*gv%H_to_Z
      call get_langmuir_number( la, g, gv, us, top, ustar(ii,jj), ii, jj, &
             h(ii,jj,:),override_ma=.false.,waves=cs)
      cs%La_turb(ii,jj) = la
    enddo
  enddo

  ! Output any desired quantities
  if (cs%id_surfacestokes_y>0) &
    call post_data(cs%id_surfacestokes_y, cs%us0_y, cs%diag)
  if (cs%id_surfacestokes_x>0) &
    call post_data(cs%id_surfacestokes_x, cs%us0_x, cs%diag)
  if (cs%id_3dstokes_y>0) &
    call post_data(cs%id_3dstokes_y, cs%us_y, cs%diag)
  if (cs%id_3dstokes_x>0) &
    call post_data(cs%id_3dstokes_x, cs%us_x, cs%diag)
  if (cs%id_La_turb>0) &
    call post_data(cs%id_La_turb, cs%La_turb, cs%diag)

end subroutine update_stokes_drift

!> A subroutine to fill the Stokes drift from a NetCDF file
!! using the data_override procedures.
subroutine surface_bands_by_data_override(day_center, G, GV, US, CS)
  use netcdf
  type(time_type),          intent(in) :: day_center !< Center of timestep
  type(wave_parameters_cs), pointer    :: CS         !< Wave structure
  type(ocean_grid_type), intent(inout) :: G          !< Grid structure
  type(verticalgrid_type),  intent(in) :: GV         !< Vertical grid structure
  type(unit_scale_type),    intent(in) :: US         !< A dimensional unit scaling type
  ! Local variables
  real    :: temp_x(szi_(g),szj_(g)) ! Pseudo-zonal Stokes drift of band at h-points [m s-1]
  real    :: temp_y(szi_(g),szj_(g)) ! Psuedo-meridional Stokes drift of band at h-points [m s-1]
  real    :: Top, MidPoint
  integer :: b
  integer :: i, j
  integer, dimension(4) :: start, counter, dims, dim_id
  character(len=12)  :: dim_name(4)
  character(20) :: varname, varread1, varread2
  integer :: rcode_fr, rcode_wn, ncid, varid_fr, varid_wn, id, ndims

  if (.not.dataoverrideisinitialized) then
    call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
    dataoverrideisinitialized = .true.

    ! Read in number of wavenumber bands in file to set number to be read in
    ! Hardcoded filename/variables
    varread1 = 'wavenumber' !Old method gives wavenumber
    varread2 = 'frequency'  !New method gives frequency
    rcode_wn = nf90_open(trim(surfbandfilename), nf90_nowrite, ncid)
    if (rcode_wn /= 0) then
      call mom_error(fatal,"error opening file "//trim(surfbandfilename)//&
            " in MOM_wave_interface.")
    endif

    ! Check if rcode_wn or rcode_fr is 0 (checks if input has wavenumber or frequency)
    rcode_wn = nf90_inq_varid(ncid, varread1, varid_wn)
    rcode_fr = nf90_inq_varid(ncid, varread2, varid_fr)

    if (rcode_wn /= 0 .and. rcode_fr /= 0) then
      call mom_error(fatal,"error finding variable "//trim(varread1)//&
         " or "//trim(varread2)//" in file "//trim(surfbandfilename)//" in MOM_wave_interface.")

    elseif (rcode_wn == 0) then
      ! wavenumbers found:
      partitionmode = 0
      rcode_wn = nf90_inquire_variable(ncid, varid_wn, ndims=ndims, &
           dimids=dims)
      if (rcode_wn /= 0) then
        call mom_error(fatal, &
             'error inquiring dimensions MOM_wave_interface.')
      endif
      rcode_wn = nf90_inquire_dimension(ncid, dims(1), dim_name(1), len=id)
      if (rcode_wn /= 0) then
        call mom_error(fatal,"error reading dimension 1 data for "// &
             trim(varread1)//" in file "// trim(surfbandfilename)//          &
             " in MOM_wave_interface.")
      endif
      rcode_wn = nf90_inq_varid(ncid, dim_name(1), dim_id(1))
      if (rcode_wn /= 0) then
        call mom_error(fatal,"error finding variable "//trim(dim_name(1))//&
          " in file "//trim(surfbandfilename)//" in MOM_wave_interace.")
      endif
      ! Allocating size of wavenumber bins
      allocate( cs%WaveNum_Cen(1:id) )
      cs%WaveNum_Cen(:) = 0.0
    elseif (rcode_fr == 0) then
      ! frequencies found:
      partitionmode = 1
      rcode_fr = nf90_inquire_variable(ncid, varid_fr, ndims=ndims, &
           dimids=dims)
      if (rcode_fr /= 0) then
        call mom_error(fatal,&
             'error inquiring dimensions MOM_wave_interface.')
      endif
      rcode_fr = nf90_inquire_dimension(ncid, dims(1), dim_name(1), len=id)
      if (rcode_fr /= 0) then
        call mom_error(fatal,"error reading dimension 1 data for "// &
             trim(varread2)//" in file "// trim(surfbandfilename)// &
             " in MOM_wave_interface.")
      endif
      rcode_fr = nf90_inq_varid(ncid, dim_name(1), dim_id(1))
      if (rcode_fr /= 0) then
        call mom_error(fatal,"error finding variable "//trim(dim_name(1))//&
             " in file "//trim(surfbandfilename)//" in MOM_wave_interace.")
      endif
      ! Allocating size of frequency bins
      allocate( cs%Freq_Cen(1:id) )
      cs%Freq_Cen(:) = 0.0
      ! Allocating size of wavenumber bins
      allocate( cs%WaveNum_Cen(1:id) )
      cs%WaveNum_Cen(:) = 0.0
      allocate( cs%STKx0(g%isdB:g%iedB,g%jsd:g%jed,1:id))
      cs%STKx0(:,:,:) = 0.0
      allocate( cs%STKy0(g%isd:g%ied,g%jsdB:g%jedB,1:id))
      cs%STKy0(:,:,:) = 0.0
    endif

    ! Reading wavenumber bins/Frequencies
    start(:) = 1 ! Set all start to 1
    counter(:) = 1 ! Set all counter to 1
    counter(1) = id ! Set counter(1) to id (number of frequency bins)
    if (partitionmode==0) then
      rcode_wn = nf90_get_var(ncid, dim_id(1), cs%WaveNum_Cen, start, counter)
      if (rcode_wn /= 0) then
        call mom_error(fatal,&
             "error reading dimension 1 values for var_name "// &
             trim(varread1)//",dim_name "//trim(dim_name(1))//  &
             " in file "// trim(surfbandfilename)//" in MOM_wave_interface")
      endif
      numbands = id
      do b = 1,numbands ; cs%WaveNum_Cen(b) = us%Z_to_m*cs%WaveNum_Cen(b) ; enddo
    elseif (partitionmode==1) then
      rcode_fr = nf90_get_var(ncid, dim_id(1), cs%Freq_Cen, start, counter)
      if (rcode_fr /= 0) then
        call mom_error(fatal,&
             "error reading dimension 1 values for var_name "// &
             trim(varread2)//",dim_name "//trim(dim_name(1))//  &
             " in file "// trim(surfbandfilename)//" in MOM_wave_interface")
      endif
      numbands = id
      do b = 1,numbands
        cs%WaveNum_Cen(b) = (2.*pi*cs%Freq_Cen(b)*us%T_to_s)**2 / (us%L_to_Z**2*gv%g_Earth)
      enddo
    endif

  endif

  do b = 1,numbands
    temp_x(:,:) = 0.0
    temp_y(:,:) = 0.0
    varname = '                    '
    write(varname,"(A3,I0)")'Usx',b
    call data_override('OCN',trim(varname), temp_x, day_center)
    varname = '                    '
    write(varname,'(A3,I0)')'Usy',b
    call data_override('OCN',trim(varname), temp_y, day_center)
    ! Disperse into halo on h-grid
    call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
    !Filter land values
    do j = g%jsd,g%jed
      do i = g%Isd,g%Ied
        if (abs(temp_x(i,j)) > 10. .or. abs(temp_y(i,j)) > 10.) then
          ! Assume land-mask and zero out
          temp_x(i,j) = 0.0
          temp_y(i,j) = 0.0
        endif
      enddo
    enddo

    ! Interpolate to u/v grids
    do j = g%jsc,g%jec
      do i = g%IscB,g%IecB
        cs%STKx0(i,j,b) = 0.5 * (temp_x(i,j) + temp_x(i+1,j))
      enddo
    enddo
    do j = g%JscB,g%JecB
      do i = g%isc,g%iec
        cs%STKy0(i,j,b) = 0.5 * (temp_y(i,j) + temp_y(i,j+1))
      enddo
    enddo
    ! Disperse into halo on u/v grids
    call pass_vector(cs%STKx0(:,:,b),cs%STKy0(:,:,b), g%Domain, to_all)
  enddo !Closes b-loop

end subroutine surface_bands_by_data_override

!> Interface to get Langmuir number based on options stored in wave structure
!!
!! Note this can be called with an unallocated Waves pointer, which is okay if we
!!  want the wind-speed only dependent Langmuir number.  Therefore, we need to be
!!  careful about what we try to access here.
subroutine get_langmuir_number( LA, G, GV, US, HBL, ustar, i, j, &
                                H, U_H, V_H, Override_MA, Waves )
  type(ocean_grid_type),   intent(in) :: G  !< Ocean grid structure
  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in) :: US !< A dimensional unit scaling type
  integer, intent(in) :: i      !< Meridional index of h-point
  integer, intent(in) :: j      !< Zonal index of h-point
  real, intent(in)    :: ustar  !< Friction velocity [Z T-1 ~> m s-1].
  real, intent(in)    :: HBL    !< (Positive) thickness of boundary layer [Z ~> m].
  logical, optional,       intent(in) :: Override_MA !< Override to use misalignment in LA
                                !! calculation. This can be used if diagnostic
                                !! LA outputs are desired that are different than
                                !! those used by the dynamical model.
  real, dimension(SZK_(GV)), optional, &
       intent(in)      :: H     !< Grid layer thickness [H ~> m or kg m-2]
  real, dimension(SZK_(GV)), optional, &
       intent(in)      :: U_H   !< Zonal velocity at H point [m s-1]
  real, dimension(SZK_(GV)), optional, &
       intent(in)      :: V_H   !< Meridional velocity at H point [m s-1]
  type(wave_parameters_cs), &
       pointer         :: Waves !< Surface wave control structure.

  real, intent(out)    :: LA    !< Langmuir number

!Local Variables
  real :: Top, bottom, midpoint
  real :: Dpt_LASL, ShearDirection, WaveDirection
  real :: LA_STKx, LA_STKy, LA_STK ! Stokes velocities in [m s-1]
  logical :: ContinueLoop, USE_MA
  real, dimension(SZK_(G)) :: US_H, VS_H
  real, dimension(NumBands) :: StkBand_X, StkBand_Y
  integer :: KK, BB

 ! Compute averaging depth for Stokes drift (negative)
  dpt_lasl = min(-0.1*us%m_to_Z, -la_frachbl*hbl)

  use_ma = la_misalignment
  if (present(override_ma)) use_ma = override_ma

  ! If requesting to use misalignment in the Langmuir number compute the Shear Direction
  if (use_ma) then
    if (.not.(present(h).and.present(u_h).and.present(v_h))) then
      call mom_error(fatal,'Get_LA_waves requested to consider misalignment.')
    endif
    continueloop = .true.
    bottom = 0.0
    do kk = 1,g%ke
      top = bottom
      midpoint = bottom + gv%H_to_Z*0.5*h(kk)
      bottom = bottom + gv%H_to_Z*h(kk)
      if (midpoint > dpt_lasl .and. kk > 1 .and. continueloop) then
        sheardirection = atan2(v_h(1)-v_h(kk),u_h(1)-u_h(kk))
        continueloop = .false.
      endif
    enddo
  endif

  if (wavemethod==testprof) then
    do kk = 1,g%ke
      us_h(kk) = 0.5*(waves%US_X(i,j,kk)+waves%US_X(i-1,j,kk))
      vs_h(kk) = 0.5*(waves%US_Y(i,j,kk)+waves%US_Y(i,j-1,kk))
    enddo
    call get_sl_average_prof( gv, dpt_lasl, h, us_h, la_stkx)
    call get_sl_average_prof( gv, dpt_lasl, h, vs_h, la_stky)
    la_stk = sqrt(la_stkx*la_stkx+la_stky*la_stky)
  elseif (wavemethod==surfbands) then
    do bb = 1,numbands
      stkband_x(bb) = 0.5*(waves%STKx0(i,j,bb)+waves%STKx0(i-1,j,bb))
      stkband_y(bb) = 0.5*(waves%STKy0(i,j,bb)+waves%STKy0(i,j-1,bb))
    enddo
    call get_sl_average_band(gv, dpt_lasl, numbands, waves%WaveNum_Cen, stkband_x, la_stkx )
    call get_sl_average_band(gv, dpt_lasl, numbands, waves%WaveNum_Cen, stkband_y, la_stky )
    la_stk = sqrt(la_stkx**2 + la_stky**2)
  elseif (wavemethod==dhh85) then
    ! Temporarily integrating profile rather than spectrum for simplicity
    do kk = 1,gv%ke
      us_h(kk) = 0.5*(waves%US_X(i,j,kk)+waves%US_X(i-1,j,kk))
      vs_h(kk) = 0.5*(waves%US_Y(i,j,kk)+waves%US_Y(i,j-1,kk))
    enddo
    call get_sl_average_prof( gv, dpt_lasl, h, us_h, la_stkx)
    call get_sl_average_prof( gv, dpt_lasl, h, vs_h, la_stky)
    la_stk = sqrt(la_stkx**2 + la_stky**2)
  elseif (wavemethod==lf17) then
    call get_stokessl_lifoxkemper(ustar, hbl*la_frachbl, gv, us, la_stk, la)
  elseif (wavemethod==null_wavemethod) then
    call mom_error(fatal, "Get_Langmuir_number called without defining a WaveMethod. "//&
                          "Suggest to make sure USE_LT is set/overridden to False or "//&
                          "choose a wave method (or set USE_LA_LI2016 to use statistical "//&
                          "waves.")
  endif

  if (.not.(wavemethod==lf17)) then
    ! This is an arbitrary lower bound on Langmuir number.
    ! We shouldn't expect values lower than this, but
    ! there is also no good reason to cap it here other then
    ! to prevent large enhancements in unconstrained parts of
    ! the curve fit parameterizations.
    ! Note the dimensional constant background Stokes velocity of 10^-10 m s-1.
    la = max(waves%La_min, sqrt(us%Z_to_m*us%s_to_T*ustar / (la_stk+1.e-10)))
  endif

  if (use_ma) then
    wavedirection = atan2(la_stky, la_stkx)
    la = la / sqrt(max(1.e-8, cos( wavedirection - sheardirection)))
  endif

  return
end subroutine get_langmuir_number

!> Get SL averaged Stokes drift from Li/FK 17 method
!!
!! Original description:
!! - This function returns the enhancement factor, given the 10-meter
!!   wind [m s-1], friction velocity [m s-1] and the boundary layer depth [m].
!!
!! Update (Jan/25):
!! - Converted from function to subroutine, now returns Langmuir number.
!! - Computs 10m wind internally, so only ustar and hbl need passed to
!!   subroutine.
!!
!! Qing Li, 160606
!! - BGR port from CVMix to MOM6 Jan/25/2017
!! - BGR change output to LA from Efactor
!! - BGR remove u10 input
!! - BGR note: fixed parameter values should be changed to "get_params"
subroutine get_stokessl_lifoxkemper(ustar, hbl, GV, US, UStokes_SL, LA)
  real, intent(in)  :: ustar !< water-side surface friction velocity [Z T-1 ~> m s-1].
  real, intent(in)  :: hbl   !< boundary layer depth [Z ~> m].
  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure
  type(unit_scale_type),   intent(in) :: US !< A dimensional unit scaling type
  real, intent(out) :: UStokes_SL !< Surface layer averaged Stokes drift [m s-1]
  real, intent(out) :: LA    !< Langmuir number
  ! Local variables
  ! parameters
  real, parameter :: &
       ! ratio of U19.5 to U10 (Holthuijsen, 2007)
       u19p5_to_u10 = 1.075, &
       ! ratio of mean frequency to peak frequency for
       ! Pierson-Moskowitz spectrum (Webb, 2011)
       fm_into_fp = 1.296, &
       ! ratio of surface Stokes drift to U10
       us_to_u10 = 0.0162, &
       ! loss ratio of Stokes transport
       r_loss = 0.667
  real :: UStokes, hm0, fm, fp, vstokes, kphil, kstar
  real :: z0, z0i, r1, r2, r3, r4, tmp, lasl_sqr_i
  real :: u10

  if (ustar > 0.0) then
    ! Computing u10 based on u_star and COARE 3.5 relationships
    call ust_2_u10_coare3p5(us%Z_to_m*us%s_to_T*ustar*sqrt(us%R_to_kg_m3*gv%Rho0/1.225), u10, gv, us)
    ! surface Stokes drift
    ustokes = us_to_u10*u10
    !
    ! significant wave height from Pierson-Moskowitz
    ! spectrum (Bouws, 1998)
    hm0 = 0.0246 *u10**2
    !
    ! peak frequency (PM, Bouws, 1998)
    tmp = 2.0 * pi * u19p5_to_u10 * u10
    fp = 0.877 * us%L_T_to_m_s**2*us%m_to_Z * gv%g_Earth / tmp
    !
    ! mean frequency
    fm = fm_into_fp * fp
    !
    ! total Stokes transport (a factor r_loss is applied to account
    !  for the effect of directional spreading, multidirectional waves
    !  and the use of PM peak frequency and PM significant wave height
    !  on estimating the Stokes transport)
    vstokes = 0.125 * pi * r_loss * fm * hm0**2
    !
    ! the general peak wavenumber for Phillips' spectrum
    ! (Breivik et al., 2016) with correction of directional spreading
    kphil = 0.176 * ustokes / vstokes
    !
    ! surface layer averaged Stokes dirft with Stokes drift profile
    ! estimated from Phillips' spectrum (Breivik et al., 2016)
    ! the directional spreading effect from Webb and Fox-Kemper, 2015
    ! is also included
    kstar = kphil * 2.56
    ! surface layer
    z0 = abs(us%Z_to_m*hbl)
    z0i = 1.0 / z0
    ! term 1 to 4
    r1 = ( 0.151 / kphil * z0i -0.84 ) * &
         ( 1.0 - exp(-2.0 * kphil * z0) )
    r2 = -( 0.84 + 0.0591 / kphil * z0i ) * &
         sqrt( 2.0 * pi * kphil * z0 ) * &
         erfc( sqrt( 2.0 * kphil * z0 ) )
    r3 = ( 0.0632 / kstar * z0i + 0.125 ) * &
         (1.0 - exp(-2.0 * kstar * z0) )
    r4 = ( 0.125 + 0.0946 / kstar * z0i ) * &
         sqrt( 2.0 * pi *kstar * z0) * &
         erfc( sqrt( 2.0 * kstar * z0 ) )
    ustokes_sl = ustokes * (0.715 + r1 + r2 + r3 + r4)
    la = sqrt(us%Z_to_m*us%s_to_T*ustar / ustokes_sl)
  else
    ustokes_sl = 0.0
    la=1.e8
  endif

end subroutine get_stokessl_lifoxkemper

!> Get SL Averaged Stokes drift from a Stokes drift Profile
subroutine get_sl_average_prof( GV, AvgDepth, H, Profile, Average )
  type(verticalgrid_type),  &
       intent(in)   :: GV       !< Ocean vertical grid structure
  real, intent(in)  :: AvgDepth !< Depth to average over (negative) [Z ~> m].
  real, dimension(SZK_(GV)), &
       intent(in)   :: H        !< Grid thickness [H ~> m or kg m-2]
  real, dimension(SZK_(GV)), &
       intent(in)   :: Profile  !< Profile of quantity to be averaged [arbitrary]
                                !! (used here for Stokes drift)
  real, intent(out) :: Average  !< Output quantity averaged over depth AvgDepth [arbitrary]
                                !! (used here for Stokes drift)
  !Local variables
  real :: top, midpoint, bottom ! Depths, negative downward [Z ~> m].
  real :: Sum
  integer :: kk

  ! Initializing sum
  sum = 0.0

  ! Integrate
  bottom = 0.0
  do kk = 1, gv%ke
    top = bottom
    midpoint = bottom - gv%H_to_Z * 0.5*h(kk)
    bottom = bottom - gv%H_to_Z * h(kk)
    if (avgdepth < bottom) then ! The whole cell is within H_LA
      sum = sum + profile(kk) * (gv%H_to_Z * h(kk))
    elseif (avgdepth < top) then ! A partial cell is within H_LA
      sum = sum + profile(kk) * (top-avgdepth)
      exit
    else
      exit
    endif
  enddo

  ! Divide by AvgDepth or the depth in the column, whichever is smaller.
  if (abs(avgdepth) <= abs(bottom)) then
    average = sum / abs(avgdepth)
  elseif (abs(bottom) > 0.0) then
    average = sum / abs(bottom)
  else
    average = 0.0
  endif

end subroutine get_sl_average_prof

!> Get SL averaged Stokes drift from the banded Spectrum method
subroutine get_sl_average_band( GV, AvgDepth, NB, WaveNumbers, SurfStokes, Average )
  type(verticalgrid_type),  &
       intent(in)     :: GV          !< Ocean vertical grid
  real, intent(in)    :: AvgDepth    !< Depth to average over [Z ~> m].
  integer, intent(in) :: NB          !< Number of bands used
  real, dimension(NB), &
       intent(in)     :: WaveNumbers !< Wavenumber corresponding to each band [Z-1 ~> m-1]
  real, dimension(NB), &
       intent(in)     :: SurfStokes  !< Surface Stokes drift for each band [m s-1]
  real, intent(out)   :: Average     !< Output average Stokes drift over depth AvgDepth [m s-1]

  ! Local variables
  integer :: bb

  ! Loop over bands
  average = 0.0
  do bb = 1, nb
    ! Factor includes analytical integration of e(2kz)
    !  - divided by (-H_LA) to get average from integral.
    average = average + surfstokes(bb) * &
              (1.-exp(-abs(avgdepth * 2.0 * wavenumbers(bb)))) / &
              abs(avgdepth * 2.0 * wavenumbers(bb))
  enddo

  return
end subroutine get_sl_average_band

!> Compute the Stokes drift at a given depth
!!
!! Taken from Qing Li (Brown)
!! use for comparing MOM6 simulation to his LES
!! computed at z mid point (I think) and not depth averaged.
!! Should be fine to integrate in frequency from 0.1 to sqrt(-0.2*grav*2pi/dz
subroutine dhh85_mid(GV, US, zpt, UStokes)
  type(verticalgrid_type), intent(in)  :: GV  !< Ocean vertical grid
  type(unit_scale_type),   intent(in)  :: US  !< A dimensional unit scaling type
  real, intent(in)  :: zpt   !< Depth to get Stokes drift [Z ~> m].
  real, intent(out) :: UStokes !< Stokes drift [m s-1]
  !
  real :: ann, Bnn, Snn, Cnn, Dnn
  real :: omega_peak, omega, u10, WA, domega
  real :: omega_min, omega_max, wavespec, Stokes
  real :: g_Earth ! Gravitational acceleration [m s-2]
  integer :: Nomega, OI

  wa = waveage
  u10 = wavewind
  g_earth = us%L_T_to_m_s**2*us%m_to_Z * gv%g_Earth

  !/
  omega_min = 0.1 ! Hz
  ! Cut off at 30cm for now...
  omega_max = 10. ! ~sqrt(0.2*g_Earth*2*pi/0.3)
  nomega = 1000
  domega = (omega_max-omega_min)/real(nomega)

  !
  if (waveagepeakfreq) then
    omega_peak = g_earth / (wa * u10)
  else
    omega_peak = 2. * pi * 0.13 * g_earth / u10
  endif
  !/
  ann = 0.006 * waveage**(-0.55)
  bnn = 1.0
  snn = 0.08 * (1.0 + 4.0 * waveage**3)
  cnn = 1.7
  if (wa < 1.) then
    cnn = cnn - 6.0*log10(wa)
  endif
  !/
  ustokes = 0.0
  omega = omega_min + 0.5*domega
  do oi = 1,nomega-1
    dnn = exp( -0.5 * (omega-omega_peak)**2 / (snn**2 * omega_peak**2) )
    ! wavespec units = m2s
    wavespec = (ann * g_earth**2 / (omega_peak*omega**4 ) ) * &
               exp(-bnn*(omega_peak/omega)**4)*cnn**dnn
    ! Stokes units m  (multiply by frequency range for units of m/s)
    stokes = 2.0 * wavespec * omega**3 * &
         exp( 2.0 * omega**2 * us%Z_to_m*zpt / g_earth) / g_earth
    ustokes = ustokes + stokes*domega
    omega = omega + domega
  enddo

  return
end subroutine dhh85_mid

!> Explicit solver for Stokes mixing.
!! Still in development do not use.
subroutine stokesmixing(G, GV, dt, h, u, v, Waves )
  type(ocean_grid_type), &
       intent(in)    :: G     !< Ocean grid
  type(verticalgrid_type), &
       intent(in)    :: GV    !< Ocean vertical grid
  real, intent(in)   :: dt    !< Time step of MOM6 [T ~> s] for explicit solver
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),&
       intent(in)    :: h     !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
       intent(inout) :: u     !< Velocity i-component [m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
       intent(inout) :: v     !< Velocity j-component [m s-1]
  type(wave_parameters_cs), &
       pointer       :: Waves !< Surface wave related control structure.
  ! Local variables
  real :: dTauUp, dTauDn ! Vertical momentum fluxes [Z T-1 m s-1]
  real :: h_Lay  ! The layer thickness at a velocity point [Z ~> m].
  integer :: i,j,k

! This is a template to think about down-Stokes mixing.
! This is not ready for use...

  do k = 1, g%ke
    do j = g%jsc, g%jec
      do i = g%iscB, g%iecB
        h_lay = gv%H_to_Z*0.5*(h(i,j,k)+h(i+1,j,k))
        dtauup = 0.0
        if (k > 1) &
          dtauup = 0.5*(waves%Kvs(i,j,k)+waves%Kvs(i+1,j,k)) * &
               (waves%us_x(i,j,k-1)-waves%us_x(i,j,k)) / &
               (0.5*(h_lay + gv%H_to_Z*0.5*(h(i,j,k-1)+h(i+1,j,k-1)) ))
        dtaudn = 0.0
        if (k < g%ke-1) &
          dtaudn = 0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i+1,j,k+1)) * &
               (waves%us_x(i,j,k)-waves%us_x(i,j,k+1)) / &
               (0.5*(h_lay + gv%H_to_Z*0.5*(h(i,j,k+1)+h(i+1,j,k+1)) ))
        u(i,j,k) = u(i,j,k) + dt * (dtauup-dtaudn) / h_lay
      enddo
    enddo
  enddo

  do k = 1, g%ke
    do j = g%jscB, g%jecB
      do i = g%isc, g%iec
        h_lay = gv%H_to_Z*0.5*(h(i,j,k)+h(i,j+1,k))
        dtauup = 0.
        if (k > 1) &
          dtauup = 0.5*(waves%Kvs(i,j,k)+waves%Kvs(i,j+1,k)) * &
               (waves%us_y(i,j,k-1)-waves%us_y(i,j,k)) / &
               (0.5*(h_lay + gv%H_to_Z*0.5*(h(i,j,k-1)+h(i,j+1,k-1)) ))
        dtaudn = 0.0
        if (k < g%ke-1) &
          dtaudn =0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i,j+1,k+1)) * &
               (waves%us_y(i,j,k)-waves%us_y(i,j,k+1)) / &
               (0.5*(h_lay + gv%H_to_Z*0.5*(h(i,j,k+1)+h(i,j+1,k+1)) ))
        v(i,j,k) = v(i,j,k) + dt * (dtauup-dtaudn) / h_lay
      enddo
    enddo
  enddo

end subroutine stokesmixing

!> Solver to add Coriolis-Stokes to model
!! Still in development and not meant for general use.
!! Can be activated (with code intervention) for LES comparison
!! CHECK THAT RIGHT TIMESTEP IS PASSED IF YOU USE THIS**
!!
!! Not accessed in the standard code.
subroutine coriolisstokes(G, GV, DT, h, u, v, WAVES, US)
  type(ocean_grid_type), &
       intent(in)    :: G     !< Ocean grid
  type(verticalgrid_type), &
       intent(in)   :: GV     !< Ocean vertical grid
  real, intent(in)  :: Dt     !< Time step of MOM6 [s] CHECK IF PASSING RIGHT TIMESTEP
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
       intent(in)    :: h     !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
       intent(inout) :: u     !< Velocity i-component [m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
       intent(inout) :: v     !< Velocity j-component [m s-1]
  type(wave_parameters_cs), &
       pointer       :: Waves !< Surface wave related control structure.
  type(unit_scale_type),   intent(in) :: US     !< A dimensional unit scaling type
  ! Local variables
  real :: DVel ! A rescaled velocity change [m s-1 T-1 ~> m s-2]
  integer :: i,j,k

  do k = 1, g%ke
    do j = g%jsc, g%jec
      do i = g%iscB, g%iecB
        dvel = 0.25*(waves%us_y(i,j+1,k)+waves%us_y(i-1,j+1,k))*g%CoriolisBu(i,j+1) + &
               0.25*(waves%us_y(i,j,k)+waves%us_y(i-1,j,k))*g%CoriolisBu(i,j)
        u(i,j,k) = u(i,j,k) + dvel*us%s_to_T*dt
      enddo
    enddo
  enddo

  do k = 1, g%ke
    do j = g%jscB, g%jecB
      do i = g%isc, g%iec
        dvel = 0.25*(waves%us_x(i+1,j,k)+waves%us_x(i+1,j-1,k))*g%CoriolisBu(i+1,j) + &
               0.25*(waves%us_x(i,j,k)+waves%us_x(i,j-1,k))*g%CoriolisBu(i,j)
        v(i,j,k) = v(i,j,k) - dvel*us%s_to_T*dt
      enddo
    enddo
  enddo
end subroutine coriolisstokes

!> Computes wind speed from ustar_air based on COARE 3.5 Cd relationship
!! Probably doesn't belong in this module, but it is used here to estimate
!! wind speed for wind-wave relationships.  Should be a fine way to estimate
!! the neutral wind-speed as written here.
subroutine ust_2_u10_coare3p5(USTair, U10, GV, US)
  real, intent(in)                    :: USTair !< Wind friction velocity [m s-1]
  real, intent(out)                   :: U10    !< 10-m neutral wind speed [m s-1]
  type(verticalgrid_type), intent(in) :: GV     !< vertical grid type
  type(unit_scale_type),   intent(in) :: US     !< A dimensional unit scaling type

  ! Local variables
  real, parameter :: vonkar = 0.4 ! Should access a get_param von karman
  real, parameter :: nu=1e-6 ! Should access a get_param air-viscosity
  real :: z0sm, z0, z0rough, u10a, alpha, CD
  integer :: CT

  ! Uses empirical formula for z0 to convert ustar_air to u10 based on the
  !  COARE 3.5 paper (Edson et al., 2013)
  ! alpha=m*U10+b
  ! Note in Edson et al. 2013, eq. 13 m is given as 0.017.  However,
  ! m=0.0017 reproduces the curve in their figure 6.

  z0sm = 0.11 * nu * us%m_to_Z / ustair !Compute z0smooth from ustar guess
  u10 = ustair/sqrt(0.001)  !Guess for u10
  u10a = 1000

  ct=0
  do while (abs(u10a/u10-1.) > 0.001)
    ct=ct+1
    u10a = u10
    alpha = min(0.028, 0.0017 * u10 - 0.005)
    z0rough = alpha * (us%m_s_to_L_T*ustair)**2 / gv%g_Earth ! Compute z0rough from ustar guess
    z0 = z0sm + z0rough
    cd = ( vonkar / log(10.*us%m_to_Z / z0) )**2 ! Compute CD from derived roughness
    u10 = ustair/sqrt(cd)  ! Compute new u10 from derived CD, while loop
                           ! ends and checks for convergence...CT counter
                           ! makes sure loop doesn't run away if function
                           ! doesn't converge.  This code was produced offline
                           ! and converged rapidly (e.g. 2 cycles)
                           ! for ustar=0.0001:0.0001:10.
    if (ct>20) then
      u10 = ustair/sqrt(0.0015) ! I don't expect to get here, but just
                              !  in case it will output a reasonable value.
      exit
    endif
  enddo
  return
end subroutine ust_2_u10_coare3p5

!> Clear pointers, deallocate memory
subroutine waves_end(CS)
  type(wave_parameters_cs), pointer :: CS !< Control structure

  if (allocated(cs%WaveNum_Cen)) deallocate( cs%WaveNum_Cen )
  if (allocated(cs%Freq_Cen))    deallocate( cs%Freq_Cen )
  if (allocated(cs%Us_x))        deallocate( cs%Us_x )
  if (allocated(cs%Us_y))        deallocate( cs%Us_y )
  if (allocated(cs%La_SL))       deallocate( cs%La_SL )
  if (allocated(cs%La_turb))     deallocate( cs%La_turb )
  if (allocated(cs%STKx0))       deallocate( cs%STKx0 )
  if (allocated(cs%STKy0))       deallocate( cs%STKy0 )
  if (allocated(cs%KvS))         deallocate( cs%KvS )
  if (allocated(cs%Us0_y))       deallocate( cs%Us0_y )
  if (allocated(cs%Us0_x))       deallocate( cs%Us0_x )

  deallocate( cs )

  return
end subroutine waves_end

!> \namespace  mom_wave_interface
!!
!! \author Brandon Reichl, 2018.
!!
!! This module should be moved as wave coupling progresses and
!! likely will should mirror the iceberg or sea-ice model set-up.
!!
!! This module is meant to contain the routines to read in and
!! interpret surface wave data for MOM6. In its original form, the
!! capabilities include setting the Stokes drift in the model (from a
!! variety of sources including prescribed, empirical, and input
!! files).  In short order, the plan is to also ammend the subroutine
!! to accept Stokes drift information from an external coupler.
!! Eventually, it will be necessary to break this file apart so that
!! general wave information may be stored in the control structure
!! and the Stokes drift effect can be isolated from processes such as
!! sea-state dependent momentum fluxes, gas fluxes, and other wave
!! related air-sea interaction and boundary layer phenomenon.
!!
!! The Stokes drift are stored on the C-grid with the conventional
!! protocol to interpolate to the h-grid to compute Langmuir number,
!! the primary quantity needed for Langmuir turbulence
!! parameterizations in both the ePBL and KPP approach.  This module
!! also computes full 3d Stokes drift profiles, which will be useful
!! if second-order type boundary layer parameterizations are
!! implemented (perhaps via GOTM, work in progress).

end module mom_wave_interface