Fixing MAGIC

dev_notes/todo.txt

@@ -2,6 +2,11 @@
 
 	*** High pri ***
 
+
+		Remember to update existing MAGIC setups with the new runoff scale mechanic.
+
+
+
 		If an input series is provided indexed over "subcatchment", but should be indexed over "water body" and "subcatchment" is a union member of "water body", that should be made to work. (and all similar cases) Just re-map the index.
 
 		If you load different modules using the same load name, you don't get a name clash!
@@ -40,6 +45,19 @@
 				SolveFreeFluoride can go into an infinite loop if F is small, but depends on Al.
 				Should be more careful that concentrations don't go negative.
 
+				How to fix it so that runoff is scaled correctly wrt rel_area automatically.
+				Difficult without some kind of sum_above or aggregate_above(rel_area, con) for graph connection.
+					But we could make that.
+					It is sort of like an in_flux_connection, except not for fluxes.
+					Importantly, this should not have the agg_weight applied, but that could be confusing??
+					It could be an aggregation for the edge compartment (if we make one)
+				Or make the user input a separate runoff scale.
+				Problematic to compute in cases where you have split runoff or if it is vertical, not horizontal.
+
+				This happens if you try to access rel_area[con.below] in code (not in agg. weight)
+				ERROR (model_building): In file modules/magic/magic_forest_drivers.txt line 107 column 19:
+				A 'below' indexing can only be applied to a parameter in a context where it can be guaranteed that the source and target of the flux are the same node type of the graph.
+
 
 
 

models/example_data/MAGIC/birkenes_test.dat

@@ -77,7 +77,7 @@ data_set {
 			[ 0.33 0.07 ]
 
 			par_real("Depth") 
-			[ 0.4 2 ]
+			[ 0.4 0.5 ]
 
 			par_real("Porosity") 
 			[ 0.5 0 ]
@@ -149,6 +149,9 @@ data_set {
 
 			par_real("Minimal compartment temperature") 
 			[ 0.4 0.4 ]
+
+			par_real("Runoff scale") 
+			[ 1 20 ]
 
 		}
 
@@ -217,6 +220,56 @@ data_set {
 
 			par_real("PO4 conc in precipitation") 
 			[ 0 ]
+
+			par_bool("Sea salt deposition") 
+			[ true ]
+
+		}
+
+		par_group("Overall climate parameters") {
+
+			par_real("Precipitation") 
+			[ 1 ]
+
+			par_real("Runoff") 
+			[ 1 ]
+
+			par_real("Air temperature") 
+			[ 4 ]
+
+		}
+
+		par_group("Sink parameters", ci) {
+
+			par_real("Ca sink") 
+			[ 0 0 ]
+
+			par_real("Mg sink") 
+			[ 0 0 ]
+
+			par_real("Na sink") 
+			[ 0 0 ]
+
+			par_real("K sink") 
+			[ 0 0 ]
+
+			par_real("NH4 sink") 
+			[ 0 0 ]
+
+			par_real("SO4 sink") 
+			[ 0 0 ]
+
+			par_real("Cl sink") 
+			[ 0 0 ]
+
+			par_real("NO3 sink") 
+			[ 0 0 ]
+
+			par_real("F sink") 
+			[ 0 0 ]
+
+			par_real("PO4 sink") 
+			[ 0 0 ]
 
 		}
 

models/example_data/MAGIC/birkenes_uncalibrated.dat

@@ -77,7 +77,7 @@ data_set {
 			[ 0.33 0.07 ]
 
 			par_real("Depth") 
-			[ 0.4 2 ]
+			[ 0.4 0.5 ]
 
 			par_real("Porosity") 
 			[ 0.5 0 ]
@@ -149,6 +149,9 @@ data_set {
 
 			par_real("Minimal compartment temperature") 
 			[ 0.4 0.4 ]
+
+			par_real("Runoff scale") 
+			[ 1 20 ]
 
 		}
 
@@ -217,6 +220,56 @@ data_set {
 
 			par_real("PO4 conc in precipitation") 
 			[ 0 ]
+
+			par_bool("Sea salt deposition") 
+			[ true ]
+
+		}
+
+		par_group("Overall climate parameters") {
+
+			par_real("Precipitation") 
+			[ 1 ]
+
+			par_real("Runoff") 
+			[ 1 ]
+
+			par_real("Air temperature") 
+			[ 4 ]
+
+		}
+
+		par_group("Sink parameters", ci) {
+
+			par_real("Ca sink") 
+			[ 0 0 ]
+
+			par_real("Mg sink") 
+			[ 0 0 ]
+
+			par_real("Na sink") 
+			[ 0 0 ]
+
+			par_real("K sink") 
+			[ 0 0 ]
+
+			par_real("NH4 sink") 
+			[ 0 0 ]
+
+			par_real("SO4 sink") 
+			[ 0 0 ]
+
+			par_real("Cl sink") 
+			[ 0 0 ]
+
+			par_real("NO3 sink") 
+			[ 0 0 ]
+
+			par_real("F sink") 
+			[ 0 0 ]
+
+			par_real("PO4 sink") 
+			[ 0 0 ]
 
 		}
 

models/modules/magic/magic_core_wrapper.txt

@@ -60,31 +60,31 @@ This module is for now just for testing calling the MAGIC core from Mobius2. Wil
 
 	var(comp.ca,  [m eq, m-2]) @initial {
 		scececa := {
-			(exch_ca->[])*cec*depth*bd if is_soil,
+			(exch->[])*cec*depth*bd if is_soil,
 			0                          otherwise
 		},
 		scececa + water*concn
 	}
 
 	var(comp.mg,  [m eq, m-2]) @initial {
 		scecemg := {
-			(exch_mg->[])*cec*depth*bd if is_soil,
+			(exch->[])*cec*depth*bd if is_soil,
 			0                          otherwise
 		},
 		scecemg + water*concn
 	}
 
 	var(comp.na,  [m eq, m-2]) @initial {
 		scecena := {
-			(exch_na->[])*cec*depth*bd if is_soil,
+			(exch->[])*cec*depth*bd if is_soil,
 			0                          otherwise
 		},
 		scecena + water*concn
 	}
 
 	var(comp.k,   [m eq, m-2]) @initial {
 		scecek := {
-			(exch_k->[])*cec*depth*bd if is_soil,
+			(exch->[])*cec*depth*bd if is_soil,
 			0                          otherwise
 		},
 		scecek + water*concn
@@ -94,7 +94,7 @@ This module is for now just for testing calling the MAGIC core from Mobius2. Wil
 
 	var(comp.so4, [m eq, m-2]) @initial {
 		smeso4 := {
-			(exch_so4->[])*depth*bd*so4_max  if is_soil,
+			(exch->[])*depth*bd*so4_max  if is_soil,
 			0                                otherwise
 		},
 		smeso4 + water*conc_all
@@ -111,11 +111,7 @@ This module is for now just for testing calling the MAGIC core from Mobius2. Wil
 
 
 	ph        : property("pH")
-	exch_ca   : property("Exchangeable Ca on soil as % of CEC")
-	exch_mg   : property("Exchangeable Mg on soil as % of CEC")
-	exch_na   : property("Exchangeable Na on soil as % of CEC")
-	exch_k    : property("Exchangeable K on soil as % of CEC")
-	exch_so4  : property("Exchangeable SO4 on soil as % of max cap")
+	exch      : property("Exchangeable on soil as % of capacity")
 	bss       : property("Base saturation of soil (ECa + EMg + ENa + EK)")
 
 	conc_hco3 : property("HCO3 (bicarbonate) ionic concentration")
@@ -181,11 +177,11 @@ This module is for now just for testing calling the MAGIC core from Mobius2. Wil
 	}
 
 	var(comp.ph, [])
-	var(comp.exch_ca, [perc]) @initial { init_eca }
-	var(comp.exch_mg, [perc]) @initial { init_emg }
-	var(comp.exch_na, [perc]) @initial { init_ena }
-	var(comp.exch_k, [perc])  @initial { init_ek }
-	var(comp.exch_so4, [perc])
+	var(comp.ca.exch, [perc], "Exchangeable Ca on soil as % of CEC") @initial { init_eca }
+	var(comp.mg.exch, [perc], "Exchangeable Mg on soil as % of CEC") @initial { init_emg }
+	var(comp.na.exch, [perc], "Exchangeable Na on soil as % of CEC") @initial { init_ena }
+	var(comp.k.exch, [perc], "Exchangeable K on soil as % of CEC")  @initial { init_ek }
+	var(comp.so4.exch, [perc], "Exchangeable SO4 on soil as % of max cap")
 	var(comp.bss, [perc])
 	var(comp.conc_hco3, [m mol, m-3])
 	var(comp.conc_co3, [m mol, m-3])
@@ -220,11 +216,11 @@ This module is for now just for testing calling the MAGIC core from Mobius2. Wil
 		result(comp.so4.conc_all)
 		result(comp.f.conc_all)
 		result(comp.ph)
-		result(comp.exch_ca)
-		result(comp.exch_mg)
-		result(comp.exch_na)
-		result(comp.exch_k)
-		result(comp.exch_so4)
+		result(comp.ca.exch)
+		result(comp.mg.exch)
+		result(comp.na.exch)
+		result(comp.k.exch)
+		result(comp.so4.exch)
 		result(comp.bss)
 		result(comp.conc_hco3)
 		result(comp.conc_co3)
@@ -277,7 +273,7 @@ This module is for now just for testing calling the MAGIC core from Mobius2. Wil
 		result(comp.k.s_al)
 		result(comp.ph)
 		result(comp.ionic_str)
-		result(comp.exch_so4)
+		result(comp.so4.exch)
 		comp.ca.concn
 		comp.mg.concn
 		comp.na.concn
@@ -288,10 +284,10 @@ This module is for now just for testing calling the MAGIC core from Mobius2. Wil
 		comp.no3.concn
 		comp.f.conc_all
 		comp.po4.concn
-		comp.exch_ca
-		comp.exch_mg
-		comp.exch_na
-		comp.exch_k
+		comp.ca.exch
+		comp.mg.exch
+		comp.na.exch
+		comp.k.exch
 		max_diff
 		is_soil
 		depth

models/modules/magic/magic_forest_drivers.txt

@@ -31,6 +31,7 @@ module("MAGIC-Forest drivers", version(0, 0, 2),
 
 	par_group("Climate parameters", comp) {
 		is_top: par_bool("This is a surface compartment", true, "If true it receives deposition")
+		ro_scale : par_real("Runoff scale", [], 1, 0, 1, "Scale runoff to per m-2 of this compartment")
 		min_t : par_real("Minimal compartment temperature", [deg_c], 0.4, -10, 10)
 		oa    : par_real("Organic acid concentration", [m mol, m-3], 0, 0, 200)
 		adjoa : par_bool("Adjust OA concentration", false, "Adjust based on SO4 concentration")
@@ -102,9 +103,15 @@ module("MAGIC-Forest drivers", version(0, 0, 2),
 		time.step_length_in_seconds*1[month-1] / (time.days_this_year * 86400[s, day-1])
 	}
 
+	# TODO: It is not the best solution to make the user input the ro_scale
+	# (even more complicated if you have more chained parameters)
+	# (scale is sum of rel_area of compartments before it including itself, divided by rel_area of itself)
+	# (before only counts those sideways, not those in a layer above)
 	var(comp.rel_runoff, [m, month-1]) {
-		qout / rel_area ->>
-	}
+		qout*ro_scale->>
+		#(qout / rel_area)->>   if is_finite(qout) & qout > 1e-6[m m, month-1],
+		#ro_par*ts              otherwise
+	} @no_store
 
 	tot_dep : property("Total deposition")
 

models/modules/simplyq.txt

@@ -1,6 +1,6 @@
 
 
-module("SimplyQ land", version(0, 5, 0),
+module("SimplyQ land", version(0, 5, 1),
 	soil  : compartment,
 	water : quantity,
 	flow  : property,

src/model_specific/MAGIC_core.h

@@ -634,7 +634,7 @@ MagicCore(const magic_input &Input, const magic_param &Param, magic_output &Resu
 	{
 		// Set increment for changing pH and begin the iterative solution loop
 		double dpH = (ChargeBalance < 0.0) ? -0.02 : 0.02;
-		for(int Iter = 0; Iter < 100; ++Iter)
+		for(int Iter = 0; Iter < 10000; ++Iter)
 		{
 			//NOTE: It could be possible to factor out some of the below since it repeats what is done above, but it is a little tricky since some small details are done differently.
 
@@ -729,8 +729,6 @@ MagicCore(const magic_input &Input, const magic_param &Param, magic_output &Resu
 			double ChargeSgn = ChargeBalance0 / ChargeBalance;
 			if(ChargeSgn < 0.0) dpH = -dpH/2.0;
 
-			//printf("pH: %g, dpH: %g  charge balance: %g\n", Result.pH, dpH, ChargeBalance0);
-
 			ChargeBalance = ChargeBalance0;
 
 			// TODO: Have to come up with a specific way to signal an error in an
@@ -747,7 +745,7 @@ MagicCore(const magic_input &Input, const magic_param &Param, magic_output &Resu
 	// Calculate H neutralized by Al species in acidimetric titration (meq/m3)- from H, AL3+, SO4, F and organic acid trivalent anion concens (mmol/m3)
 	double AlTitration = ComputeAlHEquivalent(Coeff, Result.conc_H, Result.conc_Al, Result.conc_SO4, Result.conc_F, Result.conc_A3M);
 
-	// Calulate weak acid alkalinity of solution (meq/m3) )limnological definition)
+	// Calulate weak acid alkalinity of solution (meq/m3) (limnological definition)
 	Result.WeakAcidAlk = Result.conc_HCO3 + 2.0*Result.conc_CO3 + Result.conc_H2AM + 3.0*Result.conc_A3M + Coeff.K_W/Result.conc_H - AlTitration - Result.conc_H;
 
 	// Calculate total aqueous Al concen (free + F-SO4-DOC-complexed) (mmol/m3) - from H, AL3+, SO4, F and organic acid trivalent anion concens (mmol/m3)
@@ -899,9 +897,8 @@ MagicCoreInitial(const magic_init_input &Input, const magic_param &Param, magic_
 	{
 		// Set increment for changing pH and begin the iterative solution loop
 
-		// TODO: dpH should ideally decrease if it detects that we swap back and forth
 		double dpH = (ChargeBalance < 0.0) ? -0.02 : 0.02;
-		for(int Iter = 0; Iter < 100; ++Iter)
+		for(int Iter = 0; Iter < 10000; ++Iter)
 		{
 			Result.pH += dpH;
 
@@ -933,7 +930,8 @@ MagicCoreInitial(const magic_init_input &Input, const magic_param &Param, magic_
 
 			ChargeBalance = ChargeBalance0;
 
-			// TODO: Error if it reached end of iteration!
+			//if(Iter == 999)
+			//	log_print("WARNING: Initial solution failed to converge.\n");
 		}
 	}
 
@@ -942,10 +940,7 @@ MagicCoreInitial(const magic_init_input &Input, const magic_param &Param, magic_
 
 		double exchangeable_Al = 1.0 - Input.exchangeable_Ca - Input.exchangeable_Mg - Input.exchangeable_Na - Input.exchangeable_K;
 
-		//if(exchangeable_Al < 0.0)
-		//	fatal_error(Mobius_Error::model_specific, "Initial exchangeable Ca, Mg, Na, and K sum to more than 100%\n");
-
-		double Ratio =std:: log10(exchangeable_Al / conc_Al);
+		double Ratio = std::log10(exchangeable_Al / conc_Al);
 
 		Result.Log10CaAlSelectCoeff = 2.0*Ratio + 3.0*std::log10(Input.conc_Ca / Input.exchangeable_Ca) - 6.0;
 		Result.Log10MgAlSelectCoeff = 2.0*Ratio + 3.0*std::log10(Input.conc_Mg / Input.exchangeable_Mg) - 6.0;