Improve notation for computing turbulent fluxes (#99)

* Improve notation for computing turbulent fluxes

* Minor cleanup

* Add some comments

* Another comment

* Update src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl

Co-authored-by: Simone Silvestri <[email protected]>

* Resolve dependencies

* Add using ClimaOcean to runtests setup

* Fix bug in extending -

---------

Co-authored-by: Simone Silvestri <[email protected]>

Manifest.toml

@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.4"
+julia_version = "1.10.0"
 manifest_format = "2.0"
 project_hash = "086b5c52c3f2a61f133762647db63128f8714c56"
 
@@ -255,7 +255,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.1.1+0"
+version = "1.0.5+1"
 
 [[deps.CompositionsBase]]
 git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
@@ -906,7 +906,7 @@ weakdeps = ["Adapt"]
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+4"
+version = "0.3.23+2"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]

src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl

@@ -58,7 +58,8 @@ function compute_atmosphere_ocean_fluxes!(coupled_model)
     # kernel parameters that compute fluxes in 0:Nx+1 and 0:Ny+1
     kernel_parameters = KernelParameters(kernel_size, (-1, -1))
 
-    launch!(arch, grid, kernel_parameters, _compute_atmosphere_ocean_similarity_theory_fluxes!,
+    launch!(arch, grid, kernel_parameters,
+            _compute_atmosphere_ocean_similarity_theory_fluxes!,
             similarity_theory,
             grid,
             clock,
@@ -73,7 +74,8 @@ function compute_atmosphere_ocean_fluxes!(coupled_model)
             atmosphere.boundary_layer_height,
             atmosphere.thermodynamics_parameters)   
 
-    launch!(arch, grid, kernel_parameters, _assemble_atmosphere_ocean_fluxes!,
+    launch!(arch, grid, kernel_parameters,
+            _assemble_atmosphere_ocean_fluxes!,
             centered_velocity_fluxes,
             net_tracer_fluxes,
             grid,
@@ -182,19 +184,13 @@ limit_fluxes_over_sea_ice!(args...) = nothing
                                                similarity_theory.water_vapor_saturation,
                                                surface_type)
 
-    # Thermodynamic and dynamic surface state
+    # Thermodynamic and dynamic (ocean) surface state
     𝒬₀ = thermodynamic_surface_state = AtmosphericThermodynamics.PhaseEquil_pTq(ℂₐ, pₐ, Tₒ, qₒ)
-
     h₀ = zero(grid) # surface height
     Uₒ = SVector(uₒ, vₒ)
     𝒰₀ = dynamic_ocean_state = SurfaceFluxes.StateValues(h₀, Uₒ, 𝒬₀)
 
-    Qv = similarity_theory.fields.latent_heat
-    Qc = similarity_theory.fields.sensible_heat
-    Fv = similarity_theory.fields.water_vapor
-    τx = similarity_theory.fields.x_momentum
-    τy = similarity_theory.fields.y_momentum
-
+    # Some parameters
     g = default_gravitational_acceleration
     ϰ = similarity_theory.von_karman_constant
 
@@ -209,16 +205,24 @@ limit_fluxes_over_sea_ice!(args...) = nothing
 
     # Convert back from a zonal - meridional flux to the frame of 
     # reference of the native ocean grid
-    τˣ, τʸ = convert_to_extrinsic_reference_frame(i, j, kᴺ, grid, turbulent_fluxes.x_momentum, 
-                                                                  turbulent_fluxes.y_momentum)
+    ρτxⁱʲ, ρτyⁱʲ = convert_to_extrinsic_reference_frame(i, j, kᴺ, grid,
+                                                        turbulent_fluxes.x_momentum, 
+                                                        turbulent_fluxes.y_momentum)
+
+    # Store fluxes
+    Qv = similarity_theory.fields.latent_heat
+    Qc = similarity_theory.fields.sensible_heat
+    Fv = similarity_theory.fields.water_vapor
+    ρτx = similarity_theory.fields.x_momentum
+    ρτy = similarity_theory.fields.y_momentum
 
     @inbounds begin
         # +0: cooling, -0: heating
-        Qv[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.latent_heat)
-        Qc[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.sensible_heat)
-        Fv[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.water_vapor)
-        τx[i, j, 1] = ifelse(inactive, 0, τˣ)
-        τy[i, j, 1] = ifelse(inactive, 0, τʸ)
+        Qv[i, j, 1]  = ifelse(inactive, 0, turbulent_fluxes.latent_heat)
+        Qc[i, j, 1]  = ifelse(inactive, 0, turbulent_fluxes.sensible_heat)
+        Fv[i, j, 1]  = ifelse(inactive, 0, turbulent_fluxes.water_vapor)
+        ρτx[i, j, 1] = ifelse(inactive, 0, ρτxⁱʲ)
+        ρτy[i, j, 1] = ifelse(inactive, 0, ρτyⁱʲ)
     end
 end
 
@@ -262,11 +266,11 @@ end
         Mp = interp_atmos_time_series(prescribed_freshwater_flux, X, time, atmos_args...)
         Mr = get_runoff_flux(X, time, runoff_args) 
 
-        Qc = similarity_theory_fields.sensible_heat[i, j, 1] # sensible or "conductive" heat flux
-        Qv = similarity_theory_fields.latent_heat[i, j, 1]   # latent heat flux
-        Mv = similarity_theory_fields.water_vapor[i, j, 1]   # mass flux of water vapor
-        τx = similarity_theory_fields.x_momentum[i, j, 1]    # zonal momentum flux
-        τy = similarity_theory_fields.y_momentum[i, j, 1]    # meridional momentum flux
+        Qc  = similarity_theory_fields.sensible_heat[i, j, 1] # sensible or "conductive" heat flux
+        Qv  = similarity_theory_fields.latent_heat[i, j, 1]   # latent heat flux
+        Mv  = similarity_theory_fields.water_vapor[i, j, 1]   # mass flux of water vapor
+        ρτx = similarity_theory_fields.x_momentum[i, j, 1]    # zonal momentum flux
+        ρτy = similarity_theory_fields.y_momentum[i, j, 1]    # meridional momentum flux
     end
 
     # Compute heat fluxes, bulk flux first
@@ -286,25 +290,25 @@ end
     ΣF += Fv
 
     # Compute fluxes for u, v, T, S from momentum, heat, and freshwater fluxes
-    Jᵘ = centered_velocity_fluxes.u
-    Jᵛ = centered_velocity_fluxes.v
+    τx = centered_velocity_fluxes.u
+    τy = centered_velocity_fluxes.v
     Jᵀ = net_tracer_fluxes.T
     Jˢ = net_tracer_fluxes.S
 
     ρₒ = ocean_reference_density
     cₒ = ocean_heat_capacity
 
-    atmos_ocean_Jᵘ = τx / ρₒ
-    atmos_ocean_Jᵛ = τy / ρₒ
+    atmos_ocean_τx = ρτx / ρₒ
+    atmos_ocean_τy = ρτy / ρₒ
     atmos_ocean_Jᵀ = ΣQ / (ρₒ * cₒ)
     atmos_ocean_Jˢ = - Sₒ * ΣF
 
     # Mask fluxes over land for convenience
     inactive = inactive_node(i, j, kᴺ, grid, c, c, c)
 
     @inbounds begin
-        Jᵘ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jᵘ)
-        Jᵛ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jᵛ)
+        τx[i, j, 1] = ifelse(inactive, 0, atmos_ocean_τx)
+        τy[i, j, 1] = ifelse(inactive, 0, atmos_ocean_τy)
         Jᵀ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jᵀ)
         Jˢ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jˢ)
     end

src/OceanSeaIceModels/CrossRealmFluxes/ocean_sea_ice_surface_fluxes.jl

@@ -29,6 +29,8 @@ using KernelAbstractions: @kernel, @index
 struct OceanSeaIceSurfaceFluxes{T, P, C, R, PI, PC, FT, UN}
     turbulent :: T
     prescribed :: P
+    # Add `components` which will also store components of the total fluxes
+    # (eg latent, sensible heat flux)
     total :: C
     radiation :: R
     previous_ice_thickness :: PI

src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl

@@ -30,20 +30,20 @@ import SurfaceFluxes.Parameters:
 #####
 
 struct SimilarityTheoryTurbulentFluxes{FT, UF, TP, S, W, R, B, V, F}
-    gravitational_acceleration :: FT
-    von_karman_constant :: FT
-    turbulent_prandtl_number :: FT
-    gustiness_parameter :: FT
-    stability_functions :: UF
-    thermodynamics_parameters :: TP
-    water_vapor_saturation :: S
-    water_mole_fraction :: W
-    roughness_lengths :: R
-    bulk_coefficients :: B
-    bulk_velocity :: V
-    tolerance :: FT
-    maxiter :: Int
-    fields :: F
+    gravitational_acceleration :: FT # parameter
+    von_karman_constant :: FT        # parameter
+    turbulent_prandtl_number :: FT   # parameter
+    gustiness_parameter :: FT        # bulk velocity parameter
+    stability_functions :: UF        # functions for turbulent fluxes
+    thermodynamics_parameters :: TP  # parameter group
+    water_vapor_saturation :: S      # model for computing the saturation water vapor mass
+    water_mole_fraction :: W         # mole fraction of H₂O in seawater
+    roughness_lengths :: R           # parameterization for turbulent fluxes
+    bulk_coefficients :: B           # ?
+    bulk_velocity :: V               # bulk velocity scale for turbulent fluxes
+    tolerance :: FT                  # solver option
+    maxiter :: Int                   # solver option
+    fields :: F                      # fields that store turbulent fluxes
 end
 
 const STTF = SimilarityTheoryTurbulentFluxes
@@ -229,23 +229,27 @@ end
     Σ★ = SimilarityScales(u★, u★, u★) 
 
     # The inital velocity scale assumes that
-    # the gustiness velocity `uᴳ` is equal to 0.5 ms⁻¹. 
+    # the gustiness velocity `Uᴳ` is equal to 0.5 ms⁻¹. 
     # That will be refined later on.
-    ΔUᴳ = sqrt(Δu^2 + Δv^2 + convert(eltype(Δh), 0.25))
+    FT = eltype(Δh)
+    Uᴳᵢ = convert(FT, 0.5^2)
+    ΔU = sqrt(Δu^2 + Δv^2 + Uᴳᵢ)
 
     # Initialize the solver
     iteration = 0
 
     while iterating(Σ★ - Σ₀, iteration, maxiter, similarity_theory)
         Σ₀ = Σ★
-        Σ★, ΔUᴳ = refine_characteristic_scales(Σ★, ΔUᴳ, 
-                                               similarity_theory,
-                                               surface_state,
-                                               differences,
-                                               atmos_boundary_layer_height,
-                                               thermodynamics_parameters,
-                                               gravitational_acceleration,
-                                               von_karman_constant)
+        # Refine both the characteristic scale and the effective
+        # velocity difference ΔU, including gustiness.
+        Σ★, ΔU = refine_similarity_variables(Σ★, ΔU, 
+                                             similarity_theory,
+                                             surface_state,
+                                             differences,
+                                             atmos_boundary_layer_height,
+                                             thermodynamics_parameters,
+                                             gravitational_acceleration,
+                                             von_karman_constant)
         iteration += 1
     end
 
@@ -257,9 +261,9 @@ end
     q★ = q★ / similarity_theory.turbulent_prandtl_number
 
     # `u★² ≡ sqrt(τx² + τy²)`
-    # We remove the gustiness by dividing by `ΔUᴳ`
-    τx = - u★^2 * Δu / ΔUᴳ
-    τy = - u★^2 * Δv / ΔUᴳ
+    # We remove the gustiness by dividing by `ΔU`
+    τx = - u★^2 * Δu / ΔU
+    τy = - u★^2 * Δv / ΔU
 
     𝒬ₐ = atmos_state.ts
     ρₐ = AtmosphericThermodynamics.air_density(ℂₐ, 𝒬ₐ)
@@ -326,15 +330,15 @@ end
     return Δh, Δu, Δv, Δθ, Δq
 end
 
-@inline function refine_characteristic_scales(estimated_characteristic_scales, 
-                                              velocity_scale,
-                                              similarity_theory,
-                                              surface_state,
-                                              differences,
-                                              atmos_boundary_layer_height,
-                                              thermodynamics_parameters,
-                                              gravitational_acceleration,
-                                              von_karman_constant)
+@inline function refine_similarity_variables(estimated_characteristic_scales, 
+                                             velocity_scale,
+                                             similarity_theory,
+                                             surface_state,
+                                             differences,
+                                             atmos_boundary_layer_height,
+                                             thermodynamics_parameters,
+                                             gravitational_acceleration,
+                                             von_karman_constant)
 
     # "initial" scales because we will recompute them
     u★ = estimated_characteristic_scales.momentum
@@ -358,7 +362,7 @@ end
     ℂ  = thermodynamics_parameters
     g  = gravitational_acceleration
     𝒬ₒ = surface_state.ts # thermodynamic state
-    zᵢ = atmos_boundary_layer_height
+    hᵢ = atmos_boundary_layer_height
 
     # Compute Monin-Obukhov length scale depending on a `buoyancy flux`
     b★ = buoyancy_scale(θ★, q★, 𝒬ₒ, ℂ, g)
@@ -388,10 +392,11 @@ end
 
     # Buoyancy flux characteristic scale for gustiness (Edson 2013)
     Jᵇ = - u★ * b★
-    uᴳ = β * cbrt(Jᵇ * zᵢ)
+    Uᴳ = β * cbrt(Jᵇ * hᵢ)
 
     # New velocity difference accounting for gustiness
-    ΔUᴳ = sqrt(Δu^2 + Δv^2 + uᴳ^2)
+    ΔU = sqrt(Δu^2 + Δv^2 + Uᴳ^2)
+
+    return SimilarityScales(u★, θ★, q★), ΔU
+end
 
-    return SimilarityScales(u★, θ★, q★), ΔUᴳ
-end

src/OceanSeaIceModels/CrossRealmFluxes/stability_functions.jl

@@ -1,5 +1,5 @@
-import Base: -
 import Statistics
+import Base: -
 
 #####
 ##### Struct that represents a 3-tuple of momentum, heat, and water vapor
@@ -11,18 +11,19 @@ struct SimilarityScales{U, T, Q}
     water_vapor :: Q
 end
 
--(a::SimilarityScales, b::SimilarityScales) = SimilarityScales(a.momentum - b.momentum, 
-                                                               a.temperature - b.temperature,
-                                                               a.water_vapor - b.water_vapor)
+function -(a::SimilarityScales, b::SimilarityScales)
+    Δu = a.momentum - b.momentum
+    Δθ = a.temperature - b.temperature
+    Δq = a.water_vapor - b.water_vapor
+    return SimilarityScales(Δu, Δθ, Δq)
+end
 
 Statistics.norm(a::SimilarityScales) = norm(a.momentum) + norm(a.temperature) + norm(a.water_vapor)
 
+# Edson et al. (2013)
 function edson_stability_functions(FT = Float64)
-
-    # Edson et al. (2013)
     ψu = MomentumStabilityFunction()
     ψc = ScalarStabilityFunction()
-
     return SimilarityScales(ψu, ψc, ψc)
 end
 

test/runtests_setup.jl

@@ -6,4 +6,6 @@ using Test
 using Oceananigans.Architectures: architecture
 using Oceananigans.OutputReaders: interpolate!
 
+using ClimaOcean
+
 test_architectures = [CPU()]