Trying to get the lake version of NIVAFjord to work

dev_notes/todo.txt

@@ -1,6 +1,12 @@
 
 
 	*** High pri ***
+
+
+		There is a pretty bad order-of-module load error in nivafjord_lake. The model can be wrong or not load if airsea_fpv (or airsea) is loaded in the outer and not the inner model!!!!!
+			The energy loss in the simplytox version seems unrelated. What causes it to happen here and not in the other version though?? Seems unrelated to the river too.
+
+
 
 		Add forums on github page? (github Discussions)
 

models/modules/airseaRes.txt

@@ -0,0 +1,169 @@
+
+module("AirSeaRes", version(0, 1, 0),
+	air  : compartment,
+	surf : compartment,
+
+	ice        : quantity,
+
+	evap_mm    : property,
+
+	temp      : property,
+	precip    : property,
+	wind      : property,
+	g_rad     : property,
+	pressure  : property,
+	rho       : property,
+	a_hum     : property,
+	lwd       : property,
+	sw        : property,
+	attn      : property,
+	area	  : loc,
+
+	top_water_heat    : loc, # TODO: System for composing locs also? E.g. top_water.heat
+	top_water_heat_sw : loc, # Needed since it could be on a separate connection e.g. in NIVAFjord. This is a bit awkward though..
+	top_water_temp    : loc
+) {
+"""
+Air-sea/lake heat fluxes are based off of
+[Air-Sea bulk transfer coefficients in diabatic conditions, Junsei Kondo, 1975, Boundary-Layer Meteorology 9(1), 91-112](https://link.springer.com/article/10.1007/BF00232256)
+The implementation used here is influenced by the implementation in [GOTM](https://github.com/gotm-model).
+
+The ice module is influenced by the ice module in MyLake, and uses Stefan's law for ice accumulation.
+"""
+
+	load("stdlib/seawater.txt", library("Air-sea"))
+	load("stdlib/atmospheric.txt", library("Meteorology"), library("Radiation"))
+	load("stdlib/physiochemistry.txt", library("Thermodynamics"), library("Water utils"))
+
+	rho_ice    : constant("Ice density", [k g, m-3], 917)
+	l_freeze   : constant("Latent heat of freezing", [J, k g-1], 333500)
+	lambda_ice : constant("Ice heat conduction coefficient", [W, m-1, K-1], 2.1)
+	water_alb  : constant("Water albedo", [], 0.045)  # Could be a module parameter
+
+	par_group("Ice", surf) {
+		init_ice  : par_real("Initial ice thickness", [m], 0, 0, 10)
+		ice_alb   : par_real("Ice albedo", [], 0.4, 0, 1)
+		ice_att_c : par_real("Ice attenuation coefficient", [m-1], 5)
+		th_frazil : par_real("Frazil threshold", [m], 0.05, 1e-5, 0.1)
+		freeze_temp  : par_real("Ice formation temperature", [deg_c], 0, -2, 2)
+		always_cover : par_bool("Always cover", true)
+		hum_scalor : par_real("Humidity scalor", [], 1, 0.5, 1.5)
+	}
+
+	stab   : property("Stability")
+	ced    : property("Transfer coefficient for latent heat flux")
+	chd    : property("Transfer coefficent for sensible heat flux")
+	s_hum  : property("Saturation specific humidity")
+	lwu    : property("Emitted longwave radiation")
+	albedo : property("Albedo")
+
+	energy    : property("Freeze energy")
+	indicator : property("Ice indicator")
+
+	var(air.wind, [m, s-1], "Wind speed")
+
+	var(surf.stab, []) {   surface_stability(air.wind, top_water_temp, air.temp)   }
+
+	var(surf.ced, [])  {   tc_latent_heat(air.wind, stab)   }
+
+	var(surf.chd, [])  {   tc_sensible_heat(air.wind, stab)   }
+
+	var(surf.s_hum, []) {
+		svap := saturation_vapor_pressure(top_water_temp),
+		specific_humidity_from_pressure(air.pressure, svap)
+	}
+
+	var(surf.lwu, [W, m-2]) {
+		emissivity := 0.98,
+		emissivity * black_body_radiation(top_water_temp->[K])
+	}
+
+	var(surf.albedo, []) {
+		ice_alb        if ice.indicator,
+		water_alb      otherwise
+	}
+
+	var(surf.sw, [W, m-2]) {   (1 - albedo) * (1 - ice.attn) * air.g_rad   }
+
+	flux(nowhere, top_water_heat_sw, [J, day-1], "Net shortwave") {  area * surf.sw ->>  }
+
+	flux(nowhere, top_water_heat,    [J, day-1], "Net longwave") {
+		net_rad := (1 - surf.albedo)*air.lwd - surf.lwu,    # W/m2       # TODO: Should really albedo be detracted from longwave?
+		0                              if surf.ice.indicator,
+		area * net_rad ->>           otherwise
+	}
+
+	flux(nowhere, top_water_heat, [J, day-1], "Freeze heating") {  area * surf.ice.energy ->>  }    # Energy used to melt ice instead of heating the water body (or other way around with freezing)
+
+	flux(nowhere, top_water_heat, [J, day-1], "Latent heat flux") {
+		l_vap := latent_heat_of_vaporization(top_water_temp),
+		0                                                                               if surf.ice.indicator,
+		area * (surf.ced * l_vap * air.rho * air.wind * (air.a_hum - surf.s_hum)) ->> otherwise
+	}
+
+	flux(nowhere, top_water_heat, [J, day-1], "Sensible heat flux") {
+		0                                                                                            if surf.ice.indicator,
+		area * (surf.chd * C_air * air.rho * air.wind * (air.temp->[K] - top_water_temp->[K])) ->> otherwise
+	}
+
+	var(surf.evap_mm, [m m, day-1]) {
+		rho_ref := 1025 [k g, m-3],
+		0                                                                          if surf.ice.indicator,
+		-(air.rho / rho_ref) * chd * air.wind * (air.a_hum - s_hum) * hum_scalor ->>            otherwise
+	}
+
+	flux(nowhere, top_water_heat, [J, day-1], "Precipitation heat") {
+		precip_t := max(0.4[deg_c], air.temp),
+		#0 if surf.ice.indicator,
+		#{
+		#precip_t := air.temp,
+		V := area*air.precip -> [m 3, day-1],   
+		water_temp_to_heat(V => [m 3], precip_t) => [J, day-1] # TODO: Not sure how to best express this. It is a bit annoying that the 1/day doesn't "pass through" the function
+		#} otherwise
+	}
+
+	var(surf.ice, [m], "Ice thickness") @initial { init_ice }
+
+	var(surf.ice.energy, [W, m-2]) {
+		0[W, m-2] if always_cover,
+		{
+			z_surf := 1[m],   #  Thickness of water layer that ice formation can draw energy from. #TODO: Should probably be module parameter?
+			K_ice := 200[W, m-3, deg_c-1],   # freezing/melting heat transfer coefficient ... should not affect magnitude of temperature or ice cover, only how fast they "converge"
+			e := (freeze_temp - top_water_temp)*z_surf*K_ice,
+
+			0        if (ice < 1e-6[m]) & (e < 0),    # No melting when there is no ice
+			e        otherwise
+		} otherwise
+	}
+
+	var(surf.ice.indicator, [], "Ice indicator") {   always_cover | (ice > th_frazil)   }
+
+	var(surf.ice.attn, []) {
+		0                       if !indicator,
+		1 - exp(-ice_att_c * ice) otherwise
+	}
+
+	var(surf.ice.temp, [deg_c], "Ice surface temperature") {
+		#alfa := 1 / (10[m-1] * ice),
+		#(alfa * freeze_temp + air.temp) / (1 + alfa)   if indicator,
+		#0                                              otherwise
+		0 =>>
+	}
+
+	flux(nowhere, surf.ice, [m, day-1], "Ice change") {
+		#NOTE: Stefan's law in ODE form:
+		#d(h_ice)/dt = lambda_ice*(ice_bottom_temperature - ice_surface_temperature)/(h_ice*rho_ice*l_freeze)
+		0[m, day-1] if always_cover,
+		{
+			(energy + {
+				lambda_ice * (freeze_temp->[K] - temp->[K]) / ice  if (temp <= freeze_temp) & indicator,  # TODO: Should the condition be on air.temp or ice.temp being less than freeze_temp?
+				-(1-albedo)*air.g_rad*attn                         if indicator,  # Hmm, MyLake includes all fluxes here, not just shortwave. Not sure if it matters since the thawing is so fast.
+				0                                                  otherwise
+			}) / (rho_ice * l_freeze) ->>
+			+ {
+				air.precip   if indicator & (air.temp <= freeze_temp),     # if there is ice, and temperature is below freezing, precipitation is added as ice
+				0            otherwise
+			} ->>
+		} otherwise
+	}
+}

models/modules/nivafjord/basin.txt

@@ -99,7 +99,6 @@ Authors: Magnus D. Norling
 	var(layer.dz, [m], "Layer thickness") {
 		# NOTE: This is a bit simplified since if the level changes, A would in reality also change.
 		Aavg := 0.5*(A + A[vert.below]),
-		#Aavg := A,
 		water / Aavg
 	} @initial {    dz0   }
 

models/modules/nivafjord/lake_basin.txt

@@ -0,0 +1,38 @@
+
+
+module("NIVAFjord lake", version(0, 0, 1),
+	river : compartment,
+	basin : compartment,
+	layer : compartment,
+	water : quantity,
+	flow  : property,
+	dz    : property,
+	river_lake : connection,
+	vert  : connection,
+	dims  : preamble
+) {
+"""
+This is a module that can be used if a NIVAFjord basin is to be treated as a lake with river inputs and discharge.
+(in development)
+"""
+
+	par_group("Lake discharge", basin) {
+		rl : par_real("Rating curve linear component", [m 2, s-1], 1, 0.1, 1000) 
+	}
+
+	# Right now it can only have one lake at the end. Should maybe configure so that you can have multiple lakes..
+
+	# Note: only discharge to top layer. We could go with the more detailed NIVAFjord approach, but it is probably unnecessary?
+
+	flux(river.water, layer.water[river_lake.below, vert.top], [m 3, s-1], "River discharge to lake") {
+		river.water.flow
+	}
+
+	# TODO: more detailed rating curve?
+
+	flux(layer.water[vert.top], out, [m 3, s-1], "Lake outlet discharge") {
+		excess := max(0, dz[vert.top] - dz0[vert.top]),
+		r1 * excess
+	}
+
+}

models/modules/simplytox.txt

@@ -17,15 +17,15 @@ Parameters for specifying various physiochemical properties of a contaminant.
 
 	par_group("Physiochemistry", tox) {
 		molmass     : par_real("Contaminant molar mass", [g, mol-1], 50, 0, 1000)
-		molvol      : par_real("Contaminant molar volume", [c m 3, mol-1], 20, 0, 1000)
+		molvol      : par_real("Contaminant molar volume", [c m 3, mol-1], 20, 0, 1000, "At surface pressure")
 		dUaw        : par_real("Enthalpy of phase transfer between air and water", [k J, mol-1], 0, -100, 100)
 		dUow        : par_real("Enthalpy of phase transfer between octanol and water", [k J, mol-1], 0, -100, 100)
 		#dUoa        : par_real("Enthalpy of phase transfer between octanol and air", [k J, mol-1], 0, -100, 100)
 		kH25        : par_real("Henry's constant at 25°C", [Pa, m 3, mol-1], 0, 0, 100)  #TODO: better default?
 		log10Kow25  : par_real("log10 Octanol-water partitioning coefficient at 25°C", [], 0, -10, 10)
 		#log10Koa25  : par_real("log10 Octanol-air partitioning coefficient at 25°C", [], 0, -10, 10)
-		rwoc        : par_real("Water-OC rate coefficient", [day-1],  0.5, 0, 3)
-		rfs         : par_real("Fast-Slow rate coefficient", [day-1], 0.05, 0, 3)
+		rwoc        : par_real("Water-OC rate coefficient", [day-1],  0.5, 0, 3, "Speed to reach partitioning equilibrium between water and organic carbon")
+		rfs         : par_real("Fast-Slow rate coefficient", [day-1], 0.05, 0, 3, "Speed to reach partitioning equilibrium between fast and slow soil organic carbon")
 		raw         : par_real("Water-air rate coefficient", [m, day-1], 0.05, 0, 3, "Inside soil pores only.")
 		raa         : par_real("Soil atmosphere exchange rate", [m, day-1], 0.05, 0, 10, "Between soil pore air and atmosphere")
 
@@ -234,6 +234,8 @@ Contaminant exchange between a surface water and the atmosphere.
 
 	var(comp.water.tox.moldiff, [m 2, s-1])
 
+	var(air.wind, [m, s-1]) # In case it is not declared anywhere else.
+
 	# Maybe put this one in the library also?
 	waterside_vel_low_turbulence : function(T : [K], moldiff : [m 2, s-1], wind : [m, s-1]) {
 		# This one is for freshwater only.

models/nivafjord_lake_model.txt

@@ -0,0 +1,92 @@
+
+model("NIVAFjord lake") {
+
+"""
+Note, this model is incomplete, it must be extended by one that contains an airsea module.
+
+It is almost possible to unify a lot of this one with nivafjord_model.txt . The only difficulty is the all_basins index set in the other one. We should figure out a way to do it.
+"""
+
+	basin_idx : index_set("Basin index")
+	layer_idx : index_set("Layer index") @sub(basin_idx)
+
+	air       : compartment("Atmosphere")
+	basin     : compartment("Basin",          basin_idx)
+	layer     : compartment("Fjord layer",    basin_idx, layer_idx)
+	sediment  : compartment("Fjord sediment", basin_idx, layer_idx)
+
+	water : quantity("Water")
+	ice   : quantity("Ice")
+	salt  : quantity("Salt")
+	heat  : quantity("Heat energy")
+
+	sed : quantity("Sediments")
+	oc : quantity("Organic carbon")
+	on : quantity("Organic nitrogen")
+	op : quantity("Organic phosphorous")
+	din : quantity("Nitrate")
+	phos : quantity("Phosphate")
+
+	temp     : property("Temperature")
+	salinity : property("Salinity")
+	precip   : property("Precipitation")
+	wind     : property("Wind speed")
+	g_rad    : property("Global radiation")
+	pressure : property("Pressure")
+	a_hum    : property("Actual specific humidity")
+	rho      : property("Density")
+	lwd      : property("Downwelling longwave radiation")
+	sw       : property("Shortwave radiation")
+	cos_z    : property("Cosine of the solar zenith angle") 
+	attn     : property("Attenuation")
+	freeze_t : property("Ice formation temperature")
+	indicator : property("Ice indicator")
+	a_vap    : property("Actual vapor pressure")
+	s_vap    : property("Saturation vapor pressure")
+
+	z    : property("Depth")
+	dz   : property("Thickness")
+	h    : property("Sea level")
+	evap : property("Evaporation")
+	area : property("Area")
+
+	vert    : connection("Fjord vertical") @grid1d(layer, layer_idx)
+	sw_vert : connection("Shortwave vertical") @grid1d(layer, layer_idx)
+
+	load("modules/atmospheric.txt",
+		module("Atmospheric", air, temp, wind, g_rad, pressure, a_hum, rho, lwd, cos_z, a_vap, s_vap))
+
+	#fpv : compartment("Floating photovoltaics", basin_idx)
+	#load("modules/airsea_fpv.txt",
+	#	module("AirSea FPV", air, basin, fpv, ice, heat, temp, precip, wind, g_rad, pressure, rho, a_hum, lwd, sw, attn, indicator,
+	#		evap, cos_z, loc(basin.freeze_t), loc(basin.area), loc(layer.water[vert.top]), loc(layer.water.heat[sw_vert.top])))
+
+	load("modules/nivafjord/basin.txt",
+		dims : preamble("NIVAFjord dimensions", basin_idx, layer_idx, layer),
+		module("NIVAFjord basin", air, basin, layer, water, salt, heat, temp, salinity, pressure, wind, g_rad, rho, attn, z, dz, h, area, freeze_t, sw, loc(basin.ice.indicator), vert, sw_vert, loc(sediment.heat), dims))
+
+	fjord_phyt : quantity("Fjord phytoplankton")
+	zoo        : quantity("Zooplankton")
+	o2         : quantity("O₂")
+	chl_a      : property("Chlorophyll-a")
+
+	load("modules/nivafjord/fjordchem.txt",
+		chem_par : preamble("NIVAFjord chemistry rates", fjord_phyt),
+		module("NIVAFjord chemistry", air, basin, layer, sediment, water, o2, ice, oc, din, on, phos, op, sed, fjord_phyt, zoo, chl_a, temp, salinity, wind, z, dz, indicator, attn, precip, area, cos_z, sw, vert, chem_par, dims))
+
+	load("modules/nivafjord/sediment.txt",
+		module("NIVAFjord sediments", layer, sediment, water, o2, sed, oc, din, on, phos, op, heat, area, temp, dz, vert, chem_par, dims))
+
+
+
+	par_group("Solver step") {
+		solver_h    : par_real("Initial solver step", [hr], 1)
+		solver_hmin : par_real("Solver relative min step", [], 0.01)
+	}
+
+	sol : solver("NIVAFjord solver", inca_dascru, solver_h, solver_hmin)
+	solve(sol, layer.water, basin.ice)
+
+	solve(sol, sediment.sed, sediment.heat)
+	solve(sol, air.cos_z, air.g_rad)
+}

models/nivafjord_lake_simplytox_model.txt

@@ -0,0 +1,40 @@
+
+
+model("NIVAFjord lake SimplyTox") {
+"""
+TODO: This one is incomplete.
+
+Maaybe we want a simpler biochemistry model for this one, not the entire NIVAFjord biochemistry...
+"""
+
+
+	extend("nivafjord_lake_model.txt")
+
+	extend("simplytox_model.txt") @exclude(
+		downstream : connection,
+		air : compartment,
+		water : quantity,
+		heat : quantity,
+		sed : quantity,
+		oc : quantity,
+		temp : property,
+		precip : property,
+		wind : property,
+		sol : solver
+	)
+
+	downstream : connection("Downstream") @directed_graph { river+ } @no_cycles
+
+	river_lake : connection("River-lake") @directed_graph { river+ layer } @no_cycles
+
+	load("modules/airsea.txt",
+		module("AirSea", "AirSea Fjord", air, basin, ice, heat, temp, precip, wind, g_rad, pressure, rho, a_hum, lwd, sw, attn, indicator,
+							evap, cos_z, loc(basin.freeze_t), loc(basin.area), loc(layer.water[vert.top]), loc(layer.water.heat[sw_vert.top]))
+	)
+
+	load("modules/nivafjord/lake_basin.txt",
+		module("NIVAFjord lake", river, basin, layer, water, flow, dz, river_lake, vert, dims)
+	)
+
+
+}