diff --git a/.gitignore b/.gitignore old mode 100644 new mode 100755 index e0978408..39c479d1 --- a/.gitignore +++ b/.gitignore @@ -41,3 +41,6 @@ Makefile # Runtime artifacts Raven Raven_errors.txt + +# Directories +test/ diff --git a/src/EnergyTransport.cpp b/src/EnergyTransport.cpp old mode 100644 new mode 100755 index 9c78c7d0..0dc25513 --- a/src/EnergyTransport.cpp +++ b/src/EnergyTransport.cpp @@ -270,47 +270,57 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub double A_new = pRes->GetSurfaceArea(); double A_old = pRes->GetOldSurfaceArea(); double A_avg = 0.5 * (A_new + A_old); + double V_h_new = pRes->GetHypolimnionStorage(); + double V_h_old = pRes->GetOldHypolimnionStorage(); + double V_e_new = V_new-V_h_new; + double V_e_old = V_old-V_h_old; + double A_h_new = pRes->GetMixingArea(); + double A_h_old = pRes->GetOldMixingArea(); + double A_h_avg = 0.5 * (A_h_new + A_h_old); CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); double Acorr=1.0; double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); - double hstar(0),ksed(0),Vsed=0.001; - double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_new); - double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p] / V_old); - double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / Vsed ); - double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p] / Vsed ); + + double hstar(0),ksed(0); //Vsed=0.001; + double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); + double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p] / V_e_old); + double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new ); + double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p] / V_h_old ); + + if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area - temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] - SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] - not using canopy correction! - LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] + temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] + SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] - not using canopy correction! + LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] LW =-STEFAN_BOLTZ*EMISS_WATER*0.5*(pow(T_new+ZERO_CELSIUS,4)+pow(T_old+ZERO_CELSIUS,4)); - AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] //*pHRU->GetArea()/A_avg + AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] //*pHRU->GetArea()/A_avg - hstar = pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] - Vsed = pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; - ksed = pRes->GetLakebedConductivity() / 0.5 / pRes->GetLakebedThickness();//[MJ/m2/d/K] + hstar = pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] + // Vsed = pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; + ksed = pRes->GetLakebedConductivity() / 0.5 / pRes->GetLakebedThickness(); //[MJ/m2/d/K] } - Q_sens =hstar* A_avg * (temp_air -0.5*(T_new+T_old)); - Q_cond =ksed * A_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); - Q_sw_in =(SW )*A_avg; - Q_lw_in =(LW_in )*A_avg; - Q_lw_out=(LW )*A_avg; - Q_lat =-(AET * DENSITY_WATER * LH_VAPOR)*A_avg; - Q_rain =_aMresRain[p]; + Q_sens =hstar* A_avg * (temp_air -0.5*(T_new+T_old)); + Q_cond =ksed * A_h_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); + Q_sw_in =(SW )*A_avg; + Q_lw_in =(LW_in )*A_avg; + Q_lw_out=(LW )*A_avg; + Q_lat =-(AET * DENSITY_WATER * LH_VAPOR)*A_avg; + Q_rain =_aMresRain[p]; return -(Q_sens + Q_cond + Q_sw_in + Q_lw_in + Q_lw_out + Q_lat + Q_rain) * tstep; //[MJ] } ////////////////////////////////////////////////////////////////// -/// \brief Calculates Res_mass (total reservoir energ) and ResSedMass (total sediment energy) in basin with reservoir at end of time step +/// \brief Calculates Res_mass (total reservoir energy) and ResSedMass (total sediment energy) in basin with reservoir at end of time step /// /// \param p [in] subbasin index /// \param aMout_new[] [in] Array of energy outflows at downstream end of each segment at end of current timestep [MJ/d] [size: nsegs] @@ -320,7 +330,7 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub /// \param tt [in] Time structure // void CEnthalpyModel::RouteMassInReservoir (const int p, // SB index - const double *aMout_new, // [MJ/d][size: nsegs ] + const double *aMout_new, // [MJ/d][size: nsegs] double &Res_mass, // [MJ] double &ResSedMass, // [MJ] const optStruct &Options, @@ -330,48 +340,172 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // CReservoir* pRes = _pModel->GetSubBasin(p)->GetReservoir(); int nSegments = _pModel->GetSubBasin(p)->GetNumSegments(); - double V_new=pRes->GetStorage (); - double V_old=pRes->GetOldStorage (); - double Q_new=pRes->GetOutflowRate ()*SEC_PER_DAY; - double Q_old=pRes->GetOldOutflowRate()*SEC_PER_DAY; - double A_new=pRes->GetSurfaceArea(); - double A_old=pRes->GetOldSurfaceArea(); - double A_avg=0.5*(A_new+A_old); - double tstep=Options.timestep; - + double V_new = pRes->GetStorage(); + double V_old = pRes->GetOldStorage(); + if((V_old<=0.0) || (V_new<=0.0)) { Res_mass=ResSedMass=0.0; return;} //handles dried out reservoir/lake + + double V_h_new = pRes->GetHypolimnionStorage(); + double V_h_old = pRes->GetOldHypolimnionStorage(); + double V_e_new = V_new-V_h_new; + double V_e_old = V_old-V_h_old; + double Q_new = pRes->GetOutflowRate() * SEC_PER_DAY; + double Q_old = pRes->GetOldOutflowRate() * SEC_PER_DAY; + double A_new = pRes->GetSurfaceArea(); + double A_old = pRes->GetOldSurfaceArea(); + double A_avg = 0.5 * (A_new + A_old); + double A_h_new = pRes->GetMixingArea(); + double A_h_old = pRes->GetOldMixingArea(); + double A_h_avg = 0.5 * (A_h_new + A_h_old); + + double tstep=Options.timestep; + + + double Q_dn_new(0),Q_up_new(0); + double Q_dn_old = 0.0; + double Q_up_old = 0.0; + + double Q_vert = (V_h_new - V_h_old)/tstep; + if (Q_vert > 0){ + Q_dn_new = Q_vert; Q_up_new = 0.0; + } + else + { + Q_dn_new = 0.0; Q_up_new = -Q_vert; + } + CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); - double Acorr=1.0; - - double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p]/V_old); + double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p]/V_e_old); + double Ts_old =ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/V_h_old); double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); - double hstar(0), ksed(0), Vsed=0.001; + double hstar(0); +//Vsed=0.001; + double u2(0); + if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input - Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area - - temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] - SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] - LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] + Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area + temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] + SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] + LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] + LW =-STEFAN_BOLTZ*EMISS_WATER*pow(T_old+ZERO_CELSIUS,4); //[MJ/m2/d] //TMP DEBUG -time-lagged - should include in the N-R formulation + AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] +// Vsed =pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; + hstar =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] +// ksed =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K] + } - LW =-STEFAN_BOLTZ*EMISS_WATER*pow(T_old+ZERO_CELSIUS,4); //[MJ/m2/d] //TMP DEBUG -time-lagged - should include in the N-R formulation - AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] - Vsed =pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; - hstar =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] - ksed =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K] +///////////// Option 1 to calculate kdiff +//// water density in each layer as a function of water temperature (Chapra, 2008) +//double dens_e,dens_h,vdiff,kdiff; +//double a0= 999.842594; +//double a1= 6.793952e-2; +//double a2=-9.095290e-3; +//double a3= 1.001685e-4; +//double a4=-1.120083e-6; +//double a5= 6.536332e-9; +//dens_e = a0+a1*T_old + a2*pow(T_old,2) + a3*pow(T_old,3) +a4*pow(T_old,4)+a5*pow(T_old,5); +//dens_h = a0+a1*Ts_old + a2*pow(Ts_old,2) + a3*pow(Ts_old,3) +a4*pow(Ts_old,4)+a5*pow(Ts_old,5); +//// Brunt-Väisälä frequency at thermocline depth (assume dz = 1.0m) (Jennings et al. 2012) +//double N2=1.0e-12; // value when dens_e >= dens_h +//if(dens_e kdiff = md*(Ke+Ked+Km); + u2 = pHRU->GetForcingFunctions()->wind_vel; + + double kstar = 0.0; + double wstar = 0.0; //Surface friction velocity [m/s] + double Ri = 0.0; //Richardson number + double N2 = 5.0E-12; //Brunt-Väisälä frequency when dens_e >= dens_h [1/s2] + double dens_e, dens_h; //Water density in each layer [kg/m3] + double Ke(0.0),Ked(0.0),Km(0.0),kdiff(0.0); //Diffusion coefficient [MJ/m2/K/d] + double md = 1.0; //Diffusion scaling factor in each layer + + //Water density in each layer as a function of water temperature. Reference: (Hostetler and Bartlein, 1990) + double a0 = -1.954e-5; + dens_e = DENSITY_WATER*(1.0+a0*pow(abs((T_old + ZERO_CELSIUS)-277.0),1.68)); + dens_h = DENSITY_WATER*(1.0+a0*pow(abs((Ts_old + ZERO_CELSIUS)-277.0),1.68)); + + //Estimate layer thicknesses + double Ze = V_e_old/A_old; //Average epilimnion thickness [m] + double Zh = V_h_old/A_h_old; //Average hppolimnion thickness [m] + double Zc = (Ze+Zh)/2.0; //Characteristic length [m] + + //Molecular diffusion coefficient, Km can be calculated as per Community Land Model v5 (Lawrence et al, 2020) + Km = TC_WATER/HCP_WATER; //[MJ/m2/K/d] + + //Estimate density gradient dp/dz over length Zc; ===> Brunt-Väisälä frequency across thermocline depth (Jennings et al. 2012) + if(dens_e FREEZING_TEMP){ + wstar = 0.0012*u2; //Surfaec friction velocity [m/s] + kstar = 6.6*pow(u2, -1.84) * sqrt(abs(sin(pHRU->GetLatRad()))); + Ri = 40.0*N2*pow(VON_KARMAN,2.0)*pow(Zc,2.0); + Ri /= pow(wstar,2.0); + Ri /= exp(-2*VON_KARMAN*Zc); + Ri += 1.0; + Ri = (-1 + sqrt(Ri))/20.0; + Ke = VON_KARMAN*wstar*Zc/(1.0 + 37.0*pow(Ri,2.0)) * exp(-VON_KARMAN*Zc); //[MJ/m2/K/s] + Ke *= SEC_PER_DAY; //[MJ/m2/K/d] + } + + //Enhanced diffusivity intended to represent unresolved mixing processes; ====> Community Land Model v5 (Lawrence et al, 2020) + if (N2 >= 7.5e-5){ + Ked = 1.04e-8 * pow(N2,-0.43)*SEC_PER_DAY; //[MJ/m2/K/d] + } + + //Increase the overall diffusivity for large lakes, intended to represent + //3-dimensional mixing processes such as cased by horizontal temperature; Gradient can be calculated as per Community Land Model v5 (Lawrence et al, 2020) + if ((Ze+Zh) >= 25.0){ + md = 10.0; //Diffusion scaling factor in each layer + } + + //Final diffusion coefficient + kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] + + + //TEMPORARY OUPUT ====> To test which parameter is abnormal + //cout<<"Ze = "< tolerance) { - T_guess = ConvertVolumetricEnthalpyToTemperature(E_guess/V_new); - Ts_guess= ConvertVolumetricEnthalpyToTemperature(Es_guess/Vsed); + T_guess = ConvertVolumetricEnthalpyToTemperature(E_guess/V_e_new); + Ts_guess= ConvertVolumetricEnthalpyToTemperature(Es_guess/V_h_new); - A[0][0] = 1.0 + 0.5 * tstep * Q_new/V_new; - A[1][1] = 1.0; - A[1][0] = 0.0; - A[0][1] = 0.0; + A[0][0] = 1.0 + 0.5*tstep*(Q_new + Q_dn_new)/V_e_new; + A[1][1] = 1.0 + 0.5*tstep*Q_up_new/V_h_new; + A[1][0]-= 0.5*tstep*Q_dn_new/V_e_new; + A[0][1]-= 0.5*tstep*Q_up_new/V_h_new; if (E_guess!=0.0){ - A[0][0]+= 0.5*tstep*hstar*A_new*T_guess/E_guess; - A[0][0]+= 0.5*tstep*ksed *A_avg*T_guess/E_guess; //[MJ/m2/d/K]*[m2]/[MJ/m3/K] = [m3/d] - A[1][0]-= 0.5*tstep*ksed *A_avg*T_guess/E_guess; + A[0][0]+= 0.5*tstep*hstar*A_avg* T_guess/E_guess; + A[0][0]+= 0.5*tstep*kdiff*A_h_avg*T_guess/E_guess; //[MJ/m2/d/K]*[m2]/[MJ/m3/K] = [m3/d] + A[1][0]-= 0.5*tstep*kdiff*A_h_avg*T_guess/E_guess; } if (Es_guess!=0.0){ - A[0][1]-= 0.5*tstep*ksed *A_avg*Ts_guess/Es_guess; - A[1][1]+= 0.5*tstep*ksed *A_avg*Ts_guess/Es_guess; + A[0][1]-= 0.5*tstep*kdiff*A_h_avg*Ts_guess/Es_guess; + A[1][1]+= 0.5*tstep*kdiff*A_h_avg*Ts_guess/Es_guess; } //J_ij=A_ij+E*dA_ij/dE = Jacobian - J[0][0]=A[0][0]+0.5*tstep*(hstar*A_new+ksed*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_new)/ V_new); - J[0][1]=A[0][1]-0.5*tstep*( ksed*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/Vsed )/ Vsed ); - J[1][0]=A[1][0]-0.5*tstep*( ksed*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_new)/ V_new); - J[1][1]=A[1][1]+0.5*tstep*( ksed*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/Vsed )/ Vsed ); + J[0][0]=A[0][0]+0.5*tstep*(hstar*A_avg+kdiff*A_h_avg) * (TemperatureEnthalpyDerivative(E_guess /V_e_new)/V_e_new); + J[0][1]=A[0][1]-0.5*tstep*( kdiff*A_h_avg) * (TemperatureEnthalpyDerivative(Es_guess/V_h_new)/V_h_new ); + J[1][0]=A[1][0]-0.5*tstep*( kdiff*A_h_avg) * (TemperatureEnthalpyDerivative(E_guess /V_e_new)/V_e_new); + J[1][1]=A[1][1]+0.5*tstep*( kdiff*A_h_avg) * (TemperatureEnthalpyDerivative(Es_guess/V_h_new)/V_h_new ); - R[0] =-A[0][0]*E_guess-A[0][1]*Es_guess+B[0]; - R[1] =-A[1][0]*E_guess-A[1][1]*Es_guess+B[1]; + R[0] = -A[0][0]*E_guess-A[0][1]*Es_guess+B[0]; + R[1] = -A[1][0]*E_guess-A[1][1]*Es_guess+B[1]; //invert 2x2 matrix analytically double den= (J[0][1]*J[1][0]-J[0][0]*J[1][1]); dE = (J[0][1]*R[1]-J[1][1]*R[0])/den; dEs = (J[1][0]*R[0]-J[0][0]*R[1])/den; - change =sqrt(dE*dE+dEs*dEs)/ V_new / HCP_WATER; //convert to approx temp difference (for >0C water) + + change =sqrt(dE*dE+dEs*dEs)/V_new/HCP_WATER; //convert to approx temp difference (for >0C water) E_guess +=dE; Es_guess +=dEs; @@ -454,6 +590,7 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // delete [] R; for (int i=0;i<2;i++){delete [] A[i]; delete[] J[i]; } delete[] A; delete[] J; } + ////////////////////////////////////////////////////////////////// /// \brief returns total energy lost from subbasin reach over current time step [MJ] /// \param p [in] subbasin index @@ -1173,8 +1310,10 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct if ((pRes != NULL) && (pRes->GetHRUIndex() != DOESNT_EXIST)) { double HRUarea= _pModel->GetHydroUnit(pRes->GetHRUIndex())->GetArea()*M2_PER_KM2; - double V_sed = pRes->GetLakebedThickness() * HRUarea; + //double V_sed = pRes->GetLakebedThickness() * HRUarea; double V_new = pRes->GetStorage(); + double V_h_new = pRes->GetHypolimnionStorage(); + double V_e_new = V_new-V_h_new; mult = 1.0 / HRUarea; @@ -1185,20 +1324,20 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct GetEnergyLossesFromLake(p,Q_sens,Q_cond,Q_lat,Q_sw_in,Q_lw_in,Q_lw_out,Q_rain); - double lakeTemp =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_new); - double sedTemp =ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_sed); + double lakeTemp =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); + double sedTemp =ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new); double inflowTemp=ConvertVolumetricEnthalpyToTemperature( Ein / Qin ); - double pctFroz =ConvertVolumetricEnthalpyToIceContent (_aMres[p] / V_new); + double pctFroz =ConvertVolumetricEnthalpyToIceContent (_aMres[p] / V_e_new); _LAKEOUT << mult * Ein << "," << mult * Eout << ","; _LAKEOUT << mult * Q_rain << "," << mult * Q_sens << ","; _LAKEOUT << mult * Q_cond << "," << mult * Q_lat << ","; _LAKEOUT << mult * Q_sw_in << "," << mult * Q_lw_in << "," << mult * Q_lw_out << ","; _LAKEOUT << mult * _aMres[p] << ","; - if (V_new!=0){_LAKEOUT << lakeTemp << ",";}else{_LAKEOUT<<",";} - if (V_sed!=0){_LAKEOUT << sedTemp << ",";}else{_LAKEOUT<<",";} - if (Qin !=0){_LAKEOUT << inflowTemp << ",";}else{_LAKEOUT<<",";} - if (V_new!=0){_LAKEOUT << pctFroz << ",";}else{_LAKEOUT<<",";} + if (V_e_new!=0){_LAKEOUT << lakeTemp << ",";}else{_LAKEOUT<<",";} + if (V_h_new!=0){_LAKEOUT << sedTemp << ",";}else{_LAKEOUT<<",";} + if (Qin !=0){_LAKEOUT << inflowTemp << ",";}else{_LAKEOUT<<",";} + if (V_e_new!=0){_LAKEOUT << pctFroz << ",";}else{_LAKEOUT<<",";} } } } diff --git a/src/Makefile b/src/Makefile old mode 100644 new mode 100755 index 455c92a6..ee09ffe1 --- a/src/Makefile +++ b/src/Makefile @@ -25,12 +25,12 @@ LDFLAGS := CXXFLAGS += -std=c++11 -fPIC # OPTION 1) include netcdf - uncomment following two commands (assumes netCDF path = /usr/local): -#CXXFLAGS += -Dnetcdf -#LDLIBS += -L/usr/local -lnetcdf +CXXFLAGS += -Dnetcdf +LDLIBS += -L/usr/local -lnetcdf # OPTION 2) include lp_solve for water management optimization- uncomment following two commands (assumes liblpsolve55.a in ../lib/lp_solve_unix/ folder): -CXXFLAGS += -D_LPSOLVE_ -LDLIBS += -L../lib/lp_solve_unix -llpsolve55 +#CXXFLAGS += -D_LPSOLVE_ +#LDLIBS += -L../lib/lp_solve_unix -llpsolve55 # OPTION 3) if you use a OSX/BSD system, uncomment the LDFLAGS line below # this is to allow for use a 1Gb stack, see http://linuxtoosx.blogspot.ca/2010/10/stack-overflow-increasing-stack-limit.html diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp old mode 100644 new mode 100755 index 3043dae3..abfb133e --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -18,7 +18,7 @@ void CReservoir::BaseConstructor(const string Name,const long SubID) _SBID=SubID; _pDownSB=NULL; - + _mixing_depth=5.0; _lakebed_thick=1.0; _lakebed_cond =0.0; _lake_convcoeff=2.0; @@ -403,6 +403,14 @@ long CReservoir::GetSubbasinID () const { return _SBID; } // double CReservoir::GetStorage () const { return GetVolume(_stage); } +/// \returns reservoir hypoliminion storage [m3] +// +double CReservoir::GetHypolimnionStorage () const { return GetVolume(_stage-_mixing_depth); } + +/// \returns reservoir old hypoliminion storage [m3] +// +double CReservoir::GetOldHypolimnionStorage () const { return GetVolume(_stage_last-_mixing_depth); } + ////////////////////////////////////////////////////////////////// /// \returns current stage [m] // @@ -449,6 +457,16 @@ double CReservoir::GetSurfaceArea () const { return GetArea(_stage); } // double CReservoir::GetOldSurfaceArea () const { return GetArea(_stage_last); } +////////////////////////////////////////////////////////////////// +/// \returns current mixing area [m2] +// +double CReservoir::GetMixingArea () const { return GetArea(_stage-_mixing_depth); } + +////////////////////////////////////////////////////////////////// +/// \returns current old mixing area [m2] +// +double CReservoir::GetOldMixingArea () const { return GetArea(_stage_last-_mixing_depth); } + ////////////////////////////////////////////////////////////////// /// \returns lakebed thickness [m] // diff --git a/src/Reservoir.h b/src/Reservoir.h old mode 100644 new mode 100755 index 03894f57..45fb8528 --- a/src/Reservoir.h +++ b/src/Reservoir.h @@ -55,6 +55,7 @@ class CReservoir string _name; ///< reservoir name long _SBID; ///< subbasin ID double _max_capacity; ///< maximum reservoir capacity [m3] + double _mixing_depth; double _lakebed_thick; ///< lakebed thickness [m] double _lakebed_cond; ///< lakebed thermal conductivity [MJ/m/K/d] @@ -178,6 +179,8 @@ class CReservoir double GetStorage () const; //[m3] double GetOldStorage () const; //[m3] + double GetHypolimnionStorage () const; //[m3] + double GetOldHypolimnionStorage () const; //[m3] double GetOutflowRate () const; //[m3/s] double GetOldOutflowRate () const; //[m3/s] double GetIntegratedOutflow (const double &tstep) const; //[m3] @@ -191,6 +194,8 @@ class CReservoir double GetOldStage () const; //[m] double GetSurfaceArea () const; //[m2] double GetOldSurfaceArea () const; //[m2] + double GetMixingArea () const; //[m2] + double GetOldMixingArea () const; //[m2] double GetLakebedThickness () const; //[m] double GetLakebedConductivity () const; //[MJ/m/K/d] double GetLakeConvectionCoeff () const; //[MJ/m2/d/K] diff --git a/src/Transport.h b/src/Transport.h index 4d0514a0..4b7c0178 100644 --- a/src/Transport.h +++ b/src/Transport.h @@ -217,6 +217,11 @@ class CConstituentModel double *_aMout_res; ///< array storing reservoir mass outflow [mg/d] or heat loss [MJ/d] [size: nSubBasins] double *_aMout_res_last; ///< array storing reservoir mass outflow [mg/d] or enthalpy outflow [MJ/d] at start of timestep [size: nSubBasins] double *_aMresRain; ///< array storing reservoir rain inputs [mg/d] or enthalpy input [MJ/d] [size: nSubBasins] + + double *_aMres_epo; ///< array used for storing epolimnion layer masses [mg] or enthalpy [MJ] [size: nSubBasins] + double *_aMres_epo_last; ///< array used for storing epolimnion layer masses [mg] at start of timestep [size: nSubBasins] + double *_aMres_hyp; ///< array used for storing hypolimnion layer masses [mg] or enthalpy [MJ] [size: nSubBasins] + double *_aMres_hyp_last; ///< array used for storing hypolimnion layer masses [mg] at start of timestep [size: nSubBasins] // Mass balance tracking variables double *_channel_storage; ///< array storing channel storage [mg] or [MJ] [size: nSubBasins]