Skip to content

Commit 37b761b

Browse files
analytophileanalytophile
committed
v482 - reversion to original lateral flow treatment in Thermal transport [BENCHMARKED]
Separated in catchment and in-channel routing (Transport.h; Solvers.cpp; MassRouting.cppRepo) Reporting of lateral flow energy budget component Disabled the q_lateral treatment (EnergyTransport.cpp) for loop lower index bug in InCatchment routing
1 parent 7adc224 commit 37b761b

File tree

6 files changed

+87
-42
lines changed

6 files changed

+87
-42
lines changed

‎src/ConstituentModel.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -730,7 +730,7 @@ void CConstituentModel::WriteOutputFileHeaders(const optStruct &Options)
730730
_OUTPUT<<", sensible (r) "<<kg<<", conductive (r) " << kg << ", latent (r) "<<kg<<", GW mixing (r) "<<kg<<", radiant (r) "<<kg<<", friction (r)"<<kg<< endl;//TMP DEBUG
731731
}
732732

733-
//Pollutograph / stream temperatures file
733+
// Pollutograph / stream temperatures file
734734
//--------------------------------------------------------------------
735735
if(!_is_passive) {
736736
if(_type!=ENTHALPY) { filename=_name+"Pollutographs.csv"; }
@@ -1262,14 +1262,14 @@ void CConstituentModel::WriteMinorOutput(const optStruct &Options,const time_str
12621262

12631263
double latent_flux=0;
12641264
double Q_sens(0.0),Q_cond(0.0),Q_lat(0.0),Q_GW(0.0),Q_rad(0.0),Q_fric(0.0);
1265-
double Qs,Qc,Ql,Qg,Qr,Qlw,Qlw_in,Qf,Tave;
1265+
double Qs,Qc,Ql,Qg,Qr,Qlw,Qlw_in,Qltrl,Qf,Tave;
12661266
if(_type==ENTHALPY)
12671267
{
12681268
pEnthalpyModel=(CEnthalpyModel*)(this);
12691269
latent_flux =pEnthalpyModel->GetAvgLatentHeatFlux()*(area*M2_PER_KM2)*Options.timestep; // [MJ] (loss term)t)
12701270

12711271
for(int p=0;p<_pModel->GetNumSubBasins();p++) {
1272-
pEnthalpyModel->GetEnergyLossesFromReach(p,Qs,Qc,Ql,Qg,Qr,Qlw_in,Qlw,Qf,Tave);
1272+
pEnthalpyModel->GetEnergyLossesFromReach(p,Qs,Qc,Ql,Qg,Qr,Qlw_in,Qlw,Qltrl,Qf,Tave);
12731273
Q_sens+=Qs;Q_cond+=Qc;Q_lat+=Ql;Q_GW+=Qg;Q_rad+=Qr+Qlw+Qlw_in; Q_fric+=Qf;
12741274
}
12751275
}

‎src/EnergyTransport.cpp‎

Lines changed: 48 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,8 @@ CEnthalpyModel::CEnthalpyModel(CModel *pMod,CTransportModel *pTMod,string name,c
2828
_anyGaugedLakes=false;
2929

3030
_lateral_via_convol=true;
31+
32+
_lateral_via_convol=false; //TMP DEBUG
3133
}
3234

3335
//////////////////////////////////////////////////////////////////
@@ -292,7 +294,7 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub
292294

293295
LW =-STEFAN_BOLTZ*EMISS_WATER*0.5*(pow(T_new+ZERO_CELSIUS,4)+pow(T_old+ZERO_CELSIUS,4));
294296

295-
AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] //*pHRU->GetArea()/A_avg
297+
AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] //*pHRU->GetArea()/A_avg
296298

297299
hstar = pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K]
298300
Vsed = pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2;
@@ -461,10 +463,10 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, //
461463
//
462464
double CEnthalpyModel::GetNetReachLosses(const int p)
463465
{
464-
double Q_sens,Q_cond,Q_lat,Q_GW,Q_rad_in,Q_lw_in,Q_lw_out, Q_fric,Qtot,Tave;
466+
double Q_sens,Q_cond,Q_lat,Q_GW,Q_rad_in,Q_lw_in,Q_lw_out, Q_lateral,Q_fric,Qtot,Tave;
465467
double Q=_pModel->GetSubBasin(p)->GetOutflowArray()[_pModel->GetSubBasin(p)->GetNumSegments()-1];
466468
if (Q<REAL_SMALL){return 0;}
467-
Qtot=GetEnergyLossesFromReach(p,Q_sens,Q_cond,Q_lat,Q_GW,Q_rad_in,Q_lw_in,Q_lw_out,Q_fric,Tave);
469+
Qtot=GetEnergyLossesFromReach(p,Q_sens,Q_cond,Q_lat,Q_GW,Q_rad_in,Q_lw_in,Q_lw_out,Q_lateral,Q_fric,Tave);
468470
_aTave_reach[p]=Tave;
469471
return Qtot;
470472
}
@@ -528,8 +530,8 @@ void CEnthalpyModel::Initialize(const optStruct& Options)
528530
L=_pModel->GetSubBasin(p)->GetReachLength();
529531
A=_pModel->GetSubBasin(p)->GetReferenceXSectArea();
530532
tr=L*A/Q/SEC_PER_DAY;
531-
if (tr<Options.timestep){_aMinResTime [p]=tr;}
532-
else {_aMinResTime [p]=1e-6; }
533+
if (tr<0.33*Options.timestep){_aMinResTime [p]=tr;}
534+
else {_aMinResTime [p]=0.333*Options.timestep; }
533535
}
534536

535537
// initialize stream temperatures if init_stream_temp is given
@@ -591,18 +593,22 @@ void CEnthalpyModel::Initialize(const optStruct& Options)
591593
//--------------------------------------------------------------------
592594
for(int p=0;p<nSB;p++) {
593595
PrepareForRouting(p);
596+
PrepareForInCatchmentRouting(p);
594597
}
595598

596599
if (Options.noisy){cout<<" ...Done initializing Energy Transport"<<endl;}
597600
}
601+
void CEnthalpyModel::PrepareForInCatchmentRouting(const int p)
602+
{
603+
UpdateCatchmentEnergySourceTerms(p);
604+
}
598605
//////////////////////////////////////////////////////////////////
599606
/// \brief Updates source terms for energy balance on subbasin reaches each time step
600607
/// \param p subbasin index
601608
//
602609
void CEnthalpyModel::PrepareForRouting(const int p)
603610
{
604611
UpdateReachEnergySourceTerms(p);
605-
UpdateCatchmentEnergySourceTerms(p);
606612
}
607613
//////////////////////////////////////////////////////////////////
608614
/// \brief Updates source terms for energy balance on subbasin reaches each time step
@@ -661,29 +667,31 @@ void CEnthalpyModel::UpdateReachEnergySourceTerms(const int p)
661667
double dbar =pBasin->GetRiverDepth(); //ensured to be >0
662668
double Ax =pBasin->GetXSectArea(); //ensured to be >0
663669
double L =max(pBasin->GetReachLength(),1.0);
664-
double qlat =pBasin->GetIntegratedLocalOutflow(tstep)/L; //total
665-
double qhlat =0.5*(_aMlocal[p] + _aMlocLast[p]) / L; //q_lat*h_lat
670+
double qlat =pBasin->GetIntegratedLocalOutflow(tstep)/L/tstep; //total [m3/d/m]
671+
double qhlat =0.5*(_aMlocal[p] + _aMlocLast[p]) / L; //q_lat*h_lat [MJ/d/m]
666672

667673
double kbed =_aKbed[p];//*bed_ratio? //[MJ/m2/d/K]
668674
double klin =4.0*STEFAN_BOLTZ*EMISS_WATER*pow(temp_lin,3.0);
669675
double kprime =qmix*bed_ratio*HCP_WATER; //[MJ/m2/d/K]
670676

671677
double Qf =GetReachFrictionHeat(pBasin->GetOutflowArray()[pBasin->GetNumSegments()-1],pBasin->GetBedslope(),pBasin->GetTopWidth());//[MJ/m2/d]
672678

679+
if (!_lateral_via_convol) {qhlat=0.0;qlat=0;}
680+
673681
double S(0.0); //source term [MJ/m3/d]
674682
S+=(SW+LW_in); //net incoming energy term
675683
S-=AET*DENSITY_WATER*LH_VAPOR; //latent heat flux term
676684
S+=Qf; //friction-generated heating term
677685
S+=hstar*temp_air; //convection with atmosphere
678686
S+=kbed* temp_bed; //conductive exchange with bed
679687
S-=klin*(ZERO_CELSIUS-0.75*temp_lin); //linearized outgoing radiation
680-
S+=(dbar/Ax)*qhlat; //lateral inflow - must add in update
688+
S+=(dbar/Ax)*qhlat; //lateral inflow - must add in update [MJ/m2/d]
681689
S+=kprime*temp_GW; //hyporheic mixing
682690

683691
S/=dbar;
684692

685693
//cout<<S<<" Qf:"<<Qf<<" h*:"<<hstar<<" Tlin:"<<temp_lin<<" qhlat:"<<qhlat<<" radin:"<<(SW+LW_in)<<" AET: "<<AET<<" kprime:"<<kprime<<" TGW:"<<temp_GW<<" Tbed:"<<temp_bed<<" Tai:"<<temp_air<<" dAx:"<<(dbar/Ax)<<endl;
686-
_aEnthalpyBeta[p]=(hstar + kbed + klin + qlat*dbar/Ax*HCP_WATER + kprime)/dbar/HCP_WATER;
694+
_aEnthalpyBeta[p]=(hstar + kbed + klin + qlat*(dbar/Ax)*HCP_WATER + kprime)/dbar/HCP_WATER;
687695

688696
for(int n=_nMinHist[p] - 1; n>0; n--) {
689697
_aEnthalpySource[p][n]=_aEnthalpySource[p][n-1];
@@ -735,19 +743,20 @@ double CEnthalpyModel::FunkyTemperatureIntegral(const int k, const int m, const
735743
/// \param Q_lat [out] energy gain from latent heat exchange [MJ/d]
736744
/// \param Q_GW [out] energy gain from groundwater mixing [MJ/d]
737745
/// \param Q_rad_in [out] energy gain from SW/LW radiation at surface [MJ/d]
738-
/// \param Q_lw_out [out] energy gain (negative) from surface losses [MJ/d]
746+
/// \param Q_lw_out [out] energy gain (always negative) from LW surface losses [MJ/d]
747+
/// \param Q_lateral[out] energy gain from lateral exchange term [MJ/d]
739748
/// \param Q_fric [out] energy gain from friction [MJ/d]
740749
/// \returns total energy lost from reach over current time step [MJ]
741750
//
742-
double CEnthalpyModel::GetEnergyLossesFromReach(const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_GW,double &Q_rad_in,double &Q_lw_in,double &Q_lw_out, double &Q_fric, double &Tave) const
751+
double CEnthalpyModel::GetEnergyLossesFromReach(const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_GW,double &Q_rad_in,double &Q_lw_in,double &Q_lw_out, double &Q_lateral, double &Q_fric, double &Tave) const
743752
{
744753
double tstep=_pModel->GetOptStruct()->timestep;
745754

746755
const CSubBasin *pBasin=_pModel->GetSubBasin(p);
747756

748757
if(pBasin->IsHeadwater())
749758
{
750-
Q_sens=Q_cond=Q_lat=Q_GW=Q_rad_in=Q_lw_in=Q_lw_out=Q_fric=0.0;
759+
Q_sens=Q_cond=Q_lat=Q_GW=Q_rad_in=Q_lw_in=Q_lw_out=Q_fric=Q_lateral=0.0;
751760
Tave=0.0;
752761
return 0.0;
753762
}
@@ -773,6 +782,12 @@ double CEnthalpyModel::GetEnergyLossesFromReach(const int p,double &Q_sens,doubl
773782
double dbar =pBasin->GetRiverDepth(); //averaged depth [m]
774783
double Ax =pBasin->GetXSectArea(); //[m2]
775784
double As =pBasin->GetTopWidth()*max(pBasin->GetReachLength(),1.0); //[m2]
785+
double L =max(pBasin->GetReachLength(),1.0);
786+
double qlat =pBasin->GetIntegratedLocalOutflow(tstep)/L; //total
787+
double qhlat =0.5*(_aMlocal[p] + _aMlocLast[p]) / L; //q_lat*h_lat
788+
789+
double temp_lat =qhlat/qlat/HCP_WATER;
790+
if (qlat <= 0) {temp_lat=0;}
776791

777792
double kbed =_aKbed[p]; //[MJ/m2/d/K]
778793
double klin =4.0*STEFAN_BOLTZ*EMISS_WATER*pow(temp_lin,3.0); //[MJ/m2/d/K]
@@ -791,7 +806,7 @@ double CEnthalpyModel::GetEnergyLossesFromReach(const int p,double &Q_sens,doubl
791806
if( aQin[m]<PRETTY_SMALL) { hin_hist[m]=0.0; }
792807
}
793808

794-
Q_sens = Q_cond = Q_lat = Q_GW = Q_rad_in = Q_lw_out = Q_fric =0.0;
809+
Q_sens = Q_cond = Q_lat = Q_GW = Q_rad_in = Q_lw_out = Q_lateral = Q_fric =0.0;
795810

796811
//Lagrangian terms (depend upon flowpath temperature):
797812
double dA,Tbar_km;
@@ -819,11 +834,15 @@ double CEnthalpyModel::GetEnergyLossesFromReach(const int p,double &Q_sens,doubl
819834
Q_sens +=dA*hstar *(temp_air-Tbar_km); //[m2]*[MJ/m2/d/K]*[K]=[MJ/d]
820835
Q_cond +=dA*kbed *(temp_bed-Tbar_km);
821836
Q_lw_out +=dA*klin *(0.75*temp_lin-(Tbar_km+ZERO_CELSIUS)); //linearized - works except at T~0
837+
Q_lateral+=dA*qlat *dbar/Ax*HCP_WATER*(temp_lat-Tbar_km);
838+
822839
temp_average+=dA/As*Tbar_km;
823840
if (As == 0.0) {temp_average=0.0;}
824841
}
825842
}
826-
843+
if (!_lateral_via_convol) {
844+
Q_lateral=0.0;
845+
}
827846
//Eulerian terms:
828847
Q_rad_in=SW*As; //[MJ/d]
829848
Q_lw_in =LW_in*As; //[MJ/d]
@@ -834,7 +853,7 @@ double CEnthalpyModel::GetEnergyLossesFromReach(const int p,double &Q_sens,doubl
834853

835854
delete[] hin_hist;
836855

837-
return -(Q_sens+Q_cond+Q_GW+Q_rad_in+Q_lw_in+Q_lw_out+Q_lat+Q_fric)*tstep;
856+
return -(Q_sens+Q_cond+Q_GW+Q_rad_in+Q_lw_in+Q_lw_out+Q_lat+Q_lateral+Q_fric)*tstep;
838857
}
839858
//////////////////////////////////////////////////////////////////
840859
/// \brief Calculates individual energy gain terms during in-catchment transit over current time step
@@ -913,13 +932,13 @@ double CEnthalpyModel::ApplyInCatchmentRouting (const int p,
913932
double beta=_aInCatch_a[p]+_aInCatch_b[p];
914933
beta=max(beta,1e-9); //to avoid divide by zero error
915934

916-
double Mout_new_test=0;
935+
//double Mout_new_test=0;
917936

918937
double Mout_new=0.0;
919938
double dt=tstep;
920939
double term1,term2,term3;
921940
term3=(1.0-exp(-beta*dt))/beta; //has limit dt as beta->0
922-
for(i=1;i<nMlatHist;i++)
941+
for(i=0;i<nMlatHist;i++)
923942
{
924943
term1=aMlatHist[i]*exp(-beta*i*dt);
925944
term2=0.0;
@@ -928,7 +947,7 @@ double CEnthalpyModel::ApplyInCatchmentRouting (const int p,
928947
}
929948
term2*=(aQlatHist[i]*SEC_PER_DAY);
930949
Mout_new+=aUnitHydro[i]*(term1+term2);
931-
Mout_new_test+=aUnitHydro[i]*aMlatHist[i];
950+
//Mout_new_test+=aUnitHydro[i]*aMlatHist[i];
932951
}
933952
//cout<<"IN CATCHMENT: "<<_aInCatch_b[p]<<" "<<Mout_new<<" "<<Mout_new_test<<" "<<(Mout_new-Mout_new_test)/Mout_new_test*100<<"%"<<endl;
934953

@@ -971,7 +990,11 @@ void CEnthalpyModel::ApplyConvolutionRouting(const int p,const double *aRoute
971990
term2+=_aEnthalpySource[p][i-j]*exp(-beta*(i-j)*dt)*term3;
972991
}
973992
term2*=(aQinHist[i]*SEC_PER_DAY);
993+
//needs to be aQinHist[i]+sum_0..i Qlathist[i]
974994
aMout_new[nSegments-1]+=aRouteHydro[i]*(term1+term2);
995+
//PROBLEM IS HERE - if aQinHist is zero, then the source term goes to zero
996+
//cout << "* " << i << " " << aRouteHydro[i] * (term1 + term2) << " " << term2
997+
// << " " << term3 << " "<<_aEnthalpySource[p][0]<< endl;
975998
}
976999

9771000
double dbed = _pModel->GetSubBasin(p)->GetRiverbedThickness();
@@ -1039,11 +1062,12 @@ void CEnthalpyModel::WriteOutputFileHeaders(const optStruct& Options)
10391062
_STREAMOUT<<name<<" Eout[MJ/m2/d],";
10401063
_STREAMOUT<<name<<" Q_sens[MJ/m2/d],";
10411064
_STREAMOUT<<name<<" Q_cond[MJ/m2/d],";
1042-
_STREAMOUT<<name<<" Q_lat[MJ/m2/d],";
1065+
_STREAMOUT<<name<<" Q_latent[MJ/m2/d],";
10431066
_STREAMOUT<<name<<" Q_GW[MJ/m2/d],";
10441067
_STREAMOUT<<name<<" Q_sw_in[MJ/m2/d],";
10451068
_STREAMOUT<<name<<" Q_lw_in[MJ/m2/d],";
10461069
_STREAMOUT<<name<<" Q_lw_out[MJ/m2/d],";
1070+
_STREAMOUT<<name<<" Q_lateral[MJ/m2/d],";
10471071
_STREAMOUT<<name<<" Q_fric[MJ/m2/d],";
10481072
_STREAMOUT<<name<<" channel storage[MJ/m2],";
10491073
}
@@ -1120,7 +1144,7 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct
11201144
// StreamReachEnergyBalances.csv
11211145
//--------------------------------------------------------------------
11221146
int p;
1123-
double Q_sens,Q_cond,Q_lat,Q_GW,Q_sw_in,Q_lw_in,Q_lw_out,Q_fric,Q_rain,Tave;
1147+
double Q_sens,Q_cond,Q_lat,Q_GW,Q_sw_in,Q_lw_in,Q_lw_out,Q_lateral,Q_fric,Q_rain,Tave;
11241148
double Ein,Eout,mult;
11251149
CSubBasin *pSB;
11261150

@@ -1141,15 +1165,16 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct
11411165

11421166
Ein = 0.5 * (_aMinHist [p][0] + _aMinHist[p][1]);
11431167
Eout = 0.5 * (_aMout_last[p] + _aMout [p][pSB->GetNumSegments() - 1]);
1144-
GetEnergyLossesFromReach(p, Q_sens, Q_cond, Q_lat, Q_GW, Q_sw_in,Q_lw_in,Q_lw_out, Q_fric, Tave);
1168+
GetEnergyLossesFromReach(p, Q_sens, Q_cond, Q_lat, Q_GW, Q_sw_in,Q_lw_in,Q_lw_out, Q_lateral,Q_fric, Tave);
11451169

11461170
if ((pSB->GetTopWidth() < REAL_SMALL) || ( pSB->IsHeadwater())){//running dry or no reach
1147-
_STREAMOUT << ",,,,,,,,,,,";
1171+
_STREAMOUT << ",,,,,,,,,,,,";
11481172
}
11491173
else {
11501174
_STREAMOUT << mult * Ein << "," << mult * Eout << ",";
11511175
_STREAMOUT << mult * Q_sens << "," << mult * Q_cond << "," << mult * Q_lat << ",";
11521176
_STREAMOUT << mult * Q_GW << "," << mult * Q_sw_in << "," << mult * Q_lw_in << "," << mult * Q_lw_out << ",";
1177+
_STREAMOUT << mult * Q_lateral<<",";
11531178
_STREAMOUT << mult * Q_fric << "," << mult * _channel_storage[p] << ",";
11541179
}
11551180
}

‎src/EnergyTransport.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,7 @@ class CEnthalpyModel :public CConstituentModel
5959
double GetWaterTemperature (const double *state_vars, const int iWater) const;
6060

6161
double GetEnergyLossesInTransit(const int p,double &Q_sens,double &Q_GW) const;
62-
double GetEnergyLossesFromReach(const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_GW,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out,double &Q_fric, double &Tave) const;
62+
double GetEnergyLossesFromReach(const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_GW,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out,double &Q_lateral, double &Q_fric, double &Tave) const;
6363
double GetEnergyLossesFromLake (const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out, double &Q_rain) const;
6464

6565
double GetOutflowIceFraction (const int p) const;
@@ -73,6 +73,8 @@ class CEnthalpyModel :public CConstituentModel
7373

7474
//Manipulators (inherited from CConstitModel)
7575
void Initialize (const optStruct& Options);
76+
77+
void PrepareForInCatchmentRouting(const int p);
7678
void PrepareForRouting (const int p);
7779
double ApplyInCatchmentRouting (const int p,
7880
const double *aUnitHydro,

‎src/MassRouting.cpp‎

Lines changed: 29 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -218,7 +218,21 @@ void CConstituentModel::SetLateralInfluxes(const int p,const double Mlat)
218218
void CConstituentModel::PrepareForRouting(const int p) {
219219
//does nothing - interesting in child classes
220220
}
221-
//////////////////////////////////////////////////////////////////
221+
void CConstituentModel::PrepareForInCatchmentRouting(const int p) {
222+
//does nothing - interesting in child classes
223+
}
224+
void CConstituentModel::InCatchmentRoute(const int p, double &Mlat_new, const optStruct &Options)
225+
{
226+
PrepareForInCatchmentRouting(p);
227+
228+
const double * aUnitHydro =_pModel->GetSubBasin(p)->GetUnitHydrograph();
229+
const double * aQlatHist =_pModel->GetSubBasin(p)->GetLatHistory();
230+
Mlat_new=ApplyInCatchmentRouting(p,aUnitHydro,aQlatHist,_aMlatHist[p],_nMlatHist[p], Options.timestep);
231+
232+
_aMlocLast[p] = _aMlocal[p];
233+
_aMlocal [p] = Mlat_new;
234+
}
235+
//////////////////////////////////////////////////////////////////
222236
/// \brief Creates aMoutnew [m^3/s], an array of point measurements for outflow at downstream end of each river segment
223237
/// \details Array represents measurements at the end of the current timestep. if (catchment_routing==distributed),
224238
/// lateral flow Qlat (e.g., runoff, interflow, or baseflow) is routed as part of channel routing routine otherwise,
@@ -240,22 +254,22 @@ void CConstituentModel::RouteMass(const int p, // SB index
240254
const optStruct &Options,
241255
const time_struct &tt) const
242256
{
243-
const double * aUnitHydro =_pModel->GetSubBasin(p)->GetUnitHydrograph();
244-
const double * aRouteHydro=_pModel->GetSubBasin(p)->GetRoutingHydrograph();
245-
const double * aQinHist =_pModel->GetSubBasin(p)->GetInflowHistory();
246-
const double * aQlatHist =_pModel->GetSubBasin(p)->GetLatHistory();
247-
248-
int nSegments =_pModel->GetSubBasin(p)->GetNumSegments();
249-
double seg_fraction=1.0/(double)(nSegments);
250-
251257
//==============================================================
252258
// route from catchment
253259
//==============================================================
254-
Mlat_new=ApplyInCatchmentRouting(p,aUnitHydro,aQlatHist,_aMlatHist[p],_nMlatHist[p], Options.timestep);
260+
//const double * aUnitHydro =_pModel->GetSubBasin(p)->GetUnitHydrograph();
261+
//const double * aQlatHist =_pModel->GetSubBasin(p)->GetLatHistory();
262+
//Mlat_new=ApplyInCatchmentRouting(p,aUnitHydro,aQlatHist,_aMlatHist[p],_nMlatHist[p], Options.timestep);
255263

256264
//==============================================================
257265
// route along channel
258266
//==============================================================
267+
const double * aRouteHydro=_pModel->GetSubBasin(p)->GetRoutingHydrograph();
268+
const double * aQinHist =_pModel->GetSubBasin(p)->GetInflowHistory();
269+
270+
int nSegments =_pModel->GetSubBasin(p)->GetNumSegments();
271+
double seg_fraction=1.0/(double)(nSegments);
272+
259273
if((Options.routing==ROUTE_PLUG_FLOW) || (Options.routing==ROUTE_DIFFUSIVE_WAVE))
260274
{
261275
// (sometimes fancy) convolution
@@ -270,7 +284,7 @@ void CConstituentModel::RouteMass(const int p, // SB index
270284
else {
271285
ExitGracefully("Unrecognized or unsupported constiuent routing method (:Routing command must be ROUTE_NONE, ROUTE_PLUG_FLOW, or ROUTE_DIFFUSIVE_WAVE to support transport)",STUB);
272286
}
273-
287+
274288
//all fluxes from catchment are routed directly to basin outlet (unless handled in convolution routing as source term)
275289
if ((!_lateral_via_convol) || (_pModel->GetSubBasin(p)->IsHeadwater())){
276290
aMout_new[nSegments-1]+=Mlat_new;
@@ -372,8 +386,6 @@ void CConstituentModel::UpdateMassOutflows( const int p,
372386
}
373387
MassOutflow=_aMout[p][pBasin->GetNumSegments()-1];// is now the new mass outflow from the channel
374388

375-
_aMlocLast[p] = _aMlocal[p];
376-
_aMlocal [p] = Mlat_new;
377389

378390
//Update reservoir concentrations
379391
//-----------------------------------------------------
@@ -421,8 +433,10 @@ void CConstituentModel::UpdateMassOutflows( const int p,
421433
dM-=GetNetReachLosses(p);
422434

423435
//mass change from lateral inflows
424-
dM+=0.5*(Mlat_new+_aMlat_last[p])*dt;
425-
436+
if (!_lateral_via_convol){//elsewise, in get net reach losses
437+
dM+=0.5*(Mlat_new+_aMlat_last[p])*dt;
438+
}
439+
426440
_channel_storage[p]+=dM;//[mg] or [MJ]
427441

428442
//Update rivulet storage

0 commit comments

Comments
 (0)