Skip to content

Commit

Permalink
More module factoring
Browse files Browse the repository at this point in the history
  • Loading branch information
Maginor committed Nov 1, 2024
1 parent 76e6465 commit 6e6f2ad
Show file tree
Hide file tree
Showing 11 changed files with 172 additions and 187 deletions.
15 changes: 10 additions & 5 deletions dev_notes/todo.txt
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,12 @@
*** High pri ***



# Need @else for options.


MobiView
If the plot x axis was set to nonsense since the first selected series was entirely NaN, it won't be fixed by selecting a non-NaN series.
Should not be considered as was_auto_resized if it was entirely NaN, so that the autoresize can trigger again.
If the parameter search is already selected, it doesn't react to clicks.

We could allow editing option_group parameters in MobiView (use the same recompile system as for constant parameters).

Dataset no longer checks if the correct amount of parameter values were provided?
May be intentional, but then it breaks when saving it again.
Expand Down Expand Up @@ -54,9 +52,16 @@


Specific models

Use options to combine airsea_fpv and airsea_fpv_simple
Put them in the main airsea file.

Maybe factor out common phyto module for EasyChem and FjordChem.
Also for O2 (found a difficulty with that earlier, but maybe we could resolve it?)
Also airsea_gas (found a difficulty with that earlier, but maybe we could resolve it now?)

Lots of duplication between the different airsea things in general.

DIC in nivafjord boundary (but how?)

NIVAFjord

Expand Down
3 changes: 0 additions & 3 deletions models/modules/airsea_fpv_simple.txt
Original file line number Diff line number Diff line change
Expand Up @@ -57,9 +57,6 @@ Authors: Francois Clayer, Magnus D. Norling
}

par_group("FPV", fpv) {
#FPV_alb : par_real("FPV albedo", [], 0.2, 0, 1)
#FPV_eff : par_real("FPV efficiency", [], 0.18)
#length_FPV : par_real("Length of FPV module", [m], 0.5)
swr_reduc : par_real("Reduction in SWR by FPV", [], 0.45, 0, 1) ## rest of SWR is reaching the water surface as heat
lwr_reduc : par_real("Reduction in LW back radiation by FPV", [], 0.55, 0, 1)
}
Expand Down
68 changes: 31 additions & 37 deletions models/modules/airsea_gas.txt
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@ module("AirSeaGas", version(0, 0, 0),
layer : compartment,
basin : compartment,
water : quantity,
ice : quantity,
o2 : quantity,
co2 : quantity,
ch4 : quantity,
Expand All @@ -14,66 +13,61 @@ module("AirSeaGas", version(0, 0, 0),
o2satconc : property,
ice_ind : loc,
vert : connection,
dims : preamble
dims : preamble,
compute_dic : par_bool
) {

load("stdlib/seawater.txt", library("Sea oxygen"))
load("stdlib/physiochemistry.txt", library("Chemistry"))

p_vel : property("O₂ piston velocity")
pco2_vel : property("CO₂ piston velocity")
pch4_vel : property("CH4 piston velocity")

var(basin.p_vel, [c m, hr-1]) {
o2_piston_velocity(air.wind, layer.water.temp[vert.top])
}

## TO DO: decide whether to add this in the stdlib together with p_vel for O2.
var(basin.pco2_vel, [c m, hr-1]) {
wnd := air.wind => [],
k_600 := 2.07 + 0.215* (wnd^1.7) => [],
t := layer.water.temp[vert.top] =>[],
schmidt := 1923.6-125.06*t +4.3773*t^2 -0.085681*t^3+0.00070284*t^4,
k_600* min(2.5,(schmidt/600)^(-0.666)) => [c m, hr-1]
}

## TO DO: decide whether to add this in the stdlib together with p_vel for O2.
var(basin.pch4_vel, [c m, hr-1]) {
wnd := air.wind => [],
k_600 := 2.07 + 0.215* (wnd^1.7) => [],
t := layer.water.temp[vert.top] =>[],
schmidt := 1909.4-120.78*t +4.1555*t^2 -0.080578*t^3+0.00065777*t^4,
k_600* min(2.5,(schmidt/600)^(-0.666)) => [c m, hr-1]
}

flux(out, layer.water.o2[vert.top], [k g, day-1], "Precipitation O₂") {
precip_saturation := 0.9, # Probably ok estimate? Not extremely important in any case.
cnc := precip_saturation*o2_saturation(air.temp, 0)*o2_mol_mass,
cnc := precip_saturation*o2_saturation_concentration(air.temp, 0)*o2_mol_mass,
air.precip*A[vert.top]*cnc->>
}

flux(layer.water.o2[vert.top], out, [k g, day-1], "O₂ gas exchange at surface") {
(!ice_ind) * basin.p_vel*(conc(o2[vert.top]) - o2satconc[vert.top])*A[vert.top] ->>
}

option(compute_dic) {

flux(layer.water.co2[vert.top], out, [k g, day-1], "CO₂ gas exchange at surface") {
top_t := (layer.water.temp[vert.top]->[K]) =>[],
lKh := 108.3865 + 0.01985076 * top_t -6919.53 / top_t-40.45154 * log10(top_t)+669365 /(top_t^2), ## Weyhenmeyer et al., 2012
Kh := 10 ^ (lKh=>[]),

co2_eq := 0.012[m g, l-1] *450 * (Kh * 1.013[] * 0.987[]), ## pco2 in ppm : 450 ppm, co2_eq in mgC/L
pco2_vel : property("CO₂ piston velocity")
pch4_vel : property("CH₄ piston velocity")

(!ice_ind) * basin.pco2_vel * (conc(layer.water.co2[vert.top]) - co2_eq ) *A[vert.top] ->>
}

flux(layer.water.ch4[vert.top], out, [k g, day-1], "CH4 gas exchange at surface") {
top_t := (layer.water.temp[vert.top]->[K]) =>[],
var(basin.pco2_vel, [c m, hr-1]) {
co2_piston_velocity(air.wind, layer.water.temp[vert.top])
}

Kh := exp( (-365.183 + 18103.7/top_t + 49.7554*ln(top_t) - 0.000285033*top_t)/1.9872)/55.556,
var(basin.pch4_vel, [c m, hr-1]) {
ch4_piston_velocity(air.wind, layer.water.temp[vert.top])
}

flux(layer.water.co2[vert.top], out, [k g, day-1], "CO₂ gas exchange at surface") {
#top_t := (layer.water.temp[vert.top]->[K]) =>[],
#lKh := 108.3865 + 0.01985076 * top_t -6919.53 / top_t-40.45154 * log10(top_t)+669365 /(top_t^2), ## Weyhenmeyer et al., 2012
#Kh := 10 ^ (lKh=>[]),
#co2_eq := 0.012[m g, l-1] *450 * (Kh * 1.013[] * 0.987[]), ## pco2 in ppm : 450 ppm, co2_eq in mgC/L

co2_eq := co2_saturation_concentration(layer.water.temp[vert.top]),

(!ice_ind) * basin.pco2_vel * (conc(layer.water.co2[vert.top]) - co2_eq ) *A[vert.top] ->>
}

flux(layer.water.ch4[vert.top], out, [k g, day-1], "CH₄ gas exchange at surface") {
top_t := (layer.water.temp[vert.top]->[K]) =>[],

Kh := exp( (-365.183 + 18103.7/top_t + 49.7554*ln(top_t) - 0.000285033*top_t)/1.9872)/55.556,

ch4_eq := 0.016[m g, l-1] * Kh * 2.3, ### ch4_eq in mgCH4/L, 2.3 being in ppm
(!ice_ind) * basin.pch4_vel * (conc(layer.water.ch4[vert.top]) - ch4_eq) * A[vert.top] ->>
ch4_eq := 0.016[m g, l-1] * Kh * 2.3, ### ch4_eq in mgCH4/L, 2.3 being in ppm
(!ice_ind) * basin.pch4_vel * (conc(layer.water.ch4[vert.top]) - ch4_eq) * A[vert.top] ->>
}
}
}

Expand Down
10 changes: 5 additions & 5 deletions models/modules/easychem.txt
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ Authors: Magnus D. Norling
var(river.water.o2, [k g], [m g, l-1], "River oxygen")

flux(out, river.water.o2, [k g, day-1], "River oxygen from catchment") {
catch_conc := f_o2sat*o2_saturation(temp, 0) * o2_mol_mass,
catch_conc := f_o2sat*o2_saturation_concentration(temp, 0) * o2_mol_mass,
in_flux(water)*catch_conc->>
}
}
Expand Down Expand Up @@ -125,8 +125,8 @@ Authors: François Clayer, Magnus D. Norling

f_par : constant("Fraction of PAR in SW radiation", [], 0.45) # Photosynthetically available radiation

var(epi.water.o2, [k g], [m g, l-1], "Epilimnion O₂") @initial_conc { o2_saturation(temp, 0)*init_O2*o2_mol_mass->> }
var(hyp.water.o2, [k g], [m g, l-1], "Hypolimnion O₂") @initial_conc { o2_saturation(temp, 0)*init_O2*o2_mol_mass->> }
var(epi.water.o2, [k g], [m g, l-1], "Epilimnion O₂") @initial_conc { o2_saturation_concentration(temp, 0)*init_O2*o2_mol_mass->> }
var(hyp.water.o2, [k g], [m g, l-1], "Hypolimnion O₂") @initial_conc { o2_saturation_concentration(temp, 0)*init_O2*o2_mol_mass->> }

var(epi.water.oc, [k g], [m g, l-1], "Epilimnion DOC") @initial_conc { init_c }
var(hyp.water.oc, [k g], [m g, l-1], "Hypolimnion DOC") @initial_conc { init_c }
Expand Down Expand Up @@ -167,11 +167,11 @@ Authors: François Clayer, Magnus D. Norling

o2sat : property("O₂ saturation concentration")
var(epi.water.o2sat, [m g, l-1]) {
o2_saturation(temp, 0)*o2_mol_mass -> [m g, l-1]
o2_saturation_concentration(temp, 0)*o2_mol_mass -> [m g, l-1]
}

flux(out, epi.water.o2, [k g, day-1], "Precipitation O₂") {
cnc := 0.9*o2_saturation(air.temp, 0)*o2_mol_mass -> [m g, l-1], #Not sure what the fraction should be..
cnc := 0.9*o2_saturation_concentration(air.temp, 0)*o2_mol_mass -> [m g, l-1], #Not sure what the fraction should be..
air.precip*area*cnc->>
}

Expand Down
29 changes: 15 additions & 14 deletions models/modules/nivafjord/fjordchem.txt
Original file line number Diff line number Diff line change
Expand Up @@ -59,38 +59,35 @@ preamble("NIVAFjord chemistry rates", version(0, 0, 0),

module("NIVAFjord chemistry", version(0, 0, 5),
air : compartment,
basin : compartment,
#basin : compartment,
layer : compartment,
sediment : compartment,
water : quantity,
o2 : quantity,
co2 : quantity,
ch4 : quantity,
ice : quantity,
oc : quantity,
din : quantity,
on : quantity,
phos : quantity,
op : quantity,
sed : quantity,
phyt : quantity,
zoo : quantity,
chl_a : property,
temp : property,
salin : property,
wind : property,
z : property,
dz : property,
indicator : property,
attn : property,
precip : property,
area : property,
cos_z : property,
sw : property,
o2satconc : property,
vert : connection,
chem_par : preamble,
dims : preamble
dims : preamble,

compute_dic : par_bool,
) {
#TODO: Would be nice to unify easychem and fjordchem a bit more, maybe some shared code.
# Esp. for surface O2 exchange and maybe microbes and phyto (into separate modules).
Expand Down Expand Up @@ -120,13 +117,15 @@ module("NIVAFjord chemistry", version(0, 0, 5),

# TODO: Initial states for all these:
# Main variables:
var(layer.water.o2, [k g], [m g, l-1], "Layer O₂") @initial_conc { o2_saturation(temp, salin)*init_O2*o2_mol_mass->> }
var(layer.water.o2, [k g], [m g, l-1], "Layer O₂") @initial_conc { o2_saturation_concentration(temp, salin)*init_O2*o2_mol_mass->> }

var(layer.water.oc, [k g], [m g, l-1], "Layer DOC") @initial_conc { init_DOC }

var(layer.water.co2, [k g], [m g, l-1], "Layer CO₂")

var(layer.water.ch4, [k g], [m g, l-1], "Layer CH4")
option(compute_dic) {
var(layer.water.co2, [k g], [m g, l-1], "Layer CO₂")

var(layer.water.ch4, [k g], [m g, l-1], "Layer CH₄")
}

var(layer.water.din, [k g], [m g, l-1], "Layer DIN")

Expand All @@ -149,7 +148,7 @@ module("NIVAFjord chemistry", version(0, 0, 5),
o2sat : property("O₂ saturation")

var(layer.water.o2satconc, [m g, l-1]) {
o2_saturation(temp, salin)*o2_mol_mass ->>
o2_saturation_concentration(temp, salin)*o2_mol_mass ->>
} #@no_store

var(layer.water.o2sat, []) {
Expand Down Expand Up @@ -347,8 +346,10 @@ module("NIVAFjord chemistry", version(0, 0, 5),
+ (on + sed.on)*min_rat*n_min*(o2_mol_mass/n_mol_mass)*1.5 # ON -> NH4 -> NO3 #TODO: Should the 1.5 be parametrizable?
}

flux(out, layer.water.co2, [k g, day-1], "Microbial CO2 production") {
(oc + sed.oc)*(min_rat+denit_rat) # NOTE: CO2 is in kg(C)
option(compute_dic) {
flux(out, layer.water.co2, [k g, day-1], "Microbial CO2 production") {
(oc + sed.oc)*(min_rat+denit_rat) # NOTE: CO2 is in kg(C)
}
}

var(layer.water.denit_rat, [day-1], "Denitrification rate") {
Expand Down
2 changes: 1 addition & 1 deletion models/modules/nivafjord/point_sources.txt
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ This is a module for point source inputs to NIVAFjord.
flux(out, layer.water.o2, [k g, day-1], "Point source O2 to layer") {
eff_temp := layer.water.temp,
eff_salin := 0,
o2c := o2sat*o2_mol_mass*o2_saturation(eff_temp, eff_salin),
o2c := o2sat*o2_mol_mass*o2_saturation_concentration(eff_temp, eff_salin),
o2c*eff_q/time.days_this_year->>
}

Expand Down
71 changes: 38 additions & 33 deletions models/modules/nivafjord/sediment.txt
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,8 @@ module("NIVAFjord sediments", version(0, 0, 3),
z : property,
vert : connection,
dims : preamble,
chem_par : preamble
chem_par : preamble,
compute_dic : par_bool,
) {

par_group("Sediments") {
Expand Down Expand Up @@ -79,7 +80,9 @@ module("NIVAFjord sediments", version(0, 0, 3),

var(sediment.sed.oc, [k g], [g, k g-1], "Sediment organic carbon") @initial_conc { init_c }

var(sediment.sed.ch4, [k g], [g, k g-1], "Sediment CH4")
option(compute_dic) {
var(sediment.sed.ch4, [k g], [g, k g-1], "Sediment CH4")
}

var(sediment.sed.on, [k g], [g, k g-1], "Sediment organic N") @initial_conc { init_c / cn_molar_to_mass_ratio(sed_cn) }

Expand Down Expand Up @@ -124,37 +127,39 @@ module("NIVAFjord sediments", version(0, 0, 3),
+ (sediment.sed.on)*mr*n_min*(o2_mol_mass/n_mol_mass)*1.5 # ON -> NH4 -> NO3 #TODO: Should the 1.5 be parametrizable?
}

# Simple CO2 release from the sediment = O2 sediment demand
flux(out, layer.water.co2, [k g, day-1], "CO₂ sediment microbial production") {
mr := layer.water.min_rat,
dr := layer.water.denit_rat,
(sediment.sed.oc)*(mr + dr)
}

# Simple CH4 production in the sediment: fraction scaled to CO2 production
flux(out, sediment.sed.ch4, [k g, day-1], "CH4 sediment microbial production") {
mr := layer.water.min_rat,
(sediment.sed.oc) * mr * ch4_production_scaler
}

ch4satconc : property("CH4 saturation concentration")

## Simplified model of CH4 sollubility as a function of water temp, ambient pressure and salinity: Duan and Mao, 2006, https://doi.org/10.1016/j.gca.2006.03.018
# TODO: Name the magic constants, and probably move to stdlib.
var(sediment.sed.ch4satconc, [g, k g-1]) {
ch4_mol_mass*((0.00247[mol, k g-1] - (0.00004033[mol, k g-1, deg_c-1] * sediment.temp))* (1-(0.2844-0.001775[deg_c-1]*sediment.temp) * (layer.water.salin/58.44)))* (1 + 44.6[m-1] * (layer.z[vert.above])/490)
}

## Simple CH4 bubble release from the sediment
flux(sediment.sed.ch4, out, [k g, day-1], "CH4 sediment bubbling") {
max(0, bubble_ch4*(conc(ch4) - ch4satconc)*sediment.sed) ->>
}

## Simple CH4 diffusion release from the sediment; Diffusion coefficient from
flux(sediment.sed.ch4, layer.water.ch4, [k g, day-1], "CH4 sediment diffusion release") {
# TODO: Use Mobius2 inbuilt unit conversion and use named constants (e.g. is 1000[g, l-1] rho_water?)
diff_coeff := (4.18e-11[m 2, s-1, deg_c-1]*sediment.temp + 8.06e-10[m 2, s-1])->[m 2, day-1],
diff_coeff * (conc(sediment.sed.ch4)* 1000[g, l-1] *(porosity)- conc(layer.water.ch4)) * 0.001[k g, g-1] / 0.25 [m] * sediment.area
option(compute_dic) {
# Simple CO2 release from the sediment = O2 sediment demand
flux(out, layer.water.co2, [k g, day-1], "CO₂ sediment microbial production") {
mr := layer.water.min_rat,
dr := layer.water.denit_rat,
(sediment.sed.oc)*(mr + dr)
}

# Simple CH4 production in the sediment: fraction scaled to CO2 production
flux(out, sediment.sed.ch4, [k g, day-1], "CH4 sediment microbial production") {
mr := layer.water.min_rat,
(sediment.sed.oc) * mr * ch4_production_scaler
}

ch4satconc : property("CH4 saturation concentration")

## Simplified model of CH4 sollubility as a function of water temp, ambient pressure and salinity: Duan and Mao, 2006, https://doi.org/10.1016/j.gca.2006.03.018
# TODO: Name the magic constants, and probably move to stdlib.
var(sediment.sed.ch4satconc, [g, k g-1]) {
ch4_mol_mass*((0.00247[mol, k g-1] - (0.00004033[mol, k g-1, deg_c-1] * sediment.temp))* (1-(0.2844-0.001775[deg_c-1]*sediment.temp) * (layer.water.salin/58.44)))* (1 + 44.6[m-1] * (layer.z[vert.above])/490)
}

## Simple CH4 bubble release from the sediment
flux(sediment.sed.ch4, out, [k g, day-1], "CH4 sediment bubbling") {
max(0, bubble_ch4*(conc(ch4) - ch4satconc)*sediment.sed) ->>
}

## Simple CH4 diffusion release from the sediment; Diffusion coefficient from
flux(sediment.sed.ch4, layer.water.ch4, [k g, day-1], "CH4 sediment diffusion release") {
# TODO: Use Mobius2 inbuilt unit conversion and use named constants (e.g. is 1000[g, l-1] rho_water?)
diff_coeff := (4.18e-11[m 2, s-1, deg_c-1]*sediment.temp + 8.06e-10[m 2, s-1])->[m 2, day-1],
diff_coeff * (conc(sediment.sed.ch4)* 1000[g, l-1] *(porosity)- conc(layer.water.ch4)) * 0.001[k g, g-1] / 0.25 [m] * sediment.area
}
}


Expand Down
Loading

0 comments on commit 6e6f2ad

Please sign in to comment.