From 6e6f2ada5c6241c58cef67f680c8e12db9a98c40 Mon Sep 17 00:00:00 2001 From: Magnus Dahler Norling Date: Fri, 1 Nov 2024 13:58:05 +0100 Subject: [PATCH] More module factoring --- dev_notes/todo.txt | 15 ++- models/modules/airsea_fpv_simple.txt | 3 - models/modules/airsea_gas.txt | 68 ++++++------- models/modules/easychem.txt | 10 +- models/modules/nivafjord/fjordchem.txt | 29 +++--- models/modules/nivafjord/point_sources.txt | 2 +- models/modules/nivafjord/sediment.txt | 71 +++++++------- models/nivafjord_chem_base.txt | 14 ++- src/model_application.cpp | 14 ++- src/model_declaration.cpp | 25 ++--- stdlib/seawater.txt | 108 +++++++-------------- 11 files changed, 172 insertions(+), 187 deletions(-) diff --git a/dev_notes/todo.txt b/dev_notes/todo.txt index abf5aeb7..e7f4624d 100644 --- a/dev_notes/todo.txt +++ b/dev_notes/todo.txt @@ -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. @@ -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 diff --git a/models/modules/airsea_fpv_simple.txt b/models/modules/airsea_fpv_simple.txt index 9257a29d..91eb6a1b 100644 --- a/models/modules/airsea_fpv_simple.txt +++ b/models/modules/airsea_fpv_simple.txt @@ -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) } diff --git a/models/modules/airsea_gas.txt b/models/modules/airsea_gas.txt index 19486f79..3a31eb78 100644 --- a/models/modules/airsea_gas.txt +++ b/models/modules/airsea_gas.txt @@ -4,7 +4,6 @@ module("AirSeaGas", version(0, 0, 0), layer : compartment, basin : compartment, water : quantity, - ice : quantity, o2 : quantity, co2 : quantity, ch4 : quantity, @@ -14,41 +13,22 @@ 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->> } @@ -56,24 +36,38 @@ module("AirSeaGas", version(0, 0, 0), (!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] ->> + } } } diff --git a/models/modules/easychem.txt b/models/modules/easychem.txt index af3c6d8e..4f71eb84 100644 --- a/models/modules/easychem.txt +++ b/models/modules/easychem.txt @@ -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->> } } @@ -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 } @@ -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->> } diff --git a/models/modules/nivafjord/fjordchem.txt b/models/modules/nivafjord/fjordchem.txt index ceaa0351..667c92b5 100644 --- a/models/modules/nivafjord/fjordchem.txt +++ b/models/modules/nivafjord/fjordchem.txt @@ -59,14 +59,13 @@ 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, @@ -74,23 +73,21 @@ module("NIVAFjord chemistry", version(0, 0, 5), 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). @@ -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") @@ -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, []) { @@ -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") { diff --git a/models/modules/nivafjord/point_sources.txt b/models/modules/nivafjord/point_sources.txt index 555e38be..15940e4c 100644 --- a/models/modules/nivafjord/point_sources.txt +++ b/models/modules/nivafjord/point_sources.txt @@ -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->> } diff --git a/models/modules/nivafjord/sediment.txt b/models/modules/nivafjord/sediment.txt index 8960a242..d2c00569 100644 --- a/models/modules/nivafjord/sediment.txt +++ b/models/modules/nivafjord/sediment.txt @@ -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") { @@ -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) } @@ -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 + } } diff --git a/models/nivafjord_chem_base.txt b/models/nivafjord_chem_base.txt index 0082bc29..77d93513 100644 --- a/models/nivafjord_chem_base.txt +++ b/models/nivafjord_chem_base.txt @@ -8,8 +8,12 @@ This won't run as its own model, it is just a base for model extension. extend("nivafjord_base.txt") - phyt_type : index_set("Phytoplankton type") + phyt_type : index_set("Phytoplankton type") #Not currently used.. + option_group("FjordChem configurations") { + compute_dic : par_bool("Compute DIC (CO₂, CH₄)", true) + # Eventually we could have other configurations, like include zooplankton, type of phyto module, allow multiple phyto species etc. + } sed : quantity("Sediments") oc : quantity("Organic carbon") @@ -20,7 +24,7 @@ This won't run as its own model, it is just a base for model extension. # Called 'fjord_phyt' to distinguish from freshwater in the version of the model that also has EasyLake fjord_phyt : quantity("Fjord phytoplankton")#, phyt_type) - zoo : quantity("Zooplankton") + #zoo : quantity("Zooplankton") o2 : quantity("O₂") co2 : quantity("CO₂") ch4 : quantity("CH₄") @@ -30,15 +34,15 @@ This won't run as its own model, it is just a base for model extension. load("modules/nivafjord/fjordchem.txt", chem_par : preamble("NIVAFjord chemistry rates", fjord_phyt), - module("NIVAFjord chemistry", air, basin, layer, sediment, water, o2, co2, ch4, ice, oc, din, on, phos, op, sed, fjord_phyt, zoo, chl_a, temp, salinity, wind, z, dz, indicator, attn, precip, area, cos_z, sw, o2satconc, vert, chem_par, dims)) + module("NIVAFjord chemistry", air, layer, sediment, water, o2, co2, ch4, oc, din, on, phos, op, sed, fjord_phyt, chl_a, temp, salinity, z, dz, attn, area, cos_z, sw, o2satconc, vert, chem_par, dims, compute_dic)) option(default_airsea) { load("modules/airsea_gas.txt", - module("AirSeaGas", air, layer, basin, water, ice, o2, co2, ch4, temp, wind, precip, o2satconc, loc(basin.ice.indicator), vert, dims)) + module("AirSeaGas", air, layer, basin, water, o2, co2, ch4, temp, wind, precip, o2satconc, loc(basin.ice.indicator), vert, dims, compute_dic)) } load("modules/nivafjord/sediment.txt", - module("NIVAFjord sediments", layer, sediment, water, o2, co2, ch4, salinity, sed, oc, din, on, phos, op, heat, area, temp, dz, z, vert, chem_par, dims)) + module("NIVAFjord sediments", layer, sediment, water, o2, co2, ch4, salinity, sed, oc, din, on, phos, op, heat, area, temp, dz, z, vert, chem_par, dims, compute_dic)) option(has_boundaries) { load("modules/nivafjord/boudary_chem.txt", diff --git a/src/model_application.cpp b/src/model_application.cpp index 04dc894d..9fbf7e4d 100644 --- a/src/model_application.cpp +++ b/src/model_application.cpp @@ -431,6 +431,12 @@ process_par_group_index_sets(Mobius_Model *model, Data_Set *data_set, Entity_Id } auto par_group = model->par_groups[group_id]; + + if(par_group->decl_type != par_group_data->decl_type) { + par_group_data->source_loc.print_error_header(); + fatal_error("The parameter group \"", par_group->name, "\" should be of type '", name(par_group->decl_type), "'."); + } + if(par_group->max_index_sets.empty()) { if(!par_group_data->index_sets.empty()) { par_group_data->source_loc.print_error_header(); @@ -1132,10 +1138,16 @@ Model_Application::save_to_data_set(Model_Data *save_from) { auto scope_id_data = map_id(model, data_set, par_group->scope_id); scope = data_set->get_scope(scope_id_data); } - par_group_data_id = data_set->par_groups.create_internal(scope, "", par_group->name, Decl_Type::par_group); + par_group_data_id = data_set->par_groups.create_internal(scope, "", par_group->name, par_group->decl_type); data_set->par_groups[par_group_data_id]->scope.parent_id = par_group_data_id; } auto par_group_data = data_set->par_groups[par_group_data_id]; + + if(par_group_data->decl_type != par_group->decl_type) { + par_group_data->source_loc.print_log_header(); + log_print("This parameter group is of the wrong type. Expected it to be '", name(par_group->decl_type), "'.\n"); + } + par_group_data->index_sets.clear(); bool index_sets_resolved = false; diff --git a/src/model_declaration.cpp b/src/model_declaration.cpp index 4c80b278..b5feefb4 100644 --- a/src/model_declaration.cpp +++ b/src/model_declaration.cpp @@ -1060,10 +1060,10 @@ should_exclude_decl(Model_Extension &extend, Decl_AST *decl) { } bool -should_use_option(Mobius_Model *model, Model_Options *options, Argument_AST *arg) { +should_use_option(Mobius_Model *model, Decl_Scope *scope, Model_Options *options, Argument_AST *arg) { std::string symbol = arg->chain[0].string_value; - auto find = model->top_scope[symbol]; + auto find = (*scope)[symbol]; if(!find) { arg->source_loc().print_error_header(); fatal_error("The identifier '", symbol, "' has not been declared."); @@ -1080,6 +1080,7 @@ should_use_option(Mobius_Model *model, Model_Options *options, Argument_AST *arg arg->source_loc().print_error_header(); fatal_error("Expected only one symbol in the argument (constant identifier)."); } + return constant->value.val_boolean; } @@ -1143,10 +1144,10 @@ should_use_option(Mobius_Model *model, Model_Options *options, Argument_AST *arg } void -process_option_decl(Mobius_Model *model, Model_Options *options, Decl_AST *decl, Model_Extension *extend, std::vector> &all_decls); +process_option_decl(Mobius_Model *model, Decl_Scope *scope, Model_Options *options, Decl_AST *decl, Model_Extension *extend, std::vector> &all_decls); void -process_option_body(Mobius_Model *model, Model_Options *options, Decl_Body_AST *body, Model_Extension *extend, std::vector> &all_decls) { +process_option_body(Mobius_Model *model, Decl_Scope *scope, Model_Options *options, Decl_Body_AST *body, Model_Extension *extend, std::vector> &all_decls) { for(Decl_AST *child : body->child_decls) { if(extend && should_exclude_decl(*extend, child)) continue; @@ -1161,14 +1162,14 @@ process_option_body(Mobius_Model *model, Model_Options *options, Decl_Body_AST * // It is unlikely to be a problem in practice... But we should give a proper error if it happens unlike now where // it would just say that it can't find the identifier. - process_option_decl(model, options, child, extend, all_decls); + process_option_decl(model, scope, options, child, extend, all_decls); } else all_decls.push_back({child, extend}); } } void -process_option_decl(Mobius_Model *model, Model_Options *options, Decl_AST *decl, Model_Extension *extend, std::vector> &all_decls) { +process_option_decl(Mobius_Model *model, Decl_Scope *scope, Model_Options *options, Decl_AST *decl, Model_Extension *extend, std::vector> &all_decls) { match_declaration(decl, {{Arg_Pattern::any}}, false, 1, true); // Have to use 'any' since this isn't quite covered by the matching system. bool correct = true; @@ -1183,12 +1184,12 @@ process_option_decl(Mobius_Model *model, Model_Options *options, Decl_AST *decl, fatal_error("Expected a parameter value."); } - bool use_option = should_use_option(model, options, arg); + bool use_option = should_use_option(model, scope, options, arg); if(use_option) { auto body = static_cast(decl->body); - process_option_body(model, options, body, extend, all_decls); + process_option_body(model, scope, options, body, extend, all_decls); } for(auto note : decl->notes) { @@ -1198,7 +1199,7 @@ process_option_decl(Mobius_Model *model, Model_Options *options, Decl_AST *decl, if(!use_option) { auto body = static_cast(note->body); - process_option_body(model, options, body, extend, all_decls); + process_option_body(model, scope, options, body, extend, all_decls); } } else { note->decl.print_error_header(); @@ -1338,7 +1339,7 @@ process_module_load(Mobius_Model *model, Model_Options *options, Token *load_nam } - // TODO: Factor out some of the option processing + // TODO: Factor out some of the option processing between modules and models std::set allowed_first_pass = { Decl_Type::option_group, Decl_Type::constant, Decl_Type::unit }; for(Decl_AST *child : body->child_decls) { @@ -1367,7 +1368,7 @@ process_module_load(Mobius_Model *model, Model_Options *options, Token *load_nam if(child->type == Decl_Type::option_group || child->type == Decl_Type::constant) continue; // Already processed. if(child->type == Decl_Type::option) - process_option_decl(model, options, child, nullptr, all_decls); + process_option_decl(model, &module->scope, options, child, nullptr, all_decls); else all_decls.push_back({child, nullptr}); } @@ -2049,7 +2050,7 @@ load_model(String_View file_name, Mobius_Config *config, Model_Options *options) if(child->type == Decl_Type::option_group || child->type == Decl_Type::constant) continue; // Already processed. if(child->type == Decl_Type::option) - process_option_decl(model, options, child, &extend, all_decls); + process_option_decl(model, &model->top_scope, options, child, &extend, all_decls); else all_decls.push_back({child, &extend}); } diff --git a/stdlib/seawater.txt b/stdlib/seawater.txt index f907704a..1ac8bcb4 100644 --- a/stdlib/seawater.txt +++ b/stdlib/seawater.txt @@ -115,80 +115,20 @@ Riley, J. P. & Skirrow, G. Chemical oceanography. 2 edn, Vol. 1 606 (Academic Pr # ref_diff - the diffusivity of the substance in water at temperature ref_T and salinity ref_S. ref_diff * (dynamic_viscosity_sea_water(ref_T, ref_S) / dynamic_viscosity_sea_water(T, S)) * ((T -> [K]) / (ref_T -> [K])) } - - /* - seawater_secant_bulk_modulus : function(T : [deg_c], S : [], P : [d bar]) { - T68 := T * 1.00024 => [] - h3 := -5.77905e-7 - h2 := 1.16092e-4 - h1 := 1.43713e-3 - h0 := 3.239908 - AW := h0 + (h1 + (h2 + h3*T68)*T68)*T68 - k2 := 5.2787e-8 - k1 := -6.12293e-6 - k0 := 8.50935e-5 - BW := k0 + (k1 + k2*T68)*T68 - e4 := -5.155288e-5 - e3 := 1.360477e-2 - e2 := -2.327105 - e1 := 148.4206 - e0 := 19652.21 - KW := e0 + (e1 + (e2 + (e3 + e4*T68)*T68)*T68)*T68 - j0 := 1.91075e-4 - i2 := -1.6078e-6 - i1 := -1.0981e-5 - i0 := 2.2838e-3 - SR := sqrt(S) - A := AW + (i0 + (i1 + i2*T68)*T68 + j0*SR)*S - m2 := 9.1697e-10 - m1 := 2.0816e-8 - m0 := -9.9348e-7 - B := (BW + (m0 + (m1 + m2*T68)*T68)*S) => [bar-1] - f3 := -6.1670e-5 - f2 := 1.09987e-2 - f1 := -0.603459 - f0 := 54.6746 - g2 := -5.3009e-4 - g1 := +1.6483e-2 - g0 := +7.944e-2 - K0 := (KW + (f0 + (f1 + (f2 + f3*T68)*T68)*T68 + (g0 + (g1 + g2*T68)*T68)*SR)*S) => [bar] - p := P->[bar] - K0 + (A + B*p)*p - } - - seawater_density : function(T : [deg_c], S : [], P : [h bar]) { - densP0 := seawater_dens_p0(T, S) - K := seawater_secant_bulk_modulus(T, S, P); - P := P -> [bar] - densP0/(1.0 - p/K) - } - */ } library("Sea oxygen") { """ -This contains formulas for O2 saturation and surface exchange in sea water. Based on +This contains formulas for O₂, CO₂ and CH₄ saturation and surface exchange in sea water. Based on R.F. Weiss, The solubility of nitrogen, oxygen and argon in water and seawater, Deep Sea Research and Oceanographic Abstracts, Volume 17, Issue 4, 1970, 721-735, [https://doi.org/10.1016/0011-7471(70)90037-9](https://doi.org/10.1016/0011-7471(70)90037-9). The implementation is influenced by the one in [SELMA](https://github.com/fabm-model/fabm/tree/master/src/models/selma). + +There are other undocumented sources. This should be updated soon. """ - /* - o2_saturation_alt : function (T : [deg_c], S : []) { - tk := (T->[K]*0.01)=>[], - c := -173.4292 + 249.6339/tk + 143.3483*ln(tk) - 21.8492*tk - + S*(-0.033096 + 0.014259*tk - 0.0017*tk^2), - exp(c)*1[m l, l-1]*((1.33/32)=>[m mol, m l-1]) # This is at surface pressure. TODO: Make this work better - } - - o2_piston_velocity_alt : function(wind : [m, s-1], temp : [deg_c]) { - t := temp=>[], - schmidt := 1920.4 - 135.6*t + 5.2122*t^2 -0.10939*t^3 + 0.00093777*t^4, - (0.251*wind*2*(schmidt/660)^(-0.5)=>[c m, hr-1]) - } - */ - o2_saturation : function (T : [deg_c], S : []) { + o2_saturation_concentration : function (T : [deg_c], S : []) { # Formula from Weiss 1970 (also in Selma). T_k := (T -> [K]) => [], logsat := - 135.90205 @@ -199,17 +139,43 @@ The implementation is influenced by the one in [SELMA](https://github.com/fabm-m - S*(0.017674 - 10.754/T_k + 2140.7/T_k^2), exp(logsat) => [m mol, m-3] } + + co2_saturation_concentration : function(T : [deg_c]) { + + #TODO: Maybe factor out Henry's constant computation for CO₂ + T_k := (T->[K]) =>[], + lKh := 108.3865 + 0.01985076*T_k - 6919.53/T_k - 40.45154*log10(T_k) + 669365/(T_k^2), ## Weyhenmeyer et al., 2012 + Kh := 10^(lKh=>[]), + + #TODO: Better explain the constants in this formula + 0.012*450*(Kh*1.013*0.987) => [m g, l-1] ## pco2 in ppm : 450 ppm, co2_eq in mgC/L + } + + # TODO: Document the source of the below functions! Are they all from the same paper? + + schmidt_600 : function(wind : [m, s-1]) { + wnd := wind => [], + 2.07 + 0.215*(wnd^1.7) + } o2_piston_velocity : function(wind : [m, s-1], temp : [deg_c]) { - # TODO: Document the source! - wnd := wind => [], T := temp => [], - k_600 := 2.07 + 0.215*(wnd^1.7), + k_600 := schmidt_600(wind), schmidt := 1800.6 - 120.10*T + 3.7818*T^2 - 0.047608*T^3, k_600*(schmidt/600)^(-0.666) => [c m, hr-1] - #sc := (1450 + (1.1*T - 71)*T), - #max(0.05, 5.9*(5.9*wnd - 49.3)/sqrt(sc)) if wind >= 13[m, s-1], - #max(0.05, 1.003*wnd / ((sc)^0.66)) if wind <= 3.6[m, s-1], - #max(0.05, 5.9*(2.85*wnd - 9.65) / sqrt(sc) otherwise + } + + co2_piston_velocity : function(wind : [m, s-1], temp : [deg_c]) { + T := temp => [], + k_600 := schmidt_600(wind), + 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] + } + + ch4_piston_velocity : function(wind : [m, s-1], temp : [deg_c]) { + T := temp => [], + k_600 := schmidt_600(wind), + 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] } } \ No newline at end of file