Merge branch 'main' of https://github.com/NIVANorge/Mobius2

models/data/EasyReservoir/easyreservoir_multibasin_Magat.dat

@@ -12,19 +12,19 @@ par_group("System") {
 
 par_group("Lake data") [ "Lake" ] {
 	par_real("Initial lake surface area")
-	[ 25200000 5000000 ]
+	[ 15200000 15000000 ]
 
 	par_real("Precipitation scalor")
 	[ 0.65 0.65 ]
 
 	par_real("Humidity scalor")
-	[ 3.06 3.06 ]
+	[ 3.51 3.51 ]
 
 	par_real("SWR scalor")
 	[ 0.9 0.9 ]
 
 	par_bool("Always cover (FPV)")
-	[ false true ]
+	[ false false ]
 
 	par_real("Degree of FPV coverage")
 	[ 1 0.9 ]
@@ -66,7 +66,7 @@ module("EasyReservoir", version(0, 1, 0)) {
 		[ 28 28 ]
 
 		par_real("Bottom temperature")
-		[ 26 26 ]
+		[ 27 27 ]
 
 		par_real("Epilimnion thickening rate")
 		[ 0.01 0.01 ]

models/easyReservoir_multibasin_EasyChem.txt

@@ -20,6 +20,17 @@ model("EasyReservoir-multibasin-EasyChem") {
 	nitro    : quantity("Nitrogen")
 	phos     : quantity("Phosphorus")
 	resp     : property("Bacterial respiration")
+	dic      : quantity("Dissolved inorganic carbon")
+	ch4      : quantity("Methane")
+
+	ocf       : quantity("Organic Carbon FPV")
+	phytf     : quantity("Phytoplankton FPV")
+	o2f       : quantity("Oxygen FPV")
+	nitrof    : quantity("Nitrogen FPV")
+	phosf     : quantity("Phosphorus FPV")
+	respf     : property("Bacterial respiration FPV")
+	dicf      : quantity("Dissolved inorganic carbon FPV")
+	ch4f      : quantity("Methane FPV")
 	###
 
 	wind     : property("Wind speed")
@@ -29,20 +40,32 @@ model("EasyReservoir-multibasin-EasyChem") {
 	rho      : property("Density")
 	lwd      : property("Downwelling longwave radiation")
 	sw       : property("Shortwave radiation (Light)")
+	sw_heat  : property("SWR as heat")
 	attn     : property("Attenuation coefficient")
+	deg_cov  : property("Degree of coverage")
+	sol_elec : property("Solar energy production")
+	hydro_elec: property("Hydropower energy production")
 	z_e      : property("Thickness")
     abstrac  : property("Lake abstraction")
+	ruleCurve  : property("Rule Curve")
     area     : property("Lake surface area")
 	sed_area : property("Sediment area")
 	temp     : property("Temperature")
+	tempf     : property("Temperature FPV")
 	precip   : property("Precipitation")
 	level	 : property("Water level")
+	wl_Wsaved : property("Water level including saved water")
 	epi_inflow : property("Inflow to reservoir")
 	in_Temp  : property("Inflow Temp")
+	in_Tempf  : property("Inflow Temp FPV")
+	depth    : property("Depth")
+	saved_water : property("Water saved")
 
 
 	heat : quantity("Heat energy")
+	heatf : quantity("Heat energy FPV")
 	water : quantity("Water")
+	water_FPV : quantity("WaterFPV")
 
 	evap_mm  : property("Evaporation per area")
 
@@ -52,25 +75,27 @@ model("EasyReservoir-multibasin-EasyChem") {
 		evap_scalor : par_real("Humidity scalor", [], 1, 0.5, 1.5)
 		swr_scalor : par_real("SWR scalor", [], 1, 0.5, 1.5)
 		always_cover : par_bool("Always cover (FPV)", false)
-		deg_cov : par_real("Degree of FPV coverage", [], 1, 0.2, 1)
+		Area_cov : par_real("Surface Area covered with FPV", [m 2], 1, 0.2, 1)
 	}
 
 	par_group("Lake ensemble data") {
 		A_total : par_real("Total initial surface area", [m 2], 107183, 0, 371e9)
+		constant_din_flux  : par_real("Epilimnion DIN point source", [k g, day-1], 0, 0, 1000000)
+		constant_tdp_flux  : par_real("Epilimnion TDP point source", [k g, day-1], 0, 0, 1000000)
 	}
 
 
 	#load("modules/simplysoiltemp.txt", module("Simply soil temperature", air, soil, snow_box, snow, temp))
 	#load("modules/rivertemp.txt",      module("RiverTemperature", air, soil, river, water, heat, temp))
 	load("modules/atmospheric.txt",    module("Atmospheric", air, temp, wind, g_rad, pressure, a_hum, rho, lwd))
-	load("modules/airseaRes_bis.txt",         module("AirSeaRes_bis", air, epi, ice, FPV, evap_mm, temp, precip, wind, g_rad, pressure, rho, a_hum, lwd, sw, attn, precip_scalor, evap_scalor, swr_scalor, always_cover, deg_cov,
-											loc(epi.area), ht : loc(epi.water.heat), ht, loc(epi.water.temp)))
-	load("modules/easyreservoir.txt",       module("EasyReservoir", air, epi, hyp, water, ice, heat, temp, precip, z_e, evap_mm, abstrac, area, sed_area, A_surf, A_total, precip_scalor, always_cover, level, epi_inflow, in_Temp)) #, loc(downstream), downstream
-	load("modules/easyChem.txt",       module("EasyChem", air, epi, hyp, surf, water, ice, heat, oc, phyt, o2, nitro , phos, temp, precip, wind, z_e, evap_mm, abstrac, area, sed_area, sw, epi_inflow, in_Temp, always_cover, deg_cov, resp, horz, horz_hyp, loc(epi.water.temp)))
+	load("modules/airseaRes_bis.txt",         module("AirSeaRes_bis", air, epi, ice, FPV, evap_mm, temp, precip, wind, g_rad, pressure, rho, a_hum, lwd, sw, sw_heat, attn, deg_cov, sol_elec, depth, saved_water, precip_scalor, evap_scalor, swr_scalor, always_cover, Area_cov,
+											loc(epi.area), ht : loc(epi.water.heat), ht, loc(epi.water.temp), loc(epi.water_FPV.heatf), loc(epi.water_FPV.tempf)))
+	load("modules/easyreservoir.txt",       module("EasyReservoir", air, epi, hyp, surf, water, water_FPV, ice, heat, temp, heatf, tempf, precip, z_e, evap_mm, abstrac, ruleCurve, area, sed_area, sw_heat, hydro_elec, depth, saved_water, A_surf, A_total, precip_scalor, always_cover, level, wl_Wsaved, epi_inflow, in_Temp, in_Tempf)) #, loc(downstream), downstream
+	load("modules/easyChem.txt",       module("EasyChem", air, epi, hyp, surf, water, water_FPV, ice, heat, oc, phyt, o2, dic, ch4, nitro , phos, temp,ocf, phytf, o2f, dicf, ch4f, nitrof , phosf, tempf, precip, wind, z_e, evap_mm, abstrac, area, sed_area, sw, epi_inflow, in_Temp, always_cover, resp, respf, horz, horz_hyp, constant_din_flux, constant_tdp_flux, A_surf, A_total, loc(epi.water.temp)))
 
 	simply_solver : solver("Simply solver", "INCADascru", 0.1, 1e-3)   # NOTE: The small relative minimal step is needed by SimplyTox.
 
-	solve(simply_solver, epi.water, hyp.water, epi.ice)
+	solve(simply_solver, epi.water, hyp.water, epi.water_FPV, hyp.water_FPV, epi.ice)
 
 
 
@@ -93,17 +118,31 @@ model("EasyReservoir-multibasin-EasyChem") {
 
 		flux(epi.water, horz, [m 3, day-1], "Epi horizontal mixing") {
 			(water + water[below]) * mix
-		} 
+		} #@no_carry{heat}
 
 		flux(hyp.water, horz_hyp, [m 3, day-1], "Hyp horizontal mixing") {
 			(water + water[below]) * mix
-		}
+		} #@no_carry{heat}
 
 		flux(hyp.water, horz_hyp, [m 3, day-1], "Water level adjustment") {
 			dl := epi.level - epi.level[below],
 			dl*A_surf*10[day-1]   if dl > 0,
 			0=>>                   otherwise
 		} #@no_carry { heat }
+
+		flux(epi.water_FPV, horz, [m 3, day-1], "Epi horizontal mixing FPV") {
+			(water + water[below]) * mix
+		} #@no_carry{heat}
+
+		flux(hyp.water_FPV, horz_hyp, [m 3, day-1], "Hyp horizontal mixing FPV") {
+			(water + water[below]) * mix
+		} #@no_carry{heat}
+
+		flux(hyp.water_FPV, horz_hyp, [m 3, day-1], "Water level adjustment FPV") {
+			dl := epi.wl_Wsaved - epi.wl_Wsaved[below],
+			dl*A_surf*10[day-1]   if dl > 0,
+			0=>>                   otherwise
+		} #@no_carry { heat }
 	}
 
 

models/modules/airsea_fpv_old.txt

@@ -17,17 +17,30 @@ module("AirSea FPV", version(0, 1, 0),
 	a_hum     : property,
 	lwd       : property,
 	sw        : property,
+	sw_heat   : property,
 	attn      : property,
+	deg_cov   : property,
+	sol_elec  : property,
+	depth     : property,
+	saved_water : property,
 	precip_scalor : par_real,
 	evap_scalor : par_real,
 	swr_scalor : par_real,
 	always_cover : par_bool,
-	deg_cov : par_real,
+	Area_cov : par_real,
 	area	  : loc,
 
+<<<<<<< HEAD:models/modules/airseaRes_bis.txt
+	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,
+	top_water_heatf    : loc,
+	top_water_tempf    : loc
+=======
 	top_water_heat    : loc, 
 	top_water_heat_sw : loc,
 	top_water_temp    : loc
+>>>>>>> 17aa78c63682102f8b629229bee3a4724472c0d1:models/modules/airsea_fpv_old.txt
 ) {
 """
 This is a modification to the AirSea module where you can cover a basin with floating photovoltaics (solar panels) (FPV).
@@ -60,12 +73,14 @@ This is a modification to the AirSea module where you can cover a basin with flo
 	}
 
 	stab   : property("Stability")
+	stabf   : property("Stability FPV")
 	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")
-	sw_heat : property("SWR as heat") 
+	s_humf  : property("Saturation specific humidity FPV")
+	lwuf    : property("Emitted longwave radiation FPV")
+	albedo : property("Albedo") 
 
 	energy    : property("Freeze energy")
 	Ice_indicator : property("Ice indicator")
@@ -83,6 +98,8 @@ This is a modification to the AirSea module where you can cover a basin with flo
 
 	var(surf.stab, []) {   surface_stability(air.wind, top_water_temp, air.temp)   }
 
+	var(surf.stabf, []) {   surface_stability(air.wind, top_water_tempf, air.temp)   }
+
 	var(surf.ced, [])  {   tc_latent_heat(air.wind, stab)   }
 
 	var(surf.chd, [])  {   tc_sensible_heat(air.wind, stab)   }
@@ -97,25 +114,56 @@ This is a modification to the AirSea module where you can cover a basin with flo
 		emissivity * black_body_radiation(top_water_temp->[K])
 	}
 
+	var(surf.s_humf, []) {
+		svap := saturation_vapor_pressure(top_water_tempf),
+		specific_humidity_from_pressure(air.pressure, svap)
+	}
+
+	var(surf.lwuf, [W, m-2]) {
+		emissivity := 0.98,
+		emissivity * black_body_radiation(top_water_tempf->[K])
+	}
+
+
+
 	var(surf.albedo, []) {
 		ice_alb        if ice.Ice_indicator,
 		water_alb      otherwise
 	}
 
+	var(surf.deg_cov, []) {
+		(Area_cov / area)       if FPV.FPV_indicator,
+		0      					otherwise
+	}
+
 	var(surf.sw, [W, m-2]) { 
 		irradiance :=  (1 - albedo) * (1 - ice.attn) * air.g_rad * swr_scalor,
-		irradiance * (1- deg_cov)					if FPV.FPV_indicator,
+		irradiance * max(0,(1- deg_cov))					if FPV.FPV_indicator,
 		irradiance									otherwise
 	}
+
+	var(surf.sol_elec, [W]) { 
+		0										if FPV.FPV_indicator,
+		surf.sw * 300000[m 2] * 0.18[] *4.5[]		otherwise
+	}
+	sol_elecGWh : property("real Solar production")
+	var(surf.sol_elecGWh, [W]) {
+		sol_elec * 24[] * 365[] *0.000000001[]
+	}
+
+	var(surf.saved_water, [m 3, day-1]) {
+		sol_elec / ((surf.depth + 33[m]) *0.55[] * 9.81[m, s-2] * 998 [k g, m-3]) * (24*3600[s, day-1])
+	}
+
 
 	var(surf.sw_heat, [W, m-2]) { 
 		under_FPV := max(0,surf.FPV.A_both * (surf.FPV.FPV_temp ->[K] - surf.FPV.T_s_back ->[K])) * swr_scalor,
 		normal :=  (1 - albedo) * (1 - ice.attn) * air.g_rad * swr_scalor,
-		under_FPV * deg_cov + normal * (1- deg_cov)					if FPV.FPV_indicator,
+		(under_FPV * deg_cov) + (normal * (1- deg_cov))					if FPV.FPV_indicator,
 		normal														otherwise
 	}
 
-	flux(nowhere, top_water_heat_sw, [J, day-1], "Net shortwave") {  area * surf.sw ->>  }
+	flux(nowhere, top_water_heat_sw, [J, day-1], "Net shortwave") {  area * surf.sw_heat ->>  }
 
 	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?
@@ -141,6 +189,23 @@ This is a modification to the AirSea module where you can cover a basin with flo
 		area * (surf.chd * C_air * air.rho * air.wind * (air.temp->[K] - top_water_temp->[K])) ->> otherwise
 	}
 
+
+
+	flux(nowhere, top_water_heatf, [J, day-1], "Latent heat flux FPV") {
+		l_vap := latent_heat_of_vaporization(top_water_tempf),
+		0                                                                               				if surf.ice.Ice_indicator,
+		area * (surf.ced * l_vap * air.rho * air.wind * (air.a_hum - surf.s_humf))* (1- deg_cov) ->>     if surf.FPV.FPV_indicator,
+		#area * (surf.ced * l_vap * air.rho * (air.wind - 0.8 * air.wind * log10(max(air.wind,0.000001[m, s-1])=>[])) * (air.a_hum - surf.s_humf)) ->> otherwise
+		area * (surf.ced * l_vap * air.rho * air.wind * (air.a_hum - surf.s_humf)) ->> otherwise
+	}
+
+	flux(nowhere, top_water_heatf, [J, day-1], "Sensible heat flux FPV") {
+		0                                                                                                                   if surf.ice.Ice_indicator,
+		area * (surf.chd * C_air * air.rho * air.wind * (air.temp->[K] - top_water_tempf->[K])) * (1- deg_cov) ->>			if surf.FPV.FPV_indicator,
+		#area * (surf.chd * C_air * air.rho * (air.wind - 0.8 * air.wind * log10(max(air.wind,0.000001[m, s-1])=>[])) * (air.temp->[K] - top_water_temp->[K])) ->> otherwise
+		area * (surf.chd * C_air * air.rho * air.wind * (air.temp->[K] - top_water_tempf->[K])) ->> otherwise
+	}
+
 	var(surf.evap_mm, [m m, day-1]) {
 		rho_ref := 1025 [k g, m-3],
 		evap := -(air.rho / rho_ref) * chd * air.wind * (air.a_hum - s_hum) * evap_scalor ->> ,
@@ -200,13 +265,14 @@ This is a modification to the AirSea module where you can cover a basin with flo
 	var(surf.FPV.FPV_temp, [deg_c]) {  # deg_c
 		0                      if !FPV_indicator,
 		{
-			#T_front_final := (air.temp+273.15[deg_c])  =>[K],
-			#need to iterate to updte T_front  
-				T_front := (top_water_temp+276.15[deg_c])  =>[K],
-				T_air := (air.temp+273.15[deg_c])=>[K],
+			#Iterate to updte T_front  
+			T_front := (top_water_temp+3[deg_c])  ->[K],
+			eps := 0.001[K],
+			i:{
+				T_air := air.temp->[K],
 				T_sky := (0.0552 * (T_air=>[]^(1.5))) => [K],
-				T_water := (top_water_temp + 273.15[deg_c]) =>[K],
-				h_air := (2.8[] + 3[] * air.wind=>[]) => [W, m-2, K-1],
+				T_water := top_water_temp ->[K],
+				h_air := (2.8 + 3 * air.wind=>[]) => [W, m-2, K-1],
 				h_sky := (1/(1/0.88+1/0.80-1)) * 5.67e-8 [W, m-2, K-4]* (T_front+T_sky)*(T_front^2 + T_sky^2),
 				B_back := surf.FPV.A_both + surf.FPV.h_water,
 				C_back := surf.FPV.h_water * T_water,
@@ -215,9 +281,10 @@ This is a modification to the AirSea module where you can cover a basin with flo
 				numer :=  B_front*B_back*(surf.FPV.A_both+surf.FPV.A_both)-surf.FPV.A_both*surf.FPV.A_both*B_front - surf.FPV.A_both*surf.FPV.A_both*B_back,
 				T_cell :=((B_front*B_back* (1- FPV_alb) * (1-FPV_eff)* air.g_rad  + surf.FPV.A_both * B_back *C_front  + surf.FPV.A_both * B_front *C_back ) / numer - (273.15[K])) =>[deg_c],
 				#surf.FPV.A_both * h_water / (surf.FPV.A_both + h_water) + surf.FPV.A_both / (surf.FPV.A_both + h_air + h_sky)*(h_air + h_sky *(Tc - T_sky)/(Tc - T_air))
-				T_front_final := ((surf.FPV.A_both =>[] * (T_cell=>[]+273.15[])=>[]) + C_front=>[]) / B_front=>[],
-
-			T_cell
+				T_front_final := ((surf.FPV.A_both =>[] * (T_cell->[K])=>[]) + C_front=>[]) / B_front=>[K],
+				T_cell										if abs(T_front - T_front_final) < eps,
+				{T_front <- T_front_final, iterate i}		otherwise
+			}	
 		} 	otherwise 
 	}