Skip to content

Commit

Permalink
trying to not modify the sediment model
Browse files Browse the repository at this point in the history
  • Loading branch information
@jagoosw
jagoosw committed Sep 24, 2023
1 parent 455e885 commit 67536a9
Showing 1 changed file with 74 additions and 56 deletions.
130 changes: 74 additions & 56 deletions src/Boundaries/Sediments/simple_multi_G.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,17 +29,18 @@ end

"""
SimpleMultiG(; grid
fast_decay_rate = 2/day,
slow_decay_rate = 0.2/day,
fast_redfield = 0.1509,
slow_redfield = 0.13,
fast_fraction = 0.74,
slow_fraction = 0.26,
refactory_fraction = 0.1,
nitrate_oxidation_params = (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232),
denitrification_params = (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662),
anoxic_params = (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143),
solid_dep_params = (0.233, 0.336, 982.0, - 1.548))
fast_decay_rate::FT = 2/day,
slow_decay_rate::FT = 0.2/day,
fast_redfield::FT = 0.1509,
slow_redfield::FT = 0.13,
fast_fraction::FT = 0.74,
slow_fraction::FT = 0.26,
refactory_fraction::FT = 0.1,
nitrate_oxidation_params::P1 = (A = - 1.9785, B = 0.2261, C = -0.0615, D = -0.0289, E = - 0.36109, F = - 0.0232),
denitrification_params::P2 = (A = - 3.0790, B = 1.7509, C = 0.0593, D = - 0.1923, E = 0.0604, F = 0.0662),
anoxic_params::P3 = (A = - 3.9476, B = 2.6269, C = - 0.2426, D = -1.3349, E = 0.1826, F = - 0.0143),
depth = abs(znodes(grid, Face())[1]),
solid_dep_params::P4 = (A = 0.233, B = 0.336, C = 982, D = - 1.548, depth = depth))
Return a single-layer "multi G" sediment model (`SimpleMultiG`) on `grid`, where parameters
can be optionally specified.
Expand Down Expand Up @@ -70,17 +71,37 @@ Single-layer multi-G sediment model (Float64)
```
"""
function SimpleMultiG(; grid,
fast_decay_rate = 2/day,
slow_decay_rate = 0.2/day,
fast_redfield = 0.1509,
slow_redfield = 0.13,
fast_fraction = 0.74,
slow_fraction = 0.26,
refactory_fraction = 0.1,
nitrate_oxidation_params = (- 1.9785, 0.2261, -0.0615, -0.0289, - 0.36109, - 0.0232),
denitrification_params = (- 3.0790, 1.7509, 0.0593, - 0.1923, 0.0604, 0.0662),
anoxic_params = (- 3.9476, 2.6269, - 0.2426, -1.3349, 0.1826, - 0.0143),
solid_dep_params = (0.233, 0.336, 982.0, - 1.548))
fast_decay_rate::FT = 2/day,
slow_decay_rate::FT = 0.2/day,
fast_redfield::FT = 0.1509,
slow_redfield::FT = 0.13,
fast_fraction::FT = 0.74,
slow_fraction::FT = 0.26,
refactory_fraction::FT = 0.1,
nitrate_oxidation_params::P1 = (A = - 1.9785,
B = 0.2261,
C = -0.0615,
D = -0.0289,
E = - 0.36109,
F = - 0.0232),
denitrification_params::P2 = (A = - 3.0790,
B = 1.7509,
C = 0.0593,
D = - 0.1923,
E = 0.0604,
F = 0.0662),
anoxic_params::P3 = (A = - 3.9476,
B = 2.6269,
C = - 0.2426,
D = -1.3349,
E = 0.1826,
F = - 0.0143),
depth = abs(znodes(grid, Face())[1]),
solid_dep_params::P4 = (A = 0.233,
B = 0.336,
C = 982,
D = - 1.548,
depth = depth)) where {FT, P1, P2, P3, P4}

@warn "Sediment models are an experimental feature and have not yet been validated."
@info "This sediment model is currently only compatible with models providing NH₄, NO₃, O₂, and DIC."
Expand All @@ -91,7 +112,7 @@ function SimpleMultiG(; grid,
fields = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names))
tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)),
G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)))

bottom_indices = calculate_bottom_indices(grid)

F = typeof(fields)
Expand All @@ -111,20 +132,19 @@ function SimpleMultiG(; grid,
end

adapt_structure(to, sediment::SimpleMultiG) =
SimpleMultiG(adapt(to, sediment.fast_decay_rate),
adapt(to, sediment.slow_decay_rate),
adapt(to, sediment.fast_redfield),
adapt(to, sediment.slow_redfield),
adapt(to, sediment.fast_fraction),
adapt(to, sediment.slow_fraction),
adapt(to, sediment.refactory_fraction),
adapt(to, sediment.nitrate_oxidation_params),
adapt(to, sediment.denitrification_params),
adapt(to, sediment.anoxic_params),
adapt(to, sediment.solid_dep_params),
SimpleMultiG(sediment.fast_decay_rate,
sediment.slow_decay_rate,
sediment.fast_redfield,
sediment.slow_redfield,
sediment.fast_fraction,
sediment.slow_fraction,
sediment.refactory_fraction,
sediment.nitrate_oxidation_params,
sediment.denitrification_params,
sediment.anoxic_params,
sediment.solid_dep_params,
adapt(to, sediment.fields),
nothing,
adapt(to, sediment.bottom_indices))
adapt(to, sediment.tendencies))

sediment_tracers(::SimpleMultiG) = (:C_slow, :C_fast, :C_ref, :N_slow, :N_fast, :N_ref)
sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow,
Expand All @@ -136,7 +156,7 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow,

@inline bottom_index_array(sediment::SimpleMultiG) = sediment.bottom_indices

@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, tendencies, sediment_tendencies)
@kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, timestepper)
i, j = @index(Global, NTuple)

k = bottom_index(i, j, sediment)
Expand All @@ -159,30 +179,28 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow,
Cᵐⁱⁿ = C_min_slow + C_min_fast
Nᵐⁱⁿ = N_min_slow + N_min_fast

reactivity = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1])
k = Cᵐⁱⁿ * day / (sediment.fields.C_slow[i, j, 1] + sediment.fields.C_fast[i, j, 1])

# sediment evolution
sediment_tendencies.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow
sediment_tendencies.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast
sediment_tendencies.C_ref[i, j, 1] = sediment.refactory_fraction * carbon_deposition
sediment.tendencies.Gⁿ.C_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * carbon_deposition - C_min_slow
sediment.tendencies.Gⁿ.C_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * carbon_deposition - C_min_fast
sediment.tendencies.Gⁿ.C_ref[i, j, 1] = sediment.refactory_fraction * carbon_deposition

sediment_tendencies.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow
sediment_tendencies.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast
sediment_tendencies.N_ref[i, j, 1] = sediment.refactory_fraction * nitrogen_deposition
sediment.tendencies.Gⁿ.N_slow[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.slow_fraction * nitrogen_deposition - N_min_slow
sediment.tendencies.Gⁿ.N_fast[i, j, 1] = (1 - sediment.refactory_fraction) * sediment.fast_fraction * nitrogen_deposition - N_min_fast
sediment.tendencies.Gⁿ.N_ref[i, j, 1] = sediment.refactory_fraction * nitrogen_deposition

# efflux/influx
O₂ = tracers.O₂[i, j, k]
NO₃ = tracers.NO₃[i, j, k]
NH₄ = tracers.NH₄[i, j, k]

A, B, C, D, E, F = sediment.nitrate_oxidation_params

pₙᵢₜ = exp(A +
B * log(Cᵐⁱⁿ * day) * log(O₂) +
C * log(Cᵐⁱⁿ * day) ^ 2 +
D * log(reactivity) * log(NH₄) +
E * log(Cᵐⁱⁿ * day) +
F * log(Cᵐⁱⁿ * day) * log(NH₄)) / (Nᵐⁱⁿ * day)
O₂ = tracers.O₂[i, j, 1]
NO₃ = tracers.NO₃[i, j, 1]
NH₄ = tracers.NH₄[i, j, 1]

pₙᵢₜ = exp(sediment.nitrate_oxidation_params.A +
sediment.nitrate_oxidation_params.B * log(Cᵐⁱⁿ * day) * log(O₂) +
sediment.nitrate_oxidation_params.C * log(Cᵐⁱⁿ * day) ^ 2 +
sediment.nitrate_oxidation_params.D * log(k) * log(NH₄) +
sediment.nitrate_oxidation_params.E * log(Cᵐⁱⁿ * day) +
sediment.nitrate_oxidation_params.F * log(Cᵐⁱⁿ * day) * log(NH₄)) / (Nᵐⁱⁿ * day)

#=
pᵈᵉⁿⁱᵗ = exp(sediment.denitrification_params.A +
Expand Down

0 comments on commit 67536a9

Please sign in to comment.