Skip to content

Commit

Permalink
update boundary
Browse files Browse the repository at this point in the history
  • Loading branch information
ftomei committed Dec 2, 2024
1 parent 8807716 commit 293de8a
Showing 1 changed file with 8 additions and 6 deletions.
14 changes: 8 additions & 6 deletions agrolib/soilFluxes3D/boundary.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <algorithm>

#include "physics.h"
#include "commonConstants.h"
Expand All @@ -37,7 +39,7 @@
#include "header/water.h"
#include "header/heat.h"

#include <iostream>


void initializeBoundary(Tboundary *myBoundary, int myType, float slope, float boundaryArea)
{
Expand Down Expand Up @@ -254,15 +256,15 @@ void updateBoundaryWater (double deltaT)
double avgH = (nodeListPtr[i].H + nodeListPtr[i].oldH) * 0.5; // [m]

// Surface water available for runoff [m]
double hs = MAXVALUE(avgH - (nodeListPtr[i].z + nodeListPtr[i].pond), 0.);
double hs = std::max(avgH - (nodeListPtr[i].z + nodeListPtr[i].pond), 0.);
if (hs > EPSILON_METER)
{
double maxFlow = (hs * nodeListPtr[i].volume_area) / deltaT; // [m3 s-1] maximum flow available during the time step
// Manning equation
double v = (1. / nodeListPtr[i].Soil->Roughness) * pow(hs, 2./3.) * sqrt(nodeListPtr[i].boundary->slope);
// on the surface boundaryArea is a side [m]
double flow = nodeListPtr[i].boundary->boundaryArea * hs * v; // [m3 s-1]
nodeListPtr[i].boundary->waterFlow = -MINVALUE(flow, maxFlow);
nodeListPtr[i].boundary->waterFlow = -std::min(flow, maxFlow);
}
}
else if (nodeListPtr[i].boundary->type == BOUNDARY_FREEDRAINAGE)
Expand Down Expand Up @@ -329,7 +331,7 @@ void updateBoundaryWater (double deltaT)
evapFromSoil *= (1. - surfaceWaterFraction);
evapFromSurface *= surfaceWaterFraction;

evapFromSurface = MAXVALUE(evapFromSurface, -waterVolume / deltaT);
evapFromSurface = std::max(evapFromSurface, -waterVolume / deltaT);

if (nodeListPtr[upIndex].boundary != nullptr)
nodeListPtr[upIndex].boundary->waterFlow = evapFromSurface;
Expand All @@ -339,9 +341,9 @@ void updateBoundaryWater (double deltaT)
}

if (evapFromSoil < 0.)
evapFromSoil = MAXVALUE(evapFromSoil, -(theta_from_Se(i) - nodeListPtr[i].Soil->Theta_r) * nodeListPtr[i].volume_area / deltaT);
evapFromSoil = std::max(evapFromSoil, -(theta_from_Se(i) - nodeListPtr[i].Soil->Theta_r) * nodeListPtr[i].volume_area / deltaT);
else
evapFromSoil = MINVALUE(evapFromSoil, (nodeListPtr[i].Soil->Theta_s - nodeListPtr[i].Soil->Theta_r) * nodeListPtr[i].volume_area / deltaT);
evapFromSoil = std::min(evapFromSoil, (nodeListPtr[i].Soil->Theta_s - nodeListPtr[i].Soil->Theta_r) * nodeListPtr[i].volume_area / deltaT);

nodeListPtr[i].boundary->waterFlow = evapFromSoil;
}
Expand Down

0 comments on commit 293de8a

Please sign in to comment.