src/user/MOM_wave_interface.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		!> Interface for surface waves
6		module MOM_wave_interface
7
8		use MOM_data_override, only : data_override_init, data_override
9		use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc
10		use MOM_diag_mediator, only : diag_ctrl
11		use MOM_domains, only : pass_var, pass_vector, AGRID
12		use MOM_domains, only : To_South, To_West, To_All
13		use MOM_error_handler, only : MOM_error, FATAL, WARNING
14		use MOM_file_parser, only : get_param, log_version, param_file_type
15		use MOM_forcing_type, only : mech_forcing
16		use MOM_grid, only : ocean_grid_type
17		use MOM_hor_index, only : hor_index_type
18		use MOM_io, only : file_exists, get_var_sizes, read_variable
19		use MOM_io, only : vardesc, var_desc
20		use MOM_safe_alloc, only : safe_alloc_ptr
21		use MOM_spatial_means, only : global_area_mean
22		use MOM_time_manager, only : time_type, operator(+), operator(/)
23		use MOM_unit_scaling, only : unit_scale_type
24		use MOM_variables, only : thermo_var_ptrs, surface
25		use MOM_verticalgrid, only : verticalGrid_type
26		use MOM_restart, only : register_restart_pair, MOM_restart_CS
27
28		implicit none ; private
29
30		#include <MOM_memory.h>
31
32		public MOM_wave_interface_init ! Public interface to fully initialize the wave routines.
33		public query_wave_properties ! Public interface to obtain information from the waves control structure.
34		public Update_Surface_Waves ! Public interface to update wave information at the
35		! coupler/driver level.
36		public Update_Stokes_Drift ! Public interface to update the Stokes drift profiles
37		! called in step_mom.
38		public get_Langmuir_Number ! Public interface to compute Langmuir number called from
39		! ePBL or KPP routines.
40		public Stokes_PGF ! Public interface to compute Stokes-shear induced pressure gradient force anomaly
41		public StokesMixing ! NOT READY - Public interface to add down-Stokes gradient
42		! momentum mixing (e.g. the approach of Harcourt 2013/2015)
43		public CoriolisStokes ! NOT READY - Public interface to add Coriolis-Stokes acceleration
44		! of the mean currents, needed for comparison with LES. It is
45		! presently advised against implementing in non-1d settings without
46		! serious consideration of the full 3d wave-averaged Navier-Stokes
47		! CL2 effects.
48		public Waves_end ! public interface to deallocate and free wave related memory.
49		public get_wave_method ! public interface to obtain the wave method string
50		public waves_register_restarts ! public interface to register wave restart fields
51
52		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
53		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
54		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
55		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
56
57		!> Container for all surface wave related parameters
58		type, public :: wave_parameters_CS ; private
59
60		! Main surface wave options and publicly visible variables
61		logical, public :: UseWaves = .false. !< Flag to enable surface gravity wave feature
62		logical, public :: Stokes_VF = .false. !< True if Stokes vortex force is used
63		logical, public :: Passive_Stokes_VF = .false. !< Computes Stokes VF, but doesn't affect dynamics
64		logical, public :: Stokes_PGF = .false. !< True if Stokes shear pressure Gradient force is used
65		logical, public :: robust_Stokes_PGF = .false. !< If true, use expressions to calculate the
66		!! Stokes-induced pressure gradient anomalies that are
67		!! more accurate in the limit of thin layers.
68		logical, public :: Passive_Stokes_PGF = .false. !< Keeps Stokes_PGF on, but doesn't affect dynamics
69		logical, public :: Stokes_DDT = .false. !< Developmental:
70		!! True if Stokes d/dt is used
71		logical, public :: Passive_Stokes_DDT = .false. !< Keeps Stokes_DDT on, but doesn't affect dynamics
72		logical :: Homogenize_Surfbands !< True to homogenize surface band Stokes drift in the horizontal
73
74		real, allocatable, dimension(:,:,:), public :: &
75		Us_x !< 3d zonal Stokes drift profile [L T-1 ~> m s-1]
76		!! Horizontal -> U points
77		!! Vertical -> Mid-points
78		real, allocatable, dimension(:,:,:), public :: &
79		Us_y !< 3d meridional Stokes drift profile [L T-1 ~> m s-1]
80		!! Horizontal -> V points
81		!! Vertical -> Mid-points
82		real, allocatable, dimension(:,:,:), public :: &
83		ddt_Us_x !< 3d time tendency of zonal Stokes drift profile [L T-2 ~> m s-2]
84		!! Horizontal -> U points
85		!! Vertical -> Mid-points
86		real, allocatable, dimension(:,:,:), public :: &
87		ddt_Us_y !< 3d time tendency of meridional Stokes drift profile [L T-2 ~> m s-2]
88		!! Horizontal -> V points
89		!! Vertical -> Mid-points
90		real, allocatable, dimension(:,:,:), public :: &
91		Us_x_from_ddt !< Check of 3d zonal Stokes drift profile [L T-1 ~> m s-1]
92		!! Horizontal -> U points
93		!! Vertical -> Mid-points
94		real, allocatable, dimension(:,:,:), public :: &
95		Us_y_from_ddt !< Check of 3d meridional Stokes drift profile [L T-1 ~> m s-1]
96		!! Horizontal -> V points
97		!! Vertical -> Mid-points
98		real, allocatable, dimension(:,:,:), public :: &
99		Us_x_prev !< 3d zonal Stokes drift profile, previous dynamics call [L T-1 ~> m s-1]
100		!! Horizontal -> U points
101		!! Vertical -> Mid-points
102		real, allocatable, dimension(:,:,:), public :: &
103		Us_y_prev !< 3d meridional Stokes drift profile, previous dynamics call [L T-1 ~> m s-1]
104		!! Horizontal -> V points
105		!! Vertical -> Mid-points
106		real, allocatable, dimension(:,:,:), public :: &
107		KvS !< Viscosity for Stokes Drift shear [H Z T-1 ~> m2 s-1 or Pa s]
108		real, allocatable, dimension(:), public :: &
109		WaveNum_Cen !< Wavenumber bands for read/coupled [Z-1 ~> m-1]
110		real, allocatable, dimension(:,:,:), public :: &
111		UStk_Hb !< Surface Stokes Drift spectrum (zonal) [L T-1 ~> m s-1]
112		!! Horizontal -> H-points
113		!! 3rd dimension -> Freq/Wavenumber
114		real, allocatable, dimension(:,:,:), public :: &
115		VStk_Hb !< Surface Stokes Drift spectrum (meridional) [L T-1 ~> m s-1]
116		!! Horizontal -> H-points
117		!! 3rd dimension -> Freq/Wavenumber
118		real, allocatable, dimension(:,:), public :: &
119		Omega_w2x !< wind direction ccw from model x- axis [nondim radians]
120		integer, public :: NumBands = 0 !< Number of wavenumber/frequency partitions
121		!! Must match the number of bands provided
122		!! via either coupling or file.
123
124		! The remainder of this control structure is private
125		integer :: WaveMethod = -99 !< Options for including wave information
126		!! Valid (tested) choices are:
127		!! 0 - Test Profile
128		!! 1 - Surface Stokes Drift Bands
129		!! 2 - DHH85
130		!! 3 - LF17
131		!! -99 - No waves computed, but empirical Langmuir number used.
132		logical :: LagrangianMixing !< This feature is in development and not ready
133		!! True if Stokes drift is present and mixing
134		!! should be applied to Lagrangian current
135		!! (mean current + Stokes drift).
136		!! See Reichl et al., 2016 KPP-LT approach
137		logical :: StokesMixing !< This feature is in development and not ready.
138		!! True if vertical mixing of momentum
139		!! should be applied directly to Stokes current
140		!! (with separate mixing parameter for Eulerian
141		!! mixing contribution).
142		!! See Harcourt 2013, 2015 Second-Moment approach
143		logical :: CoriolisStokes !< This feature is in development and not ready.
144		! True if Coriolis-Stokes acceleration should be applied.
145		real :: Stokes_min_thick_avg !< A layer thickness below which the cell-center Stokes drift is
146		!! used instead of the cell average [Z ~> m]. This is only used if
147		!! WAVE_INTERFACE_ANSWER_DATE < 20230101.
148		integer :: answer_date !< The vintage of the order of arithmetic and expressions in the
149		!! surface wave calculations. Values below 20230101 recover the
150		!! answers from the end of 2022, while higher values use updated
151		!! and more robust forms of the same expressions.
152
153		! Options if WaveMethod is Surface Stokes Drift Bands (1)
154		integer :: PartitionMode !< Method for partition mode (meant to check input)
155		!! 0 - wavenumbers
156		!! 1 - frequencies
157		integer :: DataSource !< Integer that specifies where the model Looks for data
158		!! Valid choices are:
159		!! 1 - FMS DataOverride Routine
160		!! 2 - Reserved For Coupler
161		!! 3 - User input (fixed values, useful for 1d testing)
162
163		! Options if using FMS DataOverride Routine
164		character(len=40) :: SurfBandFileName !< Filename if using DataOverride
165		real :: land_speed !< A large Stokes velocity that can be used to indicate land values in
166		!! a data override file [L T-1 ~> m s-1]. Stokes drift components larger
167		!! than this are set to zero in data override calls for the Stokes drift.
168		logical :: DataOver_initialized !< Flag for DataOverride Initialization
169
170		! Options for computing Langmuir number
171		real :: LA_FracHBL !< Fraction of OSBL for averaging Langmuir number [nondim]
172		real :: LA_HBL_min !< Minimum boundary layer depth for averaging Langmuir number [Z ~> m]
173		logical :: LA_Misalignment = .false. !< Flag to use misalignment in Langmuir number
174		logical :: LA_misalign_bug = .false. !< Flag to use code with a sign error when calculating the
175		!! misalignment between the shear and waves in the Langmuir number calculation.
176		real :: g_Earth !< The gravitational acceleration, equivalent to GV%g_Earth but with
177		!! different dimensional rescaling appropriate for deep-water gravity
178		!! waves [Z T-2 ~> m s-2]
179		real :: I_g_Earth !< The inverse of the gravitational acceleration, with dimensional rescaling
180		!! appropriate for deep-water gravity waves [T2 Z-1 ~> s2 m-1]
181		! Surface Wave Dependent 1d/2d/3d vars
182		real, allocatable, dimension(:) :: &
183		Freq_Cen !< Central frequency for wave bands, including a factor of 2*pi [T-1 ~> s-1]
184		real, allocatable, dimension(:) :: &
185		PrescribedSurfStkX !< Surface Stokes drift if prescribed [L T-1 ~> m s-1]
186		real, allocatable, dimension(:) :: &
187		PrescribedSurfStkY !< Surface Stokes drift if prescribed [L T-1 ~> m s-1]
188		real, allocatable, dimension(:,:) :: &
189		La_Turb !< Aligned Turbulent Langmuir number [nondim]
190		!! Horizontal -> H points
191		real, allocatable, dimension(:,:) :: &
192		US0_x !< Surface Stokes Drift (zonal) [L T-1 ~> m s-1]
193		!! Horizontal -> U points
194		real, allocatable, dimension(:,:) :: &
195		US0_y !< Surface Stokes Drift (meridional) [L T-1 ~> m s-1]
196		!! Horizontal -> V points
197		real, allocatable, dimension(:,:,:) :: &
198		STKx0 !< Stokes Drift spectrum (zonal) [L T-1 ~> m s-1]
199		!! Horizontal -> U points
200		!! 3rd dimension -> Freq/Wavenumber
201		real, allocatable, dimension(:,:,:) :: &
202		STKy0 !< Stokes Drift spectrum (meridional) [L T-1 ~> m s-1]
203		!! Horizontal -> V points
204		!! 3rd dimension -> Freq/Wavenumber
205
206		real :: La_min !< An arbitrary lower-bound on the Langmuir number [nondim].
207		!! Langmuir number is sqrt(u_star/u_stokes). When both are small
208		!! but u_star is orders of magnitude smaller, the Langmuir number could
209		!! have unintended consequences. Since both are small it can be safely
210		!! capped to avoid such consequences.
211		real :: La_Stk_backgnd !< A small background Stokes velocity used in the denominator of
212		!! some expressions for the Langmuir number [L T-1 ~> m s-1]
213
214		! Parameters used in estimating the wind speed or wave properties from the friction velocity
215		real :: VonKar = -1.0 !< The von Karman coefficient as used in the MOM_wave_interface module [nondim]
216		real :: rho_air !< A typical density of air at sea level, as used in wave calculations [R ~> kg m-3]
217		real :: nu_air !< The viscosity of air, as used in wave calculations [Z2 T-1 ~> m2 s-1]
218		real :: rho_ocn !< A typical surface density of seawater, as used in wave calculations in
219		!! comparison with the density of air [R ~> kg m-3]. The default is RHO_0.
220		real :: SWH_from_u10sq !< A factor for converting the square of the 10 m wind speed to the
221		!! significant wave height [Z T2 L-2 ~> s2 m-1]
222		real :: Charnock_min !< The minimum value of the Charnock coefficient, which relates the square of
223		!! the air friction velocity divided by the gravitational acceleration to the
224		!! wave roughness length [nondim]
225		real :: Charnock_slope_U10 !< The partial derivative of the Charnock coefficient with the 10 m wind
226		!! speed [T L-1 ~> s m-1]. Note that in eq. 13 of the Edson et al. 2013 describing
227		!! the COARE 3.5 bulk flux algorithm, this slope is given as 0.017. However, 0.0017
228		!! reproduces the curve in their figure 6, so that is the default value used in MOM6.
229		real :: Charnock_intercept !< The intercept of the fit for the Charnock coefficient in the limit of
230		!! no wind [nondim]. Note that this can be negative because CHARNOCK_MIN will keep
231		!! the final value for the Charnock coefficient from being from being negative.
232
233		! Options used with the test profile
234		real :: TP_STKX0 !< Test profile x-stokes drift amplitude [L T-1 ~> m s-1]
235		real :: TP_STKY0 !< Test profile y-stokes drift amplitude [L T-1 ~> m s-1]
236		real :: TP_WVL !< Test profile wavelength [Z ~> m]
237
238		! Options for use with the Donelan et al., 1985 (DHH85) spectrum
239		logical :: WaveAgePeakFreq !< Flag to use wave age to determine the peak frequency with DHH85
240		logical :: StaticWaves !< Flag to disable updating DHH85 Stokes drift
241		logical :: DHH85_is_set !< The if the wave properties have been set when WaveMethod = DHH85.
242		real :: WaveAge !< The fixed wave age used with the DHH85 spectrum [nondim]
243		real :: WaveWind !< Wind speed for the DHH85 spectrum [L T-1 ~> m s-1]
244		real :: omega_min !< Minimum wave frequency with the DHH85 spectrum [T-1 ~> s-1]
245		real :: omega_max !< Maximum wave frequency with the DHH85 spectrum [T-1 ~> s-1]
246
247		type(time_type), pointer :: Time !< A pointer to the ocean model's clock.
248		type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
249		!! timing of diagnostic output.
250
251		!>@{ Diagnostic handles
252		integer, public :: id_PFu_Stokes = -1 , id_PFv_Stokes = -1
253		integer, public :: id_3dstokes_x_from_ddt = -1 , id_3dstokes_y_from_ddt = -1
254		integer :: id_P_deltaStokes_L = -1, id_P_deltaStokes_i = -1
255		integer :: id_surfacestokes_x = -1 , id_surfacestokes_y = -1
256		integer :: id_3dstokes_x = -1 , id_3dstokes_y = -1
257		integer :: id_ddt_3dstokes_x = -1 , id_ddt_3dstokes_y = -1
258		integer :: id_La_turb = -1
259		!>@}
260
261		end type wave_parameters_CS
262
263		! Switches needed in import_stokes_drift
264		!>@{ Enumeration values for the wave method
265		integer, parameter :: TESTPROF = 0, SURFBANDS = 1, DHH85 = 2, LF17 = 3, EFACTOR = 4, NULL_WaveMethod = -99
266		!>@}
267		!>@{ Enumeration values for the wave data source
268		integer, parameter :: DATAOVR = 1, COUPLER = 2, INPUT = 3
269		!>@}
270
271		! Strings for the wave method
272		character*(5), parameter :: NULL_STRING = "EMPTY" !< null wave method string
273		character*(12), parameter :: TESTPROF_STRING = "TEST_PROFILE" !< test profile string
274		character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS" !< surface bands string
275		character*(5), parameter :: DHH85_STRING = "DHH85" !< DHH85 wave method string
276		character*(4), parameter :: LF17_STRING = "LF17" !< LF17 wave method string
277		character*(7), parameter :: EFACTOR_STRING = "EFACTOR" !< EFACTOR (based on vr12-ma) wave method string
278
279		contains
280
281		!> Initializes parameters related to MOM_wave_interface
282	1	subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag)
283		type(time_type), target, intent(in) :: Time !< Model time
284		type(ocean_grid_type), intent(inout) :: G !< Grid structure
285		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
286		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
287		type(param_file_type), intent(in) :: param_file !< Input parameter structure
288		type(wave_parameters_CS), pointer :: CS !< Wave parameter control structure
289		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostic Pointer
290
291		! Local variables
292		character(len=40) :: mdl = "MOM_wave_interface" !< This module's name.
293		! This include declares and sets the variable "version".
294		# include "version_variable.h"
295		character*(13) :: TMPSTRING1, TMPSTRING2
296		character*(12), parameter :: DATAOVR_STRING = "DATAOVERRIDE"
297		character*(7), parameter :: COUPLER_STRING = "COUPLER"
298		character*(5), parameter :: INPUT_STRING = "INPUT"
299		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags
300		logical :: use_waves
301		logical :: StatisticalWaves
302
303		! Dummy Check
304	1	if (.not. associated(CS)) then
305	0	call MOM_error(FATAL, "wave_interface_init called without an associated control structure.")
306	1	return
307		endif
308
309		call get_param(param_file, mdl, "USE_WAVES", use_waves, &
310	1	"If true, enables surface wave modules.", default=.false.)
311
312		! Check if using LA_LI2016
313		call get_param(param_file, mdl, "USE_LA_LI2016", StatisticalWaves, &
314	1	do_not_log=.true.,default=.false.)
315
316	1	if (.not.(use_waves .or. StatisticalWaves)) return
317
318	1	CS%UseWaves = use_waves
319	1	CS%diag => diag
320	1	CS%Time => Time
321
322	1	CS%g_Earth = GV%g_Earth_Z_T2
323	1	CS%I_g_Earth = 1.0 / CS%g_Earth
324
325		! Add any initializations needed here
326	1	CS%DataOver_initialized = .false.
327
328	1	call log_version(param_file, mdl, version)
329
330		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
331		"This sets the default value for the various _ANSWER_DATE parameters.", &
332	1	default=99991231)
333
334		call get_param(param_file, mdl, "WAVE_INTERFACE_ANSWER_DATE", CS%answer_date, &
335		"The vintage of the order of arithmetic and expressions in the surface wave "//&
336		"calculations. Values below 20230101 recover the answers from the end of 2022, "//&
337		"while higher values use updated and more robust forms of the same expressions:\n"//&
338		"\t < 20230101 - Original answers for wave interface routines\n"//&
339		"\t >= 20230101 - More robust expressions for Update_Stokes_Drift\n"//&
340		"\t >= 20230102 - More robust expressions for get_StokesSL_LiFoxKemper\n"//&
341		"\t >= 20230103 - More robust expressions for ust_2_u10_coare3p5", &
342	1	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
343	1	if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
344
345		! Langmuir number Options
346		call get_param(param_file, mdl, "LA_DEPTH_RATIO", CS%LA_FracHBL, &
347		"The depth (normalized by BLD) to average Stokes drift over in "//&
348		"Langmuir number calculation, where La = sqrt(ust/Stokes).", &
349	1	units="nondim", default=0.04)
350		call get_param(param_file, mdl, "LA_DEPTH_MIN", CS%LA_HBL_min, &
351		"The minimum depth over which to average the Stokes drift in the Langmuir "//&
352	1	"number calculation.", units="m", default=0.1, scale=US%m_to_Z)
353
354	1	if (StatisticalWaves) then
355	1	CS%WaveMethod = LF17
356	1	call set_LF17_wave_params(param_file, mdl, GV, US, CS)
357	1	if (.not.use_waves) return
358		else
359	0	CS%WaveMethod = NULL_WaveMethod
360		endif
361
362		! Wave modified physics
363		! Presently these are all in research mode
364		call get_param(param_file, mdl, "LAGRANGIAN_MIXING", CS%LagrangianMixing, &
365		"Flag to use Lagrangian Mixing of momentum", default=.false., &
366	0	do_not_log=.not.use_waves)
367	0	if (CS%LagrangianMixing) then
368		! Force Code Intervention
369	0	call MOM_error(FATAL,"Should you be enabling Lagrangian Mixing? Code not ready.")
370		endif
371		call get_param(param_file, mdl, "STOKES_MIXING", CS%StokesMixing, &
372		"Flag to use Stokes Mixing of momentum", default=.false., &
373	0	do_not_log=.not.use_waves)
374	0	if (CS%StokesMixing) then
375		! Force Code Intervention
376	0	call MOM_error(FATAL, "Should you be enabling Stokes Mixing? Code not ready.")
377		endif
378		call get_param(param_file, mdl, "CORIOLIS_STOKES", CS%CoriolisStokes, &
379		"Flag to use Coriolis Stokes acceleration", default=.false., &
380	0	do_not_log=.not.use_waves)
381	0	if (CS%CoriolisStokes) then
382		! Force Code Intervention
383	0	call MOM_error(FATAL, "Should you be enabling Coriolis-Stokes? Code not ready.")
384		endif
385
386		call get_param(param_file, mdl, "STOKES_VF", CS%Stokes_VF, &
387		"Flag to use Stokes vortex force", &
388	0	default=.false.)
389		call get_param(param_file, mdl, "PASSIVE_STOKES_VF", CS%Passive_Stokes_VF, &
390		"Flag to make Stokes vortex force diagnostic only.", &
391	0	default=.false.)
392		call get_param(param_file, mdl, "STOKES_PGF", CS%Stokes_PGF, &
393		"Flag to use Stokes-induced pressure gradient anomaly", &
394	0	default=.false.)
395		call get_param(param_file, mdl, "ROBUST_STOKES_PGF", CS%robust_Stokes_PGF, &
396		"If true, use expressions to calculate the Stokes-induced pressure gradient "//&
397		"anomalies that are more accurate in the limit of thin layers.", &
398	0	default=.false., do_not_log=.not.CS%Stokes_PGF)
399		!### Change the default for ROBUST_STOKES_PGF to True.
400		call get_param(param_file, mdl, "PASSIVE_STOKES_PGF", CS%Passive_Stokes_PGF, &
401		"Flag to make Stokes-induced pressure gradient anomaly diagnostic only.", &
402	0	default=.false.)
403		call get_param(param_file, mdl, "STOKES_DDT", CS%Stokes_DDT, &
404		"Flag to use Stokes d/dt", &
405	0	default=.false.)
406		call get_param(param_file, mdl, "PASSIVE_STOKES_DDT", CS%Passive_Stokes_DDT, &
407		"Flag to make Stokes d/dt diagnostic only", &
408	0	default=.false.)
409
410		! Get Wave Method and write to integer WaveMethod
411		call get_param(param_file,mdl,"WAVE_METHOD",TMPSTRING1, &
412		"Choice of wave method, valid options include: \n"// &
413		" TEST_PROFILE - Prescribed from surface Stokes drift \n"// &
414		" and a decay wavelength.\n"// &
415		" SURFACE_BANDS - Computed from multiple surface values \n"// &
416		" and decay wavelengths.\n"// &
417		" DHH85 - Uses Donelan et al. 1985 empirical \n"// &
418		" wave spectrum with prescribed values. \n"// &
419		" LF17 - Infers Stokes drift profile from wind \n"// &
420		" speed following Li and Fox-Kemper 2017.\n"// &
421		" EFACTOR - Applies an enhancement factor to the KPP\n"//&
422		" turbulent velocity scale received \n"// &
423		" directly from WW3 and is based on the \n"// &
424		" surface layer and projected Langmuir \n"// &
425		" number (Li 2016)\n", &
426	0	default=NULL_STRING)
427	0	select case (TRIM(TMPSTRING1))
428		case (NULL_STRING)! No Waves
429		call MOM_error(FATAL, "wave_interface_init called with no specified "//&
430	0	"WAVE_METHOD.")
431		case (TESTPROF_STRING)! Test Profile
432	0	CS%WaveMethod = TESTPROF
433		call get_param(param_file, mdl, "TP_STKX_SURF", CS%TP_STKX0, &
434		'Surface Stokes (x) for test profile', &
435	0	units='m/s', default=0.1, scale=US%m_s_to_L_T)
436		call get_param(param_file, mdl, "TP_STKY_SURF", CS%TP_STKY0, &
437		'Surface Stokes (y) for test profile', &
438	0	units='m/s', default=0.0, scale=US%m_s_to_L_T)
439		call get_param(param_file,mdl, "TP_WVL", CS%TP_WVL, &
440		'Wavelength for test profile', &
441	0	units='m', default=50.0, scale=US%m_to_Z)
442		case (SURFBANDS_STRING)! Surface Stokes Drift Bands
443	0	CS%WaveMethod = SURFBANDS
444		call get_param(param_file, mdl, "SURFBAND_MIN_THICK_AVG", CS%Stokes_min_thick_avg, &
445		"A layer thickness below which the cell-center Stokes drift is used instead of "//&
446		"the cell average. This is only used if WAVE_INTERFACE_ANSWER_DATE < 20230101.", &
447	0	units="m", default=0.1, scale=US%m_to_Z, do_not_log=(CS%answer_date>=20230101))
448		call get_param(param_file, mdl, "HOMOGENIZE_SURFBANDS", CS%Homogenize_Surfbands, &
449		"A logical which causes the code to horizontally homogenize the surface band "//&
450		"Stokes drift, which is needed in column mode to avoid round-off differences. "//&
451		"At present it only works with DATAOVERRIDE, and is not coded for COUPLER.",&
452	0	default=.false.)
453		call get_param(param_file, mdl, "SURFBAND_SOURCE", TMPSTRING2, &
454		"Choice of SURFACE_BANDS data mode, valid options include: \n"//&
455		" DATAOVERRIDE - Read from NetCDF using FMS DataOverride. \n"//&
456		" COUPLER - Look for variables from coupler pass \n"//&
457	0	" INPUT - Testing with fixed values.", default=NULL_STRING)
458	0	select case (TRIM(TMPSTRING2))
459		case (NULL_STRING)! Default
460		call MOM_error(FATAL, "wave_interface_init called with SURFACE_BANDS "//&
461	0	"but no SURFBAND_SOURCE.")
462		case (DATAOVR_STRING)! Using Data Override
463	0	CS%DataSource = DATAOVR
464		call get_param(param_file, mdl, "SURFBAND_FILENAME", CS%SurfBandFileName, &
465	0	"Filename of surface Stokes drift input band data.", default="StkSpec.nc")
466		call get_param(param_file, mdl, "SURFBAND_OVERRIDE_LAND_SPEED", CS%land_speed, &
467		"A large Stokes velocity that can be used to indicate land values in "//&
468		"a data override file. Stokes drift components larger than this are "//&
469		"set to zero in data override calls for the Stokes drift.", &
470	0	units="m s-1", default=10.0, scale=US%m_s_to_L_T)
471		case (COUPLER_STRING)! Reserved for coupling
472	0	CS%DataSource = COUPLER
473		! This is just to make something work, but it needs to be read from the wavemodel.
474		call get_param(param_file, mdl, "STK_BAND_COUPLER",CS%NumBands, &
475		"STK_BAND_COUPLER is the number of Stokes drift bands in the coupler. "//&
476		"This has to be consistent with the number of Stokes drift bands in WW3, "//&
477	0	"or the model will fail.", default=1)
478	0	allocate( CS%WaveNum_Cen(CS%NumBands), source=0.0 )
479	0	allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NumBands), source=0.0 )
480	0	allocate( CS%STKy0(G%isd:G%ied,G%jsdB:G%jedB,CS%NumBands), source=0.0 )
481	0	allocate( CS%UStk_Hb(G%isc:G%iec,G%jsc:G%jec,CS%NumBands), source=0.0 )
482	0	allocate( CS%VStk_Hb(G%isc:G%iec,G%jsc:G%jec,CS%NumBands), source=0.0 )
483	0	allocate( CS%Omega_w2x(G%isc:G%iec,G%jsc:G%jec) , source=0.0 )
484	0	CS%PartitionMode = 0
485		call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, &
486		"Central wavenumbers for surface Stokes drift bands.", &
487	0	units='rad/m', default=0.12566, scale=US%Z_to_m)
488		case (INPUT_STRING)! A method to input the Stokes band (globally uniform)
489	0	CS%DataSource = INPUT
490		call get_param(param_file, mdl, "SURFBAND_NB", CS%NumBands, &
491		"Prescribe number of wavenumber bands for Stokes drift. "//&
492		"Make sure this is consistnet w/ WAVENUMBERS, STOKES_X, and "//&
493	0	"STOKES_Y, there are no safety checks in the code.", default=1)
494	0	allocate( CS%WaveNum_Cen(1:CS%NumBands), source=0.0 )
495	0	allocate( CS%PrescribedSurfStkX(1:CS%NumBands), source=0.0 )
496	0	allocate( CS%PrescribedSurfStkY(1:CS%NumBands), source=0.0 )
497	0	allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,1:CS%NumBands), source=0.0 )
498	0	allocate( CS%STKy0(G%isd:G%ied,G%jsdB:G%jedB,1:CS%NumBands), source=0.0 )
499
500	0	CS%PartitionMode = 0
501		call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", CS%WaveNum_Cen, &
502		"Central wavenumbers for surface Stokes drift bands.", &
503	0	units='rad/m', default=0.12566, scale=US%Z_to_m)
504		call get_param(param_file, mdl, "SURFBAND_STOKES_X", CS%PrescribedSurfStkX, &
505		"X-direction surface Stokes drift for bands.", &
506	0	units='m/s', default=0.15, scale=US%m_s_to_L_T)
507		call get_param(param_file, mdl, "SURFBAND_STOKES_Y", CS%PrescribedSurfStkY, &
508		"Y-direction surface Stokes drift for bands.", &
509	0	units='m/s', default=0.0, scale=US%m_s_to_L_T)
510		case default! No method provided
511	0	call MOM_error(FATAL,'Check WAVE_METHOD.')
512		end select
513
514		case (DHH85_STRING) !Donelan et al., 1985 spectrum
515	0	CS%WaveMethod = DHH85
516		call MOM_error(WARNING,"DHH85 only ever set-up for uniform cases w/ "//&
517	0	"Stokes drift in x-direction.")
518		call get_param(param_file, mdl, "DHH85_AGE_FP", CS%WaveAgePeakFreq, &
519	0	"Choose true to use waveage in peak frequency.", default=.false.)
520		call get_param(param_file, mdl, "DHH85_AGE", CS%WaveAge, &
521		"Wave Age for DHH85 spectrum.", &
522	0	units='nondim', default=1.2)
523		call get_param(param_file, mdl, "DHH85_WIND", CS%WaveWind, &
524		"Wind speed for DHH85 spectrum.", &
525	0	units='m s-1', default=10.0, scale=US%m_s_to_L_T)
526		call get_param(param_file, mdl, "DHH85_MIN_WAVE_FREQ", CS%omega_min, &
527		"Minimum wave frequency for the DHH85 spectrum.", &
528	0	units='s-1', default=0.1, scale=US%T_to_s)
529		call get_param(param_file, mdl, "DHH85_MAX_WAVE_FREQ", CS%omega_max, &
530		"Maximum wave frequency for the DHH85 spectrum.", &
531	0	units='s-1', default=10.0, scale=US%T_to_s) ! The default is about a 30 cm cutoff wavelength.
532		call get_param(param_file, mdl, "STATIC_DHH85", CS%StaticWaves, &
533	0	"Flag to disable updating DHH85 Stokes drift.", default=.false.)
534		case (LF17_STRING) !Li and Fox-Kemper 17 wind-sea Langmuir number
535	0	CS%WaveMethod = LF17
536	0	call set_LF17_wave_params(param_file, mdl, GV, US, CS)
537		case (EFACTOR_STRING) !Li and Fox-Kemper 16
538	0	CS%WaveMethod = EFACTOR
539		case default
540	0	call MOM_error(FATAL,'Check WAVE_METHOD.')
541		end select
542
543		! Langmuir number Options (Note that CS%LA_FracHBL is set above.)
544		call get_param(param_file, mdl, "LA_MISALIGNMENT", CS%LA_Misalignment, &
545		"Flag (logical) if using misalignment between shear and waves in LA", &
546	0	default=.false.)
547		call get_param(param_file, mdl, "LA_MISALIGNMENT_BUG", CS%LA_misalign_bug, &
548		"If true, use a code with a sign error when calculating the misalignment between "//&
549		"the shear and waves when LA_MISALIGNMENT is true.", &
550	0	default=.false., do_not_log=.not.CS%LA_Misalignment)
551		call get_param(param_file, mdl, "MIN_LANGMUIR", CS%La_min, &
552		"A minimum value for all Langmuir numbers that is not physical, "//&
553		"but is likely only encountered when the wind is very small and "//&
554		"therefore its effects should be mostly benign.", &
555	0	units="nondim", default=0.05)
556		call get_param(param_file, mdl, "LANGMUIR_STOKES_BACKGROUND", CS%La_Stk_backgnd, &
557		"A small background Stokes velocity used in the denominator of some "//&
558		"expressions for the Langmuir number.", &
559	0	units="m s-1", default=1.0e-10, scale=US%m_s_to_L_T, do_not_log=(CS%WaveMethod==LF17))
560
561		! Allocate and initialize
562		! a. Stokes driftProfiles
563	0	allocate(CS%Us_x(G%isdB:G%IedB,G%jsd:G%jed,G%ke), source=0.0)
564	0	allocate(CS%Us_y(G%isd:G%Ied,G%jsdB:G%jedB,G%ke), source=0.0)
565	0	if (CS%Stokes_DDT) then
566		!allocate(CS%Us_x_prev(G%isdB:G%IedB,G%jsd:G%jed,G%ke), source=0.0)
567		!allocate(CS%Us_y_prev(G%isd:G%Ied,G%jsdB:G%jedB,G%ke), source=0.0)
568	0	allocate(CS%ddt_Us_x(G%isdB:G%IedB,G%jsd:G%jed,G%ke), source=0.0)
569	0	allocate(CS%ddt_Us_y(G%isd:G%Ied,G%jsdB:G%jedB,G%ke), source=0.0)
570	0	allocate(CS%Us_x_from_ddt(G%isdB:G%IedB,G%jsd:G%jed,G%ke), source=0.0)
571	0	allocate(CS%Us_y_from_ddt(G%isd:G%Ied,G%jsdB:G%jedB,G%ke), source=0.0)
572		endif
573		! b. Surface Values
574	0	allocate(CS%US0_x(G%isdB:G%iedB,G%jsd:G%jed), source=0.0)
575	0	allocate(CS%US0_y(G%isd:G%ied,G%jsdB:G%jedB), source=0.0)
576		! c. Langmuir number
577	0	allocate(CS%La_turb(G%isc:G%iec,G%jsc:G%jec), source=0.0)
578		! d. Viscosity for Stokes drift
579	0	if (CS%StokesMixing) then
580	0	allocate(CS%KvS(G%isd:G%Ied,G%jsd:G%jed,GV%ke), source=0.0)
581		endif
582
583		! Initialize Wave related outputs
584		CS%id_surfacestokes_y = register_diag_field('ocean_model','surface_stokes_y', &
585	0	CS%diag%axesCv1,Time,'Surface Stokes drift (y)', 'm s-1', conversion=US%L_T_to_m_s)
586		CS%id_surfacestokes_x = register_diag_field('ocean_model','surface_stokes_x', &
587	0	CS%diag%axesCu1,Time,'Surface Stokes drift (x)', 'm s-1', conversion=US%L_T_to_m_s)
588		CS%id_3dstokes_y = register_diag_field('ocean_model','3d_stokes_y', &
589	0	CS%diag%axesCvL,Time,'3d Stokes drift (y)', 'm s-1', conversion=US%L_T_to_m_s)
590		CS%id_3dstokes_x = register_diag_field('ocean_model','3d_stokes_x', &
591	0	CS%diag%axesCuL,Time,'3d Stokes drift (x)', 'm s-1', conversion=US%L_T_to_m_s)
592	0	if (CS%Stokes_DDT) then
593		CS%id_ddt_3dstokes_y = register_diag_field('ocean_model','dvdt_Stokes', &
594	0	CS%diag%axesCvL,Time,'d/dt Stokes drift (meridional)', 'm s-2', conversion=US%L_T2_to_m_s2)
595		CS%id_ddt_3dstokes_x = register_diag_field('ocean_model','dudt_Stokes', &
596	0	CS%diag%axesCuL,Time,'d/dt Stokes drift (zonal)', 'm s-2', conversion=US%L_T2_to_m_s2)
597		CS%id_3dstokes_y_from_ddt = register_diag_field('ocean_model','3d_stokes_y_from_ddt', &
598	0	CS%diag%axesCvL,Time,'3d Stokes drift from ddt (y)', 'm s-1', conversion=US%L_T_to_m_s)
599		CS%id_3dstokes_x_from_ddt = register_diag_field('ocean_model','3d_stokes_x_from_ddt', &
600	0	CS%diag%axesCuL,Time,'3d Stokes drift from ddt (x)', 'm s-1', conversion=US%L_T_to_m_s)
601		endif
602		CS%id_PFv_Stokes = register_diag_field('ocean_model','PFv_Stokes', &
603	0	CS%diag%axesCvL,Time,'PF from Stokes drift (meridional)','m s-2',conversion=US%L_T2_to_m_s2)
604		CS%id_PFu_Stokes = register_diag_field('ocean_model','PFu_Stokes', &
605	0	CS%diag%axesCuL,Time,'PF from Stokes drift (zonal)','m s-2',conversion=US%L_T2_to_m_s2)
606		CS%id_P_deltaStokes_i = register_diag_field('ocean_model','P_deltaStokes_i', &
607		CS%diag%axesTi,Time,'Interfacial pressure anomaly from Stokes drift used in PFu_Stokes',&
608	0	'm2 s-2',conversion=US%L_T_to_m_s**2)
609		CS%id_P_deltaStokes_L = register_diag_field('ocean_model','P_deltaStokes_L', &
610		CS%diag%axesTL,Time,'Layer averaged pressure anomaly from Stokes drift used in PFu_Stokes',&
611	0	'm2 s-2',conversion=US%L_T_to_m_s**2)
612		CS%id_La_turb = register_diag_field('ocean_model','La_turbulent', &
613	0	CS%diag%axesT1,Time,'Surface (turbulent) Langmuir number','nondim')
614
615		end subroutine MOM_wave_interface_init
616
617		!> Set the parameters that are used to determine the averaged Stokes drift and Langmuir numbers
618	1	subroutine set_LF17_wave_params(param_file, mdl, GV, US, CS)
619		type(param_file_type), intent(in) :: param_file !< Input parameter structure
620		character(len=*), intent(in) :: mdl !< A module name to use in the get_param calls
621		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
622		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
623		type(wave_parameters_CS), pointer :: CS !< Wave parameter control structure
624
625		! A separate routine is used to set these parameters because there are multiple ways that the
626		! underlying parameterizations are enabled.
627
628		call get_param(param_file, mdl, "VISCOSITY_AIR", CS%nu_air, &
629		"A typical viscosity of air at sea level, as used in wave calculations", &
630	1	units="m2 s-1", default=1.0e-6, scale=US%m2_s_to_Z2_T)
631		call get_param(param_file, mdl, "VON_KARMAN_WAVES", CS%vonKar, &
632		"The value the von Karman constant as used for surface wave calculations.", &
633	1	units="nondim", default=0.40) ! The default elsewhere in MOM6 is usually 0.41.
634		call get_param(param_file, mdl, "RHO_AIR", CS%rho_air, &
635		"A typical density of air at sea level, as used in wave calculations", &
636	1	units="kg m-3", default=1.225, scale=US%kg_m3_to_R)
637		call get_param(param_file, mdl, "RHO_SFC_WAVES", CS%Rho_ocn, &
638		"A typical surface density of seawater, as used in wave calculations in "//&
639		"comparison with the density of air. The default is RHO_0.", &
640	1	units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R)
641		call get_param(param_file, mdl, "WAVE_HEIGHT_SCALE_FACTOR", CS%SWH_from_u10sq, &
642		"A factor relating the square of the 10 m wind speed to the significant "//&
643		"wave height, with a default value based on the Pierson-Moskowitz spectrum.", &
644	1	units="s2 m-1", default=0.0246, scale=US%m_to_ZUS%L_T_to_m_s*2)
645		call get_param(param_file, mdl, "CHARNOCK_MIN", CS%Charnock_min, &
646		"The minimum value of the Charnock coefficient, which relates the square of "//&
647		"the air friction velocity divided by the gravitational acceleration to the "//&
648	1	"wave roughness length.", units="nondim", default=0.028)
649		call get_param(param_file, mdl, "CHARNOCK_SLOPE_U10", CS%Charnock_slope_U10, &
650		"The partial derivative of the Charnock coefficient with the 10 m wind speed. "//&
651		"Note that in eq. 13 of the Edson et al. 2013 describing the COARE 3.5 bulk "//&
652		"flux algorithm, this slope is given as 0.017. However, 0.0017 reproduces "//&
653		"the curve in their figure 6, so that is the default value used in MOM6.", &
654	1	units="s m-1", default=0.0017, scale=US%L_T_to_m_s)
655		call get_param(param_file, mdl, "CHARNOCK_0_WIND_INTERCEPT", CS%Charnock_intercept, &
656		"The intercept of the fit for the Charnock coefficient in the limit of no wind. "//&
657		"Note that this can be negative because CHARNOCK_MIN will keep the final "//&
658		"value for the Charnock coefficient from being from being negative.", &
659	1	units="nondim", default=-0.005)
660
661	1	end subroutine set_LF17_wave_params
662
663		!> This interface provides the caller with information from the waves control structure.
664	0	subroutine query_wave_properties(CS, NumBands, WaveNumbers, US)
665		type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure
666		integer, optional, intent(out) :: NumBands !< If present, this returns the number of
667		!!< wavenumber partitions in the wave discretization
668		real, dimension(:), optional, intent(out) :: Wavenumbers !< If present this returns the characteristic
669		!! wavenumbers of the wave discretization [m-1] or [Z-1 ~> m-1]
670		type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type that is used to undo
671		!! the dimensional scaling of the output variables, if present
672		integer :: n
673
674	0	if (present(NumBands)) NumBands = CS%NumBands
675	0	if (present(Wavenumbers)) then
676	0	if (size(Wavenumbers) < CS%NumBands) call MOM_error(FATAL, "query_wave_properties called "//&
677	0	"with a Wavenumbers array that is smaller than the number of bands.")
678	0	if (present(US)) then
679	0	do n=1,CS%NumBands ; Wavenumbers(n) = US%m_to_Z * CS%WaveNum_Cen(n) ; enddo
680		else
681	0	do n=1,CS%NumBands ; Wavenumbers(n) = CS%WaveNum_Cen(n) ; enddo
682		endif
683		endif
684
685	0	end subroutine query_wave_properties
686
687		!> Subroutine that handles updating of surface wave/Stokes drift related properties
688	0	subroutine Update_Surface_Waves(G, GV, US, Time_present, dt, CS, forces)
689		type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure
690		type(ocean_grid_type), intent(inout) :: G !< Grid structure
691		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
692		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
693		type(time_type), intent(in) :: Time_present !< Model Time
694		type(time_type), intent(in) :: dt !< Time increment as a time-type
695		type(mech_forcing), intent(in), optional :: forces !< MOM_forcing_type
696		! Local variables
697		type(time_type) :: Stokes_Time
698		integer :: i, j, b
699
700	0	if (CS%WaveMethod == TESTPROF) then
701		! Do nothing
702	0	elseif (CS%WaveMethod == SURFBANDS) then
703	0	if (CS%DataSource == DATAOVR) then
704		! Updating Stokes drift time to center of time increment.
705		! This choice makes sense for the thermodynamics, but for the
706		! dynamics it may be more useful to update to the end of the
707		! time increment.
708	0	Stokes_Time = Time_present + dt/2
709	0	call Surface_Bands_by_data_override(Stokes_Time, G, GV, US, CS)
710	0	elseif (CS%DataSource == COUPLER) then
711	0	if (.not.present(FORCES)) then
712		call MOM_error(FATAL,"The option SURFBAND = COUPLER can not be used with "//&
713		"this driver. If you are using a coupled driver with a wave model then "//&
714		"check the arguments in the subroutine call to Update_Surface_Waves, "//&
715	0	"otherwise select another option for SURFBAND_SOURCE.")
716		endif
717	0	if (size(CS%WaveNum_Cen) /= size(forces%stk_wavenumbers)) then
718		call MOM_error(FATAL, "Number of wavenumber bands in WW3 does not match that in MOM6. "//&
719		"Make sure that STK_BAND_COUPLER in MOM6 input is equal to the number of bands in "//&
720	0	"ww3_grid.inp, and that your mod_def.ww3 is up to date.")
721		endif
722
723	0	do b=1,CS%NumBands
724	0	CS%WaveNum_Cen(b) = forces%stk_wavenumbers(b)
725		!Interpolate from a grid to c grid
726	0	do j=G%jsc,G%jec
727	0	do I=G%iscB,G%iecB
728	0	CS%STKx0(I,j,b) = 0.5*(forces%UStkb(i,j,b)+forces%UStkb(i+1,j,b))
729		enddo
730		enddo
731	0	do J=G%jscB,G%jecB
732	0	do i=G%isc,G%iec
733	0	CS%STKY0(i,J,b) = 0.5*(forces%VStkb(i,j,b)+forces%VStkb(i,j+1,b))
734		enddo
735		enddo
736	0	call pass_vector(CS%STKx0(:,:,b),CS%STKy0(:,:,b), G%Domain)
737		enddo
738	0	do j=G%jsc,G%jec
739	0	do i=G%isc,G%iec
740	0	CS%Omega_w2x(i,j) = forces%omega_w2x(i,j)
741	0	do b=1,CS%NumBands
742	0	CS%UStk_Hb(i,j,b) = forces%UStkb(i,j,b)
743	0	CS%VStk_Hb(i,j,b) = forces%VStkb(i,j,b)
744		enddo
745		enddo
746		enddo
747	0	elseif (CS%DataSource == INPUT) then
748	0	do b=1,CS%NumBands
749	0	do j=G%jsd,G%jed
750	0	do I=G%isdB,G%iedB
751	0	CS%STKx0(I,j,b) = CS%PrescribedSurfStkX(b)
752		enddo
753		enddo
754	0	do J=G%jsdB, G%jedB
755	0	do i=G%isd,G%ied
756	0	CS%STKY0(i,J,b) = CS%PrescribedSurfStkY(b)
757		enddo
758		enddo
759		enddo
760		endif
761		endif
762
763	0	end subroutine Update_Surface_Waves
764
765		!> Constructs the Stokes Drift profile on the model grid based on
766		!! desired coupling options
767	0	subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step)
768		type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure
769		type(ocean_grid_type), intent(inout) :: G !< Grid structure
770		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
771		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
772		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
773		intent(in) :: dz !< Thickness in height units [Z ~> m]
774		real, dimension(SZI_(G),SZJ_(G)), &
775		intent(in) :: ustar !< Wind friction velocity [Z T-1 ~> m s-1].
776		real, intent(in) :: dt !< Time-step for computing Stokes-tendency [T ~> s]
777		logical, intent(in) :: dynamics_step !< True if this call is on a dynamics step
778
779		! Local Variables
780		real :: Top, MidPoint, Bottom ! Positions within the layer [Z ~> m]
781		real :: level_thick ! The thickness of each layer [Z ~> m]
782		real :: DecayScale ! A vertical decay scale in the test profile [Z-1 ~> m-1]
783		real :: CMN_FAC ! A nondimensional factor [nondim]
784		real :: WN ! Model wavenumber [Z-1 ~> m-1]
785		real :: UStokes ! A Stokes drift velocity [L T-1 ~> m s-1]
786		real :: PI ! 3.1415926535... [nondim]
787		real :: La ! The local Langmuir number [nondim]
788		integer :: i, j, k, b
789		real :: I_dt ! The inverse of the time step [T-1 ~> s-1]
790
791	0	if (CS%WaveMethod==EFACTOR) return
792
793	0	if (allocated(CS%US_x) .and. allocated(CS%US_y)) then
794	0	call pass_vector(CS%US_x(:,:,:),CS%US_y(:,:,:), G%Domain)
795		endif
796
797		! 1. If Test Profile Option is chosen
798		! Computing mid-point value from surface value and decay wavelength
799	0	if (CS%WaveMethod==TESTPROF) then
800	0	PI = 4.0*atan(1.0)
801	0	DecayScale = 4.*PI / CS%TP_WVL !4pi
802	0	do j=G%jsc,G%jec
803	0	do I=G%iscB,G%iecB
804	0	Bottom = 0.0
805	0	MidPoint = 0.0
806	0	do k = 1,GV%ke
807	0	Top = Bottom
808	0	MidPoint = Bottom - 0.25*(dz(i,j,k)+dz(i+1,j,k))
809	0	Bottom = Bottom - 0.5*(dz(i,j,k)+dz(i+1,j,k))
810	0	CS%Us_x(I,j,k) = CS%TP_STKX0exp(MidPointDecayScale)
811		enddo
812		enddo
813		enddo
814	0	do J=G%jscB,G%jecB
815	0	do i=G%isc,G%iec
816	0	Bottom = 0.0
817	0	MidPoint = 0.0
818	0	do k = 1,GV%ke
819	0	Top = Bottom
820	0	MidPoint = Bottom - 0.25*(dz(i,j,k)+dz(i,j+1,k))
821	0	Bottom = Bottom - 0.5*(dz(i,j,k)+dz(i,j+1,k))
822	0	CS%Us_y(i,J,k) = CS%TP_STKY0exp(MidPointDecayScale)
823		enddo
824		enddo
825		enddo
826	0	call pass_vector(CS%Us_x(:,:,:),CS%Us_y(:,:,:), G%Domain, To_All)
827		! 2. If Surface Bands is chosen
828		! In wavenumber mode compute integral for layer averaged Stokes drift.
829		! In frequency mode compuate value at midpoint.
830	0	elseif (CS%WaveMethod==SURFBANDS) then
831	0	CS%Us_x(:,:,:) = 0.0
832	0	CS%Us_y(:,:,:) = 0.0
833	0	CS%Us0_x(:,:) = 0.0
834	0	CS%Us0_y(:,:) = 0.0
835		! Computing X direction Stokes drift
836	0	do j=G%jsc,G%jec
837	0	do I=G%iscB,G%iecB
838		! 1. First compute the surface Stokes drift
839		! by summing over the partitions.
840	0	do b = 1,CS%NumBands
841	0	CS%US0_x(I,j) = CS%US0_x(I,j) + CS%STKx0(I,j,b)
842		enddo
843		! 2. Second compute the level averaged Stokes drift
844	0	bottom = 0.0
845	0	do k = 1,GV%ke
846	0	Top = Bottom
847	0	level_thick = 0.5*(dz(i,j,k)+dz(i+1,j,k))
848	0	MidPoint = Top - 0.5*level_thick
849	0	Bottom = Top - level_thick
850
851	0	if (CS%answer_date >= 20230101) then
852		! Use more accurate and numerically stable expressions that work even for vanished layers.
853	0	do b = 1,CS%NumBands
854	0	if (CS%PartitionMode == 0) then
855		! Average over a layer using the bin's central wavenumber.
856	0	CMN_FAC = exp(2.CS%WaveNum_Cen(b)Top) * one_minus_exp_x(2.CS%WaveNum_Cen(b)level_thick)
857		else
858		! Use an analytic expression for the average of an exponential over a layer
859	0	WN = CS%Freq_Cen(b)*2 CS%I_g_Earth
860	0	CMN_FAC = exp(2.WNTop) * one_minus_exp_x(2.WNlevel_thick)
861		endif
862	0	CS%US_x(I,j,k) = CS%US_x(I,j,k) + CS%STKx0(I,j,b)*CMN_FAC
863		enddo
864
865	0	elseif (level_thick > CS%Stokes_min_thick_avg) then
866		! -> Stokes drift in thin layers not averaged.
867	0	do b = 1,CS%NumBands
868	0	if (CS%PartitionMode == 0) then
869		! In wavenumber we are averaging over level
870		CMN_FAC = (exp(Top2.CS%WaveNum_Cen(b))-exp(Bottom2.CS%WaveNum_Cen(b))) &
871	0	/ ((Top-Bottom)(2.CS%WaveNum_Cen(b)))
872		else
873		! Use a numerical integration and then divide by layer thickness
874	0	WN = CS%Freq_Cen(b)**2 / CS%g_Earth !bgr bug-fix missing g
875	0	CMN_FAC = (exp(2.WNTop)-exp(2.WNBottom)) / (2.WN(Top-Bottom))
876		endif
877	0	CS%US_x(I,j,k) = CS%US_x(I,j,k) + CS%STKx0(I,j,b)*CMN_FAC
878		enddo
879		else ! Take the value at the midpoint
880	0	do b = 1,CS%NumBands
881	0	if (CS%PartitionMode == 0) then
882	0	CMN_FAC = exp(MidPoint * 2. * CS%WaveNum_Cen(b))
883		else
884	0	CMN_FAC = exp(MidPoint * 2. * CS%Freq_Cen(b)**2 / CS%g_Earth)
885		endif
886	0	CS%US_x(I,j,k) = CS%US_x(I,j,k) + CS%STKx0(I,j,b)*CMN_FAC
887		enddo
888		endif
889		enddo
890		enddo
891		enddo
892
893		! Computing Y direction Stokes drift
894	0	do J=G%jscB,G%jecB
895	0	do i=G%isc,G%iec
896		! Set the surface value to that at z=0
897	0	do b = 1,CS%NumBands
898	0	CS%US0_y(i,J) = CS%US0_y(i,J) + CS%STKy0(i,J,b)
899		enddo
900		! Compute the level averages.
901	0	bottom = 0.0
902	0	do k = 1,GV%ke
903	0	Top = Bottom
904	0	level_thick = 0.5*(dz(i,j,k)+dz(i,j+1,k))
905	0	MidPoint = Top - 0.5*level_thick
906	0	Bottom = Top - level_thick
907
908	0	if (CS%answer_date >= 20230101) then
909		! Use more accurate and numerically stable expressions that work even for vanished layers.
910	0	do b = 1,CS%NumBands
911	0	if (CS%PartitionMode == 0) then
912		! Average over a layer using the bin's central wavenumber.
913	0	CMN_FAC = exp(2.CS%WaveNum_Cen(b)Top) * one_minus_exp_x(2.CS%WaveNum_Cen(b)level_thick)
914		else
915		! Use an analytic expression for the average of an exponential over a layer
916	0	WN = CS%Freq_Cen(b)*2 CS%I_g_Earth
917	0	CMN_FAC = exp(2.WNTop) * one_minus_exp_x(2.WNlevel_thick)
918		endif
919	0	CS%US_y(i,J,k) = CS%US_y(i,J,k) + CS%STKy0(i,J,b)*CMN_FAC
920		enddo
921	0	elseif (level_thick > CS%Stokes_min_thick_avg) then
922		! -> Stokes drift in thin layers not averaged.
923	0	do b = 1,CS%NumBands
924	0	if (CS%PartitionMode == 0) then
925		! In wavenumber we are averaging over level
926		CMN_FAC = (exp(Top2.CS%WaveNum_Cen(b))-exp(Bottom2.CS%WaveNum_Cen(b))) &
927	0	/ ((Top-Bottom)(2.CS%WaveNum_Cen(b)))
928		else
929		! Use a numerical integration and then divide by layer thickness
930	0	WN = CS%Freq_Cen(b)**2 / CS%g_Earth !bgr bug-fix missing g
931	0	CMN_FAC = (exp(2.WNTop)-exp(2.WNBottom)) / (2.WN(Top-Bottom))
932		endif
933	0	CS%US_y(i,J,k) = CS%US_y(i,J,k) + CS%STKy0(i,J,b)*CMN_FAC
934		enddo
935		else ! Take the value at the midpoint
936	0	do b = 1,CS%NumBands
937	0	if (CS%PartitionMode == 0) then
938	0	CMN_FAC = exp(MidPoint2.CS%WaveNum_Cen(b))
939		else
940	0	CMN_FAC = exp(MidPoint * 2. * CS%Freq_Cen(b)**2 / CS%g_Earth)
941		endif
942	0	CS%US_y(i,J,k) = CS%US_y(i,J,k) + CS%STKy0(i,J,b)*CMN_FAC
943		enddo
944		endif
945		enddo
946		enddo
947		enddo
948	0	call pass_vector(CS%Us_x(:,:,:),CS%Us_y(:,:,:), G%Domain, To_All)
949	0	call pass_vector(CS%Us0_x(:,:),CS%Us0_y(:,:), G%Domain)
950	0	elseif (CS%WaveMethod == DHH85) then
951	0	if (.not.(CS%StaticWaves .and. CS%DHH85_is_set)) then
952	0	do j=G%jsc,G%jec
953	0	do I=G%iscB,G%iecB
954	0	bottom = 0.0
955	0	do k = 1,GV%ke
956	0	Top = Bottom
957	0	MidPoint = Top - 0.25*(dz(i,j,k)+dz(i+1,j,k))
958	0	Bottom = Top - 0.5*(dz(i,j,k)+dz(i+1,j,k))
959		!bgr note that this is using a u-point I on h-point ustar
960		! this code has only been previous used for uniform
961		! grid cases. This needs fixed if DHH85 is used for non
962		! uniform cases.
963	0	call DHH85_mid(GV, US, CS, MidPoint, UStokes)
964		! Putting into x-direction (no option for direction
965	0	CS%US_x(I,j,k) = UStokes
966		enddo
967		enddo
968		enddo
969	0	do J=G%jscB,G%jecB
970	0	do i=G%isc,G%iec
971	0	Bottom = 0.0
972	0	do k = 1,GV%ke
973	0	Top = Bottom
974	0	MidPoint = Bottom - 0.25*(dz(i,j,k)+dz(i,j+1,k))
975	0	Bottom = Bottom - 0.5*(dz(i,j,k)+dz(i,j+1,k))
976		!bgr note that this is using a v-point J on h-point ustar
977		! this code has only been previous used for uniform
978		! grid cases. This needs fixed if DHH85 is used for non
979		! uniform cases.
980		! call DHH85_mid(GV, US, CS, Midpoint, UStokes)
981		! Putting into x-direction, so setting y direction to 0
982	0	CS%US_y(i,J,k) = 0.0
983		! For rotational symmetry there should be the option for this to become = UStokes
984		! bgr - see note above, but this is true
985		! if this is used for anything
986		! other than simple LES comparison
987		enddo
988		enddo
989		enddo
990	0	CS%DHH85_is_set = .true.
991		endif
992	0	call pass_vector(CS%Us_x(:,:,:),CS%Us_y(:,:,:), G%Domain)
993		else! Keep this else, fallback to 0 Stokes drift
994	0	CS%Us_x(:,:,:) = 0.
995	0	CS%Us_y(:,:,:) = 0.
996		endif
997
998		! Turbulent Langmuir number is computed here and available to use anywhere.
999		! SL Langmuir number requires mixing layer depth, and therefore is computed
1000		! in the routine it is needed by (e.g. KPP or ePBL).
1001	0	do j=G%jsc, G%jec
1002	0	do i=G%isc,G%iec
1003		call get_Langmuir_Number( La, G, GV, US, dz(i,j,1), ustar(i,j), i, j, &
1004	0	dz(i,j,:), CS, Override_MA=.false.)
1005	0	CS%La_turb(i,j) = La
1006		enddo
1007		enddo
1008
1009		! Finding tendency of Stokes drift over the time step to apply
1010		! as an acceleration to the models current.
1011	0	if ( dynamics_step .and. CS%Stokes_DDT ) then
1012	0	I_dt = 1.0 / dt
1013	0	CS%ddt_us_x(:,:,:) = (CS%US_x(:,:,:) - CS%US_x_prev(:,:,:)) * I_dt
1014	0	CS%ddt_us_y(:,:,:) = (CS%US_y(:,:,:) - CS%US_y_prev(:,:,:)) * I_dt
1015	0	CS%US_x_prev(:,:,:) = CS%US_x(:,:,:)
1016	0	CS%US_y_prev(:,:,:) = CS%US_y(:,:,:)
1017		endif
1018
1019		! Output any desired quantities
1020	0	if (CS%id_surfacestokes_y>0) &
1021	0	call post_data(CS%id_surfacestokes_y, CS%us0_y, CS%diag)
1022	0	if (CS%id_surfacestokes_x>0) &
1023	0	call post_data(CS%id_surfacestokes_x, CS%us0_x, CS%diag)
1024	0	if (CS%id_3dstokes_y>0) &
1025	0	call post_data(CS%id_3dstokes_y, CS%us_y, CS%diag)
1026	0	if (CS%id_3dstokes_x>0) &
1027	0	call post_data(CS%id_3dstokes_x, CS%us_x, CS%diag)
1028	0	if (CS%Stokes_DDT) then
1029	0	if (CS%id_ddt_3dstokes_x>0) &
1030	0	call post_data(CS%id_ddt_3dstokes_x, CS%ddt_us_x, CS%diag)
1031	0	if (CS%id_ddt_3dstokes_y>0) &
1032	0	call post_data(CS%id_ddt_3dstokes_y, CS%ddt_us_y, CS%diag)
1033	0	if (CS%id_3dstokes_x_from_ddt>0) &
1034	0	call post_data(CS%id_3dstokes_x_from_ddt, CS%us_x_from_ddt, CS%diag)
1035	0	if (CS%id_3dstokes_y_from_ddt>0) &
1036	0	call post_data(CS%id_3dstokes_y_from_ddt, CS%us_y_from_ddt, CS%diag)
1037		endif
1038	0	if (CS%id_La_turb>0) &
1039	0	call post_data(CS%id_La_turb, CS%La_turb, CS%diag)
1040
1041		end subroutine Update_Stokes_Drift
1042
1043		!> Return the value of (1 - exp(-x))/x [nondim], using an accurate expression for small values of x.
1044	626926	real function one_minus_exp_x(x)
1045		real, intent(in) :: x !< The argument of the function ((1 - exp(-x))/x) [nondim]
1046		real, parameter :: C1_6 = 1.0/6.0 ! A rational fraction [nondim]
1047	626926	if (abs(x) <= 2.0e-5) then
1048		! The Taylor series expression for exp(-x) gives a more accurate expression for 64-bit reals.
1049	0	one_minus_exp_x = 1.0 - x * (0.5 - C1_6*x)
1050		else
1051	626926	one_minus_exp_x = (1.0 - exp(-x)) / x
1052		endif
1053	626926	end function one_minus_exp_x
1054
1055		!> Return the value of (1 - exp(-x)) [nondim], using an accurate expression for small values of x.
1056	0	real function one_minus_exp(x)
1057		real, intent(in) :: x !< The argument of the function ((1 - exp(-x))/x) [nondim]
1058		real, parameter :: C1_6 = 1.0/6.0 ! A rational fraction [nondim]
1059	0	if (abs(x) <= 2.0e-5) then
1060		! The Taylor series expression for exp(-x) gives a more accurate expression for 64-bit reals.
1061	0	one_minus_exp = x * (1.0 - x * (0.5 - C1_6*x))
1062		else
1063	0	one_minus_exp = 1.0 - exp(-x)
1064		endif
1065	0	end function one_minus_exp
1066
1067		!> A subroutine to fill the Stokes drift from a NetCDF file
1068		!! using the data_override procedures.
1069	0	subroutine Surface_Bands_by_data_override(Time, G, GV, US, CS)
1070		type(time_type), intent(in) :: Time !< Time to get Stokes drift bands
1071		type(wave_parameters_CS), pointer :: CS !< Wave structure
1072		type(ocean_grid_type), intent(inout) :: G !< Grid structure
1073		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1074		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1075
1076		! Local variables
1077	0	real :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal Stokes drift of band at h-points [L T-1 ~> m s-1]
1078	0	real :: temp_y(SZI_(G),SZJ_(G)) ! Pseudo-meridional Stokes drift of band at h-points [L T-1 ~> m s-1]
1079		integer, dimension(4) :: sizes ! The sizes of the various dimensions of the variable.
1080		character(len=48) :: dim_name(4) ! The names of the dimensions of the variable.
1081		character(len=20) :: varname ! The name of an input variable for data override.
1082		real :: PI ! 3.1415926535... [nondim]
1083		real :: avgx, avgy ! The global averages of temp_x and temp_y [L T-1 ~> m s-1]
1084		logical :: wavenumber_exists
1085		integer :: ndims, b, i, j
1086
1087	0	if (.not.CS%DataOver_initialized) then
1088	0	call data_override_init(G%Domain)
1089	0	CS%DataOver_initialized = .true.
1090
1091	0	if (.not.file_exists(CS%SurfBandFileName)) &
1092	0	call MOM_error(FATAL, "MOM_wave_interface is unable to find file "//trim(CS%SurfBandFileName))
1093
1094		! Check if input has wavenumber or frequency variables.
1095
1096		! Read the number of wavenumber bands in the file, if the variable 'wavenumber' exists.
1097	0	call get_var_sizes(CS%SurfBandFileName, 'wavenumber', ndims, sizes, dim_names=dim_name)
1098	0	wavenumber_exists = (ndims > -1)
1099
1100	0	if (.not.wavenumber_exists) then
1101		! Read the number of frequency bands in the file, if the variable 'frequency' exists.
1102	0	call get_var_sizes(CS%SurfBandFileName, 'frequency', ndims, sizes, dim_names=dim_name)
1103	0	if (ndims < 0) &
1104		call MOM_error(FATAL, "error finding variable 'wavenumber' or 'frequency' in file "//&
1105	0	trim(CS%SurfBandFileName)//" in MOM_wave_interface.")
1106		endif
1107
1108	0	CS%NUMBANDS = sizes(1)
1109		! Allocate the wavenumber bins
1110	0	allocate( CS%WaveNum_Cen(CS%NUMBANDS), source=0.0 )
1111
1112	0	if (wavenumber_exists) then
1113		! Wavenumbers found, so this file uses the old method:
1114	0	CS%PartitionMode = 0
1115
1116		! Reading wavenumber bins
1117	0	call read_variable(CS%SurfBandFileName, dim_name(1), CS%WaveNum_Cen, scale=US%Z_to_m)
1118
1119		else
1120		! Frequencies found, so this file uses the newer method:
1121	0	CS%PartitionMode = 1
1122
1123		! Allocate the frequency bins
1124	0	allocate( CS%Freq_Cen(CS%NUMBANDS), source=0.0 )
1125
1126		! Reading frequencies
1127	0	PI = 4.0*atan(1.0)
1128	0	call read_variable(CS%SurfBandFileName, dim_name(1), CS%Freq_Cen, scale=2.PIUS%T_to_s)
1129
1130	0	do b = 1,CS%NumBands
1131	0	CS%WaveNum_Cen(b) = CS%Freq_Cen(b)**2 / CS%g_Earth
1132		enddo
1133		endif
1134
1135	0	if (.not.allocated(CS%STKx0)) then
1136	0	allocate( CS%STKx0(G%isdB:G%iedB,G%jsd:G%jed,CS%NUMBANDS), source=0.0 )
1137		endif
1138	0	if (.not.allocated(CS%STKy0)) then
1139	0	allocate( CS%STKy0(G%isd:G%ied,G%jsdB:G%jedB,CS%NUMBANDS), source=0.0 )
1140		endif
1141		endif
1142
1143	0	do b = 1,CS%NumBands
1144	0	temp_x(:,:) = 0.0
1145	0	temp_y(:,:) = 0.0
1146	0	varname = ' '
1147	0	write(varname, "(A3,I0)") 'Usx', b
1148	0	call data_override(G%Domain, trim(varname), temp_x, Time, scale=US%m_s_to_L_T)
1149	0	varname = ' '
1150	0	write(varname, "(A3,I0)") 'Usy', b
1151	0	call data_override(G%Domain, trim(varname), temp_y, Time, scale=US%m_s_to_L_T)
1152		! Update halo on h-grid
1153	0	call pass_vector(temp_x, temp_y, G%Domain, To_All, AGRID)
1154		! Filter land values
1155	0	do j = G%jsd,G%jed
1156	0	do i = G%Isd,G%Ied
1157	0	if ((abs(temp_x(i,j)) > CS%land_speed) .or. (abs(temp_y(i,j)) > CS%land_speed)) then
1158		! Assume land-mask and zero out
1159	0	temp_x(i,j) = 0.0
1160	0	temp_y(i,j) = 0.0
1161		endif
1162		enddo
1163		enddo
1164	0	if (CS%Homogenize_Surfbands) then
1165	0	avgx = global_area_mean(temp_x, G)
1166	0	avgy = global_area_mean(temp_y, G)
1167	0	do j = G%jsd,G%jed ; do i = G%Isd,G%Ied ; if (G%mask2dT(i,j) > 0.0) then
1168	0	temp_y(i,j) = avgy
1169	0	temp_x(i,j) = avgx
1170		endif ; enddo ; enddo
1171		endif
1172
1173		! Interpolate to u/v grids
1174	0	do j = G%jsc,G%jec
1175	0	do I = G%IscB,G%IecB
1176	0	CS%STKx0(I,j,b) = 0.5 * (temp_x(i,j) + temp_x(i+1,j))
1177		enddo
1178		enddo
1179	0	do J = G%JscB,G%JecB
1180	0	do i = G%isc,G%iec
1181	0	CS%STKy0(i,J,b) = 0.5 * (temp_y(i,j) + temp_y(i,j+1))
1182		enddo
1183		enddo
1184		enddo !Closes b-loop
1185
1186		! Update halo on u/v grids
1187	0	call pass_vector(CS%STKx0(:,:,:), CS%STKy0(:,:,:), G%Domain, To_ALL)
1188
1189	0	end subroutine Surface_Bands_by_data_override
1190
1191		!> Interface to get Langmuir number based on options stored in wave structure
1192		!!
1193		!! Note this can be called with an unallocated Waves pointer, which is okay if we
1194		!! want the wind-speed only dependent Langmuir number. Therefore, we need to be
1195		!! careful about what we try to access here.
1196	626926	subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, dz, Waves, &
1197	626926	U_H, V_H, Override_MA )
1198		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1199		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1200		real, intent(out) :: LA !< Langmuir number [nondim]
1201		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1202		real, intent(in) :: HBL !< (Positive) thickness of boundary layer [Z ~> m]
1203		real, intent(in) :: ustar !< Friction velocity [Z T-1 ~> m s-1]
1204		integer, intent(in) :: i !< Meridional index of h-point
1205		integer, intent(in) :: j !< Zonal index of h-point
1206		real, dimension(SZK_(GV)), intent(in) :: dz !< Grid layer thickness [Z ~> m]
1207		type(Wave_parameters_CS), pointer :: Waves !< Surface wave control structure.
1208		real, dimension(SZK_(GV)), &
1209		optional, intent(in) :: U_H !< Zonal velocity at H point [L T-1 ~> m s-1] or [m s-1]
1210		real, dimension(SZK_(GV)), &
1211		optional, intent(in) :: V_H !< Meridional velocity at H point [L T-1 ~> m s-1] or [m s-1]
1212		logical, optional, intent(in) :: Override_MA !< Override to use misalignment in LA
1213		!! calculation. This can be used if diagnostic
1214		!! LA outputs are desired that are different than
1215		!! those used by the dynamical model.
1216
1217
1218		!Local Variables
1219		real :: Top, Bottom, MidPoint ! Positions within each layer [Z ~> m]
1220		real :: Dpt_LASL ! Averaging depth for Stokes drift [Z ~> m]
1221		real :: ShearDirection ! Shear angular direction from atan2 [radians]
1222		real :: WaveDirection ! Wave angular direction from atan2 [radians]
1223		real :: LA_STKx, LA_STKy, LA_STK ! Stokes velocities in [L T-1 ~> m s-1]
1224		logical :: ContinueLoop, USE_MA
1225	626926	real, dimension(SZK_(GV)) :: US_H, VS_H ! Profiles of Stokes velocities [L T-1 ~> m s-1]
1226	313463	real, allocatable :: StkBand_X(:), StkBand_Y(:) ! Stokes drifts by band [L T-1 ~> m s-1]
1227		integer :: k, BB
1228
1229		! Compute averaging depth for Stokes drift (negative)
1230	313463	Dpt_LASL = -1.0max(Waves%LA_FracHBLHBL, Waves%LA_HBL_min)
1231
1232	313463	USE_MA = Waves%LA_Misalignment
1233	313463	if (present(Override_MA)) USE_MA = Override_MA
1234
1235		! If requesting to use misalignment in the Langmuir number compute the Shear Direction
1236	313463	if (USE_MA) then
1237	0	if (.not.(present(U_H).and.present(V_H))) call MOM_error(FATAL, &
1238	0	"Get_LA_waves requested to consider misalignment, but velocities were not provided.")
1239	0	ContinueLoop = .true.
1240	0	bottom = 0.0
1241	0	do k = 1,GV%ke
1242	0	Top = Bottom
1243	0	MidPoint = Bottom + 0.5*dz(k)
1244	0	Bottom = Bottom + dz(k)
1245
1246	0	if (Waves%LA_Misalign_bug) then
1247		! Given the sign convention that Dpt_LASL is negative, the next line has a bug.
1248	0	if (MidPoint > Dpt_LASL .and. k > 1 .and. ContinueLoop) then
1249	0	ShearDirection = atan2(V_H(1)-V_H(k), U_H(1)-U_H(k))
1250	0	ContinueLoop = .false.
1251		endif
1252		else ! This version avoids the bug in the version above.
1253	0	if (MidPoint > abs(Dpt_LASL) .and. (k > 1) .and. ContinueLoop) then
1254	0	ShearDirection = atan2(V_H(1)-V_H(k), U_H(1)-U_H(k))
1255	0	ContinueLoop = .false.
1256		endif
1257		endif
1258		enddo
1259		endif
1260
1261	313463	if (Waves%WaveMethod==TESTPROF) then
1262	0	do k = 1,GV%ke
1263	0	US_H(k) = 0.5*(Waves%US_X(I,j,k)+Waves%US_X(I-1,j,k))
1264	0	VS_H(k) = 0.5*(Waves%US_Y(i,J,k)+Waves%US_Y(i,J-1,k))
1265		enddo
1266	0	call Get_SL_Average_Prof( GV, Dpt_LASL, dz, US_H, LA_STKx)
1267	0	call Get_SL_Average_Prof( GV, Dpt_LASL, dz, VS_H, LA_STKy)
1268	0	LA_STK = sqrt((LA_STKXLA_STKX) + (LA_STKYLA_STKY))
1269	313463	elseif (Waves%WaveMethod==SURFBANDS) then
1270	0	allocate(StkBand_X(Waves%NumBands), StkBand_Y(Waves%NumBands))
1271	0	do bb = 1,Waves%NumBands
1272	0	StkBand_X(bb) = 0.5*(Waves%STKx0(I,j,bb)+Waves%STKx0(I-1,j,bb))
1273	0	StkBand_Y(bb) = 0.5*(Waves%STKy0(i,J,bb)+Waves%STKy0(i,J-1,bb))
1274		enddo
1275	0	call Get_SL_Average_Band(GV, Dpt_LASL, Waves%NumBands, Waves%WaveNum_Cen, StkBand_X, LA_STKx )
1276	0	call Get_SL_Average_Band(GV, Dpt_LASL, Waves%NumBands, Waves%WaveNum_Cen, StkBand_Y, LA_STKy )
1277	0	LA_STK = sqrt((LA_STKX2) + (LA_STKY2))
1278	0	deallocate(StkBand_X, StkBand_Y)
1279	313463	elseif (Waves%WaveMethod==DHH85) then
1280		! Temporarily integrating profile rather than spectrum for simplicity
1281	0	do k = 1,GV%ke
1282	0	US_H(k) = 0.5*(Waves%US_X(I,j,k)+Waves%US_X(I-1,j,k))
1283	0	VS_H(k) = 0.5*(Waves%US_Y(i,J,k)+Waves%US_Y(i,J-1,k))
1284		enddo
1285	0	call Get_SL_Average_Prof( GV, Dpt_LASL, dz, US_H, LA_STKx)
1286	0	call Get_SL_Average_Prof( GV, Dpt_LASL, dz, VS_H, LA_STKy)
1287	0	LA_STK = sqrt((LA_STKX2) + (LA_STKY2))
1288	313463	elseif (Waves%WaveMethod==LF17) then
1289	313463	call get_StokesSL_LiFoxKemper(ustar, HBL*Waves%LA_FracHBL, GV, US, Waves, LA_STK, LA)
1290	0	elseif (Waves%WaveMethod==Null_WaveMethod) then
1291		call MOM_error(FATAL, "Get_Langmuir_number called without defining a WaveMethod. "//&
1292		"Suggest to make sure USE_LT is set/overridden to False or choose "//&
1293	0	"a wave method (or set USE_LA_LI2016 to use statistical waves).")
1294		endif
1295
1296	313463	if (.not.(Waves%WaveMethod==LF17)) then
1297		! This expression uses an arbitrary lower bound on Langmuir number.
1298		! We shouldn't expect values lower than this, but there is also no good reason to cap it here
1299		! other than to prevent large enhancements in unconstrained parts of the curve fit parameterizations.
1300	0	LA = max(Waves%La_min, sqrt(US%Z_to_L*ustar / (LA_STK + Waves%La_Stk_backgnd)))
1301		endif
1302
1303	313463	if (Use_MA) then
1304	0	WaveDirection = atan2(LA_STKy, LA_STKx)
1305	0	LA = LA / sqrt(max(1.e-8, cos( WaveDirection - ShearDirection)))
1306		endif
1307
1308	940389	end subroutine get_Langmuir_Number
1309
1310		!> function to return the wave method string set in the param file
1311	0	function get_wave_method(CS)
1312		character(:), allocatable :: get_wave_method
1313		type(wave_parameters_CS), pointer :: CS !< Control structure
1314
1315	0	if (associated(CS)) then
1316	0	select case(CS%WaveMethod)
1317		case (Null_WaveMethod)
1318	0	get_wave_method = NULL_STRING
1319		case (TESTPROF)
1320	0	get_wave_method = TESTPROF_STRING
1321		case (SURFBANDS)
1322	0	get_wave_method = SURFBANDS_STRING
1323		case (DHH85)
1324	0	get_wave_method = DHH85_STRING
1325		case (LF17)
1326	0	get_wave_method = LF17_STRING
1327		case (EFACTOR)
1328	0	get_wave_method = EFACTOR_STRING
1329		end select
1330		else
1331	0	get_wave_method = NULL_STRING
1332		endif
1333	0	end function get_wave_method
1334
1335		!> Get SL averaged Stokes drift from Li/FK 17 method
1336		!!
1337		!! Original description:
1338		!! - This function returns the enhancement factor, given the 10-meter
1339		!! wind [m s-1], friction velocity [m s-1] and the boundary layer depth [m].
1340		!!
1341		!! Update (Jan/25):
1342		!! - Converted from function to subroutine, now returns Langmuir number.
1343		!! - Compute 10m wind internally, so only ustar and hbl need passed to
1344		!! subroutine.
1345		!!
1346		!! Qing Li, 160606
1347		!! - BGR port from CVMix to MOM6 Jan/25/2017
1348		!! - BGR change output to LA from Efactor
1349		!! - BGR remove u10 input
1350		!! - BGR note: fixed parameter values should be changed to "get_params"
1351	313463	subroutine get_StokesSL_LiFoxKemper(ustar, hbl, GV, US, CS, UStokes_SL, LA)
1352		real, intent(in) :: ustar !< water-side surface friction velocity [Z T-1 ~> m s-1].
1353		real, intent(in) :: hbl !< boundary layer depth [Z ~> m].
1354		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1355		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1356		type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure
1357		real, intent(out) :: UStokes_SL !< Surface layer averaged Stokes drift [L T-1 ~> m s-1]
1358		real, intent(out) :: LA !< Langmuir number [nondim]
1359		! Local variables
1360		! parameters
1361		real, parameter :: u19p5_to_u10 = 1.075 ! ratio of U19.5 to U10 (Holthuijsen, 2007) [nondim]
1362		real, parameter :: fm_into_fp = 1.296 ! ratio of mean frequency to peak frequency for
1363		! Pierson-Moskowitz spectrum (Webb, 2011) [nondim]
1364		real, parameter :: us_to_u10 = 0.0162 ! ratio of surface Stokes drift to U10 [nondim]
1365		real, parameter :: r_loss = 0.667 ! loss ratio of Stokes transport [nondim]
1366		real :: UStokes ! The surface Stokes drift [L T-1 ~> m s-1]
1367		real :: hm0 ! The significant wave height [Z ~> m]
1368		real :: fm ! The mean wave frequency [T-1 ~> s-1]
1369		real :: fp ! The peak wave frequency [T-1 ~> s-1]
1370		real :: kphil ! A peak wavenumber in the Phillips spectrum [Z-1 ~> m-1]
1371		real :: kstar ! A rescaled wavenumber? [Z-1 ~> m-1]
1372		real :: vstokes ! The total Stokes transport [Z L T-1 ~> m2 s-1]
1373		real :: z0 ! The boundary layer depth [Z ~> m]
1374		real :: z0i ! The inverse of the boundary layer depth [Z-1 ~> m-1]
1375		real :: r1, r2, r3, r4 ! Nondimensional ratios [nondim]
1376		real :: r5 ! A single expression that combines r2 and r4 [nondim]
1377		real :: root_2kz ! The square root of twice the peak wavenumber times the
1378		! boundary layer depth [nondim]
1379		real :: u10 ! The 10 m wind speed [L T-1 ~> m s-1]
1380		real :: PI ! 3.1415926535... [nondim]
1381
1382	313463	PI = 4.0*atan(1.0)
1383	313463	UStokes_sl = 0.0
1384	313463	LA = 1.e8
1385	313463	if (ustar > 0.0) then
1386		! This code should be revised to minimize the number of divisions and cancel out common factors.
1387
1388		! Computing u10 based on u_star and COARE 3.5 relationships
1389	313463	call ust_2_u10_coare3p5(ustar*sqrt(CS%rho_ocn/CS%rho_air), u10, GV, US, CS)
1390		! surface Stokes drift
1391	313463	UStokes = us_to_u10*u10
1392		!
1393		! significant wave height from Pierson-Moskowitz spectrum (Bouws, 1998)
1394	313463	hm0 = CS%SWH_from_u10sq * u10**2
1395		!
1396		! peak frequency (PM, Bouws, 1998)
1397	313463	fp = 0.877 * (US%L_to_ZGV%g_Earth) / (2.0 PI * u19p5_to_u10 * u10)
1398		!
1399		! mean frequency
1400	313463	fm = fm_into_fp * fp
1401		!
1402		! total Stokes transport (a factor r_loss is applied to account
1403		! for the effect of directional spreading, multidirectional waves
1404		! and the use of PM peak frequency and PM significant wave height
1405		! on estimating the Stokes transport)
1406	313463	vstokes = 0.125 * PI * r_loss * US%Z_to_L * fm * hm0**2
1407		!
1408		! the general peak wavenumber for Phillips' spectrum
1409		! (Breivik et al., 2016) with correction of directional spreading
1410	313463	kphil = 0.176 * UStokes / vstokes
1411
1412		! Combining all of the expressions above gives kPhil as the following
1413		! where the first two lines are just a constant:
1414		! kphil = ((0.176 * us_to_u10 * u19p5_to_u10) / &
1415		! (0.50.125 r_loss * fm_into_fp * 0.877 * CS%SWH_from_u10sq**2)) / &
1416		! (GV%g_Earth * u10**2)
1417
1418		! surface layer
1419	313463	z0 = abs(hbl)
1420
1421	313463	if (CS%answer_date < 20230102) then
1422	0	z0i = 1.0 / z0
1423
1424		! Surface layer averaged Stokes drift with Stokes drift profile
1425		! estimated from Phillips' spectrum (Breivik et al., 2016)
1426		! The directional spreading effect from Webb and Fox-Kemper, 2015 is also included.
1427	0	kstar = kphil * 2.56
1428
1429		! Terms 1 to 4, as written in the appendix of Li et al. (2017)
1430		r1 = ( 0.151 / kphil * z0i - 0.84 ) * &
1431	0	( 1.0 - exp(-2.0 * kphil * z0) )
1432		r2 = -( 0.84 + 0.0591 / kphil * z0i ) * &
1433		sqrt( 2.0 * PI * kphil * z0 ) * &
1434	0	erfc( sqrt( 2.0 * kphil * z0 ) )
1435		r3 = ( 0.0632 / kstar * z0i + 0.125 ) * &
1436	0	(1.0 - exp(-2.0 * kstar * z0) )
1437		r4 = ( 0.125 + 0.0946 / kstar * z0i ) * &
1438		sqrt( 2.0 * PI * kstar * z0) * &
1439	0	erfc( sqrt( 2.0 * kstar * z0 ) )
1440	0	UStokes_sl = UStokes * (0.715 + r1 + r2 + r3 + r4)
1441		else
1442		! The following is equivalent to the code above, but avoids singularities
1443	313463	r1 = ( 0.302 - 1.68(kphilz0) ) * one_minus_exp_x(2.0 * (kphil * z0))
1444	313463	r3 = ( 0.1264 + 0.64(kphilz0) ) * one_minus_exp_x(5.12 * (kphil * z0))
1445
1446	313463	root_2kz = sqrt(2.0 * kphil * z0)
1447		! r2 = -( 0.84 + 0.05912.0 / (root_2kz2) ) sqrt(PI) * root_2kz * erfc( root_2kz )
1448		! r4 = ( 0.2 + 0.0591252.0 / (root_2kz2) ) sqrt(PI) * root_2kz * erfc( 1.6 * root_2kz )
1449
1450		! r5 = r2 + r4 (with a small correction to one coefficient to avoid a singularity when z0 = 0):
1451		! The correction leads to <1% relative differences in (r2+r4) for root_2kz > 0.05, but without
1452		! it the values of r2 + r4 are qualitatively wrong (>50% errors) for root_2kz < 0.0015 .
1453		! It has been verified that these two expressions for r5 are the same to 6 decimal places for
1454		! root_2kz between 1e-10 and 1e-3, but that the first one degrades for smaller values.
1455	313463	if (root_2kz > 1e-3) then
1456		r5 = sqrt(PI) * (root_2kz * (-0.84 * erfc(root_2kz) + 0.2 * erfc(1.6*root_2kz)) + &
1457	313463	0.1182 * (erfc(1.6*root_2kz) - erfc(root_2kz)) / root_2kz)
1458		else
1459		! It is more accurate to replace erf with the first two terms of its Taylor series
1460		! erf(z) = (2/sqrt(pi)) * z * (1. - (1/3)z2 + (1/10)z*4 - (1/42)z**6 + ...)
1461		! and then cancel or combine common terms and drop negligibly small terms.
1462	0	r5 = -0.64sqrt(PI)root_2kz + (-0.14184 + 1.0839648 * root_2kz**2)
1463		endif
1464	313463	UStokes_sl = UStokes * (0.715 + ((r1 + r3) + r5))
1465		endif
1466
1467	313463	if (UStokes_sl /= 0.0) LA = sqrt(US%Z_to_L*ustar / UStokes_sl)
1468		endif
1469
1470	313463	end subroutine Get_StokesSL_LiFoxKemper
1471
1472		!> Get SL Averaged Stokes drift from a Stokes drift Profile
1473	0	subroutine Get_SL_Average_Prof( GV, AvgDepth, dz, Profile, Average )
1474		type(verticalGrid_type), &
1475		intent(in) :: GV !< Ocean vertical grid structure
1476		real, intent(in) :: AvgDepth !< Depth to average over (negative) [Z ~> m]
1477		real, dimension(SZK_(GV)), &
1478		intent(in) :: dz !< Grid thickness [Z ~> m]
1479		real, dimension(SZK_(GV)), &
1480		intent(in) :: Profile !< Profile of quantity to be averaged in arbitrary units [A]
1481		!! (used here for Stokes drift)
1482		real, intent(out) :: Average !< Output quantity averaged over depth AvgDepth [A]
1483		!! (used here for Stokes drift)
1484		!Local variables
1485		real :: Top, Bottom ! Depths, negative downward [Z ~> m]
1486		real :: Sum ! The depth weighted vertical sum of a quantity [A Z ~> A m]
1487		integer :: k
1488
1489		! Initializing sum
1490	0	Sum = 0.0
1491
1492		! Integrate
1493	0	bottom = 0.0
1494	0	do k = 1, GV%ke
1495	0	Top = Bottom
1496	0	Bottom = Bottom - dz(k)
1497	0	if (AvgDepth < Bottom) then ! The whole cell is within H_LA
1498	0	Sum = Sum + Profile(k) * dz(k)
1499	0	elseif (AvgDepth < Top) then ! A partial cell is within H_LA
1500	0	Sum = Sum + Profile(k) * (Top-AvgDepth)
1501	0	exit
1502		else
1503	0	exit
1504		endif
1505		enddo
1506
1507		! Divide by AvgDepth or the depth in the column, whichever is smaller.
1508	0	if (abs(AvgDepth) <= abs(Bottom)) then
1509	0	Average = Sum / abs(AvgDepth)
1510	0	elseif (abs(Bottom) > 0.0) then
1511	0	Average = Sum / abs(Bottom)
1512		else
1513	0	Average = 0.0
1514		endif
1515
1516	0	end subroutine Get_SL_Average_Prof
1517
1518		!> Get SL averaged Stokes drift from the banded Spectrum method
1519	0	subroutine Get_SL_Average_Band( GV, AvgDepth, NB, WaveNumbers, SurfStokes, Average )
1520		type(verticalGrid_type), &
1521		intent(in) :: GV !< Ocean vertical grid
1522		real, intent(in) :: AvgDepth !< Depth to average over [Z ~> m].
1523		integer, intent(in) :: NB !< Number of bands used
1524		real, dimension(NB), &
1525		intent(in) :: WaveNumbers !< Wavenumber corresponding to each band [Z-1 ~> m-1]
1526		real, dimension(NB), &
1527		intent(in) :: SurfStokes !< Surface Stokes drift for each band [L T-1 ~> m s-1]
1528		real, intent(out) :: Average !< Output average Stokes drift over depth AvgDepth [L T-1 ~> m s-1]
1529
1530		! Local variables
1531		integer :: bb
1532
1533		! Loop over bands
1534	0	Average = 0.0
1535	0	do bb = 1, NB
1536		! Factor includes analytical integration of e(2kz)
1537		! - divided by (-H_LA) to get average from integral.
1538		Average = Average + SurfStokes(BB) * &
1539		(1.-EXP(-abs(AvgDepth * 2.0 * WaveNumbers(BB)))) / &
1540	0	abs(AvgDepth * 2.0 * WaveNumbers(BB))
1541
1542		! For accuracy when AvgDepth is small change the above to:
1543		! Average = Average + SurfStokes(BB) * one_minus_exp_x(abs(AvgDepth * 2.0 * WaveNumbers(BB)))
1544		enddo
1545
1546	0	end subroutine Get_SL_Average_Band
1547
1548		!> Compute the Stokes drift at a given depth
1549		!!
1550		!! Taken from Qing Li (Brown)
1551		!! use for comparing MOM6 simulation to his LES
1552		!! computed at z mid point (I think) and not depth averaged.
1553		!! Should be fine to integrate in frequency from 0.1 to sqrt(-0.2grav2pi/dz
1554	0	subroutine DHH85_mid(GV, US, CS, zpt, UStokes)
1555		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
1556		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1557		type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure
1558		real, intent(in) :: zpt !< Depth to get Stokes drift [Z ~> m].
1559		real, intent(out) :: UStokes !< Stokes drift [L T-1 ~> m s-1]
1560		!
1561		real :: ann, Bnn, Snn, Cnn, Dnn ! Nondimensional factors [nondim]
1562		real :: omega_peak ! The peak wave frequency [T-1 ~> s-1]
1563		real :: omega ! The average frequency in the band [T-1 ~> s-1]
1564		real :: domega ! The width in frequency of the band [T-1 ~> s-1]
1565		real :: u10 ! The wind speed for this spectrum [Z T-1 ~> m s-1]
1566		real :: wavespec ! The wave spectrum [L Z T ~> m2 s]
1567		real :: Stokes ! The Stokes displacement per cycle [L ~> m]
1568		real :: PI ! 3.1415926535... [nondim]
1569		integer :: Nomega ! The number of wavenumber bands
1570		integer :: OI
1571
1572	0	u10 = CS%WaveWind*US%L_to_Z
1573
1574		!/
1575	0	NOmega = 1000
1576	0	domega = (CS%omega_max - CS%omega_min) / real(NOmega)
1577
1578		!
1579	0	if (CS%WaveAgePeakFreq) then
1580	0	omega_peak = CS%g_Earth / (CS%WaveAge * u10)
1581		else
1582	0	PI = 4.0*atan(1.0)
1583	0	omega_peak = 2. * PI * 0.13 * CS%g_Earth / u10
1584		endif
1585		!/
1586	0	Ann = 0.006 * CS%WaveAge**(-0.55)
1587	0	Bnn = 1.0
1588	0	Snn = 0.08 * (1.0 + 4.0 * CS%WaveAge**3)
1589	0	Cnn = 1.7
1590	0	if (CS%WaveAge < 1.) then
1591	0	Cnn = Cnn - 6.0*log10(CS%WaveAge)
1592		endif
1593		!/
1594	0	UStokes = 0.0
1595	0	omega = CS%omega_min + 0.5*domega
1596	0	do oi = 1,nomega-1
1597	0	Dnn = exp ( -0.5 * (omega-omega_peak)2 / (Snn2 * omega_peak**2) )
1598		! wavespec units [L Z T ~> m2 s]
1599		wavespec = US%Z_to_L * (Ann * CS%g_Earth*2 / (omega_peakomega*4 ) ) &
1600	0	exp(-bnn(omega_peak/omega)4)Cnn**Dnn
1601		! Stokes units [L ~> m] (multiply by frequency range for units of [L T-1 ~> m s-1])
1602		Stokes = 2.0 * wavespec * omega*3 &
1603	0	exp( 2.0 * omega*2 zpt / CS%g_Earth) / CS%g_Earth
1604	0	UStokes = UStokes + Stokes*domega
1605	0	omega = omega + domega
1606		enddo
1607
1608	0	end subroutine DHH85_mid
1609
1610		!> Explicit solver for Stokes mixing.
1611		!! Still in development do not use.
1612	0	subroutine StokesMixing(G, GV, dt, h, dz, u, v, Waves )
1613		type(ocean_grid_type), &
1614		intent(in) :: G !< Ocean grid
1615		type(verticalGrid_type), &
1616		intent(in) :: GV !< Ocean vertical grid
1617		real, intent(in) :: dt !< Time step of MOM6 [T ~> s] for explicit solver
1618		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1619		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1620		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1621		intent(in) :: dz !< Vertical distance between interfaces around a layer [Z ~> m]
1622		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1623		intent(inout) :: u !< Velocity i-component [L T-1 ~> m s-1]
1624		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1625		intent(inout) :: v !< Velocity j-component [L T-1 ~> m s-1]
1626		type(Wave_parameters_CS), &
1627		pointer :: Waves !< Surface wave related control structure.
1628		! Local variables
1629		real :: dTauUp, dTauDn ! Vertical momentum fluxes [H L T-2 ~> m2 s-2 or Pa]
1630		real :: h_lay ! The layer thickness at a velocity point [H ~> m or kg m-2]
1631		real :: dz_lay ! The distance between interfaces at a velocity point [Z ~> m]
1632		integer :: i, j, k
1633
1634		! This is a template to think about down-Stokes mixing.
1635		! This is not ready for use...
1636
1637	0	do k = 1, GV%ke
1638	0	do j = G%jsc, G%jec
1639	0	do I = G%iscB, G%iecB
1640	0	h_lay = 0.5*(h(i,j,k)+h(i+1,j,k))
1641	0	dz_lay = 0.5*(dz(i,j,k)+dz(i+1,j,k))
1642	0	dTauUp = 0.0
1643	0	if (k > 1) &
1644		dTauUp = (0.5(waves%Kvs(i,j,k)+waves%Kvs(i+1,j,k))) &
1645		(waves%us_x(i,j,k-1)-waves%us_x(i,j,k)) / &
1646	0	(0.5(dz_lay + 0.5(dz(i,j,k-1)+dz(i+1,j,k-1)) ))
1647	0	dTauDn = 0.0
1648	0	if (k < GV%ke-1) &
1649		dTauDn = (0.5(waves%Kvs(i,j,k+1)+waves%Kvs(i+1,j,k+1))) &
1650		(waves%us_x(i,j,k)-waves%us_x(i,j,k+1)) / &
1651	0	(0.5(dz_lay + 0.5(dz(i,j,k+1)+dz(i+1,j,k+1)) ))
1652	0	u(i,j,k) = u(i,j,k) + dt * (dTauUp-dTauDn) / h_lay
1653		enddo
1654		enddo
1655		enddo
1656
1657	0	do k = 1, GV%ke
1658	0	do J = G%jscB, G%jecB
1659	0	do i = G%isc, G%iec
1660	0	h_lay = 0.5*(h(i,j,k)+h(i,j+1,k))
1661	0	dz_lay = 0.5*(dz(i,j,k)+dz(i,j+1,k))
1662	0	dTauUp = 0.
1663	0	if (k > 1) &
1664		dTauUp = (0.5(waves%Kvs(i,j,k)+waves%Kvs(i,j+1,k))) &
1665		(waves%us_y(i,j,k-1)-waves%us_y(i,j,k)) / &
1666	0	(0.5(dz_lay + 0.5(dz(i,j,k-1)+dz(i,j+1,k-1)) ))
1667	0	dTauDn = 0.0
1668	0	if (k < GV%ke-1) &
1669		dTauDn = (0.5(waves%Kvs(i,j,k+1)+waves%Kvs(i,j+1,k+1))) &
1670		(waves%us_y(i,j,k)-waves%us_y(i,j,k+1)) / &
1671	0	(0.5(dz_lay + 0.5(dz(i,j,k+1)+dz(i,j+1,k+1)) ))
1672	0	v(i,J,k) = v(i,J,k) + dt * (dTauUp-dTauDn) / h_lay
1673		enddo
1674		enddo
1675		enddo
1676
1677	0	end subroutine StokesMixing
1678
1679		!> Solver to add Coriolis-Stokes to model
1680		!! Still in development and not meant for general use.
1681		!! Can be activated (with code intervention) for LES comparison
1682		!! CHECK THAT RIGHT TIMESTEP IS PASSED IF YOU USE THIS**
1683		!!
1684		!! Not accessed in the standard code.
1685	0	subroutine CoriolisStokes(G, GV, dt, h, u, v, Waves)
1686		type(ocean_grid_type), &
1687		intent(in) :: G !< Ocean grid
1688		type(verticalGrid_type), &
1689		intent(in) :: GV !< Ocean vertical grid
1690		real, intent(in) :: dt !< Time step of MOM6 [T ~> s]
1691		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1692		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1693		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1694		intent(inout) :: u !< Velocity i-component [L T-1 ~> m s-1]
1695		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1696		intent(inout) :: v !< Velocity j-component [L T-1 ~> m s-1]
1697		type(Wave_parameters_CS), &
1698		pointer :: Waves !< Surface wave related control structure.
1699
1700		! Local variables
1701		real :: DVel ! A rescaled velocity change [L T-2 ~> m s-2]
1702		integer :: i, j, k
1703
1704	0	do k = 1, GV%ke
1705	0	do j = G%jsc, G%jec
1706	0	do I = G%iscB, G%iecB
1707		DVel = 0.25((Waves%us_y(i,J-1,k)+Waves%us_y(i+1,J-1,k)) G%CoriolisBu(I,J-1)) + &
1708	0	0.25((Waves%us_y(i,J,k)+Waves%us_y(i+1,J,k)) G%CoriolisBu(I,J))
1709	0	u(I,j,k) = u(I,j,k) + DVEL*dt
1710		enddo
1711		enddo
1712		enddo
1713
1714	0	do k = 1, GV%ke
1715	0	do J = G%jscB, G%jecB
1716	0	do i = G%isc, G%iec
1717		DVel = 0.25((Waves%us_x(I-1,j,k)+Waves%us_x(I-1,j+1,k)) G%CoriolisBu(I-1,j)) + &
1718	0	0.25((Waves%us_x(I,j,k)+Waves%us_x(I,j+1,k)) G%CoriolisBu(I,J))
1719	0	v(i,J,k) = v(i,j,k) - DVEL*dt
1720		enddo
1721		enddo
1722		enddo
1723	0	end subroutine CoriolisStokes
1724
1725		!> Computes tendency due to Stokes pressure gradient force anomaly
1726		!! including analytical integration of Stokes shear using multiple-exponential decay
1727		!! Stokes drift profile and vertical integration of the resulting pressure
1728		!! anomaly to the total pressure gradient force
1729	0	subroutine Stokes_PGF(G, GV, US, dz, u, v, PFu_Stokes, PFv_Stokes, CS )
1730		type(ocean_grid_type), &
1731		intent(in) :: G !< Ocean grid
1732		type(verticalGrid_type), &
1733		intent(in) :: GV !< Ocean vertical grid
1734		type(unit_scale_type), &
1735		intent(in) :: US !< A dimensional unit scaling type
1736		real, dimension(SZI_(G),SZJ_(G),SZK_(G)),&
1737		intent(in) :: dz !< Layer thicknesses in height units [Z ~> m]
1738		real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1739		intent(in) :: u !< Lagrangian Velocity i-component [L T-1 ~> m s-1]
1740		real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1741		intent(in) :: v !< Lagrangian Velocity j-component [L T-1 ~> m s-1]
1742		real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1743		intent(out) :: PFu_Stokes !< PGF Stokes-shear i-component [L T-2 ~> m s-2]
1744		real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1745		intent(out) :: PFv_Stokes !< PGF Stokes-shear j-component [L T-2 ~> m s-2]
1746		type(Wave_parameters_CS), &
1747		pointer :: CS !< Surface wave related control structure.
1748
1749		! Local variables
1750	0	real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: P_deltaStokes_L ! The Stokes induced pressure anomaly,
1751		! layer averaged [L2 T-2 ~> m2 s-2]
1752	0	real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: P_deltaStokes_i ! The Stokes induced pressure anomaly
1753		! at interfaces [L2 T-2 ~> m2 s-2]
1754		real :: P_Stokes_l, P_Stokes_r ! Stokes-induced pressure anomaly over layer (left/right of point) [L2 T-2 ~> m2 s-2]
1755		real :: P_Stokes_l0, P_Stokes_r0 ! Stokes-induced pressure anomaly at interface
1756		! (left/right of point) [L2 T-2 ~> m2 s-2]
1757		real :: dP_Stokes_l_dz, dP_Stokes_r_dz ! Contribution of layer to integrated Stokes pressure anomaly for summation
1758		! (left/right of point) [Z L2 T-2 ~> m3 s-2]
1759		real :: dP_lay_Stokes_l, dP_lay_Stokes_r ! Contribution of layer to integrated Stokes pressure anomaly for summation
1760		! (left/right of point) [L2 T-2 ~> m2 s-2]
1761		real :: dP_Stokes_l, dP_Stokes_r ! Net increment of Stokes pressure anomaly across layer for summation
1762		! (left/right of point) [L2 T-2 ~> m2 s-2]
1763		real :: uE_l, uE_r, vE_l, vE_r ! Eulerian velocity components (left/right of point) [L T-1 ~> m s-1]
1764		real :: uS0_l, uS0_r, vS0_l, vS0_r ! Surface Stokes velocity components (left/right of point) [L T-1 ~> m s-1]
1765	0	real :: zi_l(SZK_(G)+1), zi_r(SZK_(G)+1) ! The height of the edges of the cells (left/right of point) [Z ~> m].
1766	0	real :: idz_l(SZK_(G)), idz_r(SZK_(G)) ! The inverse thickness of the cells (left/right of point) [Z-1 ~> m-1]
1767		real :: h_l, h_r ! The thickness of the cell (left/right of point) [Z ~> m].
1768		real :: exp_top ! The decay of the surface stokes drift to the interface atop a layer [nondim]
1769		real :: dexp2kzL, dexp4kzL, dexp2kzR, dexp4kzR ! Analytical evaluation of multi-exponential decay
1770		! contribution to Stokes pressure anomalies [nondim].
1771		real :: TwoK, FourK ! Wavenumbers multiplied by a factor [Z-1 ~> m-1]
1772		real :: iTwoK, iFourK ! Inverses of wavenumbers [Z ~> m]
1773
1774		integer :: i, j, k, l
1775
1776		!---------------------------------------------------------------
1777		! Compute the Stokes contribution to the pressure gradient force
1778		!---------------------------------------------------------------
1779		! Notes on the algorithm/code:
1780		! This code requires computing velocities at bounding h points
1781		! of the u/v points to get the pressure-gradient. In this
1782		! implementation there are several redundant calculations as the
1783		! left/right points are computed at each cell while integrating
1784		! in the vertical, requiring about twice the calculations. The
1785		! velocities at the tracer points could be precomputed and
1786		! stored, but this would require more memory and cycling through
1787		! large 3d arrays while computing the pressures. This could be
1788		! explored as a way to speed up this code.
1789		!---------------------------------------------------------------
1790
1791	0	PFu_Stokes(:,:,:) = 0.0
1792	0	PFv_Stokes(:,:,:) = 0.0
1793	0	if (CS%id_P_deltaStokes_i > 0) P_deltaStokes_i(:,:,:) = 0.0
1794	0	if (CS%id_P_deltaStokes_L > 0) P_deltaStokes_L(:,:,:) = 0.0
1795
1796		! First compute PGFu. The Stokes-induced pressure anomaly diagnostic is stored from this calculation.
1797		! > Seeking PGFx at (I,j), meaning we need to compute pressure at h-points (i,j) and (i+1,j).
1798		! UL(i,j) -> found as average of I-1 & I on j
1799		! UR(i+1,j) -> found as average of I & I+1 on j
1800		! VL(i,j) -> found on i as average of J-1 & J
1801		! VR(i+1,j) -> found on i+1 as average of J-1 & J
1802		!
1803	0	do j = G%jsc, G%jec ; do I = G%iscB, G%iecB
1804	0	if (G%mask2dCu(I,j)>0.5) then
1805	0	P_Stokes_l0 = 0.0
1806	0	P_Stokes_r0 = 0.0
1807		! We don't need to precompute the grid in physical space arrays and could have done this during
1808		! the next loop, but this gives flexibility if the loop directions (integrations) are performed
1809		! upwards instead of downwards (it seems downwards is the better approach).
1810	0	zi_l(1) = 0.0
1811	0	zi_r(1) = 0.0
1812	0	do k = 1, G%ke
1813	0	h_l = dz(i,j,k)
1814	0	h_r = dz(i+1,j,k)
1815	0	zi_l(k+1) = zi_l(k) - h_l
1816	0	zi_r(k+1) = zi_r(k) - h_r
1817	0	if (.not.CS%robust_Stokes_PGF) then
1818		! When the code is properly refactored, the following hard-coded constants are unnecessary.
1819	0	Idz_l(k) = 1./max(0.1*US%m_to_Z, h_l)
1820	0	Idz_r(k) = 1./max(0.1*US%m_to_Z, h_r)
1821		endif
1822		enddo
1823	0	do k = 1,G%ke
1824		! Computing (left/right) Eulerian velocities assuming the velocity passed to this routine is the
1825		! Lagrangian velocity. This requires the wave acceleration terms to be activated together.
1826		uE_l = 0.5((u(I-1,j,k)-CS%Us_x(I-1,j,k))G%mask2dCu(I-1,j) + &
1827	0	(u(I,j,k)-CS%Us_x(I,j,k))*G%mask2dCu(I,j))
1828		uE_r = 0.5((u(I,j,k)-CS%Us_x(I,j,k))G%mask2dCu(I,j) + &
1829	0	(u(I+1,j,k)-CS%Us_x(I+1,j,k))*G%mask2dCu(I+1,j))
1830		vE_l = 0.5((v(i,J-1,k)-CS%Us_y(i,J-1,k))G%mask2dCv(i,J-1) + &
1831	0	(v(i,J,k)-CS%Us_y(i,J,k))*G%mask2dCv(i,J))
1832		vE_r = 0.5((v(i+1,J-1,k)-CS%Us_y(i+1,J-1,k))G%mask2dCv(i+1,J-1) + &
1833	0	(v(i+1,J,k)-CS%Us_y(i+1,J,k))*G%mask2dCv(i+1,J))
1834
1835	0	dP_Stokes_l_dz = 0.0
1836	0	dP_Stokes_r_dz = 0.0
1837	0	dP_Stokes_l = 0.0
1838	0	dP_Stokes_r = 0.0
1839	0	dP_lay_Stokes_l=0.0
1840	0	dP_lay_Stokes_r=0.0
1841
1842	0	do l = 1, CS%numbands
1843
1844		! Computing (left/right) surface Stokes drift velocities at wavenumber band
1845		uS0_l = 0.5(CS%Stkx0(I-1,j,l)G%mask2dCu(I-1,j) + &
1846	0	CS%Stkx0(I,j,l)*G%mask2dCu(I,j))
1847		uS0_r = 0.5(CS%Stkx0(I,j,l)G%mask2dCu(I,j) + &
1848	0	CS%Stkx0(I+1,j,l)*G%mask2dCu(I+1,j))
1849		vS0_l = 0.5(CS%Stky0(i,J-1,l)G%mask2dCv(i,J-1) + &
1850	0	CS%Stky0(i,J,l)*G%mask2dCv(i,J))
1851		vS0_r = 0.5(CS%Stky0(i+1,J-1,l)G%mask2dCv(i+1,J-1) + &
1852	0	CS%Stky0(i+1,J,l)*G%mask2dCv(i+1,J))
1853
1854		! Wavenumber terms that are useful to simplify the pressure calculations
1855	0	TwoK = 2.*CS%WaveNum_Cen(l)
1856	0	FourK = 2.*TwoK
1857	0	if (.not.CS%robust_Stokes_PGF) then
1858	0	iTwoK = 1. / TwoK
1859	0	iFourK = 1. / FourK
1860		endif
1861
1862		! Compute Pressure at interface and integrated over layer on left/right bounding points.
1863		! These are summed over wavenumber bands.
1864	0	if (G%mask2dT(i,j)>0.5) then
1865	0	if (.not.CS%robust_Stokes_PGF) then
1866	0	dexp2kzL = exp(TwoKzi_l(k))-exp(TwoKzi_l(k+1))
1867	0	dexp4kzL = exp(FourKzi_l(k))-exp(FourKzi_l(k+1))
1868		dP_Stokes_l_dz = dP_Stokes_l_dz + &
1869	0	((uE_luS0_l+vE_lvS0_l)iTwoKdexp2kzL + 0.5(uS0_luS0_l+vS0_lvS0_l)iFourK*dexp4kzL)
1870	0	dP_Stokes_l = dP_Stokes_l + (uE_luS0_l+vE_lvS0_l)dexp2kzL + 0.5(uS0_luS0_l+vS0_lvS0_l)*dexp4kzL
1871		else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers.
1872	0	exp_top = exp(TwoK*zi_l(k))
1873		dP_lay_Stokes_l = dP_lay_Stokes_l + &
1874		((((uE_luS0_l)+(vE_lvS0_l)) * exp_top) * one_minus_exp_x(TwoK*dz(i,j,k)) + &
1875	0	(0.5((uS0_l2)+(vS0_l2)) exp_top*2) one_minus_exp_x(FourK*dz(i,j,k)) )
1876		dP_Stokes_l = dP_Stokes_l + &
1877		((((uE_luS0_l)+(vE_lvS0_l)) * exp_top) * one_minus_exp(TwoK*dz(i,j,k)) + &
1878	0	(0.5((uS0_l2)+(vS0_l2)) exp_top*2) one_minus_exp(FourK*dz(i,j,k)) )
1879		endif
1880		endif
1881	0	if (G%mask2dT(i+1,j)>0.5) then
1882	0	if (.not.CS%robust_Stokes_PGF) then
1883	0	dexp2kzR = exp(TwoKzi_r(k))-exp(TwoKzi_r(k+1))
1884	0	dexp4kzR = exp(FourKzi_r(k))-exp(FourKzi_r(k+1))
1885		dP_Stokes_r_dz = dP_Stokes_r_dz + &
1886	0	((uE_ruS0_r+vE_rvS0_r)iTwoKdexp2kzR + 0.5(uS0_ruS0_r+vS0_rvS0_r)iFourK*dexp4kzR)
1887	0	dP_Stokes_r = dP_Stokes_r + (uE_ruS0_r+vE_rvS0_r)dexp2kzR + 0.5(uS0_ruS0_r+vS0_rvS0_r)*dexp4kzR
1888		else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers.
1889	0	exp_top = exp(TwoK*zi_r(k))
1890		dP_lay_Stokes_r = dP_lay_Stokes_r + &
1891		((((uE_ruS0_r)+(vE_rvS0_r)) * exp_top) * one_minus_exp_x(TwoK*dz(i+1,j,k)) + &
1892	0	(0.5((uS0_r2)+(vS0_r2)) exp_top*2) one_minus_exp_x(FourK*dz(i+1,j,k)) )
1893		dP_Stokes_r = dP_Stokes_r + &
1894		((((uE_ruS0_r)+(vE_rvS0_r)) * exp_top) * one_minus_exp(TwoK*dz(i+1,j,k)) + &
1895	0	(0.5((uS0_r2)+(vS0_r2)) exp_top*2) one_minus_exp(FourK*dz(i+1,j,k)) )
1896		endif
1897		endif
1898		enddo
1899
1900		! Summing PF over bands
1901		! > Increment the Layer averaged pressure
1902	0	if (.not.CS%robust_Stokes_PGF) then
1903	0	P_Stokes_l = P_Stokes_l0 + dP_Stokes_l_dz*Idz_l(k)
1904	0	P_Stokes_r = P_Stokes_r0 + dP_Stokes_r_dz*Idz_r(k)
1905		else
1906	0	P_Stokes_l = P_Stokes_l0 + dP_lay_Stokes_l
1907	0	P_Stokes_r = P_Stokes_r0 + dP_lay_Stokes_r
1908		endif
1909
1910		! > Increment the Interface pressure
1911	0	P_Stokes_l0 = P_Stokes_l0 + dP_Stokes_l
1912	0	P_Stokes_r0 = P_Stokes_r0 + dP_Stokes_r
1913
1914		! Pressure force anomaly is finite difference across the cell.
1915	0	PFu_Stokes(I,j,k) = (P_Stokes_r - P_Stokes_l)*G%IdxCu(I,j)
1916
1917		! Choose to output the pressure delta on the h-points from the PFu calculation.
1918	0	if (G%mask2dT(i,j)>0.5 .and. CS%id_P_deltaStokes_L > 0) P_deltaStokes_L(i,j,k) = P_Stokes_l
1919	0	if (G%mask2dT(i,j)>0.5 .and. CS%id_P_deltaStokes_i > 0) P_deltaStokes_i(i,j,k+1) = P_Stokes_l0
1920
1921		enddo
1922		endif
1923		enddo ; enddo
1924
1925		! Next compute PGFv. The Stokes-induced pressure anomaly diagnostic is stored from this calculation.
1926		! > Seeking PGFy at (i,J), meaning we need to compute pressure at h-points (i,j) and (i,j+1).
1927		! UL(i,j) -> found as average of I-1 & I on j
1928		! UR(i,j+1) -> found as average of I-1 & I on j+1
1929		! VL(i,j) -> found on i as average of J-1 & J
1930		! VR(i,j+1) -> found on i as average of J & J+1
1931		!
1932	0	do J = G%jscB, G%jecB ; do i = G%isc, G%iec
1933	0	if (G%mask2dCv(i,J)>0.5) then
1934	0	P_Stokes_l0 = 0.0
1935	0	P_Stokes_r0 = 0.0
1936	0	zi_l(1) = 0.0
1937	0	zi_r(1) = 0.0
1938	0	do k = 1, G%ke
1939	0	h_l = dz(i,j,k)
1940	0	h_r = dz(i,j+1,k)
1941	0	zi_l(k+1) = zi_l(k) - h_l
1942	0	zi_r(k+1) = zi_r(k) - h_r
1943	0	if (.not.CS%robust_Stokes_PGF) then
1944		! When the code is properly refactored, the following hard-coded constants are unnecessary.
1945	0	Idz_l(k) = 1. / max(0.1*US%m_to_Z, h_l)
1946	0	Idz_r(k) = 1. / max(0.1*US%m_to_Z, h_r)
1947		endif
1948		enddo
1949	0	do k = 1,G%ke
1950		! Computing (left/right) Eulerian velocities assuming the velocity passed to this routine is the
1951		! Lagrangian velocity. This requires the wave acceleration terms to be activated together.
1952		uE_l = 0.5((u(I-1,j,k)-CS%Us_x(I-1,j,k))G%mask2dCu(I-1,j) + &
1953	0	(u(I,j,k)-CS%Us_x(I,j,k))*G%mask2dCu(I,j))
1954		uE_r = 0.5((u(I-1,j+1,k)-CS%Us_x(I-1,j+1,k))G%mask2dCu(I-1,j+1) + &
1955	0	(u(I,j+1,k)-CS%Us_x(I,j+1,k))*G%mask2dCu(I,j+1))
1956		vE_l = 0.5((v(i,J-1,k)-CS%Us_y(i,J-1,k))G%mask2dCv(i,J-1) + &
1957	0	(v(i,J,k)-CS%Us_y(i,J,k))*G%mask2dCv(i,J))
1958		vE_r = 0.5((v(i,J,k)-CS%Us_y(i,J,k))G%mask2dCv(i,J) + &
1959	0	(v(i,J+1,k)-CS%Us_y(i,J+1,k))*G%mask2dCv(i,J+1))
1960
1961	0	dP_Stokes_l_dz = 0.0
1962	0	dP_Stokes_r_dz = 0.0
1963	0	dP_Stokes_l = 0.0
1964	0	dP_Stokes_r = 0.0
1965	0	dP_lay_Stokes_l=0.0
1966	0	dP_lay_Stokes_r=0.0
1967
1968	0	do l = 1, CS%numbands
1969
1970		! Computing (left/right) surface Stokes drift velocities at wavenumber band
1971		uS0_l = 0.5(CS%Stkx0(I-1,j,l)G%mask2dCu(I-1,j) + &
1972	0	CS%Stkx0(I,j,l)*G%mask2dCu(I,j))
1973		uS0_r = 0.5(CS%Stkx0(I-1,j+1,l)G%mask2dCu(I-1,j+1) + &
1974	0	CS%Stkx0(I,j+1,l)*G%mask2dCu(I,j+1))
1975		vS0_l = 0.5(CS%Stky0(i,J-1,l)G%mask2dCv(i,J-1) + &
1976	0	CS%Stky0(i,J,l)*G%mask2dCv(i,J))
1977		vS0_r = 0.5(CS%Stky0(i,J,l)G%mask2dCv(i,J) + &
1978	0	CS%Stky0(i,J+1,l)*G%mask2dCv(i,J+1))
1979
1980		! Wavenumber terms that are useful to simplify the pressure calculations
1981	0	TwoK = 2.*CS%WaveNum_Cen(l)
1982	0	FourK = 2.*TwoK
1983	0	if (.not.CS%robust_Stokes_PGF) then
1984	0	iTwoK = 1. / TwoK
1985	0	iFourK = 1. / FourK
1986		endif
1987
1988		! Compute Pressure at interface and integrated over layer on left/right bounding points.
1989		! These are summed over wavenumber bands.
1990	0	if (G%mask2dT(i,j)>0.5) then
1991	0	if (.not.CS%robust_Stokes_PGF) then
1992	0	dexp2kzL = exp(TwoKzi_l(k))-exp(TwoKzi_l(k+1))
1993	0	dexp4kzL = exp(FourKzi_l(k))-exp(FourKzi_l(k+1))
1994		dP_Stokes_l_dz = dP_Stokes_l_dz + &
1995	0	((uE_luS0_l+vE_lvS0_l)iTwoKdexp2kzL + 0.5(uS0_luS0_l+vS0_lvS0_l)iFourK*dexp4kzL)
1996	0	dP_Stokes_l = dP_Stokes_l + (uE_luS0_l+vE_lvS0_l)dexp2kzL + 0.5(uS0_luS0_l+vS0_lvS0_l)*dexp4kzL
1997		else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers.
1998	0	exp_top = exp(TwoK*zi_l(k))
1999		dP_lay_Stokes_l = dP_lay_Stokes_l + &
2000		((((uE_luS0_l)+(vE_lvS0_l)) * exp_top) * one_minus_exp_x(TwoK*dz(i,j,k)) + &
2001	0	(0.5((uS0_l2)+(vS0_l2)) exp_top*2) one_minus_exp_x(FourK*dz(i,j,k)) )
2002		dP_Stokes_l = dP_Stokes_l + &
2003		((((uE_luS0_l)+(vE_lvS0_l)) * exp_top) * one_minus_exp(TwoK*dz(i,j,k)) + &
2004	0	(0.5((uS0_l2)+(vS0_l2)) exp_top*2) one_minus_exp(FourK*dz(i,j,k)) )
2005		endif
2006		endif
2007	0	if (G%mask2dT(i,j+1)>0.5) then
2008	0	if (.not.CS%robust_Stokes_PGF) then
2009	0	dexp2kzR = exp(TwoKzi_r(k))-exp(TwoKzi_r(k+1))
2010	0	dexp4kzR = exp(FourKzi_r(k))-exp(FourKzi_r(k+1))
2011		dP_Stokes_r_dz = dP_Stokes_r_dz + &
2012	0	((uE_ruS0_r+vE_rvS0_r)iTwoKdexp2kzR + 0.5(uS0_ruS0_r+vS0_rvS0_r)iFourK*dexp4kzR)
2013	0	dP_Stokes_r = dP_Stokes_r + (uE_ruS0_r+vE_rvS0_r)dexp2kzR + 0.5(uS0_ruS0_r+vS0_rvS0_r)*dexp4kzR
2014		else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers.
2015	0	exp_top = exp(TwoK*zi_r(k))
2016		dP_lay_Stokes_r = dP_lay_Stokes_r + &
2017		((((uE_ruS0_r)+(vE_rvS0_r)) * exp_top) * one_minus_exp_x(TwoK*dz(i,j+1,k)) + &
2018	0	(0.5((uS0_r2)+(vS0_r2)) exp_top*2) one_minus_exp_x(FourK*dz(i,j+1,k)) )
2019		dP_Stokes_r = dP_Stokes_r + &
2020		((((uE_ruS0_r)+(vE_rvS0_r)) * exp_top) * one_minus_exp(TwoK*dz(i,j+1,k)) + &
2021	0	(0.5((uS0_r2)+(vS0_r2)) exp_top*2) one_minus_exp(FourK*dz(i,j+1,k)) )
2022		endif
2023		endif
2024		enddo
2025
2026		! Summing PF over bands
2027		! > Increment the Layer averaged pressure
2028	0	if (.not.CS%robust_Stokes_PGF) then
2029	0	P_Stokes_l = P_Stokes_l0 + dP_Stokes_l_dz*Idz_l(k)
2030	0	P_Stokes_r = P_Stokes_r0 + dP_Stokes_r_dz*Idz_r(k)
2031		else
2032	0	P_Stokes_l = P_Stokes_l0 + dP_lay_Stokes_l
2033	0	P_Stokes_r = P_Stokes_r0 + dP_lay_Stokes_r
2034		endif
2035
2036		! > Increment the Interface pressure
2037	0	P_Stokes_l0 = P_Stokes_l0 + dP_Stokes_l
2038	0	P_Stokes_r0 = P_Stokes_r0 + dP_Stokes_r
2039
2040		! Pressure force anomaly is finite difference across the cell.
2041	0	PFv_Stokes(I,j,k) = (P_Stokes_r - P_Stokes_l)*G%IdyCv(i,J)
2042
2043		enddo
2044		endif
2045		enddo ; enddo
2046
2047	0	if (CS%id_PFv_Stokes>0) &
2048	0	call post_data(CS%id_PFv_Stokes, PFv_Stokes, CS%diag)
2049	0	if (CS%id_PFu_Stokes>0) &
2050	0	call post_data(CS%id_PFu_Stokes, PFu_Stokes, CS%diag)
2051	0	if (CS%id_P_deltaStokes_L>0) &
2052	0	call post_data(CS%id_P_deltaStokes_L, P_deltaStokes_L, CS%diag)
2053	0	if (CS%id_P_deltaStokes_i>0) &
2054	0	call post_data(CS%id_P_deltaStokes_i, P_deltaStokes_i, CS%diag)
2055
2056	0	end subroutine Stokes_PGF
2057
2058		!> Computes wind speed from ustar_air based on COARE 3.5 Cd relationship
2059		!! Probably doesn't belong in this module, but it is used here to estimate
2060		!! wind speed for wind-wave relationships. Should be a fine way to estimate
2061		!! the neutral wind-speed as written here.
2062	313463	subroutine ust_2_u10_coare3p5(USTair, U10, GV, US, CS)
2063		real, intent(in) :: USTair !< Wind friction velocity [Z T-1 ~> m s-1]
2064		real, intent(out) :: U10 !< 10-m neutral wind speed [L T-1 ~> m s-1]
2065		type(verticalGrid_type), intent(in) :: GV !< vertical grid type
2066		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2067		type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure
2068
2069		! Local variables
2070		real :: z0sm, z0, z0rough ! Roughness lengths [Z ~> m]
2071		real :: ten_m_scale ! The 10 m reference height, in rescaled units [Z ~> m]
2072		real :: I_ten_m_scale ! The inverse of the 10 m reference height, in rescaled units [Z-1 ~> m-1]
2073		real :: u10a ! The previous guess for u10 [L T-1 ~> m s-1]
2074		real :: alpha ! The Charnock coeffient relating the wind friction velocity squared to the
2075		! roughness length [nondim]
2076		real :: Cd ! The drag coefficient [nondim]
2077		real :: I_sqrtCd ! The inverse of the square root of the drag coefficient [nondim]
2078		real :: I_vonKar ! The inverse of the von Karman coefficient [nondim]
2079		integer :: CT
2080
2081		! Uses empirical formula for z0 to convert ustar_air to u10 based on the
2082		! COARE 3.5 paper (Edson et al., 2013)
2083		! alpha=m*U10+b
2084		! Note in Edson et al. 2013, eq. 13 m is given as 0.017. However,
2085		! m=0.0017 reproduces the curve in their figure 6.
2086
2087	313463	if (CS%vonKar < 0.0) call MOM_error(FATAL, &
2088	0	"ust_2_u10_coare3p5 called with a negative value of Waves%vonKar")
2089
2090	313463	z0sm = 0.11 * CS%nu_air / USTair ! Compute z0smooth from ustar guess
2091	313463	u10a = 1000.0*US%m_s_to_L_T ! An insanely large upper bound for u10.
2092
2093	313463	if (CS%answer_date < 20230103) then
2094	0	u10 = US%Z_to_L*USTair / sqrt(0.001) ! Guess for u10
2095	0	ten_m_scale = 10.0*US%m_to_Z
2096	0	CT=0
2097	0	do while (abs(u10a/u10 - 1.) > 0.001)
2098	0	CT=CT+1
2099	0	u10a = u10
2100	0	alpha = min(CS%Charnock_min, CS%Charnock_slope_U10 * u10 + CS%Charnock_intercept)
2101	0	z0rough = alpha * (US%Z_to_LUSTair)*2 / GV%g_Earth ! Compute z0rough from ustar guess
2102	0	z0 = z0sm + z0rough
2103	0	Cd = ( CS%vonKar / log(ten_m_scale / z0) )**2 ! Compute Cd from derived roughness
2104	0	u10 = US%Z_to_L*USTair/sqrt(Cd) ! Compute new u10 from derived Cd, while loop
2105		! ends and checks for convergence...CT counter
2106		! makes sure loop doesn't run away if function
2107		! doesn't converge. This code was produced offline
2108		! and converged rapidly (e.g. 2 cycles)
2109		! for ustar=0.0001:0.0001:10.
2110	0	if (CT>20) then
2111	0	u10 = US%Z_to_L*USTair/sqrt(0.0015) ! I don't expect to get here, but just
2112		! in case it will output a reasonable value.
2113	0	exit
2114		endif
2115		enddo
2116
2117		else ! Use more efficient expressions that are mathematically equivalent to those above.
2118	313463	u10 = US%Z_to_LUSTair sqrt(1000.0) ! First guess for u10.
2119		! In the line above 1000 is the inverse of a plausible first guess of the drag coefficient.
2120	313463	I_vonKar = 1.0 / CS%vonKar
2121	313463	I_ten_m_scale = 0.1*US%Z_to_m
2122
2123	1440143	do CT=1,20
2124	1440143	if (abs(u10a - u10) <= 0.001*u10) exit ! Check for convergence.
2125	1126680	u10a = u10
2126	1126680	alpha = min(CS%Charnock_min, CS%Charnock_slope_U10 * u10 + CS%Charnock_intercept)
2127	1126680	z0rough = alpha * (CS%I_g_Earth * USTair**2) ! Compute z0rough from ustar guess
2128	1126680	z0 = z0sm + z0rough
2129	1126680	I_sqrtCd = abs(log(z0 * I_ten_m_scale)) * I_vonKar ! Compute Cd from derived roughness
2130	1126680	u10 = US%Z_to_LUSTair I_sqrtCd ! Compute new u10 from the derived Cd.
2131		enddo
2132
2133		! Output a reasonable estimate of u10 if the iteration has not converged. The hard-coded
2134		! number 25.82 is 1/sqrt(0.0015) to 4 decimal places, but the exact value should not matter.
2135	313463	if (abs(u10a - u10) > 0.001u10) u10 = US%Z_to_LUSTair * 25.82
2136		endif
2137
2138	313463	end subroutine ust_2_u10_coare3p5
2139
2140		!> Clear pointers, deallocate memory
2141	0	subroutine Waves_end(CS)
2142		type(wave_parameters_CS), pointer :: CS !< Control structure
2143
2144	0	if (allocated(CS%WaveNum_Cen)) deallocate( CS%WaveNum_Cen )
2145	0	if (allocated(CS%Freq_Cen)) deallocate( CS%Freq_Cen )
2146	0	if (allocated(CS%Us_x)) deallocate( CS%Us_x )
2147	0	if (allocated(CS%Us_y)) deallocate( CS%Us_y )
2148	0	if (allocated(CS%La_turb)) deallocate( CS%La_turb )
2149	0	if (allocated(CS%STKx0)) deallocate( CS%STKx0 )
2150	0	if (allocated(CS%STKy0)) deallocate( CS%STKy0 )
2151	0	if (allocated(CS%UStk_Hb)) deallocate( CS%UStk_Hb )
2152	0	if (allocated(CS%VStk_Hb)) deallocate( CS%VStk_Hb )
2153	0	if (allocated(CS%Omega_w2x)) deallocate( CS%Omega_w2x )
2154	0	if (allocated(CS%KvS)) deallocate( CS%KvS )
2155	0	if (allocated(CS%Us0_y)) deallocate( CS%Us0_y )
2156	0	if (allocated(CS%Us0_x)) deallocate( CS%Us0_x )
2157
2158	0	deallocate( CS )
2159
2160	0	end subroutine Waves_end
2161
2162		!> Register wave restart fields. To be called before MOM_wave_interface_init
2163	1	subroutine waves_register_restarts(CS, HI, GV, US, param_file, restart_CSp)
2164		type(wave_parameters_CS), pointer :: CS !< Wave parameter Control structure
2165		type(hor_index_type), intent(inout) :: HI !< Grid structure
2166		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
2167		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2168		type(param_file_type), intent(in) :: param_file !< Input parameter structure
2169		type(MOM_restart_CS), pointer :: restart_CSp !< Restart structure, data intent(inout)
2170		! Local variables
2171	31	type(vardesc) :: vd(2)
2172		logical :: use_waves
2173		logical :: StatisticalWaves
2174		logical :: time_tendency_term
2175		character(len=40) :: mdl = "MOM_wave_interface" !< This module's name.
2176
2177	1	if (associated(CS)) then
2178	0	call MOM_error(FATAL, "waves_register_restarts: Called with initialized waves control structure")
2179		endif
2180	1	allocate(CS)
2181
2182		call get_param(param_file, mdl, "USE_WAVES", use_waves, &
2183	1	"If true, enables surface wave modules.", do_not_log=.true., default=.false.)
2184
2185		! Check if using LA_LI2016
2186		call get_param(param_file,mdl,"USE_LA_LI2016",StatisticalWaves, &
2187	1	do_not_log=.true.,default=.false.)
2188
2189	1	if (.not.(use_waves .or. StatisticalWaves)) return
2190
2191	1	call get_param(param_file, mdl, "STOKES_DDT", time_tendency_term, do_not_log=.true., default=.false.)
2192
2193	1	if (time_tendency_term) then
2194		! Allocate wave fields needed for restart file
2195	0	allocate(CS%Us_x_prev(HI%isdB:HI%IedB,HI%jsd:HI%jed,GV%ke), source=0.0)
2196	0	allocate(CS%Us_y_prev(HI%isd:HI%Ied,HI%jsdB:HI%jedB,GV%ke), source=0.0)
2197
2198		! Register to restart files. If these are not found in a restart file, they stay 0.
2199		vd(1) = var_desc("Us_x_prev", "m s-1", "3d zonal Stokes drift profile",&
2200	0	hor_grid='u', z_grid='L')
2201		vd(2) = var_desc("Us_y_prev", "m s-1", "3d meridional Stokes drift profile",&
2202	0	hor_grid='v', z_grid='L')
2203		call register_restart_pair(CS%US_x_prev, CS%US_y_prev, vd(1), vd(2), .false., &
2204	0	restart_CSp, conversion=US%L_T_to_m_s)
2205		endif
2206
2207		end subroutine waves_register_restarts
2208
2209		!> \namespace mom_wave_interface
2210		!!
2211		!! \author Brandon Reichl, 2018.
2212		!!
2213		!! This module should be moved as wave coupling progresses and
2214		!! likely will should mirror the iceberg or sea-ice model set-up.
2215		!!
2216		!! This module is meant to contain the routines to read in and
2217		!! interpret surface wave data for MOM6. In its original form, the
2218		!! capabilities include setting the Stokes drift in the model (from a
2219		!! variety of sources including prescribed, empirical, and input
2220		!! files). In short order, the plan is to also amend the subroutine
2221		!! to accept Stokes drift information from an external coupler.
2222		!! Eventually, it will be necessary to break this file apart so that
2223		!! general wave information may be stored in the control structure
2224		!! and the Stokes drift effect can be isolated from processes such as
2225		!! sea-state dependent momentum fluxes, gas fluxes, and other wave
2226		!! related air-sea interaction and boundary layer phenomenon.
2227		!!
2228		!! The Stokes drift are stored on the C-grid with the conventional
2229		!! protocol to interpolate to the h-grid to compute Langmuir number,
2230		!! the primary quantity needed for Langmuir turbulence
2231		!! parameterizations in both the ePBL and KPP approach. This module
2232		!! also computes full 3d Stokes drift profiles, which will be useful
2233		!! if second-order type boundary layer parameterizations are
2234		!! implemented (perhaps via GOTM, work in progress).
2235
2236	0	end module MOM_wave_interface