Add wind-dependent Charnock coefficient (#520)

* try it out

* wind-dependent charnock

* Update roughness_lengths.jl

* Update roughness_lengths.jl

* Update roughness_lengths.jl

* Update roughness_lengths.jl

* Update roughness_lengths.jl

* Update similarity_theory_turbulent_fluxes.jl

* just a test

* another test

* Update roughness_lengths.jl

* Apply suggestions from code review

Co-authored-by: Gregory L. Wagner <[email protected]>

* docstring defaults mirror the actual defaults

* C -> ℂ

* WindDependentGravityWaveParameter -> WindDependentWaveFormulation

* Update roughness_lengths.jl

* correct wave formulation

* gravity_waves -> wave_formulation

* fixes

* leave wave_formulation = 0.011 as default

* fix docstring

* back to gustiness_parameter = 1

---------

Co-authored-by: Navid C. Constantinou <[email protected]>
Co-authored-by: Gregory L. Wagner <[email protected]>

Author: 3 people

src/OceanSeaIceModels/InterfaceComputations/roughness_lengths.jl

@@ -1,7 +1,7 @@
-struct MomentumRoughnessLength{FT, V}
+struct MomentumRoughnessLength{FT, G, V}
     gravitational_acceleration :: FT
     air_kinematic_viscosity :: V
-    gravity_wave_parameter :: FT
+    wave_formulation :: G
     laminar_parameter :: FT
     maximum_roughness_length :: FT
 end
@@ -12,17 +12,37 @@ struct ScalarRoughnessLength{FT, V, R}
     maximum_roughness_length :: FT
 end
 
+struct WindDependentWaveFormulation{FT}
+    Umax :: FT
+    ℂ₁ :: FT
+    ℂ₂ :: FT
+end
+
+"""
+    WindDependentWaveFormulation(FT = Oceananigans.defaults.FloatType;
+                                 Umax = 19, ℂ₁ = 0.0017, ℂ₂ = -0.005)
+
+A gravity wave parameter based on the wind speed `ΔU` with the formula `ℂ₁ * max(ΔU, Umax) + ℂ₂`.
+"""
+WindDependentWaveFormulation(FT=Oceananigans.defaults.FloatType; Umax = 19, ℂ₁ = 0.0017, ℂ₂ = -0.005) =
+    WindDependentWaveFormulation(convert(FT, Umax),
+                                 convert(FT, ℂ₁),
+                                 convert(FT, ℂ₂))
+
+gravity_wave_parameter(α::Number, args...) = α
+gravity_wave_parameter(α::WindDependentWaveFormulation, ΔU) = α.ℂ₁ * max(ΔU, α.Umax) + α.ℂ₂
+
 """
     ScalarRoughnessLength(FT = Float64;
                           air_kinematic_viscosity = temperature_dependent_viscosity,
                           reynolds_number_scaling_function = empirical_scaling_function,
                           maximum_roughness_length = 1.6e-4)
 
-Constructs a `ScalarRoughnessLength` object that represents the scalar roughness length
+Construct a `ScalarRoughnessLength` object that represents the scalar roughness length
 that regulates the exchange of heat and water vapor between the ocean and the atmosphere.
 
 Keyword Arguments
-==================
+=================
 
 - `air_kinematic_viscosity::Function`: The function to compute the air kinematic viscosity.
 - `reynolds_number_scaling_function::Function`: The function to compute the Reynolds number scaling factor.
@@ -43,7 +63,7 @@ end
                             gravitational_acceleration = default_gravitational_acceleration,
                             maximum_roughness_length = 1.0,
                             air_kinematic_viscosity = TemperatureDependentAirViscosity(FT),
-                            gravity_wave_parameter = 0.011,
+                            wave_formulation = 0.011,
                             laminar_parameter = 0.11)
 
 Construct a `MomentumRoughnessLength` object that represents the momentum roughness length that
@@ -55,21 +75,21 @@ Keyword Arguments
 - `gravitational_acceleration`: The gravitational acceleration. Default: `default_gravitational_acceleration`.
 - `maximum_roughness_length`: The maximum roughness length. Default: 1.0.
 - `air_kinematic_viscosity`: The air kinematic viscosity. Default: `TemperatureDependentAirViscosity(FT)`.
-- `gravity_wave_parameter`: The wave parameter. Default: 0.011.
+- `wave_formulation`: The wave parameter. Either constant or `WindDependentWaveFormulation(FT)`. Default: 0.011.
 - `laminar_parameter`: The laminar parameter. Default: 0.11.
 """
 function MomentumRoughnessLength(FT=Oceananigans.defaults.FloatType;
                                  gravitational_acceleration = default_gravitational_acceleration,
                                  maximum_roughness_length = 1.0, # An estimate?
                                  air_kinematic_viscosity = TemperatureDependentAirViscosity(FT),
-                                 gravity_wave_parameter = 0.011,
+                                 wave_formulation = 0.011,
                                  laminar_parameter = 0.11)
 
     return MomentumRoughnessLength(convert(FT, gravitational_acceleration),
-                                          air_kinematic_viscosity,
-                                          convert(FT, gravity_wave_parameter),
-                                          convert(FT, laminar_parameter),
-                                          convert(FT, maximum_roughness_length))
+                                   air_kinematic_viscosity,
+                                   wave_formulation,
+                                   convert(FT, laminar_parameter),
+                                   convert(FT, maximum_roughness_length))
 end
 
 function default_roughness_lengths(FT=Oceananigans.defaults.FloatType)
@@ -81,42 +101,43 @@ end
 
 # Temperature-dependent viscosity law
 struct TemperatureDependentAirViscosity{FT}
-    C₀ :: FT
-    C₁ :: FT
-    C₂ :: FT
-    C₃ :: FT
+    ℂ₀ :: FT
+    ℂ₁ :: FT
+    ℂ₂ :: FT
+    ℂ₃ :: FT
 end
 
 """
     TemperatureDependentAirViscosity([FT = Oceananigans.defaults.FloatType;
-                                      C₀ = 1.326e-5,
-                                      C₁ = C₀ * 6.542e-3,
-                                      C₂ = C₀ * 8.301e-6,
-                                      C₃ = - C₀ * 4.84e-9])
+                                      ℂ₀ = 1.326e-5,
+                                      ℂ₁ = ℂ₀ * 6.542e-3,
+                                      ℂ₂ = ℂ₀ * 8.301e-6,
+                                      ℂ₃ = - ℂ₀ * 4.84e-9])
 
-Constructs a `TemperatureDependentAirViscosity` object that calculates the kinematic
+Construct a `TemperatureDependentAirViscosity` object that calculates the kinematic
 viscosity of air as
+
 ```math
-C₀ + C₁ T + C₂ T^2 + C₃ T^3.
+ℂ₀ + ℂ₁ T + ℂ₂ T^2 + ℂ₃ T^3
 ```
 """
 function TemperatureDependentAirViscosity(FT = Oceananigans.defaults.FloatType;
-                                          C₀ = 1.326e-5,
-                                          C₁ = C₀ * 6.542e-3,
-                                          C₂ = C₀ * 8.301e-6,
-                                          C₃ = - C₀ * 4.84e-9)
-
-    return TemperatureDependentAirViscosity(convert(FT, C₀),
-                                            convert(FT, C₁),
-                                            convert(FT, C₂),
-                                            convert(FT, C₃))
+                                          ℂ₀ = 1.326e-5,
+                                          ℂ₁ = ℂ₀ * 6.542e-3,
+                                          ℂ₂ = ℂ₀ * 8.301e-6,
+                                          ℂ₃ = - ℂ₀ * 4.84e-9)
+
+    return TemperatureDependentAirViscosity(convert(FT, ℂ₀),
+                                            convert(FT, ℂ₁),
+                                            convert(FT, ℂ₂),
+                                            convert(FT, ℂ₃))
 end
 
 """ Calculate the air viscosity based on the temperature θ in Celsius. """
 @inline function (ν::TemperatureDependentAirViscosity)(θ)
-    FT = eltype(ν.C₀)
+    FT = eltype(ν.ℂ₀)
     T  = convert(FT, θ - celsius_to_kelvin)
-    return ν.C₀ + ν.C₁ * T + ν.C₂ * T^2 + ν.C₃ * T^3
+    return ν.ℂ₀ + ν.ℂ₁ * T + ν.ℂ₂ * T^2 + ν.ℂ₃ * T^3
 end
 
 # Fallbacks for constant roughness length!
@@ -125,11 +146,11 @@ end
 
 # Momentum roughness length should be different from scalar roughness length.
 # Temperature and water vapor can be considered the same (Edson et al 2013)
-@inline function roughness_length(ℓ::MomentumRoughnessLength{FT}, u★, 𝒬, ℂ) where FT
+@inline function roughness_length(ℓ::MomentumRoughnessLength{FT}, ΔU, u★, 𝒬, ℂ) where FT
     g  = ℓ.gravitational_acceleration
-    α  = ℓ.gravity_wave_parameter
     β  = ℓ.laminar_parameter
     ℓm = ℓ.maximum_roughness_length
+    α  = gravity_wave_parameter(ℓ.wave_formulation, ΔU)
 
     θ₀ = AtmosphericThermodynamics.air_temperature(ℂ, 𝒬)
     ν  = ℓ.air_kinematic_viscosity(θ₀)
@@ -178,4 +199,3 @@ ReynoldsScalingFunction(FT = Oceananigans.defaults.FloatType; A = 5.85e-5, b = 0
 
     return min(ℓq, ℓm)
 end
-

src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl

@@ -181,16 +181,6 @@ function iterate_interface_fluxes(flux_formulation::SimilarityTheoryFluxes,
     ϰ = flux_formulation.von_karman_constant
     L★ = ifelse(b★ == 0, Inf, - u★^2 / (ϰ * b★))
 
-    # Compute roughness length scales
-    ℓu₀ = roughness_length(ℓu, u★, 𝒬ₛ, ℂₐ)
-    ℓq₀ = roughness_length(ℓq, ℓu₀, u★, 𝒬ₛ, ℂₐ)
-    ℓθ₀ = roughness_length(ℓθ, ℓu₀, u★, 𝒬ₛ, ℂₐ)
-
-    # Transfer coefficients at height `h`
-    form = flux_formulation.similarity_form
-    χu = ϰ / similarity_profile(form, ψu, Δh, ℓu₀, L★)
-    χθ = ϰ / similarity_profile(form, ψθ, Δh, ℓθ₀, L★)
-    χq = ϰ / similarity_profile(form, ψq, Δh, ℓq₀, L★)
 
     # Buoyancy flux characteristic scale for gustiness (Edson 2013)
     h_bℓ = atmosphere_state.h_bℓ
@@ -201,8 +191,20 @@ function iterate_interface_fluxes(flux_formulation::SimilarityTheoryFluxes,
     Δu, Δv = velocity_difference(interface_properties.velocity_formulation,
                                  atmosphere_state,
                                  approximate_interface_state)
+
     ΔU = sqrt(Δu^2 + Δv^2 + Uᴳ^2)
 
+    # Compute roughness length scales
+    ℓu₀ = roughness_length(ℓu, ΔU,  u★, 𝒬ₛ, ℂₐ)
+    ℓq₀ = roughness_length(ℓq, ℓu₀, u★, 𝒬ₛ, ℂₐ)
+    ℓθ₀ = roughness_length(ℓθ, ℓu₀, u★, 𝒬ₛ, ℂₐ)
+
+    # Transfer coefficients at height `h`
+    form = flux_formulation.similarity_form
+    χu = ϰ / similarity_profile(form, ψu, Δh, ℓu₀, L★)
+    χθ = ϰ / similarity_profile(form, ψθ, Δh, ℓθ₀, L★)
+    χq = ϰ / similarity_profile(form, ψq, Δh, ℓq₀, L★)
+
     # Recompute
     u★ = χu * ΔU
     θ★ = χθ * Δθ