MOM_wave_interface.F90

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
7
8use mom_data_override, only : data_override_init, data_override
9use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc
11use mom_domains, only : pass_var, pass_vector, agrid
12use mom_domains, only : to_south, to_west, to_all
13use mom_error_handler, only : mom_error, fatal, warning
14use mom_file_parser, only : get_param, log_version, param_file_type
16use mom_grid, only : ocean_grid_type
18use mom_io, only : file_exists, get_var_sizes, read_variable
19use mom_io, only : vardesc, var_desc
20use mom_safe_alloc, only : safe_alloc_ptr
22use mom_time_manager, only : time_type, operator(+), operator(/)
26use mom_restart, only : register_restart_pair, mom_restart_cs
27
28implicit none ; private
29
30#include <MOM_memory.h>
31
32public mom_wave_interface_init ! Public interface to fully initialize the wave routines.
33public query_wave_properties ! Public interface to obtain information from the waves control structure.
34public update_surface_waves ! Public interface to update wave information at the
35 ! coupler/driver level.
36public update_stokes_drift ! Public interface to update the Stokes drift profiles
37 ! called in step_mom.
38public get_langmuir_number ! Public interface to compute Langmuir number called from
39 ! ePBL or KPP routines.
40public stokes_pgf ! Public interface to compute Stokes-shear induced pressure gradient force anomaly
41public stokesmixing ! NOT READY - Public interface to add down-Stokes gradient
42 ! momentum mixing (e.g. the approach of Harcourt 2013/2015)
43public 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.
48public waves_end ! public interface to deallocate and free wave related memory.
49public get_wave_method ! public interface to obtain the wave method string
50public 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
58type, 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
261end type wave_parameters_cs
262
263! Switches needed in import_stokes_drift
264!>@{ Enumeration values for the wave method
265integer, parameter :: testprof = 0, surfbands = 1, dhh85 = 2, lf17 = 3, efactor = 4, null_wavemethod = -99
266!>@}
267!>@{ Enumeration values for the wave data source
268integer, parameter :: dataovr = 1, coupler = 2, input = 3
269!>@}
270
271! Strings for the wave method
272character*(5), parameter :: null_string = "EMPTY" !< null wave method string
273character*(12), parameter :: testprof_string = "TEST_PROFILE" !< test profile string
274character*(13), parameter :: surfbands_string = "SURFACE_BANDS" !< surface bands string
275character*(5), parameter :: dhh85_string = "DHH85" !< DHH85 wave method string
276character*(4), parameter :: lf17_string = "LF17" !< LF17 wave method string
277character*(7), parameter :: efactor_string = "EFACTOR" !< EFACTOR (based on vr12-ma) wave method string
278
279contains
280
281!> Initializes parameters related to MOM_wave_interface
282subroutine 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 if (.not. associated(cs)) then
305 call mom_error(fatal, "wave_interface_init called without an associated control structure.")
306 return
307 endif
308
309 call get_param(param_file, mdl, "USE_WAVES", use_waves, &
310 "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 do_not_log=.true.,default=.false.)
315
316 if (.not.(use_waves .or. statisticalwaves)) return
317
318 cs%UseWaves = use_waves
319 cs%diag => diag
320 cs%Time => time
321
322 cs%g_Earth = gv%g_Earth_Z_T2
323 cs%I_g_Earth = 1.0 / cs%g_Earth
324
325 ! Add any initializations needed here
326 cs%DataOver_initialized = .false.
327
328 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 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 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
343 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 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 "number calculation.", units="m", default=0.1, scale=us%m_to_Z)
353
354 if (statisticalwaves) then
355 cs%WaveMethod = lf17
356 call set_lf17_wave_params(param_file, mdl, gv, us, cs)
357 if (.not.use_waves) return
358 else
359 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 do_not_log=.not.use_waves)
367 if (cs%LagrangianMixing) then
368 ! Force Code Intervention
369 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 do_not_log=.not.use_waves)
374 if (cs%StokesMixing) then
375 ! Force Code Intervention
376 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 do_not_log=.not.use_waves)
381 if (cs%CoriolisStokes) then
382 ! Force Code Intervention
383 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 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 default=.false.)
392 call get_param(param_file, mdl, "STOKES_PGF", cs%Stokes_PGF, &
393 "Flag to use Stokes-induced pressure gradient anomaly", &
394 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 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 default=.false.)
403 call get_param(param_file, mdl, "STOKES_DDT", cs%Stokes_DDT, &
404 "Flag to use Stokes d/dt", &
405 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 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 default=null_string)
427 select case (trim(tmpstring1))
428 case (null_string)! No Waves
429 call mom_error(fatal, "wave_interface_init called with no specified "//&
430 "WAVE_METHOD.")
431 case (testprof_string)! Test Profile
432 cs%WaveMethod = testprof
433 call get_param(param_file, mdl, "TP_STKX_SURF", cs%TP_STKX0, &
434 'Surface Stokes (x) for test profile', &
435 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 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 units='m', default=50.0, scale=us%m_to_Z)
442 case (surfbands_string)! Surface Stokes Drift Bands
443 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 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 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 " INPUT - Testing with fixed values.", default=null_string)
458 select case (trim(tmpstring2))
459 case (null_string)! Default
460 call mom_error(fatal, "wave_interface_init called with SURFACE_BANDS "//&
461 "but no SURFBAND_SOURCE.")
462 case (dataovr_string)! Using Data Override
463 cs%DataSource = dataovr
464 call get_param(param_file, mdl, "SURFBAND_FILENAME", cs%SurfBandFileName, &
465 "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 units="m s-1", default=10.0, scale=us%m_s_to_L_T)
471 case (coupler_string)! Reserved for coupling
472 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 "or the model will fail.", default=1)
478 allocate( cs%WaveNum_Cen(cs%NumBands), source=0.0 )
479 allocate( cs%STKx0(g%isdB:g%iedB,g%jsd:g%jed,cs%NumBands), source=0.0 )
480 allocate( cs%STKy0(g%isd:g%ied,g%jsdB:g%jedB,cs%NumBands), source=0.0 )
481 allocate( cs%UStk_Hb(g%isc:g%iec,g%jsc:g%jec,cs%NumBands), source=0.0 )
482 allocate( cs%VStk_Hb(g%isc:g%iec,g%jsc:g%jec,cs%NumBands), source=0.0 )
483 allocate( cs%Omega_w2x(g%isc:g%iec,g%jsc:g%jec) , source=0.0 )
484 cs%PartitionMode = 0
485 call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", cs%WaveNum_Cen, &
486 "Central wavenumbers for surface Stokes drift bands.", &
487 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 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 "STOKES_Y, there are no safety checks in the code.", default=1)
494 allocate( cs%WaveNum_Cen(1:cs%NumBands), source=0.0 )
495 allocate( cs%PrescribedSurfStkX(1:cs%NumBands), source=0.0 )
496 allocate( cs%PrescribedSurfStkY(1:cs%NumBands), source=0.0 )
497 allocate( cs%STKx0(g%isdB:g%iedB,g%jsd:g%jed,1:cs%NumBands), source=0.0 )
498 allocate( cs%STKy0(g%isd:g%ied,g%jsdB:g%jedB,1:cs%NumBands), source=0.0 )
499
500 cs%PartitionMode = 0
501 call get_param(param_file, mdl, "SURFBAND_WAVENUMBERS", cs%WaveNum_Cen, &
502 "Central wavenumbers for surface Stokes drift bands.", &
503 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 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 units='m/s', default=0.0, scale=us%m_s_to_L_T)
510 case default! No method provided
511 call mom_error(fatal,'Check WAVE_METHOD.')
512 end select
513
514 case (dhh85_string) !Donelan et al., 1985 spectrum
515 cs%WaveMethod = dhh85
516 call mom_error(warning,"DHH85 only ever set-up for uniform cases w/ "//&
517 "Stokes drift in x-direction.")
518 call get_param(param_file, mdl, "DHH85_AGE_FP", cs%WaveAgePeakFreq, &
519 "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 units='nondim', default=1.2)
523 call get_param(param_file, mdl, "DHH85_WIND", cs%WaveWind, &
524 "Wind speed for DHH85 spectrum.", &
525 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 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 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 "Flag to disable updating DHH85 Stokes drift.", default=.false.)
534 case (lf17_string) !Li and Fox-Kemper 17 wind-sea Langmuir number
535 cs%WaveMethod = lf17
536 call set_lf17_wave_params(param_file, mdl, gv, us, cs)
537 case (efactor_string) !Li and Fox-Kemper 16
538 cs%WaveMethod = efactor
539 case default
540 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 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 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 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 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 allocate(cs%Us_x(g%isdB:g%IedB,g%jsd:g%jed,g%ke), source=0.0)
564 allocate(cs%Us_y(g%isd:g%Ied,g%jsdB:g%jedB,g%ke), source=0.0)
565 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 allocate(cs%ddt_Us_x(g%isdB:g%IedB,g%jsd:g%jed,g%ke), source=0.0)
569 allocate(cs%ddt_Us_y(g%isd:g%Ied,g%jsdB:g%jedB,g%ke), source=0.0)
570 allocate(cs%Us_x_from_ddt(g%isdB:g%IedB,g%jsd:g%jed,g%ke), source=0.0)
571 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 allocate(cs%US0_x(g%isdB:g%iedB,g%jsd:g%jed), source=0.0)
575 allocate(cs%US0_y(g%isd:g%ied,g%jsdB:g%jedB), source=0.0)
576 ! c. Langmuir number
577 allocate(cs%La_turb(g%isc:g%iec,g%jsc:g%jec), source=0.0)
578 ! d. Viscosity for Stokes drift
579 if (cs%StokesMixing) then
580 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 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 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 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 cs%diag%axesCuL,time,'3d Stokes drift (x)', 'm s-1', conversion=us%L_T_to_m_s)
592 if (cs%Stokes_DDT) then
593 cs%id_ddt_3dstokes_y = register_diag_field('ocean_model','dvdt_Stokes', &
594 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 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 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 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 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 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 '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 'm2 s-2',conversion=us%L_T_to_m_s**2)
612 cs%id_La_turb = register_diag_field('ocean_model','La_turbulent', &
613 cs%diag%axesT1,time,'Surface (turbulent) Langmuir number','nondim')
614
615end subroutine mom_wave_interface_init
616
617!> Set the parameters that are used to determine the averaged Stokes drift and Langmuir numbers
618subroutine 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 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 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 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 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 units="s2 m-1", default=0.0246, scale=us%m_to_Z*us%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 "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 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 units="nondim", default=-0.005)
660
661end subroutine set_lf17_wave_params
662
663!> This interface provides the caller with information from the waves control structure.
664subroutine 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 if (present(numbands)) numbands = cs%NumBands
675 if (present(wavenumbers)) then
676 if (size(wavenumbers) < cs%NumBands) call mom_error(fatal, "query_wave_properties called "//&
677 "with a Wavenumbers array that is smaller than the number of bands.")
678 if (present(us)) then
679 do n=1,cs%NumBands ; wavenumbers(n) = us%m_to_Z * cs%WaveNum_Cen(n) ; enddo
680 else
681 do n=1,cs%NumBands ; wavenumbers(n) = cs%WaveNum_Cen(n) ; enddo
682 endif
683 endif
684
685end subroutine query_wave_properties
686
687!> Subroutine that handles updating of surface wave/Stokes drift related properties
688subroutine 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 if (cs%WaveMethod == testprof) then
701 ! Do nothing
702 elseif (cs%WaveMethod == surfbands) then
703 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 stokes_time = time_present + dt/2
709 call surface_bands_by_data_override(stokes_time, g, gv, us, cs)
710 elseif (cs%DataSource == coupler) then
711 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 "otherwise select another option for SURFBAND_SOURCE.")
716 endif
717 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 "ww3_grid.inp, and that your mod_def.ww3 is up to date.")
721 endif
722
723 do b=1,cs%NumBands
724 cs%WaveNum_Cen(b) = forces%stk_wavenumbers(b)
725 !Interpolate from a grid to c grid
726 do j=g%jsc,g%jec
727 do i=g%iscB,g%iecB
728 cs%STKx0(i,j,b) = 0.5*(forces%UStkb(i,j,b)+forces%UStkb(i+1,j,b))
729 enddo
730 enddo
731 do j=g%jscB,g%jecB
732 do i=g%isc,g%iec
733 cs%STKY0(i,j,b) = 0.5*(forces%VStkb(i,j,b)+forces%VStkb(i,j+1,b))
734 enddo
735 enddo
736 call pass_vector(cs%STKx0(:,:,b),cs%STKy0(:,:,b), g%Domain)
737 enddo
738 do j=g%jsc,g%jec
739 do i=g%isc,g%iec
740 cs%Omega_w2x(i,j) = forces%omega_w2x(i,j)
741 do b=1,cs%NumBands
742 cs%UStk_Hb(i,j,b) = forces%UStkb(i,j,b)
743 cs%VStk_Hb(i,j,b) = forces%VStkb(i,j,b)
744 enddo
745 enddo
746 enddo
747 elseif (cs%DataSource == input) then
748 do b=1,cs%NumBands
749 do j=g%jsd,g%jed
750 do i=g%isdB,g%iedB
751 cs%STKx0(i,j,b) = cs%PrescribedSurfStkX(b)
752 enddo
753 enddo
754 do j=g%jsdB, g%jedB
755 do i=g%isd,g%ied
756 cs%STKY0(i,j,b) = cs%PrescribedSurfStkY(b)
757 enddo
758 enddo
759 enddo
760 endif
761 endif
762
763end subroutine update_surface_waves
764
765!> Constructs the Stokes Drift profile on the model grid based on
766!! desired coupling options
767subroutine 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 if (cs%WaveMethod==efactor) return
792
793 if (allocated(cs%US_x) .and. allocated(cs%US_y)) then
794 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 if (cs%WaveMethod==testprof) then
800 pi = 4.0*atan(1.0)
801 decayscale = 4.*pi / cs%TP_WVL !4pi
802 do j=g%jsc,g%jec
803 do i=g%iscB,g%iecB
804 bottom = 0.0
805 midpoint = 0.0
806 do k = 1,gv%ke
807 top = bottom
808 midpoint = bottom - 0.25*(dz(i,j,k)+dz(i+1,j,k))
809 bottom = bottom - 0.5*(dz(i,j,k)+dz(i+1,j,k))
810 cs%Us_x(i,j,k) = cs%TP_STKX0*exp(midpoint*decayscale)
811 enddo
812 enddo
813 enddo
814 do j=g%jscB,g%jecB
815 do i=g%isc,g%iec
816 bottom = 0.0
817 midpoint = 0.0
818 do k = 1,gv%ke
819 top = bottom
820 midpoint = bottom - 0.25*(dz(i,j,k)+dz(i,j+1,k))
821 bottom = bottom - 0.5*(dz(i,j,k)+dz(i,j+1,k))
822 cs%Us_y(i,j,k) = cs%TP_STKY0*exp(midpoint*decayscale)
823 enddo
824 enddo
825 enddo
826 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 elseif (cs%WaveMethod==surfbands) then
831 cs%Us_x(:,:,:) = 0.0
832 cs%Us_y(:,:,:) = 0.0
833 cs%Us0_x(:,:) = 0.0
834 cs%Us0_y(:,:) = 0.0
835 ! Computing X direction Stokes drift
836 do j=g%jsc,g%jec
837 do i=g%iscB,g%iecB
838 ! 1. First compute the surface Stokes drift
839 ! by summing over the partitions.
840 do b = 1,cs%NumBands
841 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 bottom = 0.0
845 do k = 1,gv%ke
846 top = bottom
847 level_thick = 0.5*(dz(i,j,k)+dz(i+1,j,k))
848 midpoint = top - 0.5*level_thick
849 bottom = top - level_thick
850
851 if (cs%answer_date >= 20230101) then
852 ! Use more accurate and numerically stable expressions that work even for vanished layers.
853 do b = 1,cs%NumBands
854 if (cs%PartitionMode == 0) then
855 ! Average over a layer using the bin's central wavenumber.
856 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 wn = cs%Freq_Cen(b)**2 * cs%I_g_Earth
860 cmn_fac = exp(2.*wn*top) * one_minus_exp_x(2.*wn*level_thick)
861 endif
862 cs%US_x(i,j,k) = cs%US_x(i,j,k) + cs%STKx0(i,j,b)*cmn_fac
863 enddo
864
865 elseif (level_thick > cs%Stokes_min_thick_avg) then
866 ! -> Stokes drift in thin layers not averaged.
867 do b = 1,cs%NumBands
868 if (cs%PartitionMode == 0) then
869 ! In wavenumber we are averaging over level
870 cmn_fac = (exp(top*2.*cs%WaveNum_Cen(b))-exp(bottom*2.*cs%WaveNum_Cen(b))) &
871 / ((top-bottom)*(2.*cs%WaveNum_Cen(b)))
872 else
873 ! Use a numerical integration and then divide by layer thickness
874 wn = cs%Freq_Cen(b)**2 / cs%g_Earth !bgr bug-fix missing g
875 cmn_fac = (exp(2.*wn*top)-exp(2.*wn*bottom)) / (2.*wn*(top-bottom))
876 endif
877 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 do b = 1,cs%NumBands
881 if (cs%PartitionMode == 0) then
882 cmn_fac = exp(midpoint * 2. * cs%WaveNum_Cen(b))
883 else
884 cmn_fac = exp(midpoint * 2. * cs%Freq_Cen(b)**2 / cs%g_Earth)
885 endif
886 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 do j=g%jscB,g%jecB
895 do i=g%isc,g%iec
896 ! Set the surface value to that at z=0
897 do b = 1,cs%NumBands
898 cs%US0_y(i,j) = cs%US0_y(i,j) + cs%STKy0(i,j,b)
899 enddo
900 ! Compute the level averages.
901 bottom = 0.0
902 do k = 1,gv%ke
903 top = bottom
904 level_thick = 0.5*(dz(i,j,k)+dz(i,j+1,k))
905 midpoint = top - 0.5*level_thick
906 bottom = top - level_thick
907
908 if (cs%answer_date >= 20230101) then
909 ! Use more accurate and numerically stable expressions that work even for vanished layers.
910 do b = 1,cs%NumBands
911 if (cs%PartitionMode == 0) then
912 ! Average over a layer using the bin's central wavenumber.
913 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 wn = cs%Freq_Cen(b)**2 * cs%I_g_Earth
917 cmn_fac = exp(2.*wn*top) * one_minus_exp_x(2.*wn*level_thick)
918 endif
919 cs%US_y(i,j,k) = cs%US_y(i,j,k) + cs%STKy0(i,j,b)*cmn_fac
920 enddo
921 elseif (level_thick > cs%Stokes_min_thick_avg) then
922 ! -> Stokes drift in thin layers not averaged.
923 do b = 1,cs%NumBands
924 if (cs%PartitionMode == 0) then
925 ! In wavenumber we are averaging over level
926 cmn_fac = (exp(top*2.*cs%WaveNum_Cen(b))-exp(bottom*2.*cs%WaveNum_Cen(b))) &
927 / ((top-bottom)*(2.*cs%WaveNum_Cen(b)))
928 else
929 ! Use a numerical integration and then divide by layer thickness
930 wn = cs%Freq_Cen(b)**2 / cs%g_Earth !bgr bug-fix missing g
931 cmn_fac = (exp(2.*wn*top)-exp(2.*wn*bottom)) / (2.*wn*(top-bottom))
932 endif
933 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 do b = 1,cs%NumBands
937 if (cs%PartitionMode == 0) then
938 cmn_fac = exp(midpoint*2.*cs%WaveNum_Cen(b))
939 else
940 cmn_fac = exp(midpoint * 2. * cs%Freq_Cen(b)**2 / cs%g_Earth)
941 endif
942 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 call pass_vector(cs%Us_x(:,:,:),cs%Us_y(:,:,:), g%Domain, to_all)
949 call pass_vector(cs%Us0_x(:,:),cs%Us0_y(:,:), g%Domain)
950 elseif (cs%WaveMethod == dhh85) then
951 if (.not.(cs%StaticWaves .and. cs%DHH85_is_set)) then
952 do j=g%jsc,g%jec
953 do i=g%iscB,g%iecB
954 bottom = 0.0
955 do k = 1,gv%ke
956 top = bottom
957 midpoint = top - 0.25*(dz(i,j,k)+dz(i+1,j,k))
958 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 call dhh85_mid(gv, us, cs, midpoint, ustokes)
964 ! Putting into x-direction (no option for direction
965 cs%US_x(i,j,k) = ustokes
966 enddo
967 enddo
968 enddo
969 do j=g%jscB,g%jecB
970 do i=g%isc,g%iec
971 bottom = 0.0
972 do k = 1,gv%ke
973 top = bottom
974 midpoint = bottom - 0.25*(dz(i,j,k)+dz(i,j+1,k))
975 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 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 cs%DHH85_is_set = .true.
991 endif
992 call pass_vector(cs%Us_x(:,:,:),cs%Us_y(:,:,:), g%Domain)
993 else! Keep this else, fallback to 0 Stokes drift
994 cs%Us_x(:,:,:) = 0.
995 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 do j=g%jsc, g%jec
1002 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 dz(i,j,:), cs, override_ma=.false.)
1005 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 if ( dynamics_step .and. cs%Stokes_DDT ) then
1012 i_dt = 1.0 / dt
1013 cs%ddt_us_x(:,:,:) = (cs%US_x(:,:,:) - cs%US_x_prev(:,:,:)) * i_dt
1014 cs%ddt_us_y(:,:,:) = (cs%US_y(:,:,:) - cs%US_y_prev(:,:,:)) * i_dt
1015 cs%US_x_prev(:,:,:) = cs%US_x(:,:,:)
1016 cs%US_y_prev(:,:,:) = cs%US_y(:,:,:)
1017 endif
1018
1019 ! Output any desired quantities
1020 if (cs%id_surfacestokes_y>0) &
1021 call post_data(cs%id_surfacestokes_y, cs%us0_y, cs%diag)
1022 if (cs%id_surfacestokes_x>0) &
1023 call post_data(cs%id_surfacestokes_x, cs%us0_x, cs%diag)
1024 if (cs%id_3dstokes_y>0) &
1025 call post_data(cs%id_3dstokes_y, cs%us_y, cs%diag)
1026 if (cs%id_3dstokes_x>0) &
1027 call post_data(cs%id_3dstokes_x, cs%us_x, cs%diag)
1028 if (cs%Stokes_DDT) then
1029 if (cs%id_ddt_3dstokes_x>0) &
1030 call post_data(cs%id_ddt_3dstokes_x, cs%ddt_us_x, cs%diag)
1031 if (cs%id_ddt_3dstokes_y>0) &
1032 call post_data(cs%id_ddt_3dstokes_y, cs%ddt_us_y, cs%diag)
1033 if (cs%id_3dstokes_x_from_ddt>0) &
1034 call post_data(cs%id_3dstokes_x_from_ddt, cs%us_x_from_ddt, cs%diag)
1035 if (cs%id_3dstokes_y_from_ddt>0) &
1036 call post_data(cs%id_3dstokes_y_from_ddt, cs%us_y_from_ddt, cs%diag)
1037 endif
1038 if (cs%id_La_turb>0) &
1039 call post_data(cs%id_La_turb, cs%La_turb, cs%diag)
1040
1041end subroutine update_stokes_drift
1042
1043!> Return the value of (1 - exp(-x))/x [nondim], using an accurate expression for small values of x.
1044real 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 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 one_minus_exp_x = 1.0 - x * (0.5 - c1_6*x)
1050 else
1051 one_minus_exp_x = (1.0 - exp(-x)) / x
1052 endif
1053end function one_minus_exp_x
1054
1055!> Return the value of (1 - exp(-x)) [nondim], using an accurate expression for small values of x.
1056real 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 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 one_minus_exp = x * (1.0 - x * (0.5 - c1_6*x))
1062 else
1063 one_minus_exp = 1.0 - exp(-x)
1064 endif
1065end function one_minus_exp
1066
1067!> A subroutine to fill the Stokes drift from a NetCDF file
1068!! using the data_override procedures.
1069subroutine 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 real :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal Stokes drift of band at h-points [L T-1 ~> m s-1]
1078 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 if (.not.cs%DataOver_initialized) then
1088 call data_override_init(g%Domain)
1089 cs%DataOver_initialized = .true.
1090
1091 if (.not.file_exists(cs%SurfBandFileName)) &
1092 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 call get_var_sizes(cs%SurfBandFileName, 'wavenumber', ndims, sizes, dim_names=dim_name)
1098 wavenumber_exists = (ndims > -1)
1099
1100 if (.not.wavenumber_exists) then
1101 ! Read the number of frequency bands in the file, if the variable 'frequency' exists.
1102 call get_var_sizes(cs%SurfBandFileName, 'frequency', ndims, sizes, dim_names=dim_name)
1103 if (ndims < 0) &
1104 call mom_error(fatal, "error finding variable 'wavenumber' or 'frequency' in file "//&
1105 trim(cs%SurfBandFileName)//" in MOM_wave_interface.")
1106 endif
1107
1108 cs%NUMBANDS = sizes(1)
1109 ! Allocate the wavenumber bins
1110 allocate( cs%WaveNum_Cen(cs%NUMBANDS), source=0.0 )
1111
1112 if (wavenumber_exists) then
1113 ! Wavenumbers found, so this file uses the old method:
1114 cs%PartitionMode = 0
1115
1116 ! Reading wavenumber bins
1117 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 cs%PartitionMode = 1
1122
1123 ! Allocate the frequency bins
1124 allocate( cs%Freq_Cen(cs%NUMBANDS), source=0.0 )
1125
1126 ! Reading frequencies
1127 pi = 4.0*atan(1.0)
1128 call read_variable(cs%SurfBandFileName, dim_name(1), cs%Freq_Cen, scale=2.*pi*us%T_to_s)
1129
1130 do b = 1,cs%NumBands
1131 cs%WaveNum_Cen(b) = cs%Freq_Cen(b)**2 / cs%g_Earth
1132 enddo
1133 endif
1134
1135 if (.not.allocated(cs%STKx0)) then
1136 allocate( cs%STKx0(g%isdB:g%iedB,g%jsd:g%jed,cs%NUMBANDS), source=0.0 )
1137 endif
1138 if (.not.allocated(cs%STKy0)) then
1139 allocate( cs%STKy0(g%isd:g%ied,g%jsdB:g%jedB,cs%NUMBANDS), source=0.0 )
1140 endif
1141 endif
1142
1143 do b = 1,cs%NumBands
1144 temp_x(:,:) = 0.0
1145 temp_y(:,:) = 0.0
1146 varname = ' '
1147 write(varname, "(A3,I0)") 'Usx', b
1148 call data_override(g%Domain, trim(varname), temp_x, time, scale=us%m_s_to_L_T)
1149 varname = ' '
1150 write(varname, "(A3,I0)") 'Usy', b
1151 call data_override(g%Domain, trim(varname), temp_y, time, scale=us%m_s_to_L_T)
1152 ! Update halo on h-grid
1153 call pass_vector(temp_x, temp_y, g%Domain, to_all, agrid)
1154 ! Filter land values
1155 do j = g%jsd,g%jed
1156 do i = g%Isd,g%Ied
1157 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 temp_x(i,j) = 0.0
1160 temp_y(i,j) = 0.0
1161 endif
1162 enddo
1163 enddo
1164 if (cs%Homogenize_Surfbands) then
1165 avgx = global_area_mean(temp_x, g)
1166 avgy = global_area_mean(temp_y, g)
1167 do j = g%jsd,g%jed ; do i = g%Isd,g%Ied ; if (g%mask2dT(i,j) > 0.0) then
1168 temp_y(i,j) = avgy
1169 temp_x(i,j) = avgx
1170 endif ; enddo ; enddo
1171 endif
1172
1173 ! Interpolate to u/v grids
1174 do j = g%jsc,g%jec
1175 do i = g%IscB,g%IecB
1176 cs%STKx0(i,j,b) = 0.5 * (temp_x(i,j) + temp_x(i+1,j))
1177 enddo
1178 enddo
1179 do j = g%JscB,g%JecB
1180 do i = g%isc,g%iec
1181 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 call pass_vector(cs%STKx0(:,:,:), cs%STKy0(:,:,:), g%Domain, to_all)
1188
1189end 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.
1196subroutine get_langmuir_number( LA, G, GV, US, HBL, ustar, i, j, dz, Waves, &
1197 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 real, dimension(SZK_(GV)) :: us_h, vs_h ! Profiles of Stokes velocities [L T-1 ~> m s-1]
1226 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 dpt_lasl = -1.0*max(waves%LA_FracHBL*hbl, waves%LA_HBL_min)
1231
1232 use_ma = waves%LA_Misalignment
1233 if (present(override_ma)) use_ma = override_ma
1234
1235 ! If requesting to use misalignment in the Langmuir number compute the Shear Direction
1236 if (use_ma) then
1237 if (.not.(present(u_h).and.present(v_h))) call mom_error(fatal, &
1238 "Get_LA_waves requested to consider misalignment, but velocities were not provided.")
1239 continueloop = .true.
1240 bottom = 0.0
1241 do k = 1,gv%ke
1242 top = bottom
1243 midpoint = bottom + 0.5*dz(k)
1244 bottom = bottom + dz(k)
1245
1246 if (waves%LA_Misalign_bug) then
1247 ! Given the sign convention that Dpt_LASL is negative, the next line has a bug.
1248 if (midpoint > dpt_lasl .and. k > 1 .and. continueloop) then
1249 sheardirection = atan2(v_h(1)-v_h(k), u_h(1)-u_h(k))
1250 continueloop = .false.
1251 endif
1252 else ! This version avoids the bug in the version above.
1253 if (midpoint > abs(dpt_lasl) .and. (k > 1) .and. continueloop) then
1254 sheardirection = atan2(v_h(1)-v_h(k), u_h(1)-u_h(k))
1255 continueloop = .false.
1256 endif
1257 endif
1258 enddo
1259 endif
1260
1261 if (waves%WaveMethod==testprof) then
1262 do k = 1,gv%ke
1263 us_h(k) = 0.5*(waves%US_X(i,j,k)+waves%US_X(i-1,j,k))
1264 vs_h(k) = 0.5*(waves%US_Y(i,j,k)+waves%US_Y(i,j-1,k))
1265 enddo
1266 call get_sl_average_prof( gv, dpt_lasl, dz, us_h, la_stkx)
1267 call get_sl_average_prof( gv, dpt_lasl, dz, vs_h, la_stky)
1268 la_stk = sqrt((la_stkx*la_stkx) + (la_stky*la_stky))
1269 elseif (waves%WaveMethod==surfbands) then
1270 allocate(stkband_x(waves%NumBands), stkband_y(waves%NumBands))
1271 do bb = 1,waves%NumBands
1272 stkband_x(bb) = 0.5*(waves%STKx0(i,j,bb)+waves%STKx0(i-1,j,bb))
1273 stkband_y(bb) = 0.5*(waves%STKy0(i,j,bb)+waves%STKy0(i,j-1,bb))
1274 enddo
1275 call get_sl_average_band(gv, dpt_lasl, waves%NumBands, waves%WaveNum_Cen, stkband_x, la_stkx )
1276 call get_sl_average_band(gv, dpt_lasl, waves%NumBands, waves%WaveNum_Cen, stkband_y, la_stky )
1277 la_stk = sqrt((la_stkx**2) + (la_stky**2))
1278 deallocate(stkband_x, stkband_y)
1279 elseif (waves%WaveMethod==dhh85) then
1280 ! Temporarily integrating profile rather than spectrum for simplicity
1281 do k = 1,gv%ke
1282 us_h(k) = 0.5*(waves%US_X(i,j,k)+waves%US_X(i-1,j,k))
1283 vs_h(k) = 0.5*(waves%US_Y(i,j,k)+waves%US_Y(i,j-1,k))
1284 enddo
1285 call get_sl_average_prof( gv, dpt_lasl, dz, us_h, la_stkx)
1286 call get_sl_average_prof( gv, dpt_lasl, dz, vs_h, la_stky)
1287 la_stk = sqrt((la_stkx**2) + (la_stky**2))
1288 elseif (waves%WaveMethod==lf17) then
1289 call get_stokessl_lifoxkemper(ustar, hbl*waves%LA_FracHBL, gv, us, waves, la_stk, la)
1290 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 "a wave method (or set USE_LA_LI2016 to use statistical waves).")
1294 endif
1295
1296 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 la = max(waves%La_min, sqrt(us%Z_to_L*ustar / (la_stk + waves%La_Stk_backgnd)))
1301 endif
1302
1303 if (use_ma) then
1304 wavedirection = atan2(la_stky, la_stkx)
1305 la = la / sqrt(max(1.e-8, cos( wavedirection - sheardirection)))
1306 endif
1307
1308end subroutine get_langmuir_number
1309
1310!> function to return the wave method string set in the param file
1311function get_wave_method(CS)
1312 character(:), allocatable :: get_wave_method
1313 type(wave_parameters_cs), pointer :: cs !< Control structure
1314
1315 if (associated(cs)) then
1316 select case(cs%WaveMethod)
1317 case (null_wavemethod)
1318 get_wave_method = null_string
1319 case (testprof)
1320 get_wave_method = testprof_string
1321 case (surfbands)
1322 get_wave_method = surfbands_string
1323 case (dhh85)
1324 get_wave_method = dhh85_string
1325 case (lf17)
1326 get_wave_method = lf17_string
1327 case (efactor)
1328 get_wave_method = efactor_string
1329 end select
1330 else
1331 get_wave_method = null_string
1332 endif
1333end 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"
1351subroutine 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 pi = 4.0*atan(1.0)
1383 ustokes_sl = 0.0
1384 la = 1.e8
1385 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 call ust_2_u10_coare3p5(ustar*sqrt(cs%rho_ocn/cs%rho_air), u10, gv, us, cs)
1390 ! surface Stokes drift
1391 ustokes = us_to_u10*u10
1392 !
1393 ! significant wave height from Pierson-Moskowitz spectrum (Bouws, 1998)
1394 hm0 = cs%SWH_from_u10sq * u10**2
1395 !
1396 ! peak frequency (PM, Bouws, 1998)
1397 fp = 0.877 * (us%L_to_Z*gv%g_Earth) / (2.0 * pi * u19p5_to_u10 * u10)
1398 !
1399 ! mean frequency
1400 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 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 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.5*0.125 * r_loss * fm_into_fp * 0.877 * CS%SWH_from_u10sq**2)) / &
1416 ! (GV%g_Earth * u10**2)
1417
1418 ! surface layer
1419 z0 = abs(hbl)
1420
1421 if (cs%answer_date < 20230102) then
1422 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 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 ( 1.0 - exp(-2.0 * kphil * z0) )
1432 r2 = -( 0.84 + 0.0591 / kphil * z0i ) * &
1433 sqrt( 2.0 * pi * kphil * z0 ) * &
1434 erfc( sqrt( 2.0 * kphil * z0 ) )
1435 r3 = ( 0.0632 / kstar * z0i + 0.125 ) * &
1436 (1.0 - exp(-2.0 * kstar * z0) )
1437 r4 = ( 0.125 + 0.0946 / kstar * z0i ) * &
1438 sqrt( 2.0 * pi * kstar * z0) * &
1439 erfc( sqrt( 2.0 * kstar * z0 ) )
1440 ustokes_sl = ustokes * (0.715 + r1 + r2 + r3 + r4)
1441 else
1442 ! The following is equivalent to the code above, but avoids singularities
1443 r1 = ( 0.302 - 1.68*(kphil*z0) ) * one_minus_exp_x(2.0 * (kphil * z0))
1444 r3 = ( 0.1264 + 0.64*(kphil*z0) ) * one_minus_exp_x(5.12 * (kphil * z0))
1445
1446 root_2kz = sqrt(2.0 * kphil * z0)
1447 ! r2 = -( 0.84 + 0.0591*2.0 / (root_2kz**2) ) * sqrt(PI) * root_2kz * erfc( root_2kz )
1448 ! r4 = ( 0.2 + 0.059125*2.0 / (root_2kz**2) ) * 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 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 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)*z**2 + (1/10)*z**4 - (1/42)*z**6 + ...)
1461 ! and then cancel or combine common terms and drop negligibly small terms.
1462 r5 = -0.64*sqrt(pi)*root_2kz + (-0.14184 + 1.0839648 * root_2kz**2)
1463 endif
1464 ustokes_sl = ustokes * (0.715 + ((r1 + r3) + r5))
1465 endif
1466
1467 if (ustokes_sl /= 0.0) la = sqrt(us%Z_to_L*ustar / ustokes_sl)
1468 endif
1469
1470end subroutine get_stokessl_lifoxkemper
1471
1472!> Get SL Averaged Stokes drift from a Stokes drift Profile
1473subroutine 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 sum = 0.0
1491
1492 ! Integrate
1493 bottom = 0.0
1494 do k = 1, gv%ke
1495 top = bottom
1496 bottom = bottom - dz(k)
1497 if (avgdepth < bottom) then ! The whole cell is within H_LA
1498 sum = sum + profile(k) * dz(k)
1499 elseif (avgdepth < top) then ! A partial cell is within H_LA
1500 sum = sum + profile(k) * (top-avgdepth)
1501 exit
1502 else
1503 exit
1504 endif
1505 enddo
1506
1507 ! Divide by AvgDepth or the depth in the column, whichever is smaller.
1508 if (abs(avgdepth) <= abs(bottom)) then
1509 average = sum / abs(avgdepth)
1510 elseif (abs(bottom) > 0.0) then
1511 average = sum / abs(bottom)
1512 else
1513 average = 0.0
1514 endif
1515
1516end subroutine get_sl_average_prof
1517
1518!> Get SL averaged Stokes drift from the banded Spectrum method
1519subroutine 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 average = 0.0
1535 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 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
1546end 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.2*grav*2pi/dz
1554subroutine 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 u10 = cs%WaveWind*us%L_to_Z
1573
1574 !/
1575 nomega = 1000
1576 domega = (cs%omega_max - cs%omega_min) / real(nomega)
1577
1578 !
1579 if (cs%WaveAgePeakFreq) then
1580 omega_peak = cs%g_Earth / (cs%WaveAge * u10)
1581 else
1582 pi = 4.0*atan(1.0)
1583 omega_peak = 2. * pi * 0.13 * cs%g_Earth / u10
1584 endif
1585 !/
1586 ann = 0.006 * cs%WaveAge**(-0.55)
1587 bnn = 1.0
1588 snn = 0.08 * (1.0 + 4.0 * cs%WaveAge**3)
1589 cnn = 1.7
1590 if (cs%WaveAge < 1.) then
1591 cnn = cnn - 6.0*log10(cs%WaveAge)
1592 endif
1593 !/
1594 ustokes = 0.0
1595 omega = cs%omega_min + 0.5*domega
1596 do oi = 1,nomega-1
1597 dnn = exp( -0.5 * (omega-omega_peak)**2 / (snn**2 * omega_peak**2) )
1598 ! wavespec units [L Z T ~> m2 s]
1599 wavespec = us%Z_to_L * (ann * cs%g_Earth**2 / (omega_peak*omega**4 ) ) * &
1600 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 exp( 2.0 * omega**2 * zpt / cs%g_Earth) / cs%g_Earth
1604 ustokes = ustokes + stokes*domega
1605 omega = omega + domega
1606 enddo
1607
1608end subroutine dhh85_mid
1609
1610!> Explicit solver for Stokes mixing.
1611!! Still in development do not use.
1612subroutine 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 do k = 1, gv%ke
1638 do j = g%jsc, g%jec
1639 do i = g%iscB, g%iecB
1640 h_lay = 0.5*(h(i,j,k)+h(i+1,j,k))
1641 dz_lay = 0.5*(dz(i,j,k)+dz(i+1,j,k))
1642 dtauup = 0.0
1643 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.5*(dz_lay + 0.5*(dz(i,j,k-1)+dz(i+1,j,k-1)) ))
1647 dtaudn = 0.0
1648 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.5*(dz_lay + 0.5*(dz(i,j,k+1)+dz(i+1,j,k+1)) ))
1652 u(i,j,k) = u(i,j,k) + dt * (dtauup-dtaudn) / h_lay
1653 enddo
1654 enddo
1655 enddo
1656
1657 do k = 1, gv%ke
1658 do j = g%jscB, g%jecB
1659 do i = g%isc, g%iec
1660 h_lay = 0.5*(h(i,j,k)+h(i,j+1,k))
1661 dz_lay = 0.5*(dz(i,j,k)+dz(i,j+1,k))
1662 dtauup = 0.
1663 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.5*(dz_lay + 0.5*(dz(i,j,k-1)+dz(i,j+1,k-1)) ))
1667 dtaudn = 0.0
1668 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.5*(dz_lay + 0.5*(dz(i,j,k+1)+dz(i,j+1,k+1)) ))
1672 v(i,j,k) = v(i,j,k) + dt * (dtauup-dtaudn) / h_lay
1673 enddo
1674 enddo
1675 enddo
1676
1677end 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.
1685subroutine 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 do k = 1, gv%ke
1705 do j = g%jsc, g%jec
1706 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.25*((waves%us_y(i,j,k)+waves%us_y(i+1,j,k)) * g%CoriolisBu(i,j))
1709 u(i,j,k) = u(i,j,k) + dvel*dt
1710 enddo
1711 enddo
1712 enddo
1713
1714 do k = 1, gv%ke
1715 do j = g%jscB, g%jecB
1716 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.25*((waves%us_x(i,j,k)+waves%us_x(i,j+1,k)) * g%CoriolisBu(i,j))
1719 v(i,j,k) = v(i,j,k) - dvel*dt
1720 enddo
1721 enddo
1722 enddo
1723end 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
1729subroutine 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 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 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 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 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 pfu_stokes(:,:,:) = 0.0
1792 pfv_stokes(:,:,:) = 0.0
1793 if (cs%id_P_deltaStokes_i > 0) p_deltastokes_i(:,:,:) = 0.0
1794 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 do j = g%jsc, g%jec ; do i = g%iscB, g%iecB
1804 if (g%mask2dCu(i,j)>0.5) then
1805 p_stokes_l0 = 0.0
1806 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 zi_l(1) = 0.0
1811 zi_r(1) = 0.0
1812 do k = 1, g%ke
1813 h_l = dz(i,j,k)
1814 h_r = dz(i+1,j,k)
1815 zi_l(k+1) = zi_l(k) - h_l
1816 zi_r(k+1) = zi_r(k) - h_r
1817 if (.not.cs%robust_Stokes_PGF) then
1818 ! When the code is properly refactored, the following hard-coded constants are unnecessary.
1819 idz_l(k) = 1./max(0.1*us%m_to_Z, h_l)
1820 idz_r(k) = 1./max(0.1*us%m_to_Z, h_r)
1821 endif
1822 enddo
1823 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 (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 (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 (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 (v(i+1,j,k)-cs%Us_y(i+1,j,k))*g%mask2dCv(i+1,j))
1834
1835 dp_stokes_l_dz = 0.0
1836 dp_stokes_r_dz = 0.0
1837 dp_stokes_l = 0.0
1838 dp_stokes_r = 0.0
1839 dp_lay_stokes_l=0.0
1840 dp_lay_stokes_r=0.0
1841
1842 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 cs%Stkx0(i,j,l)*g%mask2dCu(i,j))
1847 us0_r = 0.5*(cs%Stkx0(i,j,l)*g%mask2dCu(i,j) + &
1848 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 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 cs%Stky0(i+1,j,l)*g%mask2dCv(i+1,j))
1853
1854 ! Wavenumber terms that are useful to simplify the pressure calculations
1855 twok = 2.*cs%WaveNum_Cen(l)
1856 fourk = 2.*twok
1857 if (.not.cs%robust_Stokes_PGF) then
1858 itwok = 1. / twok
1859 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 if (g%mask2dT(i,j)>0.5) then
1865 if (.not.cs%robust_Stokes_PGF) then
1866 dexp2kzl = exp(twok*zi_l(k))-exp(twok*zi_l(k+1))
1867 dexp4kzl = exp(fourk*zi_l(k))-exp(fourk*zi_l(k+1))
1868 dp_stokes_l_dz = dp_stokes_l_dz + &
1869 ((ue_l*us0_l+ve_l*vs0_l)*itwok*dexp2kzl + 0.5*(us0_l*us0_l+vs0_l*vs0_l)*ifourk*dexp4kzl)
1870 dp_stokes_l = dp_stokes_l + (ue_l*us0_l+ve_l*vs0_l)*dexp2kzl + 0.5*(us0_l*us0_l+vs0_l*vs0_l)*dexp4kzl
1871 else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers.
1872 exp_top = exp(twok*zi_l(k))
1873 dp_lay_stokes_l = dp_lay_stokes_l + &
1874 ((((ue_l*us0_l)+(ve_l*vs0_l)) * exp_top) * one_minus_exp_x(twok*dz(i,j,k)) + &
1875 (0.5*((us0_l**2)+(vs0_l**2)) * exp_top**2) * one_minus_exp_x(fourk*dz(i,j,k)) )
1876 dp_stokes_l = dp_stokes_l + &
1877 ((((ue_l*us0_l)+(ve_l*vs0_l)) * exp_top) * one_minus_exp(twok*dz(i,j,k)) + &
1878 (0.5*((us0_l**2)+(vs0_l**2)) * exp_top**2) * one_minus_exp(fourk*dz(i,j,k)) )
1879 endif
1880 endif
1881 if (g%mask2dT(i+1,j)>0.5) then
1882 if (.not.cs%robust_Stokes_PGF) then
1883 dexp2kzr = exp(twok*zi_r(k))-exp(twok*zi_r(k+1))
1884 dexp4kzr = exp(fourk*zi_r(k))-exp(fourk*zi_r(k+1))
1885 dp_stokes_r_dz = dp_stokes_r_dz + &
1886 ((ue_r*us0_r+ve_r*vs0_r)*itwok*dexp2kzr + 0.5*(us0_r*us0_r+vs0_r*vs0_r)*ifourk*dexp4kzr)
1887 dp_stokes_r = dp_stokes_r + (ue_r*us0_r+ve_r*vs0_r)*dexp2kzr + 0.5*(us0_r*us0_r+vs0_r*vs0_r)*dexp4kzr
1888 else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers.
1889 exp_top = exp(twok*zi_r(k))
1890 dp_lay_stokes_r = dp_lay_stokes_r + &
1891 ((((ue_r*us0_r)+(ve_r*vs0_r)) * exp_top) * one_minus_exp_x(twok*dz(i+1,j,k)) + &
1892 (0.5*((us0_r**2)+(vs0_r**2)) * exp_top**2) * one_minus_exp_x(fourk*dz(i+1,j,k)) )
1893 dp_stokes_r = dp_stokes_r + &
1894 ((((ue_r*us0_r)+(ve_r*vs0_r)) * exp_top) * one_minus_exp(twok*dz(i+1,j,k)) + &
1895 (0.5*((us0_r**2)+(vs0_r**2)) * 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 if (.not.cs%robust_Stokes_PGF) then
1903 p_stokes_l = p_stokes_l0 + dp_stokes_l_dz*idz_l(k)
1904 p_stokes_r = p_stokes_r0 + dp_stokes_r_dz*idz_r(k)
1905 else
1906 p_stokes_l = p_stokes_l0 + dp_lay_stokes_l
1907 p_stokes_r = p_stokes_r0 + dp_lay_stokes_r
1908 endif
1909
1910 ! > Increment the Interface pressure
1911 p_stokes_l0 = p_stokes_l0 + dp_stokes_l
1912 p_stokes_r0 = p_stokes_r0 + dp_stokes_r
1913
1914 ! Pressure force anomaly is finite difference across the cell.
1915 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 if (g%mask2dT(i,j)>0.5 .and. cs%id_P_deltaStokes_L > 0) p_deltastokes_l(i,j,k) = p_stokes_l
1919 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 do j = g%jscB, g%jecB ; do i = g%isc, g%iec
1933 if (g%mask2dCv(i,j)>0.5) then
1934 p_stokes_l0 = 0.0
1935 p_stokes_r0 = 0.0
1936 zi_l(1) = 0.0
1937 zi_r(1) = 0.0
1938 do k = 1, g%ke
1939 h_l = dz(i,j,k)
1940 h_r = dz(i,j+1,k)
1941 zi_l(k+1) = zi_l(k) - h_l
1942 zi_r(k+1) = zi_r(k) - h_r
1943 if (.not.cs%robust_Stokes_PGF) then
1944 ! When the code is properly refactored, the following hard-coded constants are unnecessary.
1945 idz_l(k) = 1. / max(0.1*us%m_to_Z, h_l)
1946 idz_r(k) = 1. / max(0.1*us%m_to_Z, h_r)
1947 endif
1948 enddo
1949 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 (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 (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 (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 (v(i,j+1,k)-cs%Us_y(i,j+1,k))*g%mask2dCv(i,j+1))
1960
1961 dp_stokes_l_dz = 0.0
1962 dp_stokes_r_dz = 0.0
1963 dp_stokes_l = 0.0
1964 dp_stokes_r = 0.0
1965 dp_lay_stokes_l=0.0
1966 dp_lay_stokes_r=0.0
1967
1968 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 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 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 cs%Stky0(i,j,l)*g%mask2dCv(i,j))
1977 vs0_r = 0.5*(cs%Stky0(i,j,l)*g%mask2dCv(i,j) + &
1978 cs%Stky0(i,j+1,l)*g%mask2dCv(i,j+1))
1979
1980 ! Wavenumber terms that are useful to simplify the pressure calculations
1981 twok = 2.*cs%WaveNum_Cen(l)
1982 fourk = 2.*twok
1983 if (.not.cs%robust_Stokes_PGF) then
1984 itwok = 1. / twok
1985 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 if (g%mask2dT(i,j)>0.5) then
1991 if (.not.cs%robust_Stokes_PGF) then
1992 dexp2kzl = exp(twok*zi_l(k))-exp(twok*zi_l(k+1))
1993 dexp4kzl = exp(fourk*zi_l(k))-exp(fourk*zi_l(k+1))
1994 dp_stokes_l_dz = dp_stokes_l_dz + &
1995 ((ue_l*us0_l+ve_l*vs0_l)*itwok*dexp2kzl + 0.5*(us0_l*us0_l+vs0_l*vs0_l)*ifourk*dexp4kzl)
1996 dp_stokes_l = dp_stokes_l + (ue_l*us0_l+ve_l*vs0_l)*dexp2kzl + 0.5*(us0_l*us0_l+vs0_l*vs0_l)*dexp4kzl
1997 else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers.
1998 exp_top = exp(twok*zi_l(k))
1999 dp_lay_stokes_l = dp_lay_stokes_l + &
2000 ((((ue_l*us0_l)+(ve_l*vs0_l)) * exp_top) * one_minus_exp_x(twok*dz(i,j,k)) + &
2001 (0.5*((us0_l**2)+(vs0_l**2)) * exp_top**2) * one_minus_exp_x(fourk*dz(i,j,k)) )
2002 dp_stokes_l = dp_stokes_l + &
2003 ((((ue_l*us0_l)+(ve_l*vs0_l)) * exp_top) * one_minus_exp(twok*dz(i,j,k)) + &
2004 (0.5*((us0_l**2)+(vs0_l**2)) * exp_top**2) * one_minus_exp(fourk*dz(i,j,k)) )
2005 endif
2006 endif
2007 if (g%mask2dT(i,j+1)>0.5) then
2008 if (.not.cs%robust_Stokes_PGF) then
2009 dexp2kzr = exp(twok*zi_r(k))-exp(twok*zi_r(k+1))
2010 dexp4kzr = exp(fourk*zi_r(k))-exp(fourk*zi_r(k+1))
2011 dp_stokes_r_dz = dp_stokes_r_dz + &
2012 ((ue_r*us0_r+ve_r*vs0_r)*itwok*dexp2kzr + 0.5*(us0_r*us0_r+vs0_r*vs0_r)*ifourk*dexp4kzr)
2013 dp_stokes_r = dp_stokes_r + (ue_r*us0_r+ve_r*vs0_r)*dexp2kzr + 0.5*(us0_r*us0_r+vs0_r*vs0_r)*dexp4kzr
2014 else ! These expressions are equivalent to those above for thick layers, but more accurate for thin layers.
2015 exp_top = exp(twok*zi_r(k))
2016 dp_lay_stokes_r = dp_lay_stokes_r + &
2017 ((((ue_r*us0_r)+(ve_r*vs0_r)) * exp_top) * one_minus_exp_x(twok*dz(i,j+1,k)) + &
2018 (0.5*((us0_r**2)+(vs0_r**2)) * exp_top**2) * one_minus_exp_x(fourk*dz(i,j+1,k)) )
2019 dp_stokes_r = dp_stokes_r + &
2020 ((((ue_r*us0_r)+(ve_r*vs0_r)) * exp_top) * one_minus_exp(twok*dz(i,j+1,k)) + &
2021 (0.5*((us0_r**2)+(vs0_r**2)) * 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 if (.not.cs%robust_Stokes_PGF) then
2029 p_stokes_l = p_stokes_l0 + dp_stokes_l_dz*idz_l(k)
2030 p_stokes_r = p_stokes_r0 + dp_stokes_r_dz*idz_r(k)
2031 else
2032 p_stokes_l = p_stokes_l0 + dp_lay_stokes_l
2033 p_stokes_r = p_stokes_r0 + dp_lay_stokes_r
2034 endif
2035
2036 ! > Increment the Interface pressure
2037 p_stokes_l0 = p_stokes_l0 + dp_stokes_l
2038 p_stokes_r0 = p_stokes_r0 + dp_stokes_r
2039
2040 ! Pressure force anomaly is finite difference across the cell.
2041 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 if (cs%id_PFv_Stokes>0) &
2048 call post_data(cs%id_PFv_Stokes, pfv_stokes, cs%diag)
2049 if (cs%id_PFu_Stokes>0) &
2050 call post_data(cs%id_PFu_Stokes, pfu_stokes, cs%diag)
2051 if (cs%id_P_deltaStokes_L>0) &
2052 call post_data(cs%id_P_deltaStokes_L, p_deltastokes_l, cs%diag)
2053 if (cs%id_P_deltaStokes_i>0) &
2054 call post_data(cs%id_P_deltaStokes_i, p_deltastokes_i, cs%diag)
2055
2056end 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.
2062subroutine 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 if (cs%vonKar < 0.0) call mom_error(fatal, &
2088 "ust_2_u10_coare3p5 called with a negative value of Waves%vonKar")
2089
2090 z0sm = 0.11 * cs%nu_air / ustair ! Compute z0smooth from ustar guess
2091 u10a = 1000.0*us%m_s_to_L_T ! An insanely large upper bound for u10.
2092
2093 if (cs%answer_date < 20230103) then
2094 u10 = us%Z_to_L*ustair / sqrt(0.001) ! Guess for u10
2095 ten_m_scale = 10.0*us%m_to_Z
2096 ct=0
2097 do while (abs(u10a/u10 - 1.) > 0.001)
2098 ct=ct+1
2099 u10a = u10
2100 alpha = min(cs%Charnock_min, cs%Charnock_slope_U10 * u10 + cs%Charnock_intercept)
2101 z0rough = alpha * (us%Z_to_L*ustair)**2 / gv%g_Earth ! Compute z0rough from ustar guess
2102 z0 = z0sm + z0rough
2103 cd = ( cs%vonKar / log(ten_m_scale / z0) )**2 ! Compute Cd from derived roughness
2104 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 if (ct>20) then
2111 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 exit
2114 endif
2115 enddo
2116
2117 else ! Use more efficient expressions that are mathematically equivalent to those above.
2118 u10 = us%Z_to_L*ustair * 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 i_vonkar = 1.0 / cs%vonKar
2121 i_ten_m_scale = 0.1*us%Z_to_m
2122
2123 do ct=1,20
2124 if (abs(u10a - u10) <= 0.001*u10) exit ! Check for convergence.
2125 u10a = u10
2126 alpha = min(cs%Charnock_min, cs%Charnock_slope_U10 * u10 + cs%Charnock_intercept)
2127 z0rough = alpha * (cs%I_g_Earth * ustair**2) ! Compute z0rough from ustar guess
2128 z0 = z0sm + z0rough
2129 i_sqrtcd = abs(log(z0 * i_ten_m_scale)) * i_vonkar ! Compute Cd from derived roughness
2130 u10 = us%Z_to_L*ustair * 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 if (abs(u10a - u10) > 0.001*u10) u10 = us%Z_to_L*ustair * 25.82
2136 endif
2137
2138end subroutine ust_2_u10_coare3p5
2139
2140!> Clear pointers, deallocate memory
2141subroutine waves_end(CS)
2142 type(wave_parameters_cs), pointer :: cs !< Control structure
2143
2144 if (allocated(cs%WaveNum_Cen)) deallocate( cs%WaveNum_Cen )
2145 if (allocated(cs%Freq_Cen)) deallocate( cs%Freq_Cen )
2146 if (allocated(cs%Us_x)) deallocate( cs%Us_x )
2147 if (allocated(cs%Us_y)) deallocate( cs%Us_y )
2148 if (allocated(cs%La_turb)) deallocate( cs%La_turb )
2149 if (allocated(cs%STKx0)) deallocate( cs%STKx0 )
2150 if (allocated(cs%STKy0)) deallocate( cs%STKy0 )
2151 if (allocated(cs%UStk_Hb)) deallocate( cs%UStk_Hb )
2152 if (allocated(cs%VStk_Hb)) deallocate( cs%VStk_Hb )
2153 if (allocated(cs%Omega_w2x)) deallocate( cs%Omega_w2x )
2154 if (allocated(cs%KvS)) deallocate( cs%KvS )
2155 if (allocated(cs%Us0_y)) deallocate( cs%Us0_y )
2156 if (allocated(cs%Us0_x)) deallocate( cs%Us0_x )
2157
2158 deallocate( cs )
2159
2160end subroutine waves_end
2161
2162!> Register wave restart fields. To be called before MOM_wave_interface_init
2163subroutine 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 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 if (associated(cs)) then
2178 call mom_error(fatal, "waves_register_restarts: Called with initialized waves control structure")
2179 endif
2180 allocate(cs)
2181
2182 call get_param(param_file, mdl, "USE_WAVES", use_waves, &
2183 "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 do_not_log=.true.,default=.false.)
2188
2189 if (.not.(use_waves .or. statisticalwaves)) return
2190
2191 call get_param(param_file, mdl, "STOKES_DDT", time_tendency_term, do_not_log=.true., default=.false.)
2192
2193 if (time_tendency_term) then
2194 ! Allocate wave fields needed for restart file
2195 allocate(cs%Us_x_prev(hi%isdB:hi%IedB,hi%jsd:hi%jed,gv%ke), source=0.0)
2196 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 hor_grid='u', z_grid='L')
2201 vd(2) = var_desc("Us_y_prev", "m s-1", "3d meridional Stokes drift profile",&
2202 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 restart_csp, conversion=us%L_T_to_m_s)
2205 endif
2206
2207end 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
2236end module mom_wave_interface