Add meridional heat transport diagnostic via ocean heat uptake method

[taimoorsohail]

Added a Meridional Heat Transport via computing the cumulative integral of the difference between the vertical, longitudinally-integrated OHC and the longitudinally-integrated surface heat flux, both in Joules.

(edited by @navidcy)

Closes #90

[codecov]

Codecov Report

❌ Patch coverage is 0% with 20 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/Diagnostics/meridional_heat_transport.jl	0.00%	20 Missing ⚠️

📢 Thoughts on this report? Let us know!

[navidcy]

Suggested change

	export MixedLayerDepthField, MixedLayerDepthOperand, Meridional_Heat_Transport, reset_meridional_heat_transport_state!
	export MixedLayerDepthField, MixedLayerDepthOperand, meridional_heat_transport, reset_meridional_heat_transport_state!

We use CamelCase for types (or aliases of them) and this_and_that format for methods.
Never combination of camel case + _, i.e., Camel_Case

[taimoorsohail]

Thanks, have fed into the AGENTS.md file so these minor formatting issues aren't repeated

[navidcy]

Did you add this in the AGENTS.md @taimoorsohail?

[taimoorsohail]

Why is there a global heat budget imbalance? Shouldn't the blue line integrate to zero? You are running the RYF correct? I think we need to make sure the diagnostics match at least globally before we merge?

[taimoorsohail]

What does this look like in usage?

Same as the other diagnostics; see, e.g.,

NumericalEarth.jl/examples/meridional_heat_transport_ecco.jl

Lines 81 to 87 in 8f1340f

	mht_OHC= Field(meridional_heat_transport(esm, OceanHeatContentTendencyMethod()))
	mht_vT = Field(meridional_heat_transport(esm, MeridionalHeatFluxMethod()))
	ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht_OHC, mht_vT);
	schedule = TimeInterval(3hours),
	filename = "ocean_one_degree_mht",
	overwrite_existing = true)

So would I have to regrid my tripolar grid beforehand? How would that interface with the esm input variable?

I think we should see first the interface of the regridder and then see?

But perhaps something like:

function meridional_heat_transport(esm::EarthSystemModel, method= MeridionalHeatFluxMethod;

                                   regridder=nothing, reference_temperature=0)



    if regridder !== nothing

        ... # regrid velocities

    else

        return met

    end

end

so we can pass the regridder as an optional kwarg to the diagnostic?

That could be the API for all diagnostics that can optionally include a regridder btw. We can discuss this when discussion in JuliaGeo/ConservativeRegridding.jl#74 converges and we have a method implemented there?

I have been thinking a bit about this.

I actually think something like

esm_regridded = regrid!(esm, Regridder)

That is called by the user prior to diagnostics might be better. Then call the diagnostics as normal with esm_regridded. Anyway let's wait to see what that looks like.

We can merge this PR as-is when we are all happy.

[navidcy]

Why is there a global heat budget imbalance? Shouldn't the blue line integrate to zero?

What does "global heat imbalance" mean? What do you see to say so? If it's that the curves should integrate to zero then does that need to hold for every snapshot? I can compute the integrals and check if they are zero; I can't do the integrals by eye. So it should be that $$\int_{-90^\circ}^{+90^\circ} \mathrm{MHT}(\phi) \mathrm{d}\phi = 0$$? It could be that the blue is closer to zero than the orange.

You are running the RYF correct? I think we need to make sure the diagnostics match at least globally before we merge?

I'm running with a prescribed atmosphere:

NumericalEarth.jl/examples/meridional_heat_transport_ecco.jl

Lines 46 to 47 in 46efae2

	atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(80),
	include_rivers_and_icebergs = false)

I'm also running on a near-global lat-lon grid from 80S to 80N:

NumericalEarth.jl/examples/meridional_heat_transport_ecco.jl

Line 19 in 46efae2

underlying_grid = LatitudeLongitudeGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z, longitude = (0, 360), latitude = (-80, 80))

Could that be the culprit?

I'm wondering: computing the heat transport directly via $v T$ could be wrong for a near-global grid but given that $\int v \mathrm{d}A =0$ by construction we may be not getting any apparent imbalance. But the for surface heat fluxes the missing Arctic matters?

[navidcy]

I have been thinking a bit about this.

I actually think something like

esm_regridded = regrid!(esm, Regridder)

But then we create a whole new esm on a lat-lon grid? And we are timestepping the original esm? How is the esm_regridded moving forward?

Btw, I do think adding the regridding option is out of scope for this PR, right?

[taimoorsohail]

I have been thinking a bit about this.

I actually think something like

esm_regridded = regrid!(esm, Regridder)

But then we create a whole new esm on a lat-lon grid? And we are timestepping the original esm? How is the esm_regridded moving forward?

Btw, I do think adding the regridding option is out of scope for this PR, right?

Yes out of scope for this PR.

The regridded ESM is purely for output writing so would need to be created every time the MHT is output. There wouldn't be any point in time stepping it.

[taimoorsohail]

Why is there a global heat budget imbalance? Shouldn't the blue line integrate to zero?

What does "global heat imbalance" mean? What do you see to say so? If it's that the curves should integrate to zero then does that need to hold for every snapshot? I can compute the integrals and check if they are zero; I can't do the integrals by eye. So it should be that $$\int_{-90^\circ}^{+90^\circ} \mathrm{MHT}(\phi) \mathrm{d}\phi = 0$$? It could be that the blue is closer to zero than the orange.

You are running the RYF correct? I think we need to make sure the diagnostics match at least globally before we merge?

I'm running with a prescribed atmosphere:

NumericalEarth.jl/examples/meridional_heat_transport_ecco.jl

Lines 46 to 47 in 46efae2

	atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(80),
	include_rivers_and_icebergs = false)

I'm also running on a near-global lat-lon grid from 80S to 80N:

NumericalEarth.jl/examples/meridional_heat_transport_ecco.jl

Line 19 in 46efae2

underlying_grid = LatitudeLongitudeGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z, longitude = (0, 360), latitude = (-80, 80))

Could that be the culprit?

I'm wondering: computing the heat transport directly via $v T$ could be wrong for a near-global grid but given that $\int v \mathrm{d}A =0$ by construction we may be not getting any apparent imbalance. But the for surface heat fluxes the missing Arctic matters?

The MHT is a cumulative integral in latitude, right? In that case, the blue line should start and end at 0 if the global heat budget is closed. That means that there is no residual between cumulatively integrated heat flux and OHC tendency. Does that make sense?

[taimoorsohail]

Just pinging @navidcy -- would be great to merge this!

[navidcy]

The MHT is a cumulative integral in latitude, right? In that case, the blue line should start and end at 0 if the global heat budget is closed. That means that there is no residual between cumulatively integrated heat flux and OHC tendency. Does that make sense?

So given that there is a discrepancy, does that imply that the residual term (see the docstring derivation) is not integrated to zero?

[simone-silvestri]

I think our example folder is becoming a bit large. I think we could either add the example to the examples that are run in CI or maybe add a validation folder like in Oceananigans. What do you guys think?

[navidcy]

#159 added the MHT diagnostic via $\int v T dx$.
This PR adds the MHT via the ocean heat update method.

[glwagner]

what is the difference between 𝒬ᵃᵒₙₑₜ and 𝒬ᵃᵒ

[navidcy]

I think it's part of the discussion why the budget doesn't close; there might be a difference or there may be a double counting of some parts.

But if there isn't any in the end, we should drop the "net".

[navidcy]

using NumericalEarth
using Oceananigans
using Oceananigans.Units
using Dates
using Statistics
using Printf

using CUDA; CUDA.device!(3)

arch = GPU()
Nx = 360
Ny = 180
Nz = 50

depth = 5000meters
z = ExponentialDiscretization(Nz, -depth, 0; scale = depth/4)

underlying_grid = TripolarGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z)
underlying_grid = LatitudeLongitudeGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z, longitude = (0, 360), latitude = (-80, 80))
bottom_height = regrid_bathymetry(underlying_grid;
                                  minimum_depth = 10,
                                  interpolation_passes = 10,
                                  major_basins = 2)
grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height);
                            active_cells_map=true)

free_surface       = SplitExplicitFreeSurface(grid; substeps=70)
momentum_advection = WENOVectorInvariant(order=5)
tracer_advection   = WENO(order=5)
vertical_mixing = NumericalEarth.Oceans.default_ocean_closure()
ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface,
                         closure=(vertical_mixing,))
sea_ice = sea_ice_simulation(grid, ocean; advection=tracer_advection)

date = DateTime(1993, 1, 1)
ecco_set = MetadataSet(:temperature, :salinity,
                       :sea_ice_thickness, :sea_ice_concentration;
                       dataset = ECCO4Monthly(), date)

set!(ocean.model,   ecco_set)   # T, S
set!(sea_ice.model, ecco_set)   # h, ℵ

atmosphere = JRA55PrescribedAtmosphere(arch)
land       = JRA55PrescribedLand(arch)
radiation  = JRA55PrescribedRadiation(arch)
esm = OceanSeaIceModel(ocean, sea_ice; atmosphere, land, radiation)

simulation = Simulation(esm; Δt=20minutes, stop_time=5*365days)

wall_time = Ref(time_ns())

function progress(sim)
    ocean = sim.model.ocean
    u, v, w = ocean.model.velocities
    T = ocean.model.tracers.T
    e = ocean.model.tracers.e
    Tmin, Tmax, Tavg = minimum(T), maximum(T), mean(view(T, :, :, ocean.model.grid.Nz))
    emax = maximum(e)
    umax = (maximum(abs, u), maximum(abs, v), maximum(abs, w))

    step_time = 1e-9 * (time_ns() - wall_time[])

    msg1 = @sprintf("time: %s, iter: %d", prettytime(sim), iteration(sim))
    msg2 = @sprintf(", max|uo|: (%.1e, %.1e, %.1e) m s⁻¹", umax...)
    msg3 = @sprintf(", max(e): %.2f m² s⁻²", emax)
    msg4 = @sprintf(", wall time: %s \n", prettytime(step_time))

    @info msg1 * msg2 * msg3 * msg4

    wall_time[] = time_ns()

     return nothing
end

# And add it as a callback to the simulation.
add_callback!(simulation, progress, IterationInterval(200))

mht = Field(meridional_heat_transport(esm))

ocean.output_writers[:mth] = JLD2Writer(ocean.model, (; mht);
                                        schedule = TimeInterval(3hours),
                                        filename = "ocean_one_degree_mht",
                                        overwrite_existing = true)

run!(simulation)

##

using Oceananigans

mht  = FieldTimeSeries("ocean_one_degree_mht.jld2", "mht"; backend = OnDisk())

times = mht.times
Nt = length(times)

grid = mht.grid
Ny = size(mht.grid, 2)

mht_mean  = deepcopy(mht[1][1, :, 1])

for iter in 1:Nt
    @info "iteration $iter out of $Nt"
    mht_mean  +=  mht[iter][1, :, 1]
end

@. mht_mean = mht_mean / Nt

using CairoMakie

fig = Figure()
ax = Axis(fig[1, 1], xlabel="latitude (deg)", ylabel="MHT (PW)")

φ = φnodes(grid, Face())

lines!(ax, φ, mht_mean[1:Ny+1]  / 1e15, linewidth=4)

save("mht.png", fig)

[navidcy]

(pasting here the example that I used to make the plots; x-ref: #277)