Turbulent momentum flux down Stokes gradient

[harcourt]

I've been working on putting in a MY2.5-variant parameterization of Langmuir turbulence (Harcourt,2013; 2015; + newer mods) that entails a component of the vertical momentum flux directed down the gradient of the surface waves' Stokes drift, i.e.:

<u'w'> = - K_M du/dz - K_W du^s/dz

where u is horizontal velocity, K_M is eddy viscosity, u^s is the Stokes drift and K_W is an eddy coefficient for the new component of momentum flux (K_W is a.k.a. K_M^S in pubs; '_W' for Wave seems better now). Several other proposals have been made over the past decade+ to add such a term for KPP; most other proposals entail K_M=K_M^S, so that the form of the momentum flux closure is just modified by replacing Eulerian velocity u with Lagrangian u^L=u+u^s. However, in both cases it is necessary to add additional forcing to the discretized equations for predicting u & v. So his topic might be relevant to other parallel efforts out there (an to other posts here as old as 2004).

I've chosen to implement the additional force explicitly, following the stencil of the (effectively similar) BODYFORCE option (for distributing surface stress over a near-surface depth interval), and gleaning what I can from the ROMS documentation and other subroutines. However, I'm pretty new to reading the ROMS code and would like to solicit comment on the changes I've made to Nonlinear/rhs3d.F in J.C. Warner's recent COAWST version, attached here & with a diff. I'm using the COAWST version because it's the only one coupled to a wave model that has the Vortex Force Wave interaction of Uchiyama et al (2010) coded in, and that I that I could seem to get my hands on. There are other differences with the current Rutgers ROMS version of rhs3d.F, so the diff file may help.

The changes are active under a new option CLVF_CLOSURE for 'Craik-Leibovich Vortex Force Closure', and either HARCOURT_2013 or HARCOURT_2015 to get you a down-Stokes gradient closure term following (the admittedly moving target of) different formulations. I'm then adding the forcing from -K_w du^s/dz to the rhs expressions for du/dt & dv/dt , after computing K_W elsewhere, in gls_corstep.F . Things run and produce reasonable-looking (to me) results, but I have a couple of specific questions I'd like to check, and a couple more vague ones to float:

Specific Questions:

1) Do I have the correct factors of (inverse) pn's and pm's in these expressions coded into rhs3d.F? These things are Greek to me.

2) Have I correctly used z_r and Hz? I feel like I'm missing a 1/Hz, but I think that's because we're computing the rhs of d(Hz u)/dt, yes?

3) This seems to work okay & runs stably (except n.b. spline-option issues below), but does anyone see a numerically superior approach?

Comments & Vague Questions:

4) I seem to have numerical problems in gls_corstep.F and elsewhere (not attached) when I apply the SPLINE option to computing the Stokes shear as well as the Eulerian shear when it is on. For now, defining SPLINE instead entails combining polynomial splines applied to just du/dz, with ordinary centered finite differences for du^s/dz. Distasteful and probably oh-so-wrong, but it works and preserves the utility of the spline interpolation for the Eulerian shear at the mixed layer bottom. Any thoughts on why applying polynomial spline interpolation to an exponential Stokes drift profile might cause numerical problems?

5) Stokes drift from SWAN in both Rutgers and COAWST coupled models uses a monochromatic approximation for either WEC_MELLOR option in either, or for WEC_VF in COAWST, placing all the energy at the mean wavelength, as suggested in Uchiyama et al. (2010), Eq. 2. This might work for swell, or for shoaling coastal waves for all I know, but for Langmuir turbulence below deep water equilibrium wind waves, this approximation crucially underestimates the surface Stokes drift by a factor of 4-5 and overestimates the Stokes decay depth by a similar ratio. My quick fix is to rescale u^s to r u^s|_{z=0} * (u^s/u^s|_{z=0})**r while preserving the sign of u^s, where roughly 4 < r < 5 . This preserves Stokes transport for deep water waves, but it is not at all a good approximation for a more general case. In the long run I'm looking to pass into ROMS the spectrum of the Stokes drift (i.e. one scalar spectrum for each Stokes drift component), or some other equivalent distillate of the full 2D directional spectrum in order to compute u^s,v^s more accurately on ROMS' vertical grid. Looking over what's currently done to pass wave parameters in 'block' form in the current couplings, getting such spectral information (2 vectors of length the number of SWAN frequencies, at each ROMS (x,y) coordinate) to pass from SWAN to ROMS looks like a daunting task. Anybody done this already? Thoughts on how to do it? Dire warnings?

-Ramsey Harcourt, UW/APL

Refs:
Harcourt, R.R., 2013, "A second-moment closuremodel of Langmuir turbulence", J. Phys. Oceanogr., 43, 673–697.
Harcourt, R.R., 2015, "An improved second-moment closure model of Langmuir turbulence", J. Phys. Oceanogr., 45, 84-103.
Uchiyama, Y., McWlliams, J. C., Shchepetkin, A. F., 2010, "Wave-current interaction in an oceanic circulation model with a vortex-force formalism: application to the surf zone", Ocean Modeling, 34, 16-35.