Refactor vertical vorticity to return operation rather than Field (#3340)

* Refactor vertical vorticity to return operation rather than Field

* Update src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity.jl

Co-authored-by: Navid C. Constantinou <navidcy@users.noreply.github.com>

* vertical_vorticity_field.jl -> vertical_vorticity.jl

* Update test

* Bugfix

---------

Co-authored-by: Navid C. Constantinou <navidcy@users.noreply.github.com>

Author: glwagner

src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl

@@ -151,6 +151,6 @@ include("slice_ensemble_model_mode.jl")
 ##### Some diagnostics
 #####
 
-include("vertical_vorticity_field.jl")
+include("vertical_vorticity.jl")
 
 end # module

src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity.jl

@@ -0,0 +1,17 @@
+using Oceananigans.Grids: Face, Face, Center
+using Oceananigans.Operators: ζ₃ᶠᶠᶜ
+using Oceananigans.AbstractOperations: KernelFunctionOperation
+
+"""
+    vertical_vorticity(model::HydrostaticFreeSurfaceModel)
+
+Return a `KernelFunctionOperation` that represents vertical vorticity for the
+`HydrostaticFreeSurfaceModel`. The kernel function is `Oceananigans.Operators.ζ₃ᶠᶠᶜ`,
+and thus computed consistently with the `VectorInvariant` momentum advection scheme
+for curvilinear grids.
+"""
+function vertical_vorticity(model::HydrostaticFreeSurfaceModel)
+    u, v, w = model.velocities
+    return KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᶜ, model.grid, u, v)
+end
+

src/Models/HydrostaticFreeSurfaceModels/vertical_vorticity_field.jl

@@ -1,4 +1,4 @@
-using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField, VectorInvariant
+using Oceananigans.Models.HydrostaticFreeSurfaceModels: vertical_vorticity, VectorInvariant
 
 @testset "VerticalVorticityField with HydrostaticFreeSurfaceModel" begin
 
@@ -16,10 +16,10 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField,
             ψᵢ(λ, φ, z) = rand()
             set!(model, u=ψᵢ, v=ψᵢ)
 
-            ζ = VerticalVorticityField(model)
-
+            ζ = vertical_vorticity(model)
+            @test ζ isa KernelFunctionOperation
+            ζ = Field(ζ)
             compute!(ζ)
-
             @test all(isfinite.(ζ.data))
         end
     end