From f8155cabc99453148df9d751515bc0e8789d6571 Mon Sep 17 00:00:00 2001
From: Nander Wever <nander.wever@slf.ch>
Date: Wed, 5 Jun 2024 18:51:46 +0200
Subject: [PATCH] Adding snow depth and mass corrections from inflate/deflate
 to smet output

---
 applications/snowpack/Main.cc | 11 +++++++++--
 snowpack/plugins/SmetIO.cc    | 22 +++++++++++++++++++++-
 snowpack/plugins/SmetIO.h     |  2 +-
 3 files changed, 31 insertions(+), 4 deletions(-)

diff --git a/applications/snowpack/Main.cc b/applications/snowpack/Main.cc
index d641dbc2..ee639bae 100644
--- a/applications/snowpack/Main.cc
+++ b/applications/snowpack/Main.cc
@@ -88,6 +88,7 @@ class Cumsum {
 
 		double precip;
 		double drift, snow, runoff, rain;
+		double dhs_corr, mass_corr; // inflate/deflate variables
 		vector<double> erosion; // Cumulated eroded mass; dumped to file as rate
 };
 
@@ -197,7 +198,7 @@ void Slope::setSlope(const unsigned int slope_sequence, vector<SnowStation>& vec
 
 Cumsum::Cumsum(const unsigned int nSlopes)
         : precip(0.),
-          drift(0.), snow(0.), runoff(0.), rain(0.),
+          drift(0.), snow(0.), runoff(0.), rain(0.), dhs_corr(0.), mass_corr(0.),
           erosion(nSlopes, 0.)
 {}
 
@@ -1302,6 +1303,9 @@ inline void real_main (int argc, char *argv[])
 						if (slope.mainStationDriftIndex)
 							cumsum.drift = 0.;
 						surfFluxes.hoar = 0.;
+						// Inflate/deflate sums
+						cumsum.dhs_corr += qr_Hdata.at(i_hz).dhs_corr;
+						cumsum.mass_corr += qr_Hdata.at(i_hz).mass_corr;
 					}
 					// New snow water equivalent (kg m-2), rain was dealt with in Watertransport.cc
 					surfFluxes.mass[SurfaceFluxes::MS_HNW] += vecXdata[slope.mainStation].hn
@@ -1362,8 +1366,11 @@ inline void real_main (int argc, char *argv[])
 					if (slope.mainStationDriftIndex)
 						i_hz0 = i_hz;
 					const double wind_trans24 = (slope.sector == slope.mainStation) ? qr_Hdata.at(i_hz0).wind_trans24 : qr_Hdata.at(i_hz).wind_trans24;
+					ProcessDat tmpHdata = qr_Hdata.at(i_hz);			// Temporary ProcessDat object for output
+					tmpHdata.dhs_corr = cumsum.dhs_corr; cumsum.dhs_corr = 0.;	// overwrite the inflate/deflate variables
+					tmpHdata.mass_corr = cumsum.mass_corr; cumsum.mass_corr = 0.;
 					snowpackio.writeTimeSeries(vecXdata[slope.sector], surfFluxes, Mdata,
-					                           qr_Hdata.at(i_hz), wind_trans24);
+					                           tmpHdata, wind_trans24);
 
 					if (avgsum_time_series) {
 						surfFluxes.reset(cumsum_mass);
diff --git a/snowpack/plugins/SmetIO.cc b/snowpack/plugins/SmetIO.cc
index 019d716c..407ca349 100644
--- a/snowpack/plugins/SmetIO.cc
+++ b/snowpack/plugins/SmetIO.cc
@@ -150,7 +150,7 @@ SmetIO::SmetIO(const SnowpackConfig& cfg, const RunInfo& run_info)
           in_dflt_TZ(0.), calculation_step_length(0.), ts_days_between(0.), min_depth_subsurf(0.),
           avgsum_time_series(false), useCanopyModel(false), useSoilLayers(false), research_mode(false), perp_to_slope(false), haz_write(true), useReferenceLayer(false),
           out_heat(false), out_lw(false), out_sw(false), out_meteo(false), out_haz(false), out_mass(false), out_t(false),
-          out_load(false), out_stab(false), out_canopy(false), out_soileb(false), useRichardsEq(false), enable_pref_flow(false), enable_ice_reservoir(false), read_dsm(false)
+          out_load(false), out_stab(false), out_canopy(false), out_soileb(false), out_inflate(false), useRichardsEq(false), enable_pref_flow(false), enable_ice_reservoir(false), read_dsm(false)
 {
 	cfg.getValue("TIME_ZONE", "Input", in_dflt_TZ);
 	cfg.getValue("CANOPY", "Snowpack", useCanopyModel);
@@ -192,6 +192,7 @@ SmetIO::SmetIO(const SnowpackConfig& cfg, const RunInfo& run_info)
 	cfg.getValue("OUT_MASS", "Output", out_mass);
 	cfg.getValue("OUT_METEO", "Output", out_meteo);
 	cfg.getValue("OUT_SOILEB", "Output", out_soileb);
+	cfg.getValue("INFLATE_ALLOW", "Snowpack", out_inflate, IOUtils::nothrow);
 	cfg.getValue("OUT_STAB", "Output", out_stab);
 	cfg.getValue("OUT_SW", "Output", out_sw);
 	cfg.getValue("OUT_T", "Output", out_t);
@@ -259,6 +260,7 @@ SmetIO& SmetIO::operator=(const SmetIO& source) {
 		out_stab = source.out_stab;
 		out_canopy = source.out_canopy;
 		out_soileb = source.out_soileb;
+		out_inflate = source.out_inflate;
 		useRichardsEq = source.useRichardsEq;
 		enable_pref_flow = source.enable_pref_flow;
 		enable_ice_reservoir = source.enable_ice_reservoir;
@@ -965,6 +967,8 @@ std::string SmetIO::getFieldsHeader(const SnowStation& Xdata) const
 			if (useSoilLayers) os << "Lateral_flow_soil" << " ";
 		}
 	}
+	if (out_inflate)
+		os << "dHS_corr dMass_corr" << " "; //snow depth (cm) and mass correction (kg m-2) from inflate/deflate
 	if (out_load)
 		os << "load "; //Solute load at ground surface
 	if (out_t && !fixedPositions.empty()) {
@@ -1091,6 +1095,16 @@ void SmetIO::writeTimeSeriesHeader(const SnowStation& Xdata, const double& tz, s
 			}
 		}
 	}
+	if (out_inflate) {
+		//"dHS_corr dMass_corr"
+		plot_description << "snow_depth_correction_inflate_deflate  mass_correction_inflate_deflate" << " "; //snow depth (cm) and mass correction (kg m-2) from inflate/deflate
+		plot_units << "cm kg m-2" << " ";
+		units_offset << "0 0" << " ";
+		units_multiplier << "1 1" << " ";
+		plot_color << "0xa503fc 0x03bafc" << " ";
+		plot_min << "" << " ";
+		plot_max << "" << " ";
+	}
 	if (out_load) {
 		//"load"
 		plot_description << "solute_load" << " ";
@@ -1289,6 +1303,12 @@ void SmetIO::writeTimeSeriesData(const SnowStation& Xdata, const SurfaceFluxes&
 		}
 	}
 
+	if (out_inflate) {
+		// snow depth (cm) and mass correction (kg m-2) from inflate/deflate
+		data.push_back( M_TO_CM(Hdata.dhs_corr) ); vec_precision.push_back(dflt_precision); vec_width.push_back(dflt_width);
+		data.push_back( Hdata.mass_corr ); vec_precision.push_back(dflt_precision); vec_width.push_back(dflt_width);
+	}
+
 	if (out_load) {
 		if (!Sdata.load.empty())
 			data.push_back( Sdata.load[0] );
diff --git a/snowpack/plugins/SmetIO.h b/snowpack/plugins/SmetIO.h
index 9a2f9890..8e3c4ef2 100644
--- a/snowpack/plugins/SmetIO.h
+++ b/snowpack/plugins/SmetIO.h
@@ -91,7 +91,7 @@ class SmetIO : public SnowpackIOInterface {
 		double min_depth_subsurf;
 		bool avgsum_time_series, useCanopyModel, useSoilLayers, research_mode, perp_to_slope, haz_write;
 		bool useReferenceLayer;		//Whether or not the output should be referenced to the marked reference layer (i.e., the layer with int(mk/1000)==9).
-		bool out_heat, out_lw, out_sw, out_meteo, out_haz, out_mass, out_t, out_load, out_stab, out_canopy, out_soileb;
+		bool out_heat, out_lw, out_sw, out_meteo, out_haz, out_mass, out_t, out_load, out_stab, out_canopy, out_soileb, out_inflate;
 		bool useRichardsEq;
 		bool enable_pref_flow;
 		bool enable_ice_reservoir;
-- 
GitLab