Merge pull request #10 from simone-silvestri/ss/available-potential-e…

…nergy

Integral Available Potential Energy calculation

src/Diagnostics/integrated_diagnostics.jl

@@ -12,6 +12,33 @@ function integral_kinetic_energy(u::FieldTimeSeries, v::FieldTimeSeries; stride
     return energy
 end
 
+function integral_available_potential_energy(b::FieldTimeSeries; stride = 1, start_time = 1, end_time = length(u.times))
+    energy = Float64[]
+    vol = VolumeField(u.grid)
+
+    for i in start_time:stride:end_time
+        @info "integrating index $i of $end_time"
+        αe = compute_ape_density(b[i])
+        push!(energy, sum(compute!(Field(αe * vol))))
+    end
+
+    return energy
+end
+
+function compute_ape_density(b::Field)
+    ze = calculate_z★_diagnostics(b)
+
+    αe = Field{Center, Center, Center}(ze.grid)
+
+    zfield = HeightField(ze.grid)
+
+    @info "computing resting and available potential energy density..."
+    ρ = DensityOperation(b)
+    set!(αe, (zfield - ze) * ρ)
+
+    return αe
+end
+
 function ACC_transport(u::FieldTimeSeries; stride = 1, start_time = 1, end_time = length(u.times))
 
     transport = Float64[]

src/Diagnostics/spurious_mixing.jl

@@ -20,6 +20,17 @@ function calculate_z★_diagnostics(b::FieldTimeSeries)
     return z★
 end
 
+function calculate_z★_diagnostics(b::Field)
+
+    vol = VolumeField(b.grid)
+    z★  = Field{Center, Center, Center}(b.grid)
+
+    total_area = sum(AreaField(b.grid))
+    calculate_z★!(z★, b, vol, total_area)
+
+    return z★
+end
+
 function calculate_z★!(z★::Field, b::Field, vol, total_area)
     grid = b.grid
     arch = architecture(grid)