Model reorganization

dev_notes/todo.txt

@@ -56,6 +56,8 @@
 
 			NIVAFjord
 
+				Find a way to unify nivafjord_base and nivafjord_lake_base into one model declaration (using options)
+
 				Need a complete rethink of horizontal fluxes. They cause too much numerical instability when there is large in/outflow for small layers.
 
 					There should really be an error if there is an opening to a neighboring basin on a depth below the max depth of that basin.

models/magic_model_test.txt → models/magic_model_simple.txt

@@ -1,6 +1,6 @@
 
 
-model("MAGIC (test)") {
+model("MAGIC (simple)") {
 
 	ci : index_set("Compartment index")
 
@@ -27,7 +27,6 @@ model("MAGIC (test)") {
 	f   : quantity("F")
 	po4 : quantity("PO4")
 
-
 	ts         : property("Time step size")
 	concn      : property("Ionic concentration")
 	conc_all   : property("Total concentration")
@@ -46,7 +45,10 @@ model("MAGIC (test)") {
 	load("modules/magic/magic_forest_cnp.txt",
 		module("MAGIC-Forest CNP", comp, nh4, no3, po4, ts, ext_in, bal))
 
-	# TODO: parametrize:
-	sol : solver("MAGIC solver", euler, [1/10, month], 0.01)
+	par_group("Solver resolution") {
+		step : par_real("Solver step size", [month], 0.1, 0.0001, 1)
+	}
+
+	sol : solver("MAGIC solver", euler, step)
 	solve(sol, comp.ca, comp.mg, comp.na, comp.k, comp.nh4, comp.so4, comp.cl, comp.no3, comp.f, comp.po4)
 }

models/modules/microplastic/mp_soil_column.txt

@@ -18,11 +18,20 @@ This is a module for transport of microplastic in a single *experimentally contr
 		sph : par_real("Soil grain sphericity", [], 1, 0, 1)
 		dp  : par_real("Soil grain diameter", [m m], 0.5, 0.001, 10)
 
-		grad : par_real("Hydraulic gradient", [], 0.05, 0, 1, "When the box is full")
+		grad : par_real("Hydraulic gradient", [], 0.05, 0, 1, "When the box is full") #TODO: Do this differently
+
+		roots : par_bool("Roots are present", false) #TODO: May also need the fraction of volume this takes up?
 	}
 
-	par_group("General transport") {
+	par_group("General soil properties") {
+		fcs  : par_real("Field capacity saturation", [], 0.2, 0.1, 0.4, "0.1=sand, 0.4=clay")
+		pwps : par_real("Permanent wilting point saturation", [], 0.09, 0.04, 0.22, "0.04=sand, 0.22=clay")
 
+		# Tabular values:
+		# https://extension.okstate.edu/fact-sheets/understanding-soil-water-content-and-thresholds-for-irrigation-management.html
+	}
+
+	par_group("General transport") {
 		transp_f_base : par_real("Base transport factor", [], 1e-5, 0, 1)
 		bio_down : par_real("Bioturbation down rate", [day-1], 0.01, 0, 1)
 		bio_up   : par_real("Bioturbation up rate", [day-1], 0.001, 0, 1)
@@ -40,11 +49,33 @@ This is a module for transport of microplastic in a single *experimentally contr
 	load("stdlib/basic_math.txt", library("Response"), library("Basic"))
 
 	precip : property("Precipitation")
+	pet    : property("Potential evapotranspiration")
 	theta  : property("Porosity")
 	perm   : property("Soil permeability")
 	cond   : property("Soil hydraulic conductivity")
+	etp    : property("Evapotranspiration")
+	etp_above : property("Evapotranspiration above")
+
+	var(air.precip, [m m, day-1])
+
+	var(air.pet, [m m, day-1])
+
+	# TODO: May not be that realistic that etp works from top down since roots will always take some??
+	var(layer.etp, [m m, day-1]) {
+		can_etp := is_at[vert.top] | roots,
+		{
+			maxcap := dz*theta,
+			pot := max(0, air.pet-etp_above), # Remaining unsatisfied pet
+			s_response(water, maxcap*pwps, maxcap*fcs, 0, pot)  # Turn ET off as saturation goes below field capacity and approaches PWP.
+		} if can_etp,
+		0 otherwise
+	}
+
+	var(layer.etp_above, [m m, day-1]) {
+		etp[vert.above] + etp_above[vert.above]
+	} @no_store
 
-	var(air.precip, [m m, day-1], "Precipitation")
+	flux(layer.water, out, [m m, day-1], "Evapotranspiration flux") { etp } @no_store
 
 	var(layer.theta, [], "Layer porosity") { theta0 } # Eventually allow change over time
 
@@ -62,12 +93,13 @@ This is a module for transport of microplastic in a single *experimentally contr
 	perc : property("Percolation")
 
 	# Hmm this is done a bit awkward..
-	water_tot : property("Water above")
-	var(layer.water_tot, [m m]) { water + water_tot[vert.above] } @no_store
+	water_tot : property("Water above") #Not actually above, but above including this.
+	var(layer.water_tot, [m m]) { max(0, water - dz*theta*fcs) + water_tot[vert.above] } @no_store
 	z : property("Effective depth")
 	var(layer.z, [m m]) { dz*theta + z[vert.above] } @no_store
 
 	var(layer.water.side, [m m, day-1]) {
+		# TODO: This is wrong... Should compute gradient based on z and, z of outside and distance.
 		agrad := grad*(water_tot/z),
 		agrad*cond->>
 	}
@@ -77,21 +109,30 @@ This is a module for transport of microplastic in a single *experimentally contr
 		maxcap_b := dz[vert.below] * theta[vert.below],
 
 		maxrate := cond->>,
+		mincap := maxcap*fcs,
+
+		maxperc := s_response(water, mincap, min(1.1*mincap, maxcap), 0, maxrate),
+		#s_response(water[vert.below], maxcap_b, 0.9*maxcap_b, 0, maxperc)
+		#1.1*mincap, mincap, 0, maxrate),      # Cut flow off when water saturation reaches field capacity
+		s_response(water[vert.below], 0.9*maxcap_b, maxcap_b, maxperc, 0)  # Cut flow off when water saturation below fills up
 
-		# TODO: This is just a test.
-		maxperc := s_response(water, 0, 0.1*maxcap, 0, maxrate),
-		s_response(water[vert.below], 0.8*maxcap_b, maxcap_b, maxperc, 0)
+		# TODO: Capillary action when water below is above fc and this layer is below fc?
+		# We could also compute forces on the water, but probably not necessary
+		# https://onlinelibrary.wiley.com/doi/epdf/10.1155/2020/6671479
 	}
 
-	var(layer.water, [m m], "Soil layer pore water") @initial { 0[m m] } # Assume it is dry.
+	var(layer.water, [m m], "Soil layer pore water") @initial { dz*theta*fcs } 
 
 	flux(out, layer.water[vert.top], [m m, day-1], "Precipitation to soil") {  air.precip  }
 
 	flux(layer.water, out, [m m, day-1], "Layer water lateral flow") { side }
 
 	flux(layer.water, vert, [m m, day-1], "Layer water percolation") { perc }
 
-	# TODO: Capillarity, retention of water in a layer, evapotranspiration.
+
+
+	#    Need, min saturation (=field capacity), max saturation (=1?)
+	# Also permanent wilting point (etp cutoff).
 
 
 

models/modules/nivafjord/basin.txt

@@ -6,7 +6,7 @@ preamble("NIVAFjord dimensions", version(0, 0, 0),
 
 	par_group("Layer thickness", all_layer) {
 		dz0    : par_real("Stable layer thickness", [m], 1, 0.1, 5)
-	} @fully_distribute
+	} @fully_distribute # This is because an external_computation uses this parameter, and it can't change its indexing scheme.
 
 	par_group("Layer area", layer) {
 		A      : par_real("Layer area", [m 2], 10000, 0, 1e6)

models/modules/nivafjord/point_sources.txt

@@ -3,8 +3,10 @@ module("NIVAFjord point sources", version(0, 0, 0),
 	layer : compartment,
 	water : quantity,
 	din   : quantity,
+	on    : quantity,
 	phos  : quantity,
 	sed   : quantity,
+	oc    : quantity,
 	o2    : quantity,
 	heat  : quantity,
 	temp  : property
@@ -14,36 +16,42 @@ This is a module for point source inputs to NIVAFjord.
 """
 
 	par_group("Point sources", layer) {
-		eff_q  : par_real("Water point source", [m 3, day-1], 0)
-		eff_n  : par_real("DIN point source", [k g, day-1], 0)
-		eff_p  : par_real("DIP point source", [k g, day-1], 0)
-		eff_ss : par_real("SS point source", [k g, day-1], 0)
+		eff_q  : par_real("Water point source", [m 3, year-1], 0)
+		eff_in : par_real("DIN point source", [M g, year-1], 0)    #Maybe we should introduce the shorthand ton for M g.
+		eff_on : par_real("DON point source", [M g, year-1], 0)
+		eff_p  : par_real("DIP point source", [M g, year-1], 0)
+		eff_ss : par_real("SS point source", [M g, year-1], 0)
+		ss_c   : par_real("SS point source C fraction", [g, k g-1], 0)
 		o2sat  : par_real("Point source water O2 saturation", [], 1)
 	}
 
 	load("stdlib/seawater.txt", library("Sea oxygen"))
 	load("stdlib/physiochemistry.txt", library("Chemistry"))
 
-	flux(out, layer.water, [m 3, day-1], "Point source water to layer") {  eff_q  }
+	flux(out, layer.water, [m 3, day-1], "Point source water to layer") {  eff_q/time.days_this_year  }
 
 	flux(out, layer.water.heat, [J, day-1], "Point source heat to layer") {
-		conc(layer.water.heat)*eff_q  # Assume same temperature in effluent and bay. TODO: This is not always the case.
+		conc(layer.water.heat)*eff_q/time.days_this_year  # Assume same temperature in effluent and bay. TODO: This is not always the case.
 	}
 
 	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*eff_q->>
+		o2c*eff_q/time.days_this_year->>
 	}
 
 	# TODO: CO2, CH4 content of effluent water?
 
-	flux(out, layer.water.din, [k g, day-1], "Point source DIN to layer") {  eff_n  }
+	flux(out, layer.water.din, [k g, day-1], "Point source DIN to layer") {  eff_in/time.days_this_year->>  }
 
-	flux(out, layer.water.phos, [k g, day-1], "Point source DIP to layer") {  eff_p  }
+	flux(out, layer.water.on, [k g, day-1], "Point source DON to layer") {  eff_on/time.days_this_year->>  }
 
-	flux(out, layer.water.sed, [k g, day-1], "Point source suspended solids to layer") {  eff_ss  }
+	flux(out, layer.water.phos, [k g, day-1], "Point source DIP to layer") {  eff_p/time.days_this_year->>  }
+
+	flux(out, layer.water.sed, [k g, day-1], "Point source suspended solids to layer") {  eff_ss/time.days_this_year->>  }
+
+	flux(out, layer.water.sed.oc, [k g, day-1], "Point source POC to layer") {  eff_ss*ss_c/time.days_this_year->> }
 
 	# TODO: What about carbon content of the effluent suspended solids
 }

models/mp_soil_model.txt

@@ -21,7 +21,7 @@ This is a model for simulating microplastic movement in a single experimentally
 		module("Microplastic soil column", air, layer, water, mp, vert)
 	)
 
-	sol : solver("Soil solver", inca_dascru, [10, min], 1e-2)
+	sol : solver("Soil solver", inca_dascru, [10, min], 1e-5)
 
 	solve(sol, layer.water, layer.mp)
 }

models/nivafjord_chem_base.txt

@@ -6,7 +6,7 @@ model("NIVAFjord-Chem-Base") {
 This won't run as its own model, it is just a base for model extension.
 """	
 
-	extend("nivafjord_model.txt")
+	extend("nivafjord_base.txt")
 
 	phyt_type : index_set("Phytoplankton type")
 
@@ -42,7 +42,7 @@ This won't run as its own model, it is just a base for model extension.
 		module("NIVAFjord boundary chemistry", bnd_layer, water, o2, oc, din, on, phos, op, fjord_phyt, chl_a, sed, chem_par))
 
 	load("modules/nivafjord/point_sources.txt",
-		module("NIVAFjord point sources", layer, water, din, phos, sed, o2, heat, temp))
+		module("NIVAFjord point sources", layer, water, din, on, phos, sed, oc, o2, heat, temp))
 
 	load("modules/easychem.txt",
 		module("Simple River O₂", river, water, o2, temp))

models/nivafjord_isolated_basin_fpv_model.txt

@@ -9,7 +9,7 @@ It is also currently without horizontal exchange between basins, so it only work
 
 	# TODO: FPV cover does not dynamically affect wind mixing (not sure to how large a degree that would happen).
 
-	extend("nivafjord_lake_model.txt")
+	extend("nivafjord_lake_base.txt")
 
 	tox_index : index_set("Contaminant type")
 

models/nivafjord_lake_simplytox_model.txt

@@ -2,9 +2,7 @@
 
 model("NIVAFjord lake SimplyTox") {
 """
-TODO: This one is incomplete.
-
-Maaybe we want a simpler biochemistry model for this one, not the entire NIVAFjord biochemistry...
+TODO: Maybe we want a simpler biochemistry model for this one, not the entire NIVAFjord biochemistry...
 """
 
 	extend("simplytox_model.txt") @exclude(
@@ -20,7 +18,7 @@ Maaybe we want a simpler biochemistry model for this one, not the entire NIVAFjo
 		sol : solver
 	)
 
-	extend("nivafjord_lake_model.txt")
+	extend("nivafjord_lake_base.txt")
 
 	downstream : connection("Downstream") @directed_graph { river+ } @no_cycles
 

src/model_declaration.cpp

@@ -1304,6 +1304,13 @@ Solver_Registration::process_declaration(Catalog *catalog) {
 	solver_fun = scope->resolve_argument(Reg_Type::solver_function, decl->args[1]);
 	hmin = 0.01;
 
+	// TODO: There is a weakness with the matcher that it can't know the type of an identifier, so an identifier matches any Decl_Type.
+	// Hence this hack..
+	auto h_tmp = scope->resolve_argument(Reg_Type::unrecognized, decl->args[2]);
+	if(h_tmp.reg_type == Reg_Type::parameter)
+		which += 2;
+
+
 	if(which == 0 || which == 1) {
 		h_unit = scope->resolve_argument(Reg_Type::unit, decl->args[2]);
 		if(which == 1)